TALK=f;RUN(1,1)
  DISPLAY
  Problem: The quarter of the square channel.
           Inner pressure = 4.0e8, outer pressure = 0
   1. 2D 'PLANE-STRAIN' : EpsZ = 0
   2. Comparison - a calculation on code ELCUT (www.tor.ru)
  ENDDIS  
  PHOTON USE
  p;;;;
  
  
  set prop off
  msg
  msg V1 contours
  msg
  cont DELY z 1 x 1 20 y 2 21 fil;.0001
  pause
  
  cl
  msg
  msg U1 contours
  msg
  cont DELX z 1 x 1 20 y 2 21 fil;.0001
  pause

  cl
  msg
  msg EpsX contours
  msg
  cont EPSX z 1 x 1 20 y 2 21 fil;.0001
  pause
  
  cl
  msg
  msg EpsY contours
  msg
  cont EPSY z 1 x 1 20 y 2 21 fil;.0001
  pause

  cl
  msg
  msg EpsZ contours
  msg
  cont EPSZ z 1 x 1 20 y 2 21 fil;.0001
  pause

  cl
  msg
  msg STRX contours
  msg
  cont STRX z 1 x 1 20 y 2 21 fil;.0001
  pause
    
  cl
  msg
  msg STRY contours
  msg
  cont STRY z 1 x 1 20 y 2 21 fil;.0001
  pause

  cl
  msg
  msg STXY contours
  msg
  cont STXY z 1 x 1 20 y 2 21 fil;.0001
  pause

  cl
  msg Displacements fields
  gr ou z 1
  cont DELT z 1 x 1 20 y 2 21 fil;.0001
  vec z 1 x 1 20 y 2 21 col 0


  ENDUSE

  
 ************************************************************
  Group 1. Run Title and Number
 ************************************************************
 TEXT(square channel;2D[xy])
 
  Declarations and settings
