PHOTON USE
  p;;;;;;

  *rot z ang 90
  set prop off
  do kk=1,4
  gr ou z kk
  gr ou z kk x 4 8 y 1 5
  enddo
  msg If temperature is not being computed in this case, enter /
  msg in answer to the next question about the required variable
  msg
  msg temperature
  do kk=1,4
  con tem1 z kk fi;0.001
  gr ou z kk
  gr ou z kk x 4 8 y 1 5
  enddo
  pause;con off
  cl;red
  msg y-direction strain
  do kk=1,4
  con epsy z kk fi;0.001
  gr ou z kk
  gr ou z kk x 4 8 y 1 5
  enddo
  pause;con off;red
  msg x-direction strain
  do kk=1,4
  con epsx z kk fi;0.001
  gr ou z kk
  gr ou z kk x 4 8 y 1 5
  enddo
  pause;con off;red
  msg y-direction stress
  do kk=1,4
  con stry z kk fi;0.001
  gr ou z kk
  gr ou z kk x 4 8 y 1 5
  enddo
  pause;con off;red
  msg x-direction stress
  do kk=1,4
  con strx z kk fi;0.001
  gr ou z kk
  gr ou z kk x 4 8 y 1 5
  enddo
  pause;con off;red
  msg displacement and velocity vectors
  set prop on
  * set vec ref 2
  do kk=1,4
  vec z kk sh
  set prop off
  * set vec ref 25
  vec z kk x 4 8 y 1 5 sh
  gr ou z kk
  gr ou z kk x 4 8 y 1 5
  enddo
  pause;vec cl
  red
  set prop on
  do kk=1,4
  * set vec ref 2
  vec z kk sh
  enddo
  set prop off
  do kk=1,4
  con p1 z kk 1 5 x 4 8 fil;  0.001
  enddo
  msg 'Pressure' variable in the solid
  ENDUSE
  >>>>>>>>>>>>>>>>>>>>>> Comment begins >>>>>>>>>>>>>>>>>>>>

  DISPLAY
    STRESBOX :  a  collection  of examples of PLANT features
    for stress-in-solid option.

    This 3D domain comprises four variants (  one  for  each
    Z-slab)  of  the  problem  of thermal expansion of solid
    block placed in the cold stream:

    IZ=1 contains  simple rectangular block;
    IZ=2 as above but with automatic implementation of
         boundary conditions;
    IZ=3 staircase cutted piece of the above and
    IZ=4 has got the external pressure loaded block.

    Three first three problems  have  got  the  mathematical
    statements    and    analitical   solution   for   solid
    displacements shown below.
  ENDDIS

    PLANT information :
     * Data input groups used: 13
     * Ground groups planted : 13
     * Headings used  : SORC??
     * Functions used : SISBC
     * Commands used  : None


      CASE STUDY: Thermal expansion of the uniform temperature
                  slab with WEST and SOUTH constrained walls.

                          ny =0.0
                   1       NB
                c+-------------------+
                o|///////////////////|
                n|///////////////////|
                s|///////////////////|
                t|//////// //////////|
             WB r|////// T=1 ////////|EB nx =0.0
                a|//////// //////////|
                i|///////////////////|
                n|///////////////////|
                t+-c-o-n-s-t-r-a-i-n-+
               X = 0        SB       5.0

        Solid properties
        ================
        Young's modulus, E    = 1.0
        Puasson ratio,   P    = 1/3
        Thermoexpansion
        coefficient,     alfa = 1.0

        Solid displacement equations
        ============================
      d2dxx.u + d2dyy.u - ddx.p = 0
      d2dxx.v + d2dyy.v - ddy.p = 0
                ddx.u   + ddy.v = -p/3 + 8/3 T
                              T = const
        Solid boundary conditions
        =========================
      LB : ddx.u=0.333 (4 (T + nx ) + p)
           ddx.v=-ddy.u
               p= 2 (T - nx) - 1.5 ddy.v
      HB :     u=0.0
           ddx.v=0.0
      SB :     v=0.0
           ddy.u=0.0
      NB : ddy.v=0.333 (4 (T + ny) + p)
           ddy.u=-ddx.v
               p= 2 (T - ny) - 1.5 ddx.u

        Analitical solution :
        =====================
      p = -1.3333 (ny + nx )
      u = (1.333 alfa T + 0.889 nx - 0.4444 ny ) x
      v = (1.333 alfa T + 0.889 ny - 0.4444 nx ) y

  Test case
  =========
  Uniform heating.
  Cartesian geometry.
  Uniform grid.
  Fixed stiffness    = 1.0  N/m**2
  Fixed expansivity  = 1.0
  Fixed solid temp.  = 1.0 k
  Poisson ratio      = 0.3333
  Material constrained at east x
  Zero normal stress apllied to material at north y
  begin

  Displacements of solid in x-direction should be:
  IX      3          4          5          6          7
  U1  0.000E+00  1.33333    2.66667     4.0000    5.33333
  Displacements of solid in y-direction should be:
  IY      1          2          3          4          5
  V1  1.33333    2.66667    4.00000   5.333333    6.666666
  end
  <<<<<<<<<<<<<<<<<<<<<<< Comment ends <<<<<<<<<<<<<<<<<<<<<

