TALK=f;RUN(1,1)
  DISPLAY
  Problem: one quarter of the plate.
           The plate is compressed by an upward force from below, FY
   Case 1 or 2: 'plane stress', i.e. strz=0
   Case 3 or 4: 'plane strain', i.e. epsz=0
   Case 5 or 6: plate is composed of two materials (along Y-axis) with 
                ten-fold difference in Young's modulus
   Case 7 or 8: plate is composed of two materials (along X-axis) with 
                ten-fold difference in Young's modulus

  PHOTON USE
  p;;;;
  
  

  set prop off
  cl
  msg x-displacement fields
  gr ou z 1
  cont U1 z 1 x 1 21 y 2 21 fil;.0001
  vec z 1 x 2 21 y 2 21 col 0
  pause
  cl
  msg y-displacement fields
  gr ou z 1
  cont V1 z 1 x 2 21 y 1 21 fil;.0001
  vec z 1 x 2 21 y 2 21 col 0
  pause
  cl
  msg epsX fields
  gr ou z 1
  cont EPSX z 1 x 2 21 y 2 21 fil;.0001
  pause
  cl
  msg epsY fields
  gr ou z 1
  cont EPSY z 1 x 2 21 y 2 21 fil;.0001
  pause
  cl
  msg StrX fields
  gr ou z 1
  cont STRX z 1 x 2 21 y 2 21 fil;.0001
  pause
  cl
  msg STRY fields
  gr ou z 1
  cont STRY z 1 x 2 21 y 2 21 fil;.0001
  pause


  ENDDIS


  prt0begin
  print0 1.0e-20 
  prt0end 
 
 ************************************************************
  Group 1. Run Title and Number
 ************************************************************
 TEXT(Plate; 2D[xy])
