TALK=f;RUN(1,1)
  DISPLAY
   This is a 2D (xy) Plane-Strain problem, in which a beam is bent 
   by a torqe applied at its right-hand end.
   
   The problem has an exact analytical solution which is used to 
   enable the accuracy of the computed values to be assessed by 
   being divided by the exact values, and printed as U1/T, etc , 
   which, if equal to 1.0 signify correctness. [for casenos 1-11
   only at present.]
   
   A menu is provided which enables:
   * the Poisson's ratio to be varied
   * the z-dimension to be varied
   * the beam to be made much more slender
   * the end torque to be replaced by a tension force
   * point-wise and distributed forces to be applied
   * temperature gradients to be applied in x and y directions
   * the grid-fineness and distribution to be varied.
   
  In order to ensure accuracy, the grid distribution in beam-
  thickness is such as to provide one very thin cell near each beam
  surface.  
  ENDDIS 
  PHOTON USE
  p;;;;
  
  

  set prop off
  cl
  msg x-direction-displacement field
  gr ou z 1
  cont U1 z 1 fil;.0001
  vec z 1 y 2 22 col 0
  pause
   
  cl
  msg y-direction-displacement field
  gr ou z 1
  cont V1 z 1  fil;.0001
  vec z 1  col 0
  pause;cl
  
  msg v1/v1_exact field
  gr ou z 1
  cont V1/t z 1 fi;.0001
  pause;cl
  msg u1/u1_exact field
  gr ou z 1
  cont u1/t z 1  fil;.0001
  pause;cl
  msg p1/p1_exact field
  gr ou z 1
  cont p1/t z 1  fil;.0001
   
  ENDUSE

 ************************************************************
  Group 1. Run Title and Number
 ************************************************************
 
  Declarations and settings
