|
C*******************************************************************
C**
C** v e m e 0 0 e x m 0 7
C**
C** solution of a linear structural analysis problem mesh
C** generated by the user and ISVAS3 as postprocessor
C**
C** by L. Grosz Karlsruhe, Sept. 1994
C**
C*******************************************************************
C**
PROGRAM VEMEXM
C**
C**-----------------------------------------------------------------
C**
C** Here we give an example for the use of veme00 to solve
C** a linear structural analysis problems for an arbitrary
C** 3-dimensional body tracted by surface pressure (NK=DIM=3).
C** The components of the solution vector u=(u1,u2,u3) are the
C** displacements of the body in the directions x1, x2 and
C** x3. The vector of distortions in the global Cartesian
C** coordinates is defined by
C**
C** eps1 = u1x1
C** eps2 = u2x2
C** eps3 = u3x3
C** eps12 = (u2x1+u1x2)/2
C** eps23 = (u3x2+u2x3)/2
C** eps13 = (u3x1+u1x3)/2
C**
C** where the notations of equation are used. Corresponding
C** the stress vector has the form (s1,s2,s3,s12,s23,s13). By
C** Hook's law the stress and distortion vector corresponds
C** for isotropic materials by
C**
C** | s1 | | C11 C12 C12 0 0 0 | | eps1 |
C** | s2 | | C12 C11 C12 0 0 0 | | eps2 |
C** | s3 | = | C12 C12 C11 0 0 0 | * | eps3 |
C** | s12 | | 0 0 0 2*C44 0 0 | | eps12 |
C** | s23 | | 0 0 0 0 2*C44 0 | | eps23 |
C** | s13 | | 0 0 0 0 0 2*C44 | | eps13 |
C**
C** where the C-entries are given by the two values E and NY
C** called modulus of elasticity and Poisson's number.
C**
C** If the body is tracted by surface loads (F1,F2,F3) the
C** displacement of the body is the solution of the linear
C** functional equation L(v,u)=F(v) with
C**
C** L(v,u)=volume{ v1x1 * s1 + v1x2 * s12 + v1x3 * s13 +
C** v2x1 * s12 + v2x2 * s2 + v2x3 * s23 +
C** [1] v3x1 * s13 + v3x2 * s23 + v3x3 * s33 }
C**
C** F(v)=area{ v1 * F1 + v2 * F2 + v3 * F3 }.
C**
C** Additionally u has to fulfill the prescribed Dirichlet
C** conditions, which gives the restrain conditions for the
C** displacement of the body, see equation, veme00.
C**
C** The example program solves this equation for an arbitrary
C** body with isotropic material where the material properties
C** are independent from the locality in the body. The unit of
C** forces is kN and the unit of length is cm. The mesh is given
C** the vecfem input file format, see vemu02. The values
C** for the prescribed displacements are read from the input file
C** and stored as real vector parameter. The modulus of
C** elasticity, the Poisson's number and the values of the surface
C** loads are not read from the input file but set in the
C** program. The read mesh, the computed displacements and the
C** corresponding stresses are written to the ISVAS 3 mesh and
C** result files. cyl.vem is an example of an input file.
C**
C**-----------------------------------------------------------------
C**
IMPLICIT NONE
include 'bytes.h'
C**
C**-----------------------------------------------------------------
C**
C** some parameters which may be chanced:
C**
C** FINPUT = name of the vecfem input file (input)
C** FNODE = name of the ISVAS node file (output)
C** FELEM = name of the ISVAS element file (output)
C** FDISP = name of the displacement result file (output)
C** FSTRES = name of the stress result file (output)
C** STORAGE = total storage of process in Mbytes.
C**
INTEGER STORAGE
CHARACTER*80 FINPUT,FDISP,FSTRES,FELEM,FNODE
PARAMETER (FINPUT='cyl.vem',
& FNODE='nodes.isv',
& FELEM='elements.isv',
& FDISP='displacement.isv',
& FSTRES='stress.isv',
& STORAGE=20)
C**
C**-----------------------------------------------------------------
C**
C** The length of the mesh array are set:
C** It will happen, that these lengths are to small for
C** the given mesh. Then you have to enter the correct lengths,
C** which are prescribed by VECFEM, into this declaration.
C**
INTEGER LNODN,LNEK,LRPARM,LIPARM,LDNOD
PARAMETER (LNODN =2500,
& LNEK =40000,
& LIPARM=3000,
& LRPARM=1,
& LDNOD =1500)
C**
C**-----------------------------------------------------------------
C**
C** The parameters which can be estimated:
C**
INTEGER LNOD,LU,LNOPRM,LIDPRM,LRDPRM,LIVEM,LCU,
& LBIG,LLVEM,LRVEM,LOUT
PARAMETER (LNOD=LNODN*3,
& LU=LNODN*3,
& LCU=LNODN*9,
& LNOPRM=1,
& LIDPRM=LDNOD/2,
& LRDPRM=LDNOD/2,
& LIVEM=500+LU+LDNOD/2,
& LLVEM=500,
& LRVEM=60,
& LOUT=6)
C**
C**-----------------------------------------------------------------
C**
C** RBIG should be as large as possible: the available
C** storage STORAGE is reduced by all allocated array.
C** the remaining storage is reserved for RBIG.
C**
PARAMETER ( LBIG=(STORAGE * 1 000 000)/IREAL
& - (LU+LCU+LRVEM+LNOD+LNOPRM+LRPARM+LRDPRM)
& - (LIVEM+LLVEM+LNODN+LNEK+LIPARM+LDNOD+LIDPRM)/RPI)
C**
C**-----------------------------------------------------------------
C**
C** variables and arrays :
C** --------------------
C**
DOUBLE PRECISION T,NOD(LNOD),NOPARM(LNOPRM),RPARM(LRPARM),
& RVEM(LRVEM),U(LU),RDPARM(LRDPRM),RBIG(LBIG),
& CU(LCU)
INTEGER IVEM(LIVEM),NODNUM(LNODN),NEK(LNEK),
& IPARM(LIPARM),DNOD(LDNOD),IDPARM(LIDPRM),
& IBIG(RPI*LBIG)
LOGICAL MASKL(3,3,10),MASKF(3,10),LVEM(LLVEM)
C**
C**-----------------------------------------------------------------
C**
CHARACTER*80 NAME,TEXT1,TEXT2
INTEGER MYPROC,INFO,OUTFLG,MESH,NGROUP,GINFO,GINFO1,
& G,I,J,CLASS
C**
EXTERNAL VEM630,VEM500
EXTERNAL DUMMY,USERB,USERL,USERF,DISP,STRESS
C**
C**-----------------------------------------------------------------
C**
C** The equivalence of RBIG and IBIG is very important :
C**
EQUIVALENCE (RBIG,IBIG)
C**
C**-----------------------------------------------------------------
C**
C** The COMMON block transports the material properties to
C** user routines :
C**
DOUBLE PRECISION E,NY
COMMON /PROP/E,NY
E=1.93D4
NY=.3
C**
C**-----------------------------------------------------------------
C**
C** get task ids :
C**
NAME='a.out'
CALL COMBGN(IVEM(200),MYPROC,LIVEM-203,IVEM(204),NAME,INFO)
IF (INFO.NE.0) GOTO 9999
IVEM(201)=MYPROC
IVEM(202)=0
IVEM(203)=IVEM(204)
C**
C**-----------------------------------------------------------------
C**
C** a protocol is printed only on process 1 :
C**
IF (MYPROC.EQ.1) THEN
OUTFLG=1
ELSE
OUTFLG=0
ENDIF
C**
C**-----------------------------------------------------------------
C**
C**** the parameters are copied into IVEM :
C** -----------------------------------
C**
MESH=300
IVEM(1)=MESH
IVEM(MESH+ 2)=3
IVEM(MESH+ 3)=3
C**
C**-----------------------------------------------------------------
C**
C**** read a vecfem input file :
C** ------------------------
C**
IVEM(27)=LOUT
IVEM(28)=OUTFLG
IVEM(29)=99
IF (MYPROC.EQ.1) OPEN(99,FILE=FINPUT,STATUS= 'UNKNOWN',
& FORM='FORMATTED')
CALL VEMU02 (LIVEM,IVEM,LNEK,NEK,LRPARM,RPARM,LIPARM,IPARM ,
& LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM,
& LNODN,NODNUM,LNOD,NOD,LNOPRM,NOPARM,
& LBIG,RBIG,IBIG)
IF (IVEM(2).NE.0) GOTO 9999
C**
C**-----------------------------------------------------------------
C**
C**** distribute mesh :
C** ----------------
C**
IVEM(80)=LOUT
IVEM(81)=OUTFLG
IVEM(51)=5
CALL VEMDIS (LIVEM,IVEM,LNEK,NEK,LRPARM,RPARM,LIPARM,IPARM ,
& LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM,
& LNODN,NODNUM,LNOD,NOD,LNOPRM,NOPARM,
& LBIG,RBIG,IBIG)
IF (IVEM(2).NE.0) GOTO 9999
C**
C**-----------------------------------------------------------------
C**
C**** print mesh : replace '0000' by '1111' if you want to get a
C** printing of the mesh.
C** -----------
C**
IVEM(20)=LOUT
IVEM(21)=0000*OUTFLG
IVEM(22)=2
CALL VEMU01(LIVEM,IVEM,LNEK,NEK,LRPARM,RPARM,LIPARM,IPARM,
& LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM,
& LNODN,NODNUM,LNOD,NOD,LNOPRM,NOPARM,
& LBIG,RBIG,IBIG)
IF (IVEM(2).NE.0) GOTO 9999
C**
C**-----------------------------------------------------------------
C**
C**** set masks for veme00-call :
C** --------------------------
C**
NGROUP=IVEM(MESH+4)
GINFO =IVEM(MESH+21)
GINFO1=IVEM(MESH+22)
DO 100 G=1,NGROUP
CLASS=IVEM(MESH+GINFO+GINFO1*(G-1)+3)
C****
C****** right hand side :
C****
C**** only surface loads (CLASS=2) in the right hand side
C****
IF (CLASS.EQ.2) THEN
DO 10 I=1,3
MASKF(I,G)=.TRUE.
10 CONTINUE
ELSE
DO 20 I=1,3
MASKF(I,G)=.FALSE.
20 CONTINUE
ENDIF
C****
C****** bilinear form :
C****
C**** only interaction of the test function and the solution by
C**** the tensor of elasticity in the inner elements (CLASS=3)
C****
IF (CLASS.EQ.3) THEN
DO 30 I=1,3
DO 30 J=1,3
MASKL(I,J,G)=.TRUE.
30 CONTINUE
ELSE
DO 40 I=1,3
DO 40 J=1,3
MASKL(I,J,G)=.FALSE.
40 CONTINUE
ENDIF
100 CONTINUE
C**
C**-----------------------------------------------------------------
C**
C**** call of VECFEM :
C** --------------
C**
OPEN(10,FORM='UNFORMATTED',STATUS='SCRATCH')
OPEN(11,FORM='UNFORMATTED',STATUS='SCRATCH')
LVEM(1)=.TRUE.
LVEM(5)=.FALSE.
LVEM(9)=.FALSE.
RVEM(3)=1.D-8
IVEM(3)=0
IVEM(10)=10
IVEM(11)=11
IVEM(40)=LOUT
IVEM(41)=100*OUTFLG
IVEM(45)=500
IVEM(46)=0
IVEM(52)=1
IVEM(70)=123
IVEM(71)=1
IVEM(72)=10 000
CALL VEME00 (T,LU,U,LIVEM,IVEM,LLVEM,LVEM,LRVEM,RVEM,
& LNEK, NEK ,LRPARM ,RPARM ,LIPARM ,IPARM ,
& LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM,LNODN,
& NODNUM,LNOD,NOD,LNOPRM,NOPARM,LBIG,RBIG,IBIG,
& MASKL,MASKF,USERB,USERL,USERF,
& VEM500,VEM630)
IF (IVEM(2).NE.0) GOTO 9999
C**
C**-----------------------------------------------------------------
C**
C**** write the ISVAS mesh files :
C** --------------------------
C**
IVEM(120)=LOUT
IVEM(121)=OUTFLG
IVEM(125)=99
IVEM(126)=98
TEXT1='vecfem example'
IF (MYPROC.EQ.1) OPEN(99,FILE=FNODE,FORM='FORMATTED')
IF (MYPROC.EQ.1) OPEN(98,FILE=FELEM,FORM='FORMATTED')
CALL VEMISV (TEXT1,LIVEM,IVEM,LNEK,NEK,LRPARM,RPARM,LIPARM,IPARM ,
& LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM,
& LNODN,NODNUM,LNOD,NOD,LNOPRM,NOPARM,
& LBIG,RBIG,IBIG)
IF (IVEM(2).NE.0) GOTO 9999
C**
C**-----------------------------------------------------------------
C**
C**** the displacements are computed at the geometrical nodes :
C** -------------------------------------------------------
C**
IVEM(4)=30
IVEM(30)=LOUT
IVEM(31)=OUTFLG
IVEM(32)=LNODN
IVEM(33)=3
CALL VEMU05 (T,LCU,CU,LU,U,LIVEM,IVEM,
& LNEK, NEK ,LRPARM ,RPARM ,LIPARM ,IPARM ,
& LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM,LNODN,
& NODNUM,LNOD,NOD,LNOPRM,NOPARM,LBIG,RBIG,IBIG,
& DISP)
IF (IVEM(2).NE.0) GOTO 9999
C**
C**-----------------------------------------------------------------
C**
C**** the displacements are written to the ISVAS 3 result file :
C** -------------------------------------------------------
C**
IVEM(120)=LOUT
IVEM(121)=OUTFLG
IVEM(127)=99
IVEM(128)=IVEM(32)
IVEM(129)=IVEM(33)
IVEM(130)=1
TEXT1='displacements'
TEXT2='vecfem example'
IF (MYPROC.EQ.1) OPEN(99,FILE=FDISP,FORM='FORMATTED')
CALL VEIS97 (TEXT1,TEXT2,T,LCU,CU,LIVEM,IVEM,
& LNEK, NEK ,LRPARM ,RPARM ,LIPARM ,IPARM ,
& LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM,LNODN,
& NODNUM,LNOD,NOD,LNOPRM,NOPARM,LBIG,RBIG,IBIG)
IF (IVEM(2).NE.0) GOTO 9999
C**
C**-----------------------------------------------------------------
C**
C**** the stresses are computed at the geometrical nodes :
C** ---------------------------------------------------
C**
IVEM(30)=LOUT
IVEM(31)=OUTFLG
IVEM(32)=LNODN
IVEM(33)=9
CALL VEMU05 (T,LCU,CU,LU,U,LIVEM,IVEM,
& LNEK, NEK ,LRPARM ,RPARM ,LIPARM ,IPARM ,
& LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM,LNODN,
& NODNUM,LNOD,NOD,LNOPRM,NOPARM,LBIG,RBIG,IBIG,
& STRESS)
IF (IVEM(2).NE.0) GOTO 9999
C**
C**-----------------------------------------------------------------
C**
C**** the stresses are written to the ISVAS 3 result file :
C** ---------------------------------------------------
C**
IVEM(120)=LOUT
IVEM(121)=OUTFLG
IVEM(127)=99
IVEM(128)=IVEM(32)
IVEM(129)=IVEM(33)
IVEM(130)=2
TEXT1='stresses'
TEXT2='vecfem example'
IF (MYPROC.EQ.1) OPEN(99,FILE=FSTRES,FORM='FORMATTED')
CALL VEIS97 (TEXT1,TEXT2,T,LCU,CU,LIVEM,IVEM,
& LNEK, NEK ,LRPARM ,RPARM ,LIPARM ,IPARM ,
& LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM,LNODN,
& NODNUM,LNOD,NOD,LNOPRM,NOPARM,LBIG,RBIG,IBIG)
IF (IVEM(2).NE.0) GOTO 9999
C**
C**-----------------------------------------------------------------
C**
C**** end of calculation
C** ------------------
C**
9999 CALL COMEND(IVEM(200),INFO)
E N D
SUBROUTINE STRESS(T,GROUP,LAST,NELIS,
& NRSP,RSPARM,NRVP,RVP1,RVPARM,
& NISP,ISPARM,NIVP,IVP1,IVPARM,
& L,DIM,X,NK,U,DUDX,NOP,NOPARM,DNOPDX,N,CU)
C**
C*******************************************************************
C**
C** The routine STRESS evaluates the stress vector for the given
C** displacement. The full stress tensor has to be calculated in
C** spite it is symmetrical. The succession in the output vector
C** is prescribed by ISVAS 3, see userc, vemu05, veis97.
C**
C*******************************************************************
C**
INTEGER GROUP,LAST,NELIS,L,DIM,NK,N,
& NRSP,RVP1,NRVP,NISP,IVP1,NIVP,NOP
DOUBLE PRECISION T,X(L,DIM),U(L,NK),DUDX(L,NK,DIM),
& RSPARM(NRSP),RVPARM(RVP1,NRVP),
& NOPARM(L,NOP),DNOPDX(L,NOP,DIM),CU(L,N)
INTEGER ISPARM(NISP),IVPARM(IVP1,NIVP)
C**
C**-----------------------------------------------------------------
C**
INTEGER Z
DOUBLE PRECISION EPS1,EPS2,EPS3,EPS12,EPS13,EPS23,E,NY,
& C11,C44,C12
COMMON /PROP/E,NY
C**
C**-----------------------------------------------------------------
C**
C**** start of calculation :
C** --------------------
C**
C11=E*(1.-NY)/(1+NY)/(1-2*NY)
C44=E/2./(1+NY)
C12=E*NY/(1+NY)/(1-2*NY)
DO 100 Z=1,NELIS
EPS1 =DUDX(Z,1,1)
EPS2 =DUDX(Z,2,2)
EPS3 =DUDX(Z,3,3)
EPS12=DUDX(Z,1,2)+DUDX(Z,2,1)
EPS23=DUDX(Z,2,3)+DUDX(Z,3,2)
EPS13=DUDX(Z,1,3)+DUDX(Z,3,1)
CU(Z,1)=C11*EPS1+C12*EPS2+C12*EPS3
CU(Z,2)=C44*EPS12
CU(Z,3)=C44*EPS13
CU(Z,4)=CU(Z,2)
CU(Z,5)=C12*EPS1+C11*EPS2+C12*EPS3
CU(Z,6)=C44*EPS23
CU(Z,7)=CU(Z,3)
CU(Z,8)=CU(Z,6)
CU(Z,9)=C12*EPS1+C12*EPS2+C11*EPS3
100 CONTINUE
C**
C**-----------------------------------------------------------------
C**
C**** end of calculation
C** ------------------
C**
R E T U R N
C**---end of STRESS-------------------------------------------------
E N D
SUBROUTINE DISP(T,GROUP,LAST,NELIS,
& NRSP,RSPARM,NRVP,RVP1,RVPARM,
& NISP,ISPARM,NIVP,IVP1,IVPARM,
& L,DIM,X,NK,U,DUDX,NOP,NOPARM,DNOPDX,N,CU)
C**
C*******************************************************************
C**
C** The routine DISP copies the input solution to the output
C** CU. So the calling routine vemu05 interpolates the
C** displacements from the global nodes onto the geometrical
C** nodes, see userc, vemu05.
C**
C*******************************************************************
C**
INTEGER GROUP,LAST,NELIS,L,DIM,NK,N,
& NRSP,RVP1,NRVP,NISP,IVP1,NIVP,NOP
DOUBLE PRECISION T,X(L,DIM),U(L,NK),DUDX(L,NK,DIM),
& RSPARM(NRSP),RVPARM(RVP1,NRVP),
& NOPARM(L,NOP),DNOPDX(L,NOP,DIM),CU(L,N)
INTEGER ISPARM(NISP),IVPARM(IVP1,NIVP)
C**
C**-----------------------------------------------------------------
C**
INTEGER Z,I
C**
C**-----------------------------------------------------------------
C**
C**** start of calculation :
C** --------------------
C**
DO 10 I=1,MIN(N,NK)
DO 10 Z=1,NELIS
CU(Z,I) = U(Z,I)
10 CONTINUE
C**
C**-----------------------------------------------------------------
C**
C**** end of calculation
C** ------------------
C**
R E T U R N
C**---end of DISP---------------------------------------------------
E N D
SUBROUTINE USERB(T,COMPU,RHS,
& NRSDP,RSDPRM,NRVDP,RVDP1,RVDPRM,
& NISDP,ISDPRM,NIVDP,IVDP1,IVDPRM,
& NDC,DIM,X,NOP,NOPARM,B)
C**
C*******************************************************************
C**
C** The routine USERB sets the Dirichlet boundary conditions,
C** see userb. For this application these conditions gives
C** the restrain conditions for the displacements. The
C** prescribed values are defined in the input file and set to
C** the first Dirichlet vector parameter, see vemu02. Therefore
C** this parameter is copied to the output vector B.
C**
C*******************************************************************
C**
INTEGER COMPU,RHS,NRSDP,NRVDP,RVDP1,NISDP,NIVDP,IVDP1,
& NDC,DIM,NOP
DOUBLE PRECISION T,RSDPRM(NRSDP),RVDPRM(RVDP1,NRVDP),
& X(NDC,DIM),NOPARM(NDC,NOP),B(NDC)
INTEGER ISDPRM(NISDP),IVDPRM(IVDP1,NIVDP)
C**
C**-----------------------------------------------------------------
C**
INTEGER Z
C**
C**-----------------------------------------------------------------
C**
C**** start of calculation :
C** --------------------
C**
DO 10 Z=1,NDC
B(Z)=RVDPRM(Z,1)
10 CONTINUE
C**
C**-----------------------------------------------------------------
C**
C**** end of calculation
C** ------------------
C**
R E T U R N
C**---end of USERB--------------------------------------------------
E N D
SUBROUTINE USERL(T,GROUP,CLASS,COMPV,COMPU,LAST,
& NELIS,L,DIM,X,TAU,NK,U,DUDX,
& LT,UT,DUTDX,NOP,NOPARM,DNOPDX,
& NRSP,RSPARM,NRVP,RVP1,RVPARM,
& NISP,ISPARM,NIVP,IVP1,IVPARM,
& L3,L2,L1,L0)
C**
C*******************************************************************
C**
C** The routine USERL sets the coefficients of the bilinear form,
C** see userl. The values of the coefficients L3 (L2=L1=L0=0 !)
C** is read from the definition [1] of the bilinear form:
C**
C** L(v,u)=volume{ C11*v1x1*u1x1+C44*v1x2*u1x2+C44*v1x3*u1x3+
C** C12*v1x1*u2x2+C44*v1x2*u2x1+
C** C12*v1x1*u3x3+C44*v1x3*u3x1+
C** C12*v2x2*u1x1+C44*v2x1*u1x2+
C** C44*v2x1*u2x1+C11*v2x2*u2x2+C44*v2x3*u2x3+
C** C12*v2x2*u3x3+C44*v2x3*u3x2+
C** C12*v3x3*u1x1+C44*v3x1*u1x3+
C** C12*v3x3*u2x2+C44*v3x2*u2x3+
C** C44*v3x1*u3x1+C44*v3x2*u3x2+C11*v3x3*u3x3 }
C**
C*******************************************************************
C**
INTEGER GROUP,CLASS,COMPV,COMPU,LAST,NELIS,L,LT,
& DIM,NK,NOP,NRSP,RVP1,NRVP,NISP,IVP1,NIVP
DOUBLE PRECISION T,X(L,DIM),TAU(L,DIM,CLASS),U(L,NK),UT(LT,NK),
& DUDX(L,NK,CLASS),DUTDX(LT,NK,CLASS),
& NOPARM(L,NOP),DNOPDX(L,NOP,CLASS),
& RSPARM(NRSP),RVPARM(RVP1,NRVP),
& L3(L,CLASS,CLASS),L2(L,CLASS),L1(L,CLASS),L0(L)
INTEGER ISPARM(NISP),IVPARM(IVP1,NIVP)
C**
C**-----------------------------------------------------------------
C**
INTEGER Z
DOUBLE PRECISION E,NY,C11,C44,C12
COMMON /PROP/E,NY
C**
C**-----------------------------------------------------------------
C**
C**** start of calculation :
C** ---------------------
C**
IF (CLASS.EQ.3) THEN
C11=E*(1.-NY)/(1+NY)/(1-2*NY)
C44=E/2./(1+NY)
C12=E*NY/(1+NY)/(1-2*NY)
IF ((COMPV.EQ.1).AND.(COMPU.EQ.1)) THEN
DO 10 Z=1,NELIS
L3(Z,1,1)=C11
L3(Z,2,2)=C44
L3(Z,3,3)=C44
10 CONTINUE
ENDIF
IF ((COMPV.EQ.2).AND.(COMPU.EQ.2)) THEN
DO 20 Z=1,NELIS
L3(Z,1,1)=C44
L3(Z,2,2)=C11
L3(Z,3,3)=C44
20 CONTINUE
ENDIF
IF ((COMPV.EQ.3).AND.(COMPU.EQ.3)) THEN
DO 30 Z=1,NELIS
L3(Z,1,1)=C44
L3(Z,2,2)=C44
L3(Z,3,3)=C11
30 CONTINUE
ENDIF
IF (((COMPV.EQ.1) .AND. (COMPU.EQ.2)) .OR.
& ((COMPV.EQ.2) .AND. (COMPU.EQ.1))) THEN
DO 40 Z=1,NELIS
L3(Z,COMPV,COMPU)=C12
L3(Z,COMPU,COMPV)=C44
40 CONTINUE
ENDIF
IF (((COMPV.EQ.1) .AND. (COMPU.EQ.3)) .OR.
& ((COMPV.EQ.3) .AND. (COMPU.EQ.1))) THEN
DO 50 Z=1,NELIS
L3(Z,COMPV,COMPU)=C12
L3(Z,COMPU,COMPV)=C44
50 CONTINUE
ENDIF
IF (((COMPV.EQ.2) .AND. (COMPU.EQ.3)) .OR.
& ((COMPV.EQ.3) .AND. (COMPU.EQ.2))) THEN
DO 60 Z=1,NELIS
L3(Z,COMPV,COMPU)=C12
L3(Z,COMPU,COMPV)=C44
60 CONTINUE
ENDIF
ENDIF
C**
C**-----------------------------------------------------------------
C**
C**** end of calculation
C** ------------------
C**
R E T U R N
C**---end of USERL--------------------------------------------------
E N D
SUBROUTINE USERF (T,GROUP,CLASS,COMPV,RHS,LAST,
& NELIS,L,DIM,X,TAU,NK,U,DUDX,
& LT,UT,DUTDX,NOP,NOPARM,DNOPDX,
& NRSP,RSPARM,NRVP,RVP1,RVPARM,
& NISP,ISPARM,NIVP,IVP1,IVPARM,
& F1,F0)
C**
C*******************************************************************
C**
C** The routine USERF sets the coefficients of the linear form,
C** see userf. The values of the coefficients F0 (F1=0 !) is
C** read from the definition [1] of the linear form:
C**
C** F(v)=area{ v1 * F1 + v2 * F2 + v3 * F3 }.
C**
C** The values of the node loads (F1,F2,F3) shall represent a
C** unit pressure onto the surface in direction of the surface
C** normal.
C**
C*******************************************************************
C**
INTEGER GROUP,CLASS,COMPV,RHS,LAST,NELIS,L,LT,DIM,NK,NOP,
& NRSP,RVP1,NRVP,NISP,IVP1,NIVP
DOUBLE PRECISION T,X(L,DIM),TAU(L,DIM,CLASS),U(L,NK),UT(LT,NK),
& DUDX(L,NK,CLASS),DUTDX(LT,NK,CLASS),
& NOPARM(L,NOP),DNOPDX(L,NOP,CLASS),
& RSPARM(NRSP),RVPARM(RVP1,NRVP),
& F1(L,CLASS),F0(L)
INTEGER ISPARM(NISP),IVPARM(IVP1,NIVP)
C**
C**-----------------------------------------------------------------
C**
INTEGER Z
DOUBLE PRECISION N(3),NORM
C**
C**-----------------------------------------------------------------
C**
C**** start of calculation :
C** --------------------
C**
C**
IF (CLASS.EQ.2) THEN
DO 10 Z=1,NELIS
C**
C** first the surface normal is computed:
C**
N(1)= TAU(Z,2,1)*TAU(Z,3,2)-TAU(Z,2,2)*TAU(Z,3,1)
N(2)=-(TAU(Z,1,1)*TAU(Z,3,2)-TAU(Z,1,2)*TAU(Z,3,1))
N(3)= TAU(Z,1,1)*TAU(Z,2,2)-TAU(Z,2,1)*TAU(Z,1,2)
NORM=SQRT(N(1)**2+N(2)**2+N(3)**2)
C**
C** The COMPV-th component of the surface load is set :
C**
F0(Z)=-N(COMPV)/NORM
10 CONTINUE
ENDIF
C**
C**-----------------------------------------------------------------
C**
C**** end of calculation
C** ------------------
C**
R E T U R N
C**---end of USERF--------------------------------------------------
|