TALK=f;RUN(1,1)
  DISPLAY
  A long, hollow, thick-walled cylinder, immersed in an outer fluid, 
  contains a second fluid having a different pressure.
  
  Since conditions vary with radius alone, nx and nz are set = 1 .
  
  Two constraints on its length are considered, namely:
  1. the cylinder length may extend freely; and
  2. the cylinder length may not extend at all.
     
  To enable the computed stresses and strains to be compared with
  theoretical (text-book) values, the latter have been computed.
  
  So as to make comparison easy, variables named P1/T etc have also 
  been computed. They represent the calculated values divided by
  text-book ones. They should all be, and indeed are, close to 1.0;
  but zero is printed when the text-book value is zero
  
  Parameters which may be changed by editing this file include:
  (a) Material properties:  Young's modulus and Poisson's ratio
  (b) Geometric quantities: inner and outer radius, and length.
  (c) Mechanical loadings:  inner and outer pressures
  (d) The grid fineness (via NYCYL)
  (e) The grid-fineness distribution (via POWER)

  The text-book solution may be found in section 14.4 of Benham, 
  Crawford and Armstrong, Mechanics of Engineering Materials; and
  also in Love's Mathematical Theory of Elasticity p144 .
  
  ENDDIS
 ************************************************************
  Group 1. Run Title and Number
 ************************************************************
