TALK=f;RUN(1,1)
  DISPLAY
  A cubical block is subjected to y-direction tension, under 
  three possible lateral-constraint conditions:
  1. no constraint
  2. lateral displacement in x and z directions prevented 
  3. lateral displacement prevented only in x direction
  
  Because the y-direction displacements depend only on y, NX and NZ, 
  the number of x- and z-direction intervals, are both set to 1.
  
  The boolean variable 'direct' toggles which end is fixed and which
  is pulled.
  
  The stresses and strains are calculated; and their textbook
  values also. the former divided by the latter are printed as: 
    SX/T, EX/T, etc, the values of which should be close to 1.0
  (but they will be printed as 0.0 if the theoretical value 0.0)
  
  Variables which may be varied by hand-editing include:
  * grid-fineness NY ; Poisson's Ratio POISSON ; 
  * Young's Modulus YOUNG ; applied stress APSTR
  * block dimensions XULAST, YVLAST, ZWLAST
  * grid-uniformity POWER (in GRDPWR)
  ENDDIS
 ************************************************************
  Group 1. Run Title and Number
 ************************************************************
 
 TEXT(cube in y-direction tension
integer(caseno)     
mesg(if caseno = 1 .. y-direction tension, x & z free
mesg(if caseno = 2 .. y-direction tension, x & z fixed
mesg(if caseno = 3 .. y-direction tension, x free, z fixed
caseno=3
label ask
mesg(caseno=:caseno: Enter 1, 2, 3 or blank
readvdu(caseno,int,caseno)
if(caseno.lt.1) then
 goto ask
endif  
if(caseno.gt.3) then
 goto ask
endif 
caseno
real(appstr)
appstr = 1.e9     ! applied stress 
appstr
boolean(direct)
direct=t
direct
 ************************************************************
  Group 2. Time dependence
 STEADY  =    T
 ************************************************************
  Group 3. X-Direction Grid Spacing
 CARTES  =    t
 NX      =         1
 ************************************************************
  Group 4. Y-Direction Grid Spacing
 NY      = 10
 YVLAST  = 1.0*NY
 REAL(POWER)
 POWER=0.5
 GRDPWR(Y,NY,YVLAST,POWER)
 ************************************************************
  Group 5. Z-Direction Grid Spacing
 NZ      =         1
 ZWLAST  = 1.000000E+00
 ************************************************************
  Group 7. Variables: STOREd,SOLVEd,NAMEd
 ONEPHS  =    T
 SOLVE(DISY)
 STORE(PRPS,DEN1,ENUL,DVO1,DRH1)
 STORE(STRX,SXTH,STRY,SYTH,STRZ,SZTH,EPSX,EXTH,EPSY,EYTH,EPSZ,EZTH)
 ************************************************************
  Group 8. 
 STRA    =    T  ! To activate stress analysis
 ************************************************************
  GROUP 9. PROPERTIES 
 CSG10='Q1'           ! signal use of the following properties line
                      ! which correspond to steel
  MATFLG=T;NMAT=1
  160    7800.0    0.3       473.0   43.0    0.37e-5   0.5E-11 
REAL(YOUNG, POISSON)
POISSON = 0.3          ! must conform with matflg value
YOUNG   = 1/(0.5E-11)  ! must conform with matflg value 
 ************************************************************
  GROUP 11. INITIAL VALUES
 FIINIT(PRPS)=160
 FIINIT(DISY)=0.0
 
 ************************************************************
  GROUP 13. BOUNDARY & SPECIAL SOURCES
 
 IF(DIRECT) THEN
 PATCH(SOU,SWALL,1,1,1,1,1,1,1,1)
 COVAL(SOU,DISY,1.0, 0.0)
 
 PATCH(NFORCE,NORTH,1,1,NY,NY,1,1,1,1)
 COVAL(NFORCE,DISY,FIXFLU,APPSTR)
 
 ELSE
 
 PATCH(NOR,NWALL,1,1,NY,NY,1,1,1,1)
 COVAL(NOR,DISY,1, 0.0)
 
 PATCH(NFORCE,NORT,1,1,1,1,1,1,1,1)
 COVAL(NFORCE,DISY,FIXFLU,-APPSTR)
 ENDIF
 ************************************************************
  GROUP 15. TERMINATE SWEEPS
 RESREF(V1)=0.0  
 LSWEEP  = 100
 ISG21   = LSWEEP   ! to ensure at least isg21 sweeps
 ************************************************************
  GROUP 17. RELAXATION
 ************************************************************
  GROUP 23.FIELD PRINT-OUT & PLOT CONTROL
TSTSWP=-1  
#MAXMIN
#ENDPAUSE
REAL(P1TH,EXTH,EYTH,EZTH,SXTH,SYTH,SZTH,TERM)
  Formulae for the theoretical values
  
IF(CASENO.EQ.1) THEN
SYTH=APPSTR
EYTH=APPSTR/YOUNG
EXTH=-EYTH*POISSON
EZTH=-EYTH*POISSON
SXTH=0.0
SZTH=0.0
ENDIF

IF(CASENO.EQ.2) THEN
EXTH=0.0
EZTH=0.0
TERM=APPSTR
EYTH=(APPSTR/YOUNG)*(1+POISSON)*(1-2*POISSON)/(1-POISSON)    
SYTH=TERM 
SXTH=TERM*POISSON/(1-POISSON)    
SZTH=TERM*POISSON/(1-POISSON)    
ENDIF

IF(CASENO.EQ.3) THEN
EZTH=0
SXTH=0
SYTH=APPSTR
EYTH=(SYTH/YOUNG)*(1+POISSON)*(1-POISSON)
EXTH=-EYTH*POISSON/(1-POISSON)
P1TH=EXTH+EYTH
SZTH=P1TH*YOUNG*POISSON/(1-2*POISSON)/(1+POISSON)
ENDIF

  Effect of constraints in x and z directions
IF(CASENO.EQ.2) THEN
SPEDAT(BOUNDARY,XCONST,R,1.E20)  ! prevent x displacement
SPEDAT(BOUNDARY,ZCONST,R,1.E20)  ! prevent z displacement
ENDIF
IF(CASENO.EQ.3) THEN
SPEDAT(BOUNDARY,ZCONST,R,1.E20)  ! prevent z displacement
ENDIF
#$S001
STOP