|
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-------------------------------------------------- |