TALK=f;RUN(1,1)
  DISPLAY
  The bending of a beam is calculated in accordance with the
  ordinary-differential-equations method, known as the
  Engineer's theory, for three cases: 
  1. uniform bending moment, 
  2. uniform load, and 
  3. point load.
  ENDDIS
 ************************************************************
  Group 1. Run Title and Number
 ************************************************************
TEXT(bent beam; engineer's theory;s105 
TITLE
real(Force,SV0,FY,LX,LY,LZ,POISSON,YOUNG, IMomZ,Moment,LF) 
FY=40.0e8  ! H/m^2
LX=1.20
LY=5.e-3
LZ = LY
IMomZ = 1.0     

LF = 0.01

 INTEGER(CASENO)       ! start of menu
mesg(Bent beam; ENGINEER'S THEORY
mesg(caseno 1 : Moment is constant 
mesg(caseno 2 : Uniform load
mesg(caseno 3 : Point load
caseno=1
mesgm(caseno = :caseno: Enter another if not OK
readvdu(caseno, int, 1)
caseno

if(caseno.eq.1) then
 Force = FY*LX*LZ        !  [N]
 Moment = Force * LX     !  [N*m] 
 SV0 = Moment/(2*IMomZ) 
endif
if(caseno.eq.2) then
 Force=(FY*LZ)           !  [N/m]
 SV0 = Force/(4*IMomZ) 
endif
if(caseno.eq.3) then
 Force= FY*LZ*LX         !  [N]
 Moment = Force * LX     !  [N*m] 
 SV0 = Moment/(2*IMomZ) 
endif 

YOUNG = 2.0e11
POISSON=0.0
    
INTEGER(NXBODY,NYBODY)    ! grid-fineness menu
 NXBODY = 20
 NYBODY = 1
 CARTES  =    T
 NREGX=2                             ! 2 regions
 IREGX=1;GRDPWR(X,NXBODY,LX,1.0)  
 IREGX=2;GRDPWR(X,1,LF*LX,1.0)     ! single outer fluid cell 

 ************************************************************
  Group 4. Y-Direction Grid Spacing
 NREGY=3                            ! 3 regions
 IREGY=1;GRDPWR(Y,1,LF*LY,1.0)     ! single inner fluid cell
 IREGY=2;GRDPWR(Y,NYBODY,LY,1.0)  
 IREGY=3;GRDPWR(Y,1,LF*LY,1.0)     ! single outer fluid cell 
 ************************************************************
  Group 5. Z-Direction Grid Spacing
 NZ=1
 ZWLAST  = LZ
 ************************************************************
  Group 7. Variables: STOREd,SOLVEd,NAMEd
 ONEPHS  =    T
 SOLVE(P1,V1,U1)
 STORE(PRPS)
 STORE(STRY)
 STORE(EPSY,EPSX,EPSZ)
 STORE(V1Th,V1/t)
 STORE(VISL,DVO1,DRH1)

 ************************************************************
  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
  
 CSG10='Q1'                  ! materials with differing POISSON ratios
  MATFLG=T;NMAT=1         
  161    7800.0    0.0       473.0   43.0      1.0e-5   0.5e-11 
INTEGER(IPRPS)
 IPRPS=161
 ************************************************************
  GROUP 11. INITIAL VALUES
fiinit(p1)=0.0
fiinit(u1)=0.0
fiinit(v1)=0.0
 
fiinit(V1Th)=  0.0
fiinit(V1/t)=  0.0

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

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

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

PATCH(AXESVV,WWALL,1,1,1,ny-1,1,1,1,1)  ! V1 fixed on  left
coval(axesvv,v1,1,0.0)
SPEDAT(BOUNDARY,ZCONST,R,1.0e20)

patch(beamend,east,nx-1,nx-1,1,2,1,1,1,1)
patch(beam,volume,1,nx-1,1,2,1,1,1,1)   !  Volume Source as Moment
real(Tau0,CCC)

if(caseno.eq.1) then
 (source of v1 at beam is COVAL(FIXFLU,:SV0:))
 Tau0 = -SV0*LX
 (source of v1 at beamend is COVAL(FIXFLU,:Tau0:))

 char(fV,xxg)
 CCC=-SV0/YOUNG
 xxg=(XG)
 fV=:CCC:*:XXG:^2
 (stored var v1th is :fV: with swps)
 (stored var V1/t is V1/V1th)
 

endif
if(caseno.eq.2) then
 char(SV)
 SV=:SV0:*(:LX:-XG)^2
 (source of v1 at beam is COVAL(FIXFLU,:SV:))

 Tau0 = -SV0*LX**3/3
 (source of v1 at beamend is COVAL(FIXFLU,:Tau0:))

 char(fV1,fV2,xxg)
 CCC=-Force/24.0/YOUNG/ImomZ
 xxg=(XG)
 fV1=:CCC:*(:XXG:^4 - 4*:LX:*
 fV2=:XXG:^3+6*:LX:^2*:xxg:^2)
 (stored var v1th is :fV1::fV2: with swps)
 (stored var V1/t is V1/V1th)
endif
 
if(caseno.eq.3) then
 char(SV)
 SV=:SV0:*(:LX:-XG)
 (source of v1 at beam is COVAL(FIXFLU,:SV:))

 Tau0 = -SV0*LX**2/2
 (source of v1 at beamend is COVAL(FIXFLU,:Tau0:))

 char(fV1,fV2,xxg)
 CCC=-SV0/3.0/YOUNG
 xxg=(XG)
 fV1=:CCC:*(-:XXG:^3 + 3*:LX:*
 fV2=:XXG:^2)
 (stored var v1th is :fV1::fV2: with swps)
 (stored var V1/t is V1/V1th)
 
endif

    
 ************************************************************
  GROUP 15. TERMINATE SWEEPS
 LSWEEP  =   100
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
output(drh1,n,n,n,n,n,n)
output(visl,n,n,n,n,n,n)
patch(vprof,profil,1,nx,1,1,1,1,1,1)
coval(vprof,v1,0.0,0.0)
TSTSWP = - 1   ! graphic-mode

NYPRIN=1;  ! IYPRL=NY; IYPRF=NY-5
IXMON = NX-2
IYMON = ny/2
IZMON = 1
#maxmin
#endpause
#$s003

STOP