TALK=f;RUN(1,1)
  DISPLAY
  A rectangular plate with a centrally-placed CIRCULAR hole is
  extended by uniform equal forces on its left and right
  surfaces, which are allowed to distort.
  By reason of symmetry, only one quarter of the plate is analysed.
  
  Comparison - a analytical solution (Timoshenko and Goodier).
  See also: I. Demirdzic, S. Musafija 
  "Finite Volume method for stress analysis in Complex Domains"
  Int. J. for numerical methods in engineering",
  v. 37, 3751-3756 (1994).)  

  ENDDIS
  
  PHOTON USE
  p;;;;
  
  
  set prop off
  msg
  msg PRPS contours
  msg
  cont PRPS z 1 fil;.0001
  gr Z 1
  pause
  
  cl
  msg
  msg V1 contours
  msg
  cont V1 z 1 fil;.0001
  vec z 1 col 0
  pause
  
  cl
  msg
  msg U1 contours
  msg
  cont U1 z 1 fil;.0001
  vec z 1 col 0
  pause

  cl
  msg
  msg STRX contours
  msg
  cont STRX z 1 fil;.0001
  pause
    
  cl
  msg
  msg SXTH contours
  msg
  cont SXTH z 1 fil;.0001
  pause

  cl
  msg
  msg SX/T contours
  msg
  cont SX/T z 1 fil;.0001
  pause

  cl
  msg
  msg STRY contours
  msg
  cont STRY z 1 fil;.0001
  pause

  cl
  msg
  msg SYTH contours
  msg
  cont SYTH z 1 fil;.0001
  pause

  cl
  msg
  msg SY/T contours
  msg
  cont SY/T z 1 fil;.0001
  pause

  cl
  msg
  msg STXY contours
  msg
  cont STXY z 1 fil;.0001
  pause

  cl
  msg
  msg SXYT contours
  msg
  cont SXYT z 1 fil;.0001
  pause
  
  cl
  msg
  msg XY/T contours
  msg
  cont XY/T z 1 fil;.0001
  pause
  
  ENDUSE

  
 ************************************************************
  Group 1. Run Title and Number
 ************************************************************
 TEXT(2D xy Plate with CIRCULAR hole; s202
libref=202
TITLE 


integer(caseno)
caseno = 1   

 
  Declarations and settings
REAL(FX,LX,LY,R0,POISSON,YOUNG) 
FX= 1.0e4  
LX=2.0
LY=LX 
R0=0.5
YOUNG   = 1/1.0E-7   ! Young's modulus
POISSON=0.3           ! Poisson's ratio
INTEGER(NXBDY1,NYBDY1)
INTEGER(NXBDY2,NYBDY2)

 ************************************************************
  Group 2. Time dependence
 STEADY  =    T
 ************************************************************
  Group 3. X-Direction Grid Spacing
 CARTES  =    T
 NXBDY1 = 30
 NXBDY2 = 30
 NREGX=2
 IREGX=1;GRDPWR(X,NXBDY1,1.1*R0,1)  
 IREGX=2;GRDPWR(X,NXBDY2,1.1*R0-LX,1.08)  
 ************************************************************
  Group 4. Y-Direction Grid Spacing
 NYBDY1 = 30
 NYBDY2 = 30
 NREGY=2
 IREGY=1;GRDPWR(Y,NYBDY1,1.1*R0,1)  
 IREGY=2;GRDPWR(Y,NYBDY2,1.1*R0-LY,1.08)  
 ************************************************************
  Group 5. Z-Direction Grid Spacing
 NZ=1
 ZWLAST  = 0.001
 ************************************************************
  Group 7. Variables: STOREd,SOLVEd,NAMEd
 ONEPHS  =    T
 SOLVE(P1,V1,U1)
 STORE(PRPS)
 STORE(STRX,STRY,STRZ,STXY)
 STORE(EPSY,EPSX,EPSZ)

 STORE(R7,TET7,SXTH,SX/T,SYTH,SY/T)
 STORE(SXYT,XY/T)
 
 ************************************************************
  GROUP 8. ITERATION NUMBERS ETC
RESFAC = 1.e-7    ! 
RESREF(U1)=-1 
RESREF(V1)=-1  ! to prevent premature exit from solver
RESREF(P1)=-1 

LITER(V1) = 50
LITER(U1) = 50 
LITER(P1) = 20 

 ************************************************************
  GROUP 9. PROPERTIES
  
 CSG10='Q1'                  ! material properties
  MATFLG=T;NMAT=2         
    1    7800.0    0.3       473.0   43.0      1.0e-5   1.0e-7 
  160    7800.0    0.3       473.0   43.0      1.0e-5   1.0e-7 
 
 ************************************************************
  GROUP 11. INITIAL VALUES
 fiinit(p1)=0.0
 fiinit(u1)=0.0
 fiinit(v1)=0.0

 fiinit(PRPS)=160
 (INITIAL of PRPS is 1 with IF(SQRT(XG^2+YG^2).LE.:R0:))
 SPEDAT(SET,MATERIAL,  1,L,T)

 ************************************************************
  GROUP 13. BOUNDARY & SPECIAL SOURCES



PATCH(LEFT,WWALL,1,1,1,NY,1,1,1,1)         ! Left-end - fixed 
COVAL(Left,U1,1,0.0)

PATCH(Bot,SWALL,1,NX,1,1,1,1,1,1)         ! Bottom-end - fixed 
COVAL(Bot,V1,1,0.0)

PATCH(dmn,cell,1,NX,1,ny,1,1,1,1)
COVAL(DMN,p1,FIXVAL,0)

if(caseno.eq.1) then   !   stress on boundary
 real(r02)
 r02 = R0**2
 ***** TOP : V - normal ******
 char(fV1,fV2,TV1,RV1)
 RV1=(:R02:/(:LY:^2+XG^2))
 TV1=ATAN(:LY:/XG)
 fV1=:FX:*(-:RV1:*(0.5*COS(2*:TV1:)-
 fV2=COS(4*:TV1:))-1.5*:RV1:^2*COS(4*:TV1:))
 PATCH(TopV,NORTH,1,NX,NY,NY,1,1,1,1)
 (source of V1 at TOPV is COVAL(FIXFLU,:fV1::fV2:))
 
 ***** TOP : U - tangential ******
 char(fU1,fU2,TU1,RU1)
 RU1=(:R02:/(:LY:^2+XU^2))
 TU1=ATAN(:LY:/XU)
 fU1=:FX:*(-:RU1:*(0.5*SIN(2*:TU1:)+
 fU2=SIN(4*:TU1:))+1.5*:RU1:^2*SIN(4*:TU1:))
 PATCH(TopU,NORTH,1,NX-1,NY,NY,1,1,1,1)
 (source of U1 at TOPU is COVAL(FIXFLU,:fU1::fU2:))

 ***** RIGHT : U - normal ******
 char(fU3,fU4,TU2,RU2)
 RU2=(:R02:/(:LX:^2+YG^2))
 TU2=ATAN(YG/:LX:)
 fU3=:FX:*(1-:RU2:*(1.5*COS(2*:TU2:)+
 fU4=COS(4*:TU2:))+1.5*:RU2:^2*COS(4*:TU2:))
 PATCH(RIGHTU,EAST,NX,NX,1,NY,1,1,1,1)
 (source of U1 at RIGHTU is COVAL(FIXFLU,:fU3::fU4:))

 ***** RIGHT : V - tangential ******
 char(fV3,fV4,TV2,RV2)
 RV2=(:R02:/(:LX:^2+YV^2))
 TV2=ATAN(YV/:LX:)
 fV3=:FX:*(-:RV2:*(0.5*SIN(2*:TV2:)+
 fV4=SIN(4*:TV2:))+1.5*:RV2:^2*SIN(4*:TV2:))
 PATCH(RIGHTV,EAST,NX,NX,1,NY-1,1,1,1,1)
 (source of V1 at RIGHTV is COVAL(FIXFLU,:fV3::fV4:))
 
else
 PATCH(FORY,EAST,NX,NX,1,NY,1,1,1,1)        ! Right-end - compression
 COVAL(FORY,U1,FIXFLU,FX)
endif

     ! PLANE-STRAIN,   STRZ = 0  

     SPEDAT(BOUNDARY,ZCONST,R,1.e20)
 
 ************************************************************
  GROUP 15. TERMINATE SWEEPS
 LSWEEP  =   400
 ISG21=LSWEEP
 ************************************************************
  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

 ************************************************************
  GROUP 23.FIELD PRINT-OUT & PLOT CONTROL
TSTSWP = - 1   ! graphic-mode
IXMON = NX-2
IYMON = 2
IZMON = 1
NYPRIN = 4
NXPRIN = 3
#endpause
#maxmin
#$s003

  inform7begin
char(form1,form2)
(STORED VAR R7 IS SQRT(XG^2+YG^2))
(STORED VAR TET7 IS ATAN(YG/(XG+1.e-10)))

form1=:FX:*(1.0-(:R0:/R7)^2*(1.5*COS(2*TET7)+
form2=COS(4*TET7))+1.5*(:R0:/R7)^4*COS(4*TET7))
(STORED VAR SXTH IS :form1::form2: with imat>100)
  (STORED VAR SX-T IS STRX-SXTH with imat>100)
(STORED VAR SX/T IS STRX/SXTH with imat>100)

char(form3,form4)
form3=:FX:*(-(:R0:/R7)^2*(0.5*COS(2*TET7)-
form4=COS(4*TET7))-1.5*(:R0:/R7)^4*COS(4*TET7))
(STORED VAR SYTH IS :form3::form4: with imat>100)
  (STORED VAR SY-T IS STRY-SYTH with imat>100)
(STORED VAR SY/T IS STRY/SYTH with imat>100)
  inform7end

char(ff1,ff2,ff3,Tf1,Rf1)
real(rf2)
rf2=R0**2
Rf1=(:Rf2:/(YV^2+XU^2))
Tf1=ATAN(YV/XU)
ff1=:FX:*(-:Rf1:*(0.5*SIN(2*:Tf1:)+
ff2=SIN(4*:Tf1:))+
ff3=1.5*:Rf1:^2*SIN(4*:Tf1:))
(STORED VAR SXYT IS :ff1::ff2::ff3: with imat>100)
  (STORED VAR XY-T IS STXY-SXYT with imat>100)
(STORED VAR XY/T IS STXY/SXYT with imat>100)


  
       restrt(all)
  
STOP