integer(caseno)     
mesg(caseno = 1 or 2 : z free
mesg(caseno = 3 or 4 : z fixed
mesg(caseno = 5 or 6 : z free two materials (along Y-axis)
mesg(caseno = 7 or 8 : z free two materials (along X-axis)
mesg(caseno = 9 or 10 : z free two materials (along X-axis and Y-axis)
mesg(caseno odd = Stress-Model, even = Dilatation-model

caseno=1
label ask
mesg(caseno=:caseno: Enter 1 .. 8 or blank
readvdu(caseno,int,1)
if(caseno.lt.1) then
 goto ask
endif  
if(caseno.gt.10) then
 goto ask
endif 
caseno
Integer(SolvMod)
  ! SolvMod =  1 for Dilatation-model (caseno = even) 
  !            2 for Stress-model     (caseno = odd )
  !
solvmod=1 + (caseno+1)/2-caseno/2
solvmod
IG(40) = solvmod-1 

libref=caseno 
  Declarations and settings
REAL(FY,LX,LY,POISSON,YOUNG) 
FY= 40.0e6  ! H/m^2 = 40 N/mm^2
LX=90.e-3
         lx=lx*100
LY=120.e-3 
YOUNG   = 1/0.5E-11   ! Young's modulus
POISSON=0.3           ! Poisson's ratio
INTEGER(NXBODY,NYBODY)

 ************************************************************
  Group 2. Time dependence
 STEADY  =    T
 ************************************************************
  Group 3. X-Direction Grid Spacing
 CARTES  =    T
 NXBODY = 20
 NREGX=3                             ! 3 regions
 IREGX=1;GRDPWR(X,1,0.01*LX,1.0)     ! single inner fluid cell
 IREGX=2;GRDPWR(X,NXBODY,LX,1)  
 IREGX=3;GRDPWR(X,1,0.01*LX,1.0)     ! single outer fluid cell 

 ************************************************************
  Group 4. Y-Direction Grid Spacing
 NYBODY = 20
 NREGY=3                             ! 3 regions
 IREGY=1;GRDPWR(Y,1,0.01*LY,1.0)     ! single inner fluid cell
 IREGY=2;GRDPWR(Y,NYBODY,LY,1)  
 IREGY=3;GRDPWR(Y,1,0.01*LY,1.0)     ! single outer fluid cell 
 ************************************************************
  Group 5. Z-Direction Grid Spacing
 NZ=1
 ZWLAST  = 0.001
 ************************************************************
  Group 7. Variables: STOREd,SOLVEd,NAMEd
 ONEPHS  =    T
 SOLVE(P1,V1,U1)
 SOLUTN(P1  ,Y,Y,Y,N,N,N)
 SOLUTN(U1  ,Y,Y,Y,Y,N,Y)
 SOLUTN(V1  ,Y,Y,Y,Y,N,Y)
 
 STORE(PRPS)
 STORE(STRX,STRY,STRZ,STXY)
 STORE(EPSY,EPSX,EPSZ)

 STORE(U1T,V1T,U1/T,V1/T)

 ************************************************************
  GROUP 8. ITERATION NUMBERS ETC
RESFAC=1.e-7
RESREF(V1)=0.0 
RESREF(U1)=0.0  ! to prevent premature exit
LITER(V1) = 50 ! from solver
LITER(U1) = 50 
LITER(P1) = 50 

 ************************************************************
  GROUP 9. PROPERTIES
  
 CSG10='Q1'                  ! materials with various POISSON ratios
  MATFLG=T;NMAT=2         
  160    7800.0    0.3       473.0   43.0      1.0e-5   0.5e-11 
  161    7800.0    0.3       473.0   43.0      1.0e-5   0.5e-10
 
 ************************************************************
  GROUP 11. INITIAL VALUES
fiinit(p1)=0.0
fiinit(u1)=0.0
fiinit(v1)=0.0
 
 FIINIT(PRPS)=0
 PATCH(BODY1,INIVAL,2,NX-1,2,NY-1,1,1,1,1)
 INIT(BODY1,PRPS,FIXVAL,160)
 if(caseno.eq.5.or.caseno.eq.6) then
+  PATCH(BODY2,INIVAL,2,NX-1,ny/2+1,NY-1,1,1,1,1)
+  INIT(BODY2,PRPS,FIXVAL,161)
 endif
 if(caseno.eq.7.or.caseno.eq.8) then
+  PATCH(BODY2,INIVAL,NX/2+1,NX-1,2,NY-1,1,1,1,1)
+  INIT(BODY2,PRPS,FIXVAL,161)
 endif
 if(caseno.eq.9.or.caseno.eq.10) then
+  PATCH(BODY2,INIVAL,2,NX/2+1,2,NY/2,1,1,1,1)
+  INIT(BODY2,PRPS,FIXVAL,161)
+  PATCH(BODY3,INIVAL,NX/2+1,NX-1,NY/2+1,NY-1,1,1,1,1)
+  INIT(BODY3,PRPS,FIXVAL,161)
 endif
 ************************************************************
  GROUP 13. BOUNDARY & SPECIAL SOURCES
 
PATCH(UP,north,2,NX-1,NY-1,NY-1,1,1,1,1) ! top - fixed 
COVAL(UP,V1,FIXVAL,0.0)

PATCH(FORC01,NORTH,2,NX-1,1,1,1,1,1,1) ! bottom - upward force
COVAL(FORC01,V1,FIXFLU,FY)

PATCH(AXESZZ,EAST,1,1,2,NY-1,1,1,1,1)            ! LEFT - fixed
COVAL(AXESZZ,U1,FIXVAL,0.0)

IF(CASENO.EQ.3.OR.CASENO.EQ.4) THEN
SPEDAT(BOUNDARY,ZCONST,R,1.e20)
ELSE
SPEDAT(BOUNDARY,ZCONST,R,0.0)
ENDIF 
 ************************************************************
  GROUP 15. TERMINATE SWEEPS
 LSWEEP  =      200
 ISG21=LSWEEP
 ************************************************************
  GROUP 17. RELAXATION
     #CONPROM
 RELAX(P1  ,LINRLX, 1.000000E+00)
 if(solvmod.eq.1) then
+  spedat(rlxfac,rlxu1d,r,0.5)  
+  spedat(rlxfac,rlxv1d,r,0.5)
 else
+  spedat(rlxfac,rlxu1d,r,1)  
+  spedat(rlxfac,rlxv1d,r,1)
 endif   
 
 ************************************************************
  GROUP 19. DATA TRANSMITTED TO GROUND
 STRA    =    T
 PARSOL  =    F
 ISG52   =    3   !   probe & res

 ************************************************************
  GROUP 23.FIELD PRINT-OUT & PLOT CONTROL
TSTSWP = - 1   ! graphic-mode
NYPRIN = 1
IXMON = NX-2
IYMON = 2
IZMON = 1
      #conprom 

  inform7begin
real(CEPSX,CEPSY,CEPS)
IF(CASENO.EQ.2) THEN
 CEPSX = FY/YOUNG*(1+POISSON)*POISSON
 CEPSY =  -FY/YOUNG*(1-POISSON**2)
ELSE
 CEPSX = FY/YOUNG*POISSON
 CEPSY =  -FY/YOUNG
ENDIF
CEPS = CEPSX+CEPSY
CEPSX
CEPSY 
CEPS
   **** CALCULATE analytical solution ***
(STORED VAR U1T IS :CEPSX:*(XU-0.01*:LX:) with imat>100)
(STORED VAR V1T IS :CEPSY:*(YV -1.01*:LY:) with imat>100)

(STORED VAR U1/T IS U1/(U1T+1.e-20) with imat>100)
(STORED VAR V1/T IS V1/(V1T+1.e-20) with imat>100)

  inform7end


STOP