REAL(PIN,LX,LY,R0,POISSON,YOUNG) 
PIN=4.0e8  ! H/m^2
LX=60.e-3
LY=60.e-3 
R0=30.e-3
YOUNG   = 1/0.5E-11   ! Young's modulus
POISSON=0.3           ! Poisson's ratio
INTEGER(NXC,NYC)
INTEGER(CASENO)       ! start of menu
mesg(caseno 1 or 2 : default settings
mesg(caseno 3 or 4 : Poisson's ratio = 0.0 ; should affect result
mesg(caseno 5 or 6 : zwlast = 1.0 m ; unlikely to affect result
mesg(caseno 7 or 8 : focus print-out on channel corner
mesg(caseno 9 or 10 : outside fluid pressure = inside pressure
mesga(caseno odd = simultaneous, even = Stone solver
caseno=1
mesgm(caseno = :caseno: Enter another if not OK
readvdu(caseno, int, 1)
caseno 
libref=caseno     ! useful because this appears in graphical monitor
Integer(SolvMod)
  ! SolvMod =  1 for Stone Solver        (caseno = even) 
  !            2 for simultaneous Solver (caseno = odd )
  !
solvmod=1 + (caseno+1)/2-caseno/2
solvmod

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

 ************************************************************
  Group 4. Y-Direction Grid Spacing
 NYC = 10
 NREGY=3
 IREGY=1;GRDPWR(Y,1,0.01*LY,1.0)     ! single outer fluid cell 
 IREGY=2;GRDPWR(Y,NYC,LY-R0,1.0)
 IREGy=3;GRDPWR(Y,NYC,R0,1)  
 ************************************************************
  Group 5. Z-Direction Grid Spacing
 NZ=1
 ZWLAST  = 0.001
 if(caseno.eq.5.or.caseno.eq.6) then
 ZWLAST  = 1.0
 endif
 zwlast
 ************************************************************
  Group 7. Variables: STOREd,SOLVEd,NAMEd
 ONEPHS  =    T
 SOLVE(P1,V1,U1)
 STORE(PRPS)
 STORE(STRX,STRY,STRZ)
 STORE(EPSY,EPSX,EPSZ)

 STORE(STXY)
 STORE(Delt,DelX,DelY)

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

 ************************************************************
  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.0       473.0   43.0      1.0e-5   0.5e-11 
INTEGER(IPRPS)
IPRPS=160         ! used in initial-condition setting below
IF(CASENO.EQ.3.OR.CASENO.EQ.4) THEN
 IPRPS=161
 POISSON=0.0      ! for use in exact-solution calculation
ENDIF
 
 ************************************************************
  GROUP 11. INITIAL VALUES
fiinit(p1)=0.0
fiinit(u1)=0.0
fiinit(v1)=0.0
 
fiinit(Delt)=0.0


 FIINIT(PRPS)=IPRPS
 PATCH(INNER,INIVAL,1,NXC,NYC+2,NY,1,1,1,1)
 INIT(INNER,PRPS,FIXVAL,0)
 PATCH(OUTER1,INIVAL,1,NX,1,1,1,1,1,1)
 INIT(OUTER1,PRPS,FIXVAL,0)
 PATCH(OUTER2,INIVAL,NX,NX,1,NY,1,1,1,1)
 INIT(OUTER2,PRPS,FIXVAL,0)

 ************************************************************
  GROUP 13. BOUNDARY & SPECIAL SOURCES
 
PATCH(UP,north,NXC+1,NX-1,NY-1,NY-1,1,1,1,1)            ! UP - fixed 
COVAL(UP,V1,FIXVAL,0.0)

PATCH(LEFT,WEST,1,1,2,NYC+1,1,1,1,1)            ! LEFT - fixed
COVAL(LEFT,U1,FIXVAL,0.0)

PATCH(IN,CELL,1,NXC,NYC+2,NY,1,1,1,1)       ! pressurizing inner
COVAL(IN,P1,FIXVAL,PIN)                       ! fluid

PATCH(UIN,CELL,1,NXC-1,NYC+2,NY,1,1,1,1)
COVAL(UIN,U1,FIXVAL,0.0)                 
PATCH(VIN,CELL,1,NXC,NYC+2,NY,1,1,1,1)
COVAL(VIN,V1,FIXVAL,0.0)                 



PATCH(OUT1,CELL,1,NX,1,1,1,1,1,1)
INIT(OUT1,P1,FIXVAL,0)

PATCH(OUT2,CELL,NX,NX,2,NY,1,1,1,1)
INIT(OUT2,P1,FIXVAL,0)
  
if(caseno.eq.9.or.caseno.eq.10) then
 INIT(OUT1,P1,FIXVAL,pin)
 coval(OUT1,P1,FIXVAL,pin)
 INIT(OUT2,P1,FIXVAL,pin)
 coval(OUT2,P1,FIXVAL,pin)
endif 


     ! EpsZ = 0

 SPEDAT(BOUNDARY,ZCONST,R,1.0e20)
 
 ************************************************************
  GROUP 15. TERMINATE SWEEPS
 LSWEEP  =       300
 ISG21=LSWEEP
if(solvmod.eq.2) then
 LG(40) = T  ! sisSOlver
else
 lsweep=lsweep*2
endif
 ************************************************************
  GROUP 17. RELAXATION
     #CONPROM
 RELAX(P1  ,LINRLX, 1.000000E+00)
 relax(U1,linrlx,1.0)  
 relax(V1,linrlx,1.0)  

 ************************************************************
  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
  look at bottom right-hand corner of solid
nxprin=1;nyprin=1
ixprf=nx-4;ixprl=nx
iyprf=1;iyprl=5

if(caseno.eq.7.or.caseno.eq.8) then
 look at bottom right-hand corner of fluid
ixprf=nxc-2;ixprl=nxc+2                  
iyprf=nyc-3;iyprl=nyc+1      
endif  
  inform7begin
  cell-centre displacements, for comparison with ELCUT
real(R4)
(STORED VAR DelX IS 0.5*(U1+U1[-1]) with imat>100)
(STORED VAR DelY IS 0.5*(V1+V1[&-1]) with imat>100)
(STORED VAR Delt IS SQRT(DELX^2+DELY^2) with imat>100)

  inform7end


STOP