TALK=f;RUN(1,1)
  photon use
  p;;;;
  
  set prop off
  con u1 z 1 fi;.001
  vec z 1 co 1
  pause
  con v1 z 1 fi;.001
  vec z 1 co 1
  pause
  con p1 z 1 fi;.001
  vec z 1 co 1
  pause
  con u1/t z 1 fi;.001
  vec z 1 co 1
  pause
  con v1/t z 1 fi;.001
  vec z 1 co 1
  pause
  enduse
  DISPLAY
   This file presents data-input sets which correspond to 2D (xy)
   situations capable of analytical representation in terms of 
   polynomial stress functions, as described by Timoshenko and 
   Goodier, 'Theory of Elasticity', Chapter 3, Article 18.
   
   Caseno=1: Plate with uniform shear stress but zero direct stress
             (second-degree polynomial with a2=0, b2=finite, c2=0)
   The boundary conditons are first expressed by prescribing the 
   displacements at the edges; the shear stress is compared with the 
   exact solution.
   Then the same problem is solved by prescribing the shear forces 
   at the boundaries; the computed displacements are compared with
   those of the exact solution
   
   Caseno=2: PLate with linear variation of shear stress and 
   correspondingly non-uniform direct stresses         
      Analytical solution:
   1.   StrY = B3 * Y ; StrXY = - B3 * X; StrX = StrZ = 0   
 
      U = - P/E * B3 * X * Y;
      V = B3/(2*E)* [ Y^2 - (2+P)*X^2]

   2. Boundary Condition for U:
      LEFT  : U = 0;
      RIGHT : free
      UP    : StrXY = - B3*X
      DOWN  : StrXY =   B3*X
      
   3. Boundary Condition for V:
      LEFT  : V(Y=LY/2) = 0;
      RIGHT : StrXY = - B3*LX
      UP    : StrY  =   B3*LY/2
      DOWN  : StrY  =   B3*LY/2

  ENDDIS

 ************************************************************
  Group 1. Run Title and Number
 ************************************************************
TEXT(2D Plate with polynomial stress function; s200
libref=200
 

integer(caseno)
mesg(caseno=1: uniform shear
mesg(caseno=2: linear shear
caseno=1
mesg(caseno = :caseno:. OK ?(Y/n)
readvdu(ans,char,y)
if(:ans:.eq.n) then
caseno=2
endif
caseno


  Declarations and settings
REAL(FY,LX,LY,POISSON,YOUNG) 
FY= 40.0e6  ! H/m^2 = 40 N/mm^2
LX=120.e-3
LY=20.e-3
REAL(B3)
B3 = FY/(LY/2)
 
YOUNG   = 1/0.5E-11   ! Young's modulus
POISSON=0.3           ! Poisson's ratio
INTEGER(NXBODY,NYBODY)
char(formU,formV)
REAL(DLY)
DLY = LY/2
real(CO1,CO2,CO3)
CO1 = -POISSON/YOUNG*B3
CO2 = 0.5/YOUNG*B3
CO3 = (2+POISSON) 
formU=:CO1:*XU*(YG-:DLY:)
formV=:CO2:*((YV-:DLY:)^2-:CO3:*XG^2)

 ************************************************************
  Group 2. Time dependence
 STEADY  =    T
 ************************************************************
  Group 3. X-Direction Grid Spacing
 CARTES  =    T
 NXBODY = 10
 GRDPWR(X,NXBODY,LX,1)  

 ************************************************************
  Group 4. Y-Direction Grid Spacing
 NYBODY = 10
 GRDPWR(Y,NYBODY,LY,1)  
 ************************************************************
  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,N,N,Y)
 SOLUTN(V1  ,Y,Y,Y,N,N,Y)
 
 STORE(PRPS,ENUL,DRH1)
 STORE(STRX,STRY,STRZ,STXY)
 STORE(EPSY,EPSX,EPSZ)
 STORE(U1T,V1T,U1/T,V1/T)

  inform7begin
   **** CALCULATE analytical solution ***
if(caseno.eq.1) then

else

(STORED VAR U1T IS :FORMU:  with swps)
(STORED VAR V1T IS :FORMV:  with swps)

(STORED VAR U1/T IS U1/(U1T+1.e-20))
(STORED VAR V1/T IS V1/(V1T+1.e-20))
 
endif
  inform7end

 ************************************************************
  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
 SPEDAT(SET,MATERIAL,160,L,T)
 CSG10='Q1'                  ! materials with various POISSON ratios
  MATFLG=T;NMAT=1         
  160    7800.0    0.3       473.0   43.0      1.0e-5   0.5e-11 
 
 ************************************************************
  GROUP 11. INITIAL VALUES
fiinit(p1)=0.0
fiinit(u1)=0.0
fiinit(v1)=0.0

  inform11begin
 (INITIAL PRPS  is 160)
  inform11end

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

if(caseno.eq.1) then

PATCH(LEFTV,CELL,2,2,1,NY-1,1,1,1,1)            ! LEFT - V1 fixed at x=2 
COVAL(LEFTV,V1,1.E2,0.0)

PATCH(RIGHTV,CELL,NX,NX,1,NY-1,1,1,1,1)            ! RIGHT - V1 fixed in LY/2 
COVAL(RIGHTV,V1,1.E2,1.E-5)

PATCH(TOPU,CELL,1,NX-1,NY,NY,1,1,1,1)
COVAL(TOPU,U1,1.E2,0.0)

PATCH(BOTU,CELL,1,NX-1,1,1,1,1,1,1)
COVAL(BOTU,U1,1.E2,0.0)

else
 
PATCH(LEFTU,EAST,1,1,1,NY,1,1,1,1)            ! LEFT - U1/V1 fixed from analitical    
PATCH(LEFTV,NORTH,1,1,1,NY-1,1,1,1,1)

 ****** Forces V ******** 
PATCH(UPV,NORTH,1,NX,NY,NY,1,1,1,1)
COVAL(UPV,V1,FIXFLU,FY)

PATCH(DOWNV,SOUTH,1,NX,1,1,1,1,1,1)
COVAL(DOWNV,V1,FIXFLU,FY)

PATCH(RIGHTV,EAST,NX,NX,1,NY-1,1,1,1,1)
COVAL(RIGHTV,V1,FIXFLU,-B3*LX)

 ****** Forces U ******** 
PATCH(UPU,NORTH,1,NX-1,NY,NY,1,1,1,1)

PATCH(DOWNU,SOUTH,1,NX-1,1,1,1,1,1,1)

  inform13begin
(SOURCE of U1 at LEFTU is COVAL(1.e6,:formU:))
(SOURCE of V1 at LEFTV is COVAL(1.e6,:formV:))
(SOURCE of U1 at UPU is COVAL(FIXFLU,:-B3:*XU))
(SOURCE of U1 at DOWNU is COVAL(FIXFLU,:B3:*XU))
  inform13end

endif


SPEDAT(BOUNDARY,ZCONST,R,0.0)   ! free to expand in z direction

 ************************************************************
  GROUP 15. TERMINATE SWEEPS
 LSWEEP  =      1000
 ISG21=LSWEEP
 ************************************************************
  GROUP 17. RELAXATION
 RELAX(P1  ,LINRLX, 1.000000E+00)

 ************************************************************
  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
NXPRIN = 1
IXMON = NX-2
IYMON = 2
IZMON = 1
#maxmin
 #endpause
 #$s003

STOP