TEXT(2dxy uniform heating.
  ** Declaration
BOOLEAN(CALSTR)
INTEGER(IYNORT,IXWES,IXEA,IXFIX)
REAL(WIN,POISSN,TIN,TBODY)
REAL(EXCOLI,EXCOC1,EXCOC2,STIFFN,STIFC1,STIFC2,DSTRSE,DSTRSN,DSTRSW)
READQ1=T
STRA=T;CALSTR=T

  ** Setting defaults
   * Inlet velocity, inlet and body temperatures
WIN=1.0;TIN=0.0;TBODY=1.0
   * Constant linear thermal-expansion coefficient, excoli > 0
EXCOLI=1.e-05
   * Constant Young's modulus (stiffness)
STIFFN=1.e5
   * Consant Poiison ratio
POISSN= 0.3333
   * Direct stresses on high, north and low surfaces
DSTRSE=0.0
DSTRSN=0.0
DSTRSW=0.0

  ** Grid settings
NZ=4
IYNORT=5;NY=8;IXWES=4;NX=10;IXEA=NX-2;IXFIX=IXWES-1
ZWLAST=NZ
YVLAST=NY ; XULAST=NX
+    GRDPWR(X,NX,XULAST,1.0)
+    GRDPWR(Y,NY,YVLAST,1.0)
+    GRDPWR(Z,NZ,ZWLAST,1.0)

  ** Solution data
   *  Solve for P1 by whole-field method
SOLUTN(P1,Y,Y,Y,N,N,N)
SOLVE(V1);SOLUTN(V1,Y,Y,N,N,N,Y)
SOLVE(U1);SOLUTN(U1,Y,Y,N,N,N,Y)
SOLVE(TEM1);SOLUTN(TEM1,Y,Y,Y,N,N,Y)
   *  Store variables
STORE(W1)
STORE(HPOR,PRPS,MARK)
OUTPUT(PRPS,N,N,N,N,N,N)
ISOLX=0;ISOLY=0;ISOLZ=0
   *  Terms for temperature
TERMS(TEM1,N,Y,Y,Y,Y,Y)

  **  GROUP 9.  Properties of the medium (or media).
;ENUL=0.01
PRNDTL(TEM1)=CONDFILE

  **  GROUP 11. Initialization of fields of variables
INIADD=F
FIINIT(P1)=1.0;FIINIT(U1)=0.0;FIINIT(V1)=0.0;FIINIT(W1)=0.0
FIINIT(HPOR,MARK)=0.0
   * working fluid is air
FIINIT(PRPS)=0.0
   * Initialize Temperature and density (to air density) Field
FIINIT(TEM1)=TIN
  **  GROUP 13. Boundary conditions and special sources
   * inlet boundary condition, name INLET (at WEST)
PATCH(INLET,WEST,1,1,1,IYNORT,1,NZ,1,LSTEP)
COVAL(INLET,P1,FIXFLU,1.2*WIN);COVAL(INLET,U1,ONLYMS,WIN)
COVAL(INLET,TEM1,ONLYMS,TIN)
   * outlet boundary condition, name EXIT (at EAST)
PATCH(EXIT,EAST,NX,NX,1,NY,1,NZ,1,LSTEP)
COVAL(EXIT,P1,1.0,0.0)
COVAL(EXIT,TEM1,ONLYMS,SAME)
   * wall boundary conditions on north, south, low and high walls
PATCH(NSIDE,NWALL,1,NX,NY,NY,1,NZ,1,LSTEP)
COVAL(NSIDE,U1,1.0,0.0)
PATCH(SSIDE,SWALL,IXWES,IXEA,IYNORT+1,IYNORT+1,1,NZ,1,LSTEP)
COVAL(SSIDE,U1,1.0,0.0)
+ PATCH(ESIDE,WWALL,IXEA+1,IXEA+1,1,IYNORT,1,NZ,1,LSTEP)
+ COVAL(ESIDE,V1,1.0,0.0)
+ PATCH(WSIDE,EWALL,IXFIX,IXFIX,1,IYNORT,1,NZ,1,LSTEP)
+ COVAL(WSIDE,V1,1.0,0.0)

  ** Once-for-all-slabs settings
   * Body properties are those of steel for rectangular
     blocks at IZ=1 and 2.
PATCH(BODZ12,INIVAL,IXWES,IXEA,1,IYNORT,1,2,1,1)
INIT(BODZ12,PRPS,0.0,111.0)
   * Body properties and mark values are those of steel for
     staircase blocks at IZ=3.
PATCH(BODY1,INIVAL,4,8,1,1,1,NZ,1,1)
INIT(BODY1,PRPS,0.0,111.0)
INIT(BODY1,MARK,0.0,111.0)
PATCH(BODY2,INIVAL,4,8,2,2,1,NZ,1,1)
INIT(BODY2,PRPS,0.0,111.0)
INIT(BODY2,MARK,0.0,111.0)
PATCH(BODY3,INIVAL,4,6,3,3,1,NZ,1,1)
INIT(BODY3,PRPS,0.0,111.0)
INIT(BODY3,MARK,0.0,111.0)
PATCH(BODY4,INIVAL,4,5,4,4,1,NZ,1,1)
INIT(BODY4,PRPS,0.0,111.0)
INIT(BODY4,MARK,0.0,111.0)
PATCH(BODY5,INIVAL,4,5,5,5,1,NZ,1,1)
INIT(BODY5,MARK,0.0,111.0)
INIT(BODY5,PRPS,0.0,111.0)
   * Body properties are those of steel for rectangular
     block at IZ=4.
PATCH(BODZ4,INIVAL,IXWES,IXEA,1,IYNORT,NZ,NZ,1,1)
INIT(BODZ4,PRPS,0.0,111.0)
   * Heat-source boundary conditions
PATCH(HOT12,CELL,IXWES,IXEA,1,IYNORT,1,2,1,LSTEP)
COVAL(HOT12,TEM1,FIXVAL,TBODY)
PATCH(HOT4,CELL,IXWES,IXEA,1,IYNORT,NZ,NZ,1,LSTEP)
COVAL(HOT4,TEM1,FIXVAL,TBODY)

  ** Once-for-three-slabs settings
   * fix displacement at ix=IXFIX  to zero
PATCH(FIXED,EAST,IXFIX,IXFIX,1,IYNORT,1,3,1,LSTEP)
COVAL(FIXED,U1,FIXVAL,0.0)

  ** Case 1 (IZ=1) :Thermal expansion of the block
     ---------------------------------------------
   * direct-stress condition on the north face of body
PATCH(NDRSTR,NORTH,IXWES,IXEA,IYNORT,IYNORT,1,1,1,LSTEP)
COVAL(NDRSTR,V1,FIXFLU,DSTRSN*(1.+POISSN))
   * direct-stress condition on the east face of body
PATCH(EDRSTR,EAST,IXEA,IXEA,1,IYNORT,1,1,1,LSTEP)
COVAL(EDRSTR,U1,FIXFLU,DSTRSE*(1.+POISSN))

  ** Case 2 (IZ=2) : As Case 1 but with automatic setting of
     zero-stress boundary conditions
     -------------------------------
NAMSAT=MOSG
    PLANTBEGIN
INTEGER(IYNORT,IXWES,IXEA,IXFIX)
IYNORT=5;IXWES=4;IXEA=NX-2;IXFIX=IXWES-1
   * north face of body
PATCH(NSISBC,NORTH,IXWES,NX,1,NY,2,2,1,LSTEP)
   CO =SISBC(FIXFLU)
   VAL=SISBC(0.0)
COVAL(NSISBC,V1,GRND,GRND)
   * east face of body
PATCH(ESISBC,EAST,IXWES,NX,1,NY,2,2,1,LSTEP)
   CO =SISBC(FIXFLU)
   VAL=SISBC(0.0)
COVAL(ESISBC,U1,GRND,GRND)
  >>>>>>>>>>>>>>>>>>>>>> Comment begins >>>>>>>>>>>>>>>>>>>>
    SISBC function  is  used to introduce zero-normal stress
    boundary conditions at the solid fluid interface.
  <<<<<<<<<<<<<<<<<<<<<<< Comment ends <<<<<<<<<<<<<<<<<<<<<

  ** Case 3 (IZ=3) : Thermal expansion of staircase block
     ----------------------------------------------------
   * Heat-source boundary conditions
 PATCH(SS111H,CELL,1,NX,1,NY,1,NZ,1,LSTEP)
     CO=1.e10
     VAL=:TBODY:
 COVAL(SS111H,TEM1,GRND,GRND)
   * direct-stress condition on east face of body
 PATCH(SS111NOR,NORTH,1,NX,1,NY,NZ,NZ,1,LSTEP)
   CO =SISBC(FIXFLU)
   VAL=SISBC(0.0/FIXFLU)
 COVAL(SS111NOR,V1,GRND,GRND)
   * direct-stress condition on east face of body
 PATCH(SS111EAS,EAST,1,NX,1,NY,NZ,NZ,1,LSTEP)
   CO =SISBC(FIXFLU)
   VAL=SISBC(0.0/FIXFLU)
 COVAL(SS111EAS,U1,GRND,GRND)
  >>>>>>>>>>>>>>>>>>>>>> Comment begins >>>>>>>>>>>>>>>>>>>>
    SISBC function   is   used  to  introduce  automatically
    (guided  by  PRPS  distribition  )  zero-normal   stress
    boundary conditions at the solid 111 -fluid interface as
    indicated by the PATCH name SS111??.
  <<<<<<<<<<<<<<<<<<<<<<< Comment ends <<<<<<<<<<<<<<<<<<<<<

  ** Case 4 (IZ=4) : Pressure loaded block
     -------------------------------------
   * on west face of body
PATCH(WFACE,EAST,IXWES-1,IXWES-1,1,IYNORT,NZ,NZ,1,LSTEP)
    VAL=P1/:STIFFN:
COVAL(WFACE,U1,FIXFLU,GRND)
   * on north face of body
PATCH(NFACE,NORTH,IXWES,IXEA,IYNORT,IYNORT,NZ,NZ,1,LSTEP)
    VAL=-NORTH(P1)/:STIFFN:
COVAL(NFACE,V1,FIXFLU,GRND)
   * on east face of body
PATCH(EFACE,EAST,IXEA,IXEA,1,IYNORT,NZ,NZ,1,LSTEP)
    VAL=-P1[+1,,]/:STIFFN:
COVAL(EFACE,U1,FIXFLU,GRND)
  >>>>>>>>>>>>>>>>>>>>>> Comment begins >>>>>>>>>>>>>>>>>>>>
    The external fluid pressures are used  above  as  normal
    stress conditons at the solid boundaries.
  <<<<<<<<<<<<<<<<<<<<<<< Comment ends <<<<<<<<<<<<<<<<<<<<<
      PLANTEND
   * fix displacement at the west body base  to zero
PATCH(BASE,EAST,IXWES-1,IXEA,1,1,NZ,NZ,1,LSTEP)
COVAL(BASE,U1,FIXVAL,0.)

  ** Special data
SPEDAT(SET,STRAIN,CALSTR,L,:CALSTR:)
SPEDAT(SET,STRAIN,POISSN,R,:POISSN:)
SPEDAT(SET,STRAIN,EXCOLI,R,:EXCOLI:)
SPEDAT(SET,STRAIN,STIFFN,R,:STIFFN:)

  ** Solution criteria
   *  GROUP 16. Termination criteria for inner iterations.
LITER(TEM1)=10;ENDIT(TEM1)=1.E-20
LITER(P1)=20;LITER(U1)=1;LITER(V1)=1;LITER(W1)=1
SELREF=T;RESFAC=1.E-6
LSWEEP=150 ; TSTSWP=-1

   *  GROUP 17. Under-relaxation and related devices.
RELAX(P1,LINRLX,0.8)

  ** Output data
   *  GROUP 21. Frequency and extent of field printout.
       Storage & output, for stress and strain calculation.
IXPRF=3;IXPRL=9;IYPRL=6
+    STORE(EPSY,STRY)
+    OUTPUT(EPSY,Y,N,N,N,N,N) ; OUTPUT(STRY,Y,N,N,N,N,N)
+    FIINIT(EPSY)=0.0;FIINIT(STRY)=0.0
+    STORE(EPSX,STRX)
+    OUTPUT(EPSX,Y,N,N,N,N,N) ; OUTPUT(STRX,Y,N,N,N,N,N)
+    FIINIT(EPSX)=0.0;FIINIT(STRX)=0.0
+    STORE(EPST) ; OUTPUT(EPST,Y,N,N,N,N,N);FIINIT(EPST)=0.0
OUTPUT(P1,Y,Y,Y,Y,Y,Y)
OUTPUT(PRPS,N,N,N,N,N,N)

   *  GROUP 22. Location of spot-value & frequency of
                residual printout.
   *  Assign cell-indices of spot-point monitoring location
IZMON=1;IYMON=IYNORT-1;IXMON=(IXWES+IXEA)/2
nxprin=1;nyprin=1;nzprin=1

   *  GROUP 23. Variable-by-variable field printout and plot
                and/or tabulation of spot-values and residuals.
NPRINT=LSWEEP
dmpstk=t
stop