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 percentage errors on the x- and y-direction displacements to be
computed and printed.

A menu is provided which enables:
* the standard and new (simultaneous) solvers to be compared
* 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.

ENDDIS

PHOTON USE
p;;;;

set prop off
cl
msg x-direction-displacement field
gr ou z 1
cont U1 z 1 x 2 22 y 1 22 fil;.0001
vec z 1 x 2 22 y 2 22 col 0
pause

cl
msg y-direction-displacement field
gr ou z 1
cont V1 z 1 x 2 22 y 1 22 fil;.0001
vec z 1 x 2 22 y 2 22 col 0

ENDUSE

prt0begin
print0 1.0e-20
prt0end

************************************************************
Group 1. Run Title and Number
************************************************************

Declarations and settings
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
YOUNG   = 1/0.5E-11   ! Young's modulus
YOUNG
POISSON=0.3           ! Poisson's ratio
mesg(caseno 1 or 2 : default settings
mesg(caseno 3 or 4 : Poisson's ratio = 0.0
mesg(caseno 5 or 6 : zwlast = 1.0 m
mesg(caseno 7 or 8 : lx = 10.0 * ly
mesg(caseno 9 or 10 : not  bending but tension
mesg(caseno 11 or 12 : point downward load at right-hand end
mesg(caseno 13 or 14 : point load but no rotation at RH end
mesg(caseno 15 or 16 : uniform downward load; RH end fixed
mesg(caseno odd = StressModel, even = ModModel361
caseno=1
mesgm(caseno = :caseno: Enter another if not OK
caseno
libref=caseno     ! useful because this appears in graphical monitor
Integer(SolvMod)
! SolvMod =  1 for Stone Solver        (caseno = even)
!            2 for simultaneous Solver (caseno = odd )
!
solvmod=1 + (caseno+1)/2-caseno/2
solvmod
TEXT(bent beam in plane strain; s220
************************************************************
Group 2. Time dependence
************************************************************
NXBODY = 21
NYBODY = 21
REAL(POWERX,POWERY)
mesgm(in beam, nx and ny are :nxbody: and :nybody: OK? Y/n
if(:ans:.eq.n) then
mesgm(insert other nxbody
nxbody
mesgm(insert other nybody
nybody
endif
mesgm(grid intervals are uniform. OK? Y/n
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
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
powerx
endif

Group 3. X-Direction Grid Spacing

if(caseno.eq.7) then
lx=ly*10.0
endif
if(caseno.eq.8) then
lx=ly*10.0
endif
lx
CARTES  =    T
NREGX=3                             ! 3 regions
IREGX=1;GRDPWR(X,1,0.01*LX,1.0)     ! single inner fluid cell
IREGX=2;GRDPWR(X,NXBODY,LX,POWERX)
IREGX=3;GRDPWR(X,1,0.01*LX,1.0)     ! single outer fluid cell

************************************************************
Group 4. Y-Direction Grid Spacing
NREGY=3                            ! 3 regions
IREGY=1;GRDPWR(Y,1,0.01*LY,1.0)     ! single inner fluid cell
IREGY=2;GRDPWR(Y,-NYBODY,LY,POWERY)
IREGY=3;GRDPWR(Y,1,0.01*LY,1.0)     ! single outer fluid cell
************************************************************
Group 5. Z-Direction Grid Spacing
NZ=1
ZWLAST  = 0.001
if(caseno.eq.5.or.caseno.eq.6) 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,UERR,VERR,DILT)
IF(CASENO.GT.16.AND.CASENO.LT.21) THEN
SOLVE(TEM1)
STORE(EPST,DVO1)
ENDIF

************************************************************
GROUP 8. ITERATION NUMBERS ETC
RESFAC=1.e-8
RESREF(V1)=-1
RESREF(U1)=-1
RESREF(P1)=-1
IF(SOLVMOD.EQ.1) THEN   ! combined with the LSWEEP settings below
LITER(V1) = 50         ! these set up a 'fair' contest between the
LITER(U1) = 50         ! two solvers
ELSE
LITER(V1) = 50
LITER(U1) = 50
ENDIF

LITER(P1) = 4

************************************************************
GROUP 9. PROPERTIES

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.CASENO.EQ.4) THEN
IPRPS=161
POISSON=0.3      ! 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(UERR)=  0.0
fiinit(VERR)=  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.OR.CASENO.EQ.18) 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)
fx=0
ENDIF

IF(CASENO.EQ.19.OR.CASENO.EQ.20) 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)
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)
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.or.caseno.eq.10) then
form=:fx:*:LY:                        ! simple tension
TEXT(beam in simple tension; s220
endif

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

if(caseno.eq.15.or.caseno.eq.16) 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.and.caseno.lt.17) 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(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
COVAL(AXESVV,V1,FIXVAL,0.0)

! PLANE-STRAIN,   EPSZ = 0

SPEDAT(BOUNDARY,ZCONST,R,1.0e20)

if(caseno.eq.13.or.caseno.eq.14) then
TEXT(cantilever; no end rotation: s220
endif
if(caseno.eq.15.or.caseno.eq.16) then
TEXT(uniform load; no end rotation: s220
endif
title
************************************************************
GROUP 15. TERMINATE SWEEPS
IF(SOLVMOD.EQ.1) THEN
LG(40)  = F
spedat(rlxfac,rlxu1d,r,0.5)
spedat(rlxfac,rlxv1d,r,0.5)
LSWEEP  =   1000
IF(CASENO.EQ.8) THEN
LSWEEP=2000
ENDIF
ELSE
LG(40)  = T
LSWEEP  =   400
IF(CASENO.EQ.7) THEN
LSWEEP= 600    ! 6
ENDIF
ENDIF
lsweep
ISG21=LSWEEP
************************************************************
GROUP 17. RELAXATION
#CONPROM

RELAX(P1  ,LINRLX, 1.000000E+00)

spedat(rlxfac,rlxu1d,r,1)
spedat(rlxfac,rlxv1d,r,1)

************************************************************
GROUP 19. DATA TRANSMITTED TO GROUND
STRA    =    T
PARSOL  =    F
ISG52   =    3   !   probe & res

************************************************************
GROUP 23.FIELD PRINT-OUT & PLOT CONTROL
TSTSWP = - 1   ! graphic-mode
NXPRIN=1; IXPRL=NX-1; IXPRF=NX-5
NYPRIN=1; IYPRL=NY-1; IYPRF=NY-5
IXMON = NX-2
IYMON = ny/2
IZMON = 1

inform7begin
real(RA,RB,RC,RD)
RA = FX/YOUNG*(1-POISSON**2)
RB = -FX/YOUNG/2*(1-POISSON**2)
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:)
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.or.caseno.eq.10) then
fr3=:RA:*(XU-0.01*:LX:)*:LY:
fr1=(:rb:*:ly:*(YV-0.01*:ly:)*:rc:  ! this appears to be incorrect ( correct VIA)
fr2=*2.0)
fr4=:ly:
endif
Exact solutions not yet provided for cases 11 and above
(STORED VAR V1T IS :FR1::FR2: with imat>100)
(STORED VAR DILT IS :RD:*:FR4: with imat>100)
(STORED VAR U1T IS :FR3: with imat>100)
(STORED VAR UERR IS ABS(U1-U1T)/(ABS(U1T)+1.e-10)*100 with imat>100)
(STORED VAR VERR IS ABS(V1-V1T)/(ABS(V1T)+1.e-10)*100 with imat>100)
inform7end
STOP