*
#$b001
  display
  *
 
     This example Q1 demonstrates the use of 'Set internal
     points of a blok' option for BFC-type grid. This option
     is used to perform a three-dimensional interpolation
     for the internal grid coordinates of a sub-domain; it
     is used in three-dimensional cases when internal
     points of a certain shape need smoothing; that
     shape is then declared as a sub-domain over which
     trans-finite or Laplace solver is used to perform
     interpolation.
     In this example 'Set internal points of a block'
     option is used to distribute grid within a cube
     that has been deformed by a half-cylinder.
CT1=TEXT(Sub-domain_smoothing                B542
#$b002
  *
     Set the grid type to BFC
CT1=BFC=T
#$b002
  *
REAL(DIAMETER,SIDE,RT1)
 
     Enter cylinder diameter (default 0.5m)
READVDU(diameter,real,0.5)
LOOP
LOK=T
     Enter length of cube side (default 1.0m)
READVDU(side,real,1.0)
RT1=DIAMETER+(SIDE/3)
if(side.lt.rt1) then
MESG(Please, cube side must be greater than  :rt1:
LOK=F
endif
if(lok) exit
ENDLOOP
rt1=side/3
INTEGER(NAB,NBC)
 
     Input data used in this example:
MESG(    cylinder diameter = :diameter:
MESG(    cube side       = :side:
  *
 
      -- Create half-cylinder --
     Define points A,B,C and D
MESG(Point A positioned at (0,0,:rt1:)
MESG(GSET(P,A,0,0,:rt1:)
GSET(P,A,0,0,rt1)
MESG(Point B positioned at (0,0,:rt1:+diameter)
MESG(GSET(P,B,0,0,:rt1:+diameter)
GSET(P,B,0,0,RT1+DIAMETER)
MESG(Point C positioned at (side,0,:rt1:+diameter))
MESG(GSET(P,C,side,0,:rt1:+diameter)
GSET(P,C,side,0,rt1+diameter)
MESG(Point D positioned at (side,0,:rt1:)
CT1=GSET(P,D,SIDE,0,:rt1:)
#$b003
  *
     Draw two arc lines(AB & CD) and two straight
     lines (BC & DA) to create a half-cylinder.
     Lines AB&CD must have equal number of cells, as
     well as lines BC&DA.
 
     Enter number of cells for AB(&CD) lines (default 10)
readvdu(nab,int,10)
     Enter power distribution (for symmetric power-law
     precede the number with S. (default 1.0)
readvdu(pwr1,char,1.0)
MESG(Arc AB (A-B) goes through point (0,diameter/2,:rt1:+diameter/2)
MESG(  GSET(L,AB,A,B,nab,pwr1,ARC,0,diameter/2,:rt1:+(diameter/2))
GSET(L,AB,A,B,nab,pwr1,ARC,0,diameter/2,rt1+(diameter/2))
 
     Enter number of cells for BC(&DA) lines (default 10)
readvdu(nbc,int,10)
     Enter power distribution (for symmetric power-law
     precede the number with S. (default 1.0)
readvdu(pwr2,char,1.0)
     GSET(L,BC,B,C,nbc,pwr2)
GSET(L,BC,B,C,nbc,pwr2)
     Arc CD (C-D) goes through point
MESG((side,diameter/2,:rt1:+diameter/2)
     GSET(L,CD,C,D,nab,pwr1,ARC,side,
MESG(diameter/2,:rt1:+(diameter/2))
GSET(L,CD,C,D,nab,pwr1,ARC,side,diameter/2,rt1+(diameter/2))
CT1=GSET(L,DA,D,A,nbc,pwr2)
#$B004
  *
     Define frame ABCD (half-cylinder) which has corner points
     A,B,C,D
CT1=GSET(F,ABCD,A,-,B,-,C,-,D,-)
#$b003
  *
nx=nbc
     Enter number of cells for J direction (default 10)
readvdu(ny,int,10)
it1=nab+4
MESG(Enter number of cells for K direction (default :it1:)
readvdu(nz,int,it1)
if(nz.le.nab+3) then
nz=nab+4
endif
     Set the BFC grid dimension as nx x ny x nz cells with
     cube dimensions (reference length) side x side x side.
CT1=GSET(D,nx,ny,nz,side,side,side)
#$B004
  *
 
     Match grid plane J1 on frame ABCD using the
     trans-finite method
it1=(nz-nab)/2
CT1=GSET(M,ABCD,+K+I,1,1,:it1:+1,TRANS)
#$B004
  *
     Cube grid lines now pass through the cylinder. In
     order to fit the grid to the new shape, the whole
     domain is defined as 'sub-domain block' over which
     Laplace solver is used to make internal grid
     as orthogonal as possible.
CT1=GSET(B,1,nx,1,ny,1,nz,LAP15.TTNNTT)
#$B004
  enddis