C******************************************************************* C** C** v e m e 0 0 e x m 0 9 C** C** solution of a 2-dimensional, linear structural analysis C** problem with element depending material coefficients. The mesh C** is generated by this program. C** C** by L. Grosz Karlsruhe, June 1995 C** C******************************************************************* C** C** Here we give an example for the use of veme00 to solve C** a 2D linear structural analysis problems for a curved beam C** tracted by a surface load (plane stress, NK=DIM=2). The C** components of the solution vector u=(u1,u2) are the C** displacements of the body in the directions x1, x2. C** The vector of distortions in the global Cartesian C** coordinates is defined by C** C** eps1 = u1x1 C** eps2 = u2x2 C** eps12 = (u2x1+u1x2)/2 C** C** where the notations of equation are used. Corresponding C** the stress vector has the form (s1,s2,s12). By Hook's law C** the stress and distortion vector corresponds for isotropic C** materials by C** C** | s1 | | C11 C12 0 | | eps1 | C** | s2 | | C12 C11 0 | | eps2 | C** | s12 | | 0 0 2*C44 | | eps12 | 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 a suface load p=(p1,p2) 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)=area{ v1x1 * s1 + v1x2 * s12 + C** v2x1 * s12 + v2x2 * s2 } C** [1] C** C** F(v)=line{ v1 * p1 + v2 * p2 }. C** C** Additionally u has to fulfill the prescribed Dirichlet C** conditions, which gives the restraint conditions for the C** displacement of the body, see equation, veme00. C** C** The domain is a curved 2D beam of length L and height H. The C** boundaries along the beam are parabolic curve: C** C** C** p=(p1,p2)=(0,-1) surface load C** ||||||||||||||| C** C** / -- \ C** / \ C** / --- \ + C** + | / \ | | H0 C** H | | / \ | | C** + / \ + C** C** +-------------+ C** L C** C** The example program generates a subdivision of the beam into C** quadrilateral element and the upper edge into line elements, C** where the surface load p=(p1,p2) attacks. The modulus of C** elasticity and the Poisson's number depend on the element. C** They are the average values on this element. The values C** will be set by a random sample, which is used in the C** simulation of wooden beams. The generated mesh and the C** computed displacements and the stresses are written to the C** I-DEAS universal files. C** PROGRAM VEMEXM C** C**----------------------------------------------------------------- C** IMPLICIT NONE include 'bytes.h' C** C**----------------------------------------------------------------- C** C** Set the dimensions of the beam: C** DOUBLE PRECISION L,H,H0 PARAMETER (L=10,H=1,H0=.6) C** C**----------------------------------------------------------------- C** C** some parameters which may be chanced: C** C** ELEM1 = number of elements along the beam. C** ELEM2 = number of elements along the beam height. C** STORE = total storage of process in Mbytes. C** INTEGER ELEM1,ELEM2,STORE PARAMETER (ELEM2=4, & ELEM1=ELEM2*(L/H), & STORE=5) C** C** ELEM1 is selected, so that the elements are nearly C** square. C** C**----------------------------------------------------------------- C** C** N1 = number of nodes along the beam C** N2 = number of nodes in direction beam height C** C** other parameters are explained in mesh. C** INTEGER N1,N2,NK,DIM,NN,LOUT PARAMETER (N1=2*ELEM1+1, & N2=2*ELEM2+1, & NK=2, & DIM=2, & NN=N1*N2, & LOUT=6) C** C**----------------------------------------------------------------- C** C** the length of the array for the mesh are set: C** INTEGER LU,LDISP,LSTRES,LNODN,LNOD,LNOPRM,LNEK,LRPARM, & LIPARM,LDNOD,LIDPRM,LRDPRM,LIVEM,LRVEM,LLVEM,LBIG PARAMETER (LU =NN*NK, & LDISP =NN*3, & LSTRES=NN*6, & LNODN =NN, & LNOD =NN*DIM, & LNOPRM=1, & LNEK=(8*(ELEM1*ELEM2+1)+3*(ELEM1+1))*2, & LIPARM=ELEM1*ELEM2+ELEM1+2, & LRPARM=2*(ELEM1*ELEM2+1), & LDNOD =2*NK*2, & LIDPRM=1, & LRDPRM=1, & LIVEM =700+LU+LDNOD/2, & LLVEM =500, & LRVEM =60) C** C**----------------------------------------------------------------- C** C** RBIG should be as large as possible: the available C** storage STORE is reduced by all allocated array. C** the remaining storage is reserved for RBIG. C** PARAMETER ( LBIG=(STORE * 1 000 000)/IREAL & - (LU+LDISP+LSTRES+LNOD+LNOPRM+LRPARM+LRDPRM) & - (LIVEM+LNODN+LNEK+LIPARM+LDNOD+LIDPRM)/RPI ) C** C**----------------------------------------------------------------- C** C** variables and arrays : C** -------------------- C** DOUBLE PRECISION T,NOD(LNOD),NOPARM(LNOPRM),RPARM(LRPARM), & RDPARM(LRDPRM),RBIG(LBIG),U(LU),RVEM(LRVEM), & GDISP(LDISP),ESTRES(LSTRES) INTEGER IVEM(LIVEM),NODNUM(LNODN),NEK(LNEK), & IPARM(LIPARM),DNOD(LDNOD),IDPARM(LIDPRM), & IBIG(RPI*LBIG) LOGICAL MASKL(NK,NK,2),MASKF(NK,2),LVEM(LLVEM) C*** INTEGER INFO,S,NE1,ADNEK1,NE2,ADNEK2,ADIVP2,ADIVP1,NDC1, & NDC2,ADDCG1,ADDCG2,MYPROC,Z1,Z2,NGROUP,MESH, & GINFO,GINFO1,DINFO,DINFO1,ELNUM,ADRVP1 DOUBLE PRECISION E0,NY0,X1,X2 CHARACTER*80 TEXT1,TEXT2,NAME C*** EXTERNAL VEM630,VEM500 EXTERNAL DUMMY,USERB,USERL,USERF,USERC,DISP,STRESS C** C**----------------------------------------------------------------- C** C** The equivalence of RBIG and IBIG is very important : C** EQUIVALENCE (RBIG,IBIG) C** C**----------------------------------------------------------------- C** C** set the average value of E-modulus and the Poisson's number: C** E0=1.93D4 NY0=.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) IF (IVEM(200).NE.1) THEN PRINT*,'This program does not run on a parallel computer !' GOTO 9999 ENDIF C** C**----------------------------------------------------------------- C** C**** set some VECFEM parameters for mesh management: C** ---------------------------------------------- C** MESH =210 NGROUP=2 GINFO =30 GINFO1=23+2*NK DINFO =GINFO+GINFO1*NGROUP DINFO1=17 C** C**----------------------------------------------------------------- C** C**** the parameters are copied into IVEM : C** ----------------------------------- C** IVEM(1)=MESH IVEM(MESH+ 1)=N1*N2 IVEM(MESH+ 2)=NK IVEM(MESH+ 3)=DIM IVEM(MESH+ 4)=NGROUP IVEM(MESH+ 5)=NN IVEM(MESH+13)=NN IVEM(MESH+14)=0 IVEM(MESH+15)=0 IVEM(MESH+18)=0 IVEM(MESH+21)=GINFO IVEM(MESH+22)=GINFO1 IVEM(MESH+23)=DINFO IVEM(MESH+24)=DINFO1 C** C**----------------------------------------------------------------- C** C**** the generation of the geometrical nodes : C** --------------------------------------- C** DO 10 Z2=1,N2 DO 10 Z1=1,N1 X1=DBLE(Z1-1)/DBLE(N1-1) X2=DBLE(Z2-1)/DBLE(N2-1) NOD(Z1+N1*(Z2-1) )= L * (X1-0.5) NOD(Z1+N1*(Z2-1)+NN)= 4. * (1.-X1)*X1 * H0 + H*X2 NODNUM(Z1+N1*(Z2-1))=Z1+N1*(Z2-1) 10 CONTINUE C** C**----------------------------------------------------------------- C** C**** the generation of the elements : C** ------------------------------- C** C** The domain is covered by quadrilateral elements of order 2 C** and the upper boundary, where the surface load attacks, C** is described by line elments of order 2. The following C** picture illustrates the construction of the quadrilateral C** element with lower node S: C** C** S+2*N1 S+2*N1+1 S+2*N1+2 C** 4-----7-----3 C** | | C** S+N1 8 6 S+N1+2 C** | | C** 1-----5-----4 C** S S+1 S+2 C** C** ADNEK1 defines the start address of the quadrilateral C** elements in NEK and ADIVP1 defines the start address of C** the element id number assigned to the elements. The element C** id number is unique. Since the modulus of elasticity and the C** Poisson's number varry with the elements, they are assigned C** to the elements as real vector parameter (NRVP=2). ADRVP1 C** sets the address of the first parameter. NE1 is the number C** of quadrilateral elements. ELNUM counts the elements in C** group 1: C** ADNEK1=1 ADIVP1=1 ADRVP1=1 NE1=ELEM1*ELEM2 DO 20 Z2=1,ELEM2 DO 20 Z1=1,ELEM1 ELNUM=Z1+ELEM1*(Z2-1) S =2*(Z1-1) + 2*(Z2-1) * N1 + 1 IPARM(ADIVP1-1+ELNUM)=ELNUM NEK(ADNEK1-1+ELNUM )=S NEK(ADNEK1-1+ELNUM+ NE1)=S + 2 NEK(ADNEK1-1+ELNUM+2*NE1)=S + 2 *N1 + 2 NEK(ADNEK1-1+ELNUM+3*NE1)=S + 2 *N1 NEK(ADNEK1-1+ELNUM+4*NE1)=S + 1 NEK(ADNEK1-1+ELNUM+5*NE1)=S + N1 + 2 NEK(ADNEK1-1+ELNUM+6*NE1)=S + 2 *N1 + 1 NEK(ADNEK1-1+ELNUM+7*NE1)=S + N1 C** C** The modulus of elasticty and the Poisson's number on C** the element is set by a random sample, so that the C** average value will be E0 and NY0 with a dispersion of C** of 20%. Actually we have only realized an idea of C** of this model. C** RPARM(ADRVP1-1+ELNUM )= E0 * (1.+.2*SIN(FLOAT(Z1**2+Z2))) RPARM(ADRVP1-1+ELNUM+NE1)= NY0 * (1.+.2*SIN(FLOAT(Z1+Z2**2))) 20 CONTINUE C** C** ADNEK2 defines the start address of the line elements C** in NEK and ADIVP2 defines the start address of the C** element id number assigned to the elements. The entries 1 to C** 8*NE1 in NEK and 1 to NE1 in IPARM are already used by C** the elements in group 1. NE2 is the number of line elements C** on the upper boundary. C** C** S S+1 S+2 C** 2-----3-----1 C** ADNEK2=ADNEK1+8*NE1 ADIVP2=ADIVP1+NE1 NE2=ELEM1 DO 34 Z1=1,ELEM1 S =2*(Z1-1) + (N2-1) * N1 + 1 IPARM(ADIVP2-1+Z1)=Z1+ELEM1*ELEM2 NEK(ADNEK2-1+Z1 )=S + 2 NEK(ADNEK2-1+Z1+ NE2)=S NEK(ADNEK2-1+Z1+2*NE2)=S + 1 34 CONTINUE C** C**----------------------------------------------------------------- C** C** the start addresses, etc are written to IVEM: C** -------------------------------------------- C** C** group 1: quadrilateral elements C** IVEM(MESH+GINFO ) = NE1 IVEM(MESH+GINFO+ 1) = 8 IVEM(MESH+GINFO+ 2) = 4 IVEM(MESH+GINFO+ 3) = 2 IVEM(MESH+GINFO+ 4) = ADNEK1 IVEM(MESH+GINFO+ 5) = NE1 IVEM(MESH+GINFO+ 8) = 0 IVEM(MESH+GINFO+ 9) = ADRVP1 IVEM(MESH+GINFO+10) = NE1 IVEM(MESH+GINFO+11) = 2 IVEM(MESH+GINFO+13) = 0 IVEM(MESH+GINFO+14) = ADIVP1 IVEM(MESH+GINFO+15) = NE1 IVEM(MESH+GINFO+16) = 1 IVEM(MESH+GINFO+20) = ADNEK1 IVEM(MESH+GINFO+21) = NE1 IVEM(MESH+GINFO+23) = 8 C** C** group 2: line elements C** IVEM(MESH+GINFO+GINFO1 ) = NE2 IVEM(MESH+GINFO+GINFO1+ 1) = 3 IVEM(MESH+GINFO+GINFO1+ 2) = 2 IVEM(MESH+GINFO+GINFO1+ 3) = 1 IVEM(MESH+GINFO+GINFO1+ 4) = ADNEK2 IVEM(MESH+GINFO+GINFO1+ 5) = NE2 IVEM(MESH+GINFO+GINFO1+ 8) = 0 IVEM(MESH+GINFO+GINFO1+11) = 0 IVEM(MESH+GINFO+GINFO1+13) = 0 IVEM(MESH+GINFO+GINFO1+14) = ADIVP2 IVEM(MESH+GINFO+GINFO1+15) = NE2 IVEM(MESH+GINFO+GINFO1+16) = 1 IVEM(MESH+GINFO+GINFO1+20) = ADNEK2 IVEM(MESH+GINFO+GINFO1+21) = NE2 IVEM(MESH+GINFO+GINFO1+23) = 3 C** C**----------------------------------------------------------------- C** C** generation of the nodes with Dirichlet conditions : C** ------------------------------------------------- C** C** NDC1/NDC2 counts the Dirichlet condition generated for C** component 1/2. ADDCG1/ADDCG2 gives the start address of C** node ids in DNOD. C** C** for component 1 the left and the right contact point of the C** get Dirichlet conditions: C** NDC1=1 ADDCG1=1 DNOD(ADDCG1 )=1 C** C** for component 2 the left contact point of the beam get a C** Dirichlet condition: C** NDC2=2 ADDCG2=NDC1+ADDCG1 DNOD(ADDCG2 )=1 DNOD(ADDCG2+1)=N1 C** C**----------------------------------------------------------------- C** C** the start addresses,etc are written to IVEM: C** C** component 1: C** IVEM(MESH+DINFO ) = NDC1 IVEM(MESH+DINFO+ 1) = ADDCG1 IVEM(MESH+DINFO+ 2) = ADDCG1 IVEM(MESH+DINFO+ 4) = 0 IVEM(MESH+DINFO+ 7) = 0 IVEM(MESH+DINFO+ 9) = 0 IVEM(MESH+DINFO+12) = 0 C** C** component 2: C** IVEM(MESH+DINFO+DINFO1 ) = NDC2 IVEM(MESH+DINFO+DINFO1+ 1) = ADDCG2 IVEM(MESH+DINFO+DINFO1+ 2) = ADDCG2 IVEM(MESH+DINFO+DINFO1+ 4) = 0 IVEM(MESH+DINFO+DINFO1+ 7) = 0 IVEM(MESH+DINFO+DINFO1+ 9) = 0 IVEM(MESH+DINFO+DINFO1+12) = 0 C** C**----------------------------------------------------------------- C** C**** print mesh : C** ---------- C** IVEM(20)=LOUT IVEM(21)=1111 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**** distribute mesh : C** ---------------- C** IVEM(80)=LOUT IVEM(81)=1 IVEM(51)=2 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**** 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-4 IVEM(3)=0 IVEM(10)=10 IVEM(11)=11 IVEM(40)=LOUT IVEM(41)=100 IVEM(45)=500 IVEM(46)=0 IVEM(52)=1 IVEM(70)=123 IVEM(71)=1 IVEM(72)=10 000 MASKL(1,1,1)=.TRUE. MASKL(1,2,1)=.TRUE. MASKL(2,1,1)=.TRUE. MASKL(2,2,1)=.TRUE. MASKL(1,1,2)=.FALSE. MASKL(1,2,2)=.FALSE. MASKL(2,1,2)=.FALSE. MASKL(2,2,2)=.FALSE. MASKF(1,1)=.FALSE. MASKF(2,1)=.FALSE. MASKF(1,2)=.FALSE. MASKF(2,2)=.TRUE. 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**** the displacements are computed at the geometrical nodes : C** ------------------------------------------------------- C** IVEM(4)=30 IVEM(30)=LOUT IVEM(31)=1 IVEM(32)=NN IVEM(33)=2 CALL VEMU05 (T,LDISP,GDISP,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 stresses are computed at the center points of elements : C** ---------------------------------------------------------- C** IVEM(30)=LOUT IVEM(31)=1 IVEM(33)=6 CALL VEMU03 (T,LSTRES,ESTRES,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 displacements are written to the I-DEAS result file : C** ------------------------------------------------------- C** IVEM(120)=LOUT IVEM(121)=1 IVEM(127)=99 IVEM(128)=NN IVEM(129)=2 IVEM(130)=1 IVEM(131)=2 IVEM(132)=8 IVEM(137)=1 IVEM(138)=-1 TEXT1='displacements' TEXT2='veme00exm09' OPEN(99,FILE='disp.unv',STATUS= 'UNKNOWN',FORM='FORMATTED') CALL VEID97 (TEXT1,TEXT2,T,LDISP,GDISP, & LIVEM,IVEM,LNEK,NEK,LRPARM,RPARM,LIPARM,IPARM, & LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM,LNODN, & NODNUM,LNOD,NOD,LNOPRM,NOPARM, & LBIG,RBIG,IBIG) CLOSE(99) IF (IVEM(2).NE.0) GOTO 9999 C** C**----------------------------------------------------------------- C** C**** the stresses are written to the I-DEAS element result file : C** ---------------------------------------------------------- C** IVEM(120)=LOUT IVEM(121)=1 IVEM(133)=99 IVEM(134)=1 IVEM(135)=4 IVEM(136)=2 IVEM(137)=2 IVEM(138)=-1 TEXT1='stresses' TEXT2='veme00exm09' OPEN(99,FILE='stress.unv',STATUS= 'UNKNOWN',FORM='FORMATTED') CALL VEID99 (TEXT1,TEXT2,T,LSTRES,ESTRES, & LIVEM,IVEM,LNEK,NEK,LRPARM,RPARM,LIPARM,IPARM, & LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM,LNODN, & NODNUM,LNOD,NOD,LNOPRM,NOPARM, & LBIG,RBIG,IBIG) CLOSE(99) IF (IVEM(2).NE.0) GOTO 9999 C** C**----------------------------------------------------------------- C** C**** the mesh arrays written to the I-DEAS universal file : C** ---------------------------------------------------- C** IVEM(120)=LOUT IVEM(121)=1 IVEM(124)=0 IVEM(125)=99 TEXT1='veme00exm09' OPEN(99,FILE='mesh.unv',STATUS= 'UNKNOWN',FORM='FORMATTED') CALL VEMIDE (TEXT1, & LIVEM,IVEM,LNEK,NEK,LRPARM,RPARM,LIPARM,IPARM, & LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM, & LNODN,NODNUM,LNOD,NOD,LNOPRM,NOPARM, & LBIG,RBIG,IBIG) CLOSE(99) 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 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) IMPLICIT NONE 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 C** C**----------------------------------------------------------------- C** C**** start of calculation : C** -------------------- C** DO 10 Z=1,NELIS CU(Z,1) = U(Z,1) CU(Z,2) = U(Z,2) 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 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) IMPLICIT NONE C** C******************************************************************* C** C** The routine STRESS evaluates the stress vector for the given C** displacement. The succession in the output vector is C** prescribed by I-DEAS, see userc, vemu03, veid99. 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,EPS12,E,NY,C11,C44,C12 C** C**----------------------------------------------------------------- C** C**** start of calculation : C** -------------------- C** DO 100 Z=1,NELIS E =RVPARM(LAST+Z,1) NY=RVPARM(LAST+Z,2) C11=E*(1-NY)/(1+NY) C44=E/2/(1+NY) C12=E*NY/(1+NY)/(1-NY) EPS1 =DUDX(Z,1,1) EPS2 =DUDX(Z,2,2) EPS12=DUDX(Z,1,2)+DUDX(Z,2,1) CU(Z,1)=C11*EPS1+C12*EPS2 CU(Z,2)=C44*EPS12 CU(Z,3)=0 CU(Z,4)=C12*EPS1+C11*EPS2 CU(Z,5)=0 CU(Z,6)=0 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 USERB(T,COMPU,RHS, & NRSDP,RSDPRM,NRVDP,RVDP1,RVDPRM, & NISDP,ISDPRM,NIVDP,IVDP1,IVDPRM, & NDC,DIM,X,NOP,NOPARM,B) IMPLICIT NONE 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, which are equal C** to zero for both directions. Therefore therfore zeros are the C** values of the output vector B for all components. 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** IF (COMPU.EQ.1) THEN DO 10 Z=1,NDC B(Z)=0 10 CONTINUE ENDIF IF (COMPU.EQ.2) THEN DO 20 Z=1,NDC B(Z)=0 20 CONTINUE ENDIF 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) IMPLICIT NONE 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)=area{ C11*v1x1*u1x1+C44*v1x2*u1x2+ C** C12*v1x1*u2x2+C44*v1x2*u2x1+ C** C12*v2x2*u1x1+C44*v2x1*u1x2+ C** C44*v2x1*u2x1+C11*v2x2*u2x2 } 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 C** C**----------------------------------------------------------------- C** C**** start of calculation : C** --------------------- C** IF (CLASS.EQ.2) THEN IF ((COMPV.EQ.1).AND.(COMPU.EQ.1)) THEN DO 10 Z=1,NELIS E =RVPARM(LAST+Z,1) NY=RVPARM(LAST+Z,2) C11=E*(1-NY)/(1+NY) C44=E/2/(1+NY) L3(Z,1,1)=C11 L3(Z,2,2)=C44 10 CONTINUE ENDIF IF ((COMPV.EQ.2).AND.(COMPU.EQ.2)) THEN DO 20 Z=1,NELIS E =RVPARM(LAST+Z,1) NY=RVPARM(LAST+Z,2) C11=E*(1-NY)/(1+NY) C44=E/2/(1+NY) L3(Z,1,1)=C44 L3(Z,2,2)=C11 20 CONTINUE ENDIF IF (((COMPV.EQ.1) .AND. (COMPU.EQ.2)) .OR. & ((COMPV.EQ.2) .AND. (COMPU.EQ.1))) THEN DO 40 Z=1,NELIS E =RVPARM(LAST+Z,1) NY=RVPARM(LAST+Z,2) C44=E/2/(1+NY) C12=E*NY/(1+NY)/(1-NY) L3(Z,COMPV,COMPU)=C12 L3(Z,COMPU,COMPV)=C44 40 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) IMPLICIT NONE 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)=line{ - v2 }. C** C** Therefore the output vector F0 for the second component of C** the test function is equal to -1. 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 C** C**----------------------------------------------------------------- C** C**** start of calculation : C** -------------------- C** IF (CLASS.EQ.1) THEN IF (COMPV.EQ.2) THEN DO 10 Z=1,NELIS F0(Z)=-1. 10 CONTINUE ENDIF ENDIF C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C** R E T U R N C**---end of USERF--------------------------------------------------