REAL(LX0,LY0) 
LX0=120.e-3
LY0=40.e-3
INTEGER(CASENO)       ! start of menu
mesg(caseno 1 : default settings
mesg(caseno 3 : Poisson's ratio = 0.0
mesg(caseno 5 : zwlast = 1.0 m
mesg(caseno 7 : lx = 10.0 * ly
mesg(caseno 9 : not  bending but tension
mesg(caseno 11 : point downward load at right-hand end
mesg(caseno 13 : point load but no rotation at RH end
mesg(caseno 15 : uniform downward load; RH end fixed
mesg(caseno 17 : no load; x-direction temperature gradient
mesg(caseno 19 : no load; y-direction temperature gradient
caseno=1
mesgm(caseno = :caseno: Enter another if not OK
readvdu(caseno, int, 1)
caseno 
libref=caseno     ! useful because this appears in graphical monitor
 TEXT(2D bent beam in plane strain; s211
libref=211 
TITLE
INTEGER(NXBODY,NYBODY)    ! grid-fineness menu
 NXBODY = 20
 NYBODY = 20
REAL(POWERX,POWERY) 
mesgm(in beam, nx and ny are :nxbody: and :nybody: OK? Y/n
readvdu(ans,char,y)
if(:ans:.eq.n) then
mesga(Note that your choices will  not be reflected in the
mesg(PHOTON USE file)
mesgm(insert other nxbody
readvdu(nxbody,int,:nxbody:)
nxbody
mesgm(insert other nybody, not less than 3
readvdu(nybody,int,:nybody:)
nybody
endif
mesga(grid intervals are uniform except for thin cells
mesgb(near beam surface. OK? (Y/n)
readvdu(ans,char,y)
if(:ans:.eq.n) then
mesga(To use power-law x-grid insert powerx
mesg(powerx > 1 gives close spacing at low x
mesg(powerx < 1 gives close spacing at high x
mesg(insert powerx
readvdu(powerx,real,1.0)
powerx

mesga(To use symmetric power-law y-grid insert powery
mesg(powery > 1 gives close spacing at low and high y
mesg(powery < 1 gives close spacing at intermediate y
mesg(insert powery
readvdu(powery,real,1.0)
powery
endif


  Group 3. X-Direction Grid Spacing
  
 if(caseno.eq.7) then
 lx0=ly0*10.0
 endif
 lx0
 CARTES  =    T
 NREGX=3                             ! 3 regions
 IREGX=1;GRDPWR(X,1,0.01*LX0,1.0)     ! single inner fluid cell
 IREGX=2;GRDPWR(X,NXBODY,LX0,POWERX)  
 IREGX=3;GRDPWR(X,1,0.01*LX0,1.0)     ! single outer fluid cell 

 ************************************************************
  Group 4. Y-Direction Grid Spacing
 NREGY=3                            ! 3 regions
 IREGY=1;GRDPWR(Y,1,0.01*LY0,1.0)     ! single inner fluid cell
 IREGY=2;GRDPWR(Y,-NYBODY,LY0,POWERY)  
 IREGY=3;GRDPWR(Y,1,0.01*LY0,1.0)     ! single outer fluid cell 
 ************************************************************
  Group 5. Z-Direction Grid Spacing
 NZ=1
 ZWLAST  = 0.001
 if(caseno.eq.5) then
 ZWLAST  = 1.0
 endif
 zwlast
 ************************************************************
  Group 7. Variables: STOREd,SOLVEd,NAMEd
 ONEPHS  =    T
 SOLVE(P1,V1,U1)
 STORE(PRPS)
 STORE(STRX,STRY,STRZ,STXY)
 STORE(EPSY,EPSX,EPSZ)
 STORE(U1T,V1T,DILT)
 IF(CASENO.EQ.17.OR.CASENO.EQ.19) THEN
 SOLVE(TEM1)
 STORE(EPST,DVO1)
 ENDIF
 
 STORE(VISL,DVO1,DRH1)


  inform7begin
  
REAL(FX,LX,LY,POISSON,YOUNG) 
FX=40.0e6  ! H/m^2 = 40 N/mm^2    
           ! STRX(Right) = FX/(LY/2) * (Y - LY/2) = Sigma0 * (Y - LY/2)
LX=120.e-3
LY=40.e-3

FX=FX/(LY/2)          ! force gradient  
YOUNG   = 1/0.5E-11   ! Young's modulus
YOUNG
POISSON=0.3           ! Poisson's ratio
  POISSON=0.0           ! Poisson's ratio
  
 if(caseno.eq.7) then
 lx=ly*10.0
 endif
  
real(RA,RB,RC,RD)
RA = FX/YOUNG*(1-POISSON**2)     ! reference is required for the origin
RB = -FX/YOUNG/2*(1-POISSON**2)  ! of these formulae
RC = POISSON/(1-POISSON)
RD = FX/YOUNG*(1+POISSON)*(1-2*POISSON)
RA
RB
RC
RD

   **** CALCULATE analytical solution ***
char(FR1,FR2,FR3,FR4)
   FR4=(YV-0.51*:LY:)  ! YV was incorrect
FR4=(Yg-0.51*:LY:)
FR1=:RB:*((YV-0.51*:LY:)^2*:RC:
FR2= +(XG-0.01*:LX:)^2)
FR3= :RA:*(XU-0.01*:LX:)*(YG-0.51*:LY:)  

if(caseno.eq.9) then
 fr3=:RA:*(XU-0.01*:LX:)*:LY:
 fr1=(:rb:*:ly:*(YV-0.01*:ly:)*:rc:  
 fr2=*2.0)
 fr4=:ly:
endif
  Exact solutions not yet provided for cases 11 and above
(STORED VAR V1T IS :FR1::FR2:)
(STORED VAR DILT IS :RD:*:FR4: with imat>100)
(STORED VAR U1T IS :FR3: with imat>100)
(STORED VAR u1/t IS u1/u1t with imat>100)
(STORED VAR v1/t IS v1/v1t)
(STORED VAR p1/t IS p1/dilt with imat>100)
  inform7end




 ************************************************************
  GROUP 8. ITERATION NUMBERS ETC
RESFAC=1.e-8 
RESREF(V1)=0.0 
RESREF(U1)=0.0
RESREF(P1)=0.0 
endit(v1)=0.0
liter(v1)= 100
LITER(P1) = 4 
 ************************************************************
  GROUP 9. PROPERTIES
SETPRPS(1, 0,Q1)
  
 CSG10='Q1'                  ! materials with differing POISSON ratios
  MATFLG=T;NMAT=2         
  160    7800.0    0.3       473.0   43.0      1.0e-5   0.5e-11 
  161    7800.0    0.0       473.0   43.0      1.0e-5   0.5e-11 
INTEGER(IPRPS)
IPRPS=160         ! used in initial-condition setting below
IF(CASENO.EQ.3.OR.POISSON.EQ.0.0) THEN
 IPRPS=161
 POISSON=0.0      ! for use in exact-solution calculation
ENDIF
 ************************************************************
  GROUP 11. INITIAL VALUES
fiinit(p1)=0.0
fiinit(u1)=0.0
fiinit(v1)=0.0
 
fiinit(U1T) = 0.0
fiinit(V1T)=  0.0
fiinit(DILT)=  0.0

 FIINIT(PRPS)=0
 PATCH(BODY,INIVAL,2,NX-1,2,NY-1,1,1,1,1)
 INIT(BODY,PRPS,FIXVAL,IPRPS)

 ************************************************************
  GROUP 13. BOUNDARY & SPECIAL SOURCES
 
IF(CASENO.EQ.17) THEN  ! x-wise temperature gradient
PATCH(LEFTT,CELL,2,2,2,NY-1,1,1,1,1)
COVAL(LEFTT,TEM1,FIXVAL,-10)
PATCH(RIGHTT,CELL,NX-1,NX-1,2,NY-1,1,1,1,1)
COVAL(RIGHTT,TEM1,FIXVAL,10)
PATCH(RIGHT,NORTH,NX-1,NX-1,1,1,1,1,1,1)
COVAL(RIGHT,V1,FIXVAL,0.0)
TEXT(x-wise temperature gradient; s111
fx=0
ENDIF

IF(CASENO.EQ.19) THEN  ! y-wise temperature gradient
PATCH(bottomT,CELL,2,nx-1,2,2,1,1,1,1)
COVAL(bottomT,TEM1,FIXVAL,-10)
PATCH(topT,CELL,2,NX-1,NY-1,ny-1,1,1,1,1)
COVAL(topT,TEM1,FIXVAL,10)
  PATCH(RIGHT,NORTH,NX-1,NX-1,1,1,1,1,1,1)
  COVAL(RIGHT,V1,FIXVAL,0.0)
TEXT(y-wise temperature gradient; s111
fx=0
ENDIF


PATCH(AXESUU,WEST,1,1,2,NY-1,1,1,1,1)            ! LEFT : x-fixed
COVAL(AXESUU,U1,1.e6,0.0)

  inform13begin
char(form)
form=:FX:*(YG-0.51*:LY:)  ! 0.51 because thickness is LY+2*0.01*LY 

PATCH(FORC01,EAST,NX-1,NX-1,2,NY-1,1,1,1,1)      ! Right hand end

if(caseno.eq.9) then
form=:fx:*:LY:                        ! simple tension
 TEXT(beam in simple tension; s220
endif

if(caseno.gt.10) then
(SOURCE of V1 at FORC01 is COVAL(FIXFLU,-0.1*:fx:)) ! downward force
form=0.0 
 TEXT(cantilevered beam with point load; s111
ENDIF

if(caseno.eq.15) then
patch(topforce,north,2,nx-1,ny-1,ny-1,1,1,1,1)
coval(topforce,v1,fixflu,-0.1*:fx:/:ly:)  ! uniform downward force
patch(right,north,nx-1,nx-1,1,1,1,1,1,1)
coval(right,v1,fixval,0.0)
endif

if(caseno.gt.12) then
(SOURCE of U1 at FORC01 is COVAL(FIXVAL,0.0))    ! U1 fixed to zero
else
  (SOURCE of U1 at FORC01 is COVAL(FIXFLU,form))   ! bending or tension
endif
patch(rightend,cell,nx-1,nx-1,2,ny-1,1,1,1,1)  ! better end condition for u
if(nybody.gt.1) then
 (source of u1 at rightend is coval(1.e10,u1t))
endif
  patch(rightenv,ewall,nx-1,nx-1,1,ny-1,1,1,1,1)  ! better end condition for v
  (source of v1 at rightenv is coval(1.0,v1t))

patch(rightenv,cell,nx-1,nx-1,1,ny-1,1,1,1,1)  ! better end condition for v
(source of v1 at rightenv is coval(1.0e10,v1t))

  PATCH(AXESVV,cell,2,2,NY/2+1,NY/2+1,1,1,1,1)     ! V1 - fixed
PATCH(AXESVV,cell,2,2,1,1,1,1,1,1)     ! V1 fixed at bottom on left
PATCH(AXESVV,cell,2,2,1,ny-1,1,1,1,1)     ! V1 fixed at bottom on left
  PATCH(AXESVV,cell,2,3,1,ny-1,1,1,1,1)     ! V1 fixed at bottom on left
  COVAL(AXESVV,V1,FIXVAL,0.0)
(SOURCE OF V1 AT AXESVV IS COVAL(1.E10,V1T))
     ! PLANE-STRAIN,   EPSZ = 0  

SPEDAT(BOUNDARY,ZCONST,R,1.0e20)

if(caseno.eq.13) then
TEXT(cantilever; no end rotation: s12
endif
if(caseno.eq.15) then
TEXT(uniform load; no end rotation: s111
endif 
title

  if(caseno.ne.7) then
       spedat(rlxfac,rlxu1d,r,0.5)  
       spedat(rlxfac,rlxv1d,r,0.1)  
  else
       spedat(rlxfac,rlxu1d,r,1.5)  
       spedat(rlxfac,rlxv1d,r,1.5)  
  endif


  inform13end


 ************************************************************
  GROUP 15. TERMINATE SWEEPS
 LSWEEP  =   1600
lsweep
ISG21=LSWEEP
 ************************************************************
  GROUP 17. RELAXATION
     #CONPROM
     
 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
output(drh1,n,n,n,n,n,n)
output(visl,n,n,n,n,n,n)
TSTSWP = - 1   ! graphic-mode
  NXPRIN=1;  ! IXPRL=NX; IXPRF=NX-5
  NYPRIN=1;  ! IYPRL=NY; IYPRF=NY-5
IXMON = NX-2
IYMON = ny/2
IZMON = 1
#maxmin
#endpause
#$s003

STOP