TALK=f;RUN(1,1)
  DISPLAY
  Problem: Vertical prismatic (0.3x03x2.5 m) bar & gravity force 
           -Rho*g.
           At the upper end - a uniform normal stress equal to 
           total weight of the bar Rho*g*LZ.
   1. Comparison - a analytical solution:
      U = -P*Rho*g*x*z/E,
      V = -P*Rho*g*z*y/E,
      W = Rho*g/E/2*[ (z^2 - LZ^2) + P*(x^2+y^2)]

  ENDDIS

 
 ************************************************************
  Group 1. Run Title and Number
 ************************************************************
TEXT(3D Prismatic bar & gravity force; s301)
libref=301
  Declarations and settings
REAL(Rho,GF,FZ,LX,LY,LZ,POISSON,YOUNG)
Rho = 1.e4
GF=9.8 
FZ= Rho*GF
LX=0.15
LZ=2.5
LY=0.15 
YOUNG   = 1/0.5E-11   ! Young's modulus
POISSON=0.3           ! Poisson's ratio
INTEGER(NXBODY,NYBODY,NZBODY)

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

 ************************************************************
  Group 4. Y-Direction Grid Spacing
 NYBODY = 6
 GRDPWR(Y,NYBODY,LY,1)  
 ************************************************************
  Group 5. Z-Direction Grid Spacing
 NZBODY = 6
 GRDPWR(Z,NZBODY,LZ,1)  
 ************************************************************
  Group 7. Variables: STOREd,SOLVEd,NAMEd
 ONEPHS  =    T
 SOLVE(P1,V1,U1,W1)
 SOLUTN(P1  ,Y,Y,Y,N,N,N)
 SOLUTN(U1  ,Y,Y,Y,N,N,Y)
 SOLUTN(V1  ,Y,Y,Y,N,N,Y)
 SOLUTN(W1  ,Y,Y,Y,N,N,Y)
 
 STORE(PRPS)
 STORE(DRH1,VISL) 

     STORE(STRZ)  
 STORE(EPSY,EPSX,EPSZ)
 STORE(U1T,V1T,W1T,U1/T,V1/T,W1/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 
RESREF(W1)=0.0 
LITER(W1) = 50 

 ************************************************************
  GROUP 9. PROPERTIES
  
 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
fiinit(w1)=0.0
 
 FIINIT(PRPS)=160
 ************************************************************
  GROUP 13. BOUNDARY & SPECIAL SOURCES
 
PATCH(FORS01,HIGH,1,NX,1,NY,NZ,NZ,1,1)            ! Upper - sprain
COVAL(FORS01,W1,FIXFLU,FZ*LZ)

PATCH(FORCGRA,VOLUME,1,NX,1,NY,1,NZ,1,1)          ! Volume - gravity force
COVAL(FORCGRA,W1,FIXFLU,-FZ)

 ***** analytical solution *****
real(CC00,CC11)
CC00 = FZ/YOUNG*POISSON/2
char(FrZ1,FrZ2)
CC11 = FZ/YOUNG/2
FrZ1=:CC11:*(ZW^2-:LZ:^2)+
FrZ2=:CC00:*((XG-0.5*:LX:)^2+(YG-0.5*:LY:)^2)

real(CC22)
CC22 = -FZ/YOUNG*POISSON
char(FrX,FrY)
FrX=:CC22:*ZG*(XU-0.5*:LX:)
FrY=:CC22:*ZG*(YV-0.5*:LY:)


PATCH(ONEW,CELL,1,nx,1,ny,1,1,1,1)                           ! Fixed W1
(source of W1 at ONEW is COVAL(1.0e5,:FrZ1::FrZ2:))

PATCH(ONEU,CELL,1,nx-1,1,ny,1,1,1,1)                           ! Fixed U1/V1
(source of U1 at ONEU is COVAL(1.0e5,:FRX:))
PATCH(ONEV,CELL,1,nx,1,ny-1,1,1,1,1)                           ! Fixed U1/V1
(source of V1 at ONEV is COVAL(1.0e5,:FRY:))

 ************************************************************
  GROUP 15. TERMINATE SWEEPS
 LSWEEP  =      2000
 ISG21=LSWEEP
 ************************************************************
  GROUP 17. RELAXATION
     #CONPROM
 RELAX(P1  ,LINRLX, 1.000000E+00)
     spedat(rlxfac,rlxu1d,r,0.5)  
     spedat(rlxfac,rlxv1d,r,0.5)  
     spedat(rlxfac,rlxw1d,r,0.5)  
 ************************************************************
  GROUP 19. DATA TRANSMITTED TO GROUND
 STRA    =    T
 PARSOL  =    F
      ISG52   =    3   !   probe & res
#maxmin

 ************************************************************
  GROUP 23.FIELD PRINT-OUT & PLOT CONTROL
TSTSWP = - 1   ! graphic-mode
NYPRIN = 1
NXPRIN = 1
NZPRIN = 1
IXMON = NX-2
IYMON = NY-2
IZMON = NZ-2

output(prps,n,n,n,n,n,n) 
output(drh1,n,n,n,n,n,n) 
output(visl,n,n,n,n,n,n) 
    output(epsZ,n,n,n,n,n,n) 

      #conprom 

  inform7begin
(STORED VAR W1T IS :FrZ1::FrZ2:)

(STORED VAR U1T IS :frX:)
(STORED VAR V1T IS :frY:)

(STORED VAR U1/T IS U1/(U1T+1.e-20))
(STORED VAR V1/T IS V1/(V1T+1.e-20))
(STORED VAR W1/T IS W1/(W1T+1.e-20))
  inform7end
     restrt(all) 
STOP