readq1=t
  READQ1_BEGIN 
 TEXT(Pressurised long cylinder; 1Dy
  Declarations and settings
REAL(PI,PO,RI,RO,POISSON,YOUNG) 
PI=1.e7  ! inner pressure 
PO=1.e6 ! outer pressure
   
    
RI=0.2   ! inner radius
RO=0.8   ! outer radius
YOUNG   = 1/(0.5E-11) ! Young's modulus
POISSON=0.3           ! Poisson's ratio
INTEGER(CASENO,NYCYL) ! nycyl is y-direction cell number in solid
  
  Length-constraint settings
mesga(CASENO=1 means expansion in z allowed         
mesgb(CASENO=2 means expansion in z not allowed    
caseno=1
mesgm(caseno is :caseno: OK? (Y/n)
readvdu(ans,char,Y)
if(:ans:.eq.n) then
caseno=2  
endif
caseno
NYCYL=11
   NYCYL=31
   NYCYL=100

REAL(POWER)  ! used for creating non uniform grid
POWER=2
  READQ1_END 
 ************************************************************
  Group 2. Time dependence
 STEADY  =    T
 ************************************************************
  Group 3. X-Direction Grid Spacing
 CARTES  =    F
 XULAST  = 1.000000E-03
 XFRAC (   1) =  1.000000E+00
 ************************************************************
  Group 4. Y-Direction Grid Spacing
 NREGY=3                             ! 3 regions
 IREGY=1;GRDPWR(Y,1,RI,1.0)          ! single inner fluid cell
 IREGY=2;GRDPWR(Y,NYCYL,RO-RI,POWER) ! ncyl solid cells
 IREGY=3;GRDPWR(Y,1,0.01*RO,1.0)     ! single outer fluid cell 
 ************************************************************
  Group 5. Z-Direction Grid Spacing
 NZ=1
 ZWLAST  = 1.000000E+00
 ************************************************************
  Group 7. Variables: STOREd,SOLVEd,NAMEd
 ONEPHS  =    T
 SOLVE(DISY)
 STORE(P1)
 STORE(P1TH,PRPS,DEN1,ENUL,DVO1,DRH1)    ! TH means theoretical
 STORE(STRX,SXTH,EPSX,EXTH,STRY,SYTH)    ! provide storage for
 STORE(EPSY,EYTH,STRZ,SZTH,EPSZ,EZTH)    ! text-book values
 ************************************************************
  GROUP 8. ITERATION NUMBERS ETC

 ************************************************************
  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 
  161    7800.0    0.2       473.0   43.0      1.0e-5   0.5E-11 
  162    7800.0    0.1       473.0   43.0      1.0e-5   0.5E-11 
  163    7800.0    0.0       473.0   43.0      1.0e-5   0.5E-11 
  164    7800.0    0.4       473.0   43.0      1.0e-5   0.5E-11 
  165    7800.0    0.05      473.0   43.0      1.0e-5   0.5E-11 
  166    7800.0    0.01      473.0   43.0      1.0e-5   0.5E-11 
 
 ************************************************************
  GROUP 11. INITIAL VALUES
FIINIT(PRPS)=0
PATCH(INNER,INIVAL,1,NX,1,1,1,1,1,1)       ! inner fluid
COVAL(INNER,P1,FIXVAL,PI)

PATCH(OUTER,INIVAL,1,NX,NY,NY,1,1,1,1)     ! outer fluid
COVAL(OUTER,P1,FIXVAL,PO)

PATCH(CYLINDER,INIVAL,1,NX,2,NY-1,1,1,1,1) ! solid
INIT(CYLINDER,PRPS,FIXVAL,160)               ! material. Insert other

 ************************************************************
  GROUP 13. BOUNDARY & SPECIAL SOURCES
 
PATCH(IN,CELL,1,NX,1,1,1,1,1,1)            ! pressurizing inner
COVAL(IN,P1,FIXVAL,PI)                       ! fluid
 
PATCH(OUT,CELL,1,NX,NY,NY,1,1,1,1)         ! pressurizing outer
COVAL(OUT,P1,FIXVAL,PO)                      ! fluid             

 SPEDAT(BOUNDARY,XCONST,R,1.E20) ! because total angle is fixed
  
 IF(CASENO.EQ.1) THEN
 SPEDAT(BOUNDARY,ZCONST,R,0.0)   ! because axial motion allowed
 ENDIF
 
 IF(CASENO.EQ.2) THEN
 SPEDAT(BOUNDARY,ZCONST,R,1.E20) ! because axial motion not allowed
 ENDIF
 
 ************************************************************
  GROUP 15. TERMINATE SWEEPS
 LSWEEP  =  200
 ISG21   =   LSWEEP
 ************************************************************
  GROUP 17. RELAXATION
#CONPROM
 ************************************************************
  GROUP 19. DATA TRANSMITTED TO GROUND
 STRA    =    T
 ************************************************************
  GROUP 23.FIELD PRINT-OUT & PLOT CONTROL
output(den1,n,n,n,n,n,n)  
output(drh1,n,n,n,n,n,n)  
output(dvo1,n,n,n,n,n,n)  
output(enul,n,n,n,n,n,n)  
output(prps,n,n,n,n,n,n)  
NYPRIN=1
TSTSWP=-1  
IYPRF=1
IYPRL=NY

ISG52   =    3   !   probe & res


  Theoretical values according to Benham et al
  and Love (Mathematical Theory of Elasticity)
  to be conveyed to earth by way of InForm
  statements.
  
  inform7begin       
  ! Declarations first:
REAL(RODRI, RISQ, ROSQ, RODRISQ, CONI, CONO, TERM)
REAL( AA, BB, EE)
REAL(LL,GG)  ! Lame's constants
  ! deductions
                                  ! Lame's constants Lamda+LL, mu=GG
LL = YOUNG*POISSON/((1+POISSON)*(1-2.0*POISSON)) ! in terms of YOUNG
GG = YOUNG/(2.0*(1.0+POISSON))           ! and POISSON

RISQ    = RI * RI
ROSQ    = RO * RO
RODRISQ = ROSQ / RISQ
  ! settings from formulae in the above text books 
                                               ! all cases
CONI    = PI / (RODRISQ - 1)
CONO    = PO * RODRISQ / (RODRISQ - 1)

(STORED VAR SYTH IS :CONI:*(1.0 - :ROSQ:/ Rg^2 ) - $
:CONO:*(1.0 - :RISQ:/RG^2) with imat>100)

(STORED VAR SXTH IS :CONI:*(1.0 + :ROSQ:/ Rg^2 ) - $
:CONO:*(1.0 + :RISQ:/RG^2) with imat>100)

BB = (PI-PO) * RISQ * ROSQ / ( 2.0*GG * (ROSQ-RISQ) )   ! Love's B

AA = (PI*RISQ-PO*ROSQ) / ( 2.0*(LL+GG) * (ROSQ-RISQ) )  ! Love's A
                                                        ! first part
IF(CASENO.EQ.1) THEN ! ee and aa for no stress
EE = - LL * (PI*RISQ-PO*ROSQ) / ( GG * (3.0*LL+2.0*GG) )! Love's e
ee
EE = EE / (ROSQ-RISQ)                                   ! concluded
ee
AA = AA - EE*LL/( 2.0* (LL + GG))            ! Love's A concluded
(STORED VAR SZTH IS 0.0 WITH IMAT>100)
TERM = - LL / (GG * (3.0*LL + 2.0*GG) )
TERM = TERM * ( PI*RISQ - PO*ROSQ)
TERM = TERM / ( ROSQ - RISQ )
(STORED VAR EZTH IS :TERM: WITH IMAT>100)
ENDIF                                                 

IF(CASENO.EQ.2) THEN           ! fixed end
(STORED VAR EZTH IS 0 WITH IMAT>100)
TERM = LL/(LL+GG)  
TERM = TERM * ( PI*RISQ - PO*ROSQ)
TERM = TERM / ( ROSQ - RISQ )
(STORED VAR SZTH IS :TERM: WITH IMAT>100)
ENDIF

  ! print to screen and satlog.txt                                             
caseno
AA
BB
EE

                             ! all cases
(STORED VAR V1TH IS :AA:*RG + BB/RG)
(STORED VAR EXTH IS (SXTH - :POISSON:*(SYTH+SZTH))/:YOUNG:)
(STORED VAR EYTH IS (SYTH - :POISSON:*(SXTH+SZTH))/:YOUNG:)
(STORED VAR P1TH IS (EXTH + EYTH + EZTH) WITH IMAT>100)
(STORED VAR V1/T IS DISY/V1TH)
(STORED VAR EX/T IS EPSX/EXTH)
(STORED VAR EY/T IS EPSY/EYTH)
IF(CASENO.EQ.1) THEN
(STORED VAR EZ/T IS EPSZ/EZTH)                             
ENDIF
IF(CASENO.EQ.2) THEN         ! cases 1 and 2 distinguished  
(STORED VAR SZ/T IS STRZ/SZTH)                             
ENDIF
  inform7end
STOP