C EPSI2 SOURCE SP204843 25/07/03 21:15:05 12308 SUBROUTINE EPSI2(IPMAIL,IPMINT,MELE,IELE, & IVADEP,NBPTEL,LRE,NSTRS,LHOOK, & MFR,NDEP,IPORE,IREPS2,NBPGAU,IVAEPS,UZDPG,RYDPG,RXDPG,IIPDPG, & IDERI,ivamat,ivade2,mate,nmatt,cmate,ngra,noer,kerr) C---------------------------------------------------------------------* C C calcul des deformations C C massif, poreux, joints poreux, incompressibles C---------------------------------------------------------------------* C * C entrees : * C ________ * C * C ipmail pointeur sur un segment meleme * C ipmint pointeur sur un segment minte * C mele numero de l'element fini * C iele numero geometrique de l'element * C nbpgau nombre de point d'integration pour la rigidite * C ivadep pointeur sur le chamelem de deplacements * C nbptel nombre de points par element * C lre nombre de ddl dans la matrice de rigidite * C nstrs nombre de composante de contraintes/deformations * C pour une matrice de hooke * C lhook dimension de la matrice de hooke * C mfr numero de la formulation de l'element fini * C ndep nombre de composantes de deplacements * C ipore nombre de fonctions de forme * C iresp2 flag pour indiquer si on veut les contraintes * C de piola-kirchhoff * C uzdpg = deformation au point nsdpge support de la * C rydpf = deformation plane generalisee * C rxdpg = * C * C sorties : * C ________ * C * C ivaeps pointeur sur un segment mptval contenant les * C les melvals de deformations * C---------------------------------------------------------------------* C Pour MEMOIRE : Si MELE element incompressible alors MFR = 31 C---------------------------------------------------------------------* IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC CCREEL -INC CCHAMP -INC CCGEOME -INC SMCOORD -INC SMCHAML -INC SMCHPOI -INC SMELEME -INC SMINTE -INC TMPTVAL SEGMENT MWRK1 REAL*8 DDHOOK(NSTRS,NSTRS),XDDL(LRE),XSTRS(NSTRS) REAL*8 XE(3,NBBB),DDHOMU(NSTRS,NSTRS) REAL*8 SHPWRK(6,NBNO),BGENE(LHOOK,LRE) REAL*8 XE1(3,NBBB),XE2(3,NBBB),xstrs2(NSTRS) REAL*8 xjac(3,3),valmat(20) ENDSEGMENT SEGMENT MWRK2 REAL*8 TENS(9),tentra(9),xddls2(lre) REAL*8 BGR(NGRA,LRE),BB(2,NGRA),gradi(ngra),R(ngra),u(ngra) ENDSEGMENT SEGMENT MWRK3 REAL*8 BPSS(3,3),XEL(3,NBBB) REAL*8 XNTH(LPP,LPP),XNTB(LPP,LPP),XNTT(LPP) ENDSEGMENT SEGMENT MWRK5 REAL*8 XGENE(NSTN,LRN) ENDSEGMENT SEGMENT MTRACE REAL*8 TRACE(NBPTEL) ENDSEGMENT CHARACTER*8 CMATE DIMENSION A(4,60),BBX(3,60),UDPGE(3) DIMENSION IN(6),JN(6),ITAB(3,3),PP(4,4) real*8 valcar(12),var(3) real*8 cobma(lhook) DATA IN/1,2,3,1,1,2/ DATA JN/1,2,3,2,3,3/ DATA ITAB(1,1),ITAB(1,2),ITAB(1,3)/1,4,5/ DATA ITAB(2,1),ITAB(2,2),ITAB(2,3)/4,2,6/ DATA ITAB(3,1),ITAB(3,2),ITAB(3,3)/5,6,3/ real*8 s(2) s(1)=0.d0 s(2)=0.d0 kerr=0 MWRK1 = 0 MWRK2 = 0 MWRK3 = 0 MWRK5 = 0 MTRACE = 0 C Introduction du point autour duquel se fait le mouvement C de la section en defo plane generalisee C IIPDPG = numero du noeud/point support si defini pour le modele C NDPGE > 0 si prise en compte du point support IF (IIPDPG.GT.0) THEN IF (IFOUR.EQ.-3) THEN NDPGE=3 UDPGE(1)=UZDPG UDPGE(2)=RYDPG UDPGE(3)=RXDPG IREF=(IIPDPG-1)*(IDIM+1) XDPGE=XCOOR(IREF+1) YDPGE=XCOOR(IREF+2) ELSE IF (IFOUR.EQ.11) THEN NDPGE=2 UDPGE(1)=UZDPG UDPGE(2)=RXDPG UDPGE(3)=XZero XDPGE=XZero YDPGE=XZero ELSE IF (IFOUR.EQ. 7 .OR. IFOUR.EQ. 8 .OR. IFOUR.EQ. 9 .OR. & IFOUR.EQ.10 .OR. IFOUR.EQ.14) THEN NDPGE=1 UDPGE(1)=UZDPG UDPGE(2)=XZero UDPGE(3)=XZero XDPGE=XZero YDPGE=XZero else write(ioimp,*) 'EPSI2 : ERREUR NDPGE' call erreur(5) return ENDIF ELSE NDPGE=0 UDPGE(1)=UZDPG UDPGE(2)=XZero UDPGE(3)=XZero XDPGE=XZero YDPGE=XZero ENDIF MELEME=IPMAIL NBNN =NUM(/1) NBELEM=NUM(/2) NHRM=NIFOUR MINTE=IPMINT NBBB=NBNN C Petite verification prealable (normalement inutile) mptval = IVAEPS if (NSTRS.ne.ival(/1)) then write(ioimp,*) 'EPSI3 : incoherence NSTRS & IVAEPS' call erreur(5) return endif do icomp = 1, NSTRS melval = IVAL(ICOMP) if (melval.le.0) then write(ioimp,*) 'EPSI3 : incoherence IVAEPS ival(',icomp,')=0' call erreur(5) return endif if (NBPTEL.NE.melval.velche(/1)) then write(ioimp,*) 'EPSI3 : incoherence NSTRS & IVAEPS' call erreur(5) return endif if (NBELEM .NE. melval.velche(/2)) then write(ioimp,*) 'EPSI3 : incoherence NSTRS & IVAEPS' call erreur(5) return endif enddo C_______________________________________________________________________ C C numero des etiquettes : C etiquettes de 1 a 98 pour traitement specifique a l element C dans la zone specifique a chaque element commencant par : C 5 continue C element 5 etiquettes 1005 2005 3005 4005 ... C 44 continue C element 44 etiquettes 1044 2044 3044 4044 ... C_______________________________________________________________________ C IF(MELE.GE.1.AND.MELE.LE.100) THEN C CABL SEG2 SEG3 TRI3 TRI4 TRI6 TRI7 QUA4 QUA5 QUA8 GOTO ( 99, 99, 99, 4, 99, 4, 99, 4, 99, 4 C QUA9 RAC2 RAC3 CUB8 CU20 PRI6 PR15 LIA3 LIA4 LIA6 1 , 99, 99, 99, 4, 4, 4, 4, 99, 99, 99 C LIA8 MULT TET4 TE10 PYR5 PY13 COQ3 DKT POUT LISP 2 , 99, 99, 4, 4, 4, 4, 99, 99, 99, 99 C FAC3 FAC4 FAC6 FAC8 LTR3 LQU4 LCU8 LPR6 LTE4 LPY5 3 , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 C COQ8 TUYA TUFI COQ2 POI1 BARR RACO LSU2 COQ4 LISM 4 , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 C COF3 RES2 LSU3 LSU4 LICO COQ6 CVS2 CVS3 CVT3 CVT6 5 , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 C CVQ4 CVQ8 THP5 TH13 THP6 TH15 THC8 TH20 ICT3 ICQ4 6 , 99, 99, 99, 99, 99, 99, 99, 99, 4, 4 C ICT6 ICQ8 ICC8 ICT4 ICP6 IC20 IC10 IC15 TRIP QUAP 7 , 4, 4, 4, 4, 4, 4, 4, 4, 79, 79 C CUBP TETP PRIP TIMO JOI2 JOI3 JOT3 JOI4 JOI6 JOI8 8 , 79, 79, 79, 99, 99, 99, 99, 99, 99, 99 C LISC TRIH DST LIC4 CERC TUYO LSE2 LITU HYT3 HYQ4 9 , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99) c cccccc . ,MELE ELSEIF(MELE.GE.101.AND.MELE.LE.200) THEN C HYT4 HYP6 HYC8 TRIS QUAS POIS FOR3 JOP3 JOP6 JOP8 GOTO ( 99, 99, 99, 99, 99, 99, 99, 80, 80, 80 C POL3 POL4 POL5 POL6 POL7 POL8 POL9 PO10 PO11 PO12 1 , 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 C PO13 PO14 BAR3 BAEX LIA2 QUAH CUBH ROT3 SEF2 TRF3 2 , 4, 4, 99, 99, 99, 99, 99, 99, 99, 99 C QUF4 CUF8 PRF6 TEF4 PYF5 MSE3 MTR6 MQU9 MC27 MP18 3 , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 C MT10 MP14 SEF3 TRF7 QUF9 CF27 PF21 TF15 PF19 SEG6 4 , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 C TR21 QU36 C216 P126 TE56 PY91 TRH6 BSE2 BTR4 BQU5 5 , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 C BCU9 BPR7 BTE5 BPY6 FRO4 SEGS POJS JCT3 JCI4 JGI2 6 , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 C JGT3 JGI4 TRIQ QUAQ CUBQ TETQ PRIQ TRIR QUAR CUBR 7 , 99, 99, 173, 173, 173, 173, 173, 173, 173, 173 C TETR PRIR Q4RI Q8RI JOQ3 JOQ6 JOQ8 JOR3 JOR6 JOR8 8 , 173, 173, 4, 4, 185, 185, 185, 185, 185, 185 C T1D2 T1D3 M1D2 M1D3 LC03 LC07 LC09 LC27 LC21 LC15 9 , 99, 99, 4, 4, 99, 99, 99, 99, 99, 99) c cccccc . ,MELE-100 ELSEIF(MELE.GE.201.AND.MELE.LE.300) THEN C LC19 LS03 LS07 LS09 LS27 LS21 LS15 LS19 BS03 BS07 GOTO ( 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 C BS09 BS27 BS21 BS15 BS19 MC03 MC07 MC09 MC27 MC21 1 , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 C MC15 MC19 M103 M107 M109 M127 M121 M115 M119 MS03 2 , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 C MS07 MS09 MS27 MS21 MS15 MS19 QC03 QC07 QC09 QC27 3 , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 C QC21 QC15 QC19 Q103 Q107 Q109 Q127 Q121 Q115 Q119 4 , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 C QS03 QS07 QS09 QS27 QS21 QS15 QS19 CIFL SURE SHB8 5 , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 C CAF2 CAF3 XQ4R XC8R JOI1 ZCO2 ZCO3 ZCO4 TUY2 TUY3 6 , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 C COS2 COA2 ICY5 IC13 CU27 PR21 TE15 PY19 C20R P15R 7 , 99, 99, 4, 4, 4, 4, 4, 4, 4, 4) c cccccc . ,MELE-200 ENDIF GOTO 99 C C_______________________________________________________________________ C C elements massifs et elements incompressibles MECANIQUE C_______________________________________________________________________ C 4 CONTINUE IF (MFR.EQ.71 .OR. MFR.EQ.73) GOTO 97173 C IDERI <= 2 pour lineaire et quadratique et = 5 pour utilisateur C =============================================================== IF ( IDERI.LE.2.OR.IDERI.EQ.5 ) THEN C Elements massifs en FORMULATION 'MECANIQUE' C ------------------------------------------- NBNO=NBNN NDDD=NDEP-NDPGE C C Donnees liees a l'element de reference C SEGINI,MWRK1 IF (Ideri.eq.2) SEGINI,MTRACE C C boucle sur les elements C DO 3004 IB=1,NBELEM C C on cherche les deplacements C MPTVAL=IVADEP IE=1 DO IGAU=1,NBNN DO ICOMP=1,NDDD MELVAL=IVAL(ICOMP) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) XDDL(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO ENDDO IF (NDPGE.GT.0) THEN DO i=1,NDPGE XDDL(IE)=UDPGE(i) IE=IE+1 ENDDO ENDIF C C on cherche les coordonnees des noeuds de l element ib C CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) C on se met a mi-pas if (ideri.eq.5) then do iyu=1,nbnn i_z = (iyu-1)*nddd do i=1,idim XE(i,iyu)= xe(i,iyu) + xddl( i + i_z )*0.5D0 enddo enddo endif C C boucle sur les points de gauss C ISDJC=0 C C Calcul des coeff de modification de b-barre (elts incompres) C= NOM : ICT3, ICQ4, ICT6, ICQ8, ICC8, ICT4, ICP6, IC20, IC10, IC15 C= MELE : 69 , 70 , 71 , 72 , 73 , 74 , 75 , 76 , 77 , 78 IF (MFR.EQ.31) THEN CALL BBCALC(MELE,NBNN,MFR,IDIM,XE, 1 NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU, 2 NSTRS,LRE,IFOUR,NHRM,A,BBX,SHPTOT,SHPWRK, 3 BGENE,XDPGE,YDPGE,PP,NOER) IF (NOER.NE.0) THEN CALL ERREUR(noer) RETURN ENDIF ENDIF DO IGAU=1,NBPTEL C CALL BMATST(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU, 1 MELE,MFR,NBNN,LRE,IFOUR,NSTRS,NHRM,1.D0,XE, 2 SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) IF (DJAC.EQ.0.D0) THEN kerr=259 if (noer.eq.0) THEN INTERR(1)=IB CALL ERREUR(259) endif GOTO 9904 ENDIF IF (DJAC.LT.0.D0) ISDJC=ISDJC+1 C En cas d'elements incompressibles : BGENE selon la methode B-BARRE IF (MFR.EQ.31) THEN CALL BBAR(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU, & MELE,NBNN,LRE,IFOUR,NSTRS,XE,DJAC,A,BBX,BGENE) ENDIF C CALL BST(BGENE,XDDL,LRE,NSTRS,XSTRS) C C calcul des eps 2 C IF (Ideri.eq.2) & CALL BST2(SHPWRK,XDDL,XE,NBNO,IFOUR,XSTRS,TRACE, & IGAU,XDPGE,YDPGE,UDPGE,NHRM) C C remplissage du segment contenant les deformations C MPTVAL=IVAEPS DO ICOMP=1,NSTRS MELVAL=IVAL(ICOMP) VELCHE(IGAU,IB)=XSTRS(ICOMP) ENDDO C ENDDO C C fin de la boucle sur les points de gauss C C** IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN IF (ISDJC.NE.0.AND.ISDJC.NE.NBPTEL) THEN kerr=195 if (noer.eq.0) then INTERR(1)=IB CALL ERREUR(195) endif GOTO 9904 ENDIF C C correction sur la partie quadratique de la deformation dans le cas C des elements incompressibles C IF (Ideri.eq.2) THEN IF (MFR.EQ.31) THEN CALL BBST2(TRACE,NBPTEL,IFOUR,MELE,POIGAU,QSIGAU, & ETAGAU,DZEGAU,SHPTOT,NBNO,SHPWRK,XE,PP) MPTVAL=IVAEPS L=2 IF (IDIM.EQ.3 .OR. IFOUR.EQ.0) L=3 DO ICOMP=1,L MELVAL=IVAL(ICOMP) DO IGAU=1,NBPTEL VELCHE(IGAU,IB)=VELCHE(IGAU,IB)+TRACE(IGAU) ENDDO ENDDO ENDIF ENDIF 3004 CONTINUE C C fin de la boucle sur les elements C 9904 CONTINUE C =============================================================== C Cas de la derivee de Truesdell C =============================================================== ELSE IF (IDERI.EQ.3) THEN NBNO=NBNN NDDD=NDEP-NDPGE SEGINI,MWRK1 C IF (IREPS2.EQ.1) SEGINI,MTRACE C on cherche le max des variations des champs pour savoir s'il faut C appeler hookis plusieurs fois mptval=IVAMAT nbgmat=0 nelmat=0 DO IM=1,NMATT MELVAL=IVAL(IM) IF (MELVAL.NE.0) THEN nelmat=Max(nelmat,VELCHE(/2)) nbgmat=Max(nbgmat,VELCHE(/1)) ENDIF VALMAT(IM) = 0.D0 ENDDO C C boucle sur les elements C DO 34 IB=1,NBELEM C C on cherche les deplacements C MPTVAL=IVADEP IE=1 DO IGAU=1,NBNN DO ICOMP=1,NDDD MELVAL=IVAL(ICOMP) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) XDDL(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO ENDDO IF (NDPGE.GT.0) THEN DO i=1,NDPGE XDDL(IE)=UDPGE(i) IE=IE+1 ENDDO ENDIF C C on cherche les coordonnees des noeuds de l element ib C CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) C on ajoute aux coordonnees la moitie du deplacement pour faire C la configuration a mi-pas do iyu=1,idim i_z = (iyu-1)*nddd do i=1,nbnn XE(i,iyu)= XE(i,iyu) + xddl(i+i_z)*0.5D0 enddo enddo C C boucle sur les points de gauss C ISDJC=0 C DO 54 IGAU=1,NBPTEL CALL BMATST(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU, 1 MELE,MFR,NBNN,LRE,IFOUR,NSTRS,NHRM,1.D0,XE, 2 SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) IF (DJAC.EQ.0.D0) THEN kerr=259 if (noer.eq.0) THEN INTERR(1)=IB CALL ERREUR(259) endif GOTO 994 ELSE IF (DJAC.LT.0.D0) THEN ISDJC=ISDJC+1 ENDIF C C on cherche les matrices de Hooke C if(nbgmat+nelmat.gt.2 . or . ib+igau.eq.2) then mptval=ivamat DO IM=1,NMATT MELVAL=IVAL(IM) IF (MELVAL.NE.0) THEN IBMN=MIN(IB ,VELCHE(/2)) IGMN=MIN(IGAU,VELCHE(/1)) VALMAT(IM)=VELCHE(IGMN,IBMN) ELSE VALMAT(IM)=0.D0 ENDIF ENDDO kcas=2 CALL HOOKIS(VALMAT,VALCAR,VAR,MFR,IB,IGAU,EXCEN,EPAIST, + INAT,MELE,NPINT,IFOUR,KCAS,NBGMAT,Nelmat, + S,SECT,LHOOK,DDHOMU,DDHOOK, + COBMA,XMOB,IRETOU) endif do i=1,nstrs do iyu=1,nstrs ddhomu(iyu,i)=ddhook(iyu,i) enddo enddo CALL DBST(BGENE,DDHomu,XDDL,LRE,NSTRS,XSTRS) C xstrs contient la contrainte on va faire pica xstrs zdep05 DO INO = 1, NBNN i_z = (ino-1)*nddd DO ID=1,IDIM XE1(ID,INO)=XE(ID,INO) XE2(ID,INO)=XE(ID,INO)-xddl( id + i_z )*0.5D0 ENDDO ENDDO DO IYU=1,3 DO i=1,3 XJAC(i,iyu)=0.D0 enddo enddo CALL HPRIME(XE1,NBNN,IDIM,SHPtot,IGAU,SHPWRK,DJAC) C C CALCUL DE LA MATRICE F C DO IF=1,IDIM DO IE=1,IDIM R1 = 0.D0 DO ID=1,NBNN R1 = R1 + SHPWRK(IF+1,ID)*XE2(IE,ID) ENDDO XJAC(IE,IF) = R1 ENDDO ENDDO IF(IDIM.EQ.2) THEN XJAC(3,3)=1.D0 IF(IFOUR.EQ.0) THEN C CC CAS AXISYMETRIQUE C R1=0.D0 R2=0.D0 DO 150 ID=1,NBNN R1=R1+SHPWRK(1,ID)*XE1(1,ID) R2=R2+SHPWRK(1,ID)*XE2(1,ID) 150 CONTINUE XJAC(3,3)=R2/(R1+1.D-20) ENDIF ENDIF CC CALCUL DE DETERMINANT DE F C IF(IDIM.EQ.3) THEN DETF=XJAC(1,1)*(XJAC(2,2)*XJAC(3,3)-XJAC(3,2)*XJAC(2,3)) DETF=DETF-XJAC(2,1)*(XJAC(1,2)*XJAC(3,3)-XJAC(3,2)*XJAC(1,3)) DETF=DETF+XJAC(3,1)*(XJAC(1,2)*XJAC(2,3)-XJAC(1,3)*XJAC(2,2)) ELSE IF(IDIM.EQ.2) THEN DETF=XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1) DETF = DETF * XJAC (3,3) ENDIF DETF=1.D0/(DETF+1.D-20) C C CALCUL DES CONTRAINTES DE CAUCHY C DO ID=1,NSTRS IND=IN(ID) JND=JN(ID) R1=0.D0 DO IE=1,IDIM DO IF=1,IDIM ICO=ITAB(IE,IF) R1 = R1 + XSTRS(ICO)*XJAC(IND,IE)*XJAC(JND,IF) ENDDO ENDDO XSTRS2(ID)= R1 * DETF ENDDO C C PEGON : ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE C IF(IDIM.EQ.2) THEN xstrs2(3)=xstrs(3)*XJAC(3,3)*XJAC(3,3)*DETF ENDIF C fin du calcul de capi dans dans xstrs2 la contrainte de kirchoff C on continu en calculant epsi sur config initiale DO INO=1,NBNN i_z = (ino-1) * nddd DO ID=1,IDIM XE(ID,INO)=XE2(ID,INO)+xddl( id + i_z )*0.5D0 ENDDO ENDDO C inversion loi de hooke CALL INVALM(DDHOMU,LHOOK,LHOOK,KERRE,0.D0) DO ID=1,LHOOK R1 = 0.D0 DO J=1,LHOOK R1 = R1 + DDHOMU(ID,J)*xstrs2(J) ENDDO xstrs(ID)= R1 ENDDO C C remplissage du segment contenant les deformations C MPTVAL=IVAEPS DO ICOMP=1,NSTRS MELVAL=IVAL(ICOMP) VELCHE(IGAU,IB)=XSTRS(ICOMP) ENDDO 54 continue IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN kerr=195 if (noer.eq.0) then INTERR(1)=IB CALL ERREUR(195) GOTO 994 endif ENDIF 34 CONTINUE 994 CONTINUE C fin du truesdell C =============================================================== C debut du jaumann C =============================================================== ELSE IF (IDERI.EQ.4) THEN NBNO=NBNN C* NDDD=NDEP C* IF (IFOUR.EQ.-3) NDDD=NDEP-3 NDDD=NDEP-NDPGE C SEGINI,MWRK1,MTRACE,MWRK2 C boucle sur les elements C DO 394 IB=1,NBELEM C C on cherche les deplacements C MPTVAL=IVADEP IE=1 DO IGAU=1,NBNN DO ICOMP=1,NDDD MELVAL=IVAL(ICOMP) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) XDDL(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO ENDDO IF (NDPGE.GT.0) THEN DO i=1,NDPGE XDDL(IE)=UDPGE(i) IE=IE+1 ENDDO ENDIF C C on cherche les coordonnees des noeuds de l element ib C CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) C C on se met sur la config a mi pas (XE) xe1 est la config initiale C do iyu=1,nbnn i_z = (iyu-1)*nddd do iou=1,idim XE1(iou,iyu)= xe(iou,iyu) XE(iou,iyu)= xe(iou,iyu) + xddl( iou+ i_z )*0.5d0 enddo enddo C C boucle sur les points de gauss C ISDJC=0 C C Calcul des coeff de modification de b-barre (elts incompres) C= NOM : ICT3, ICQ4, ICT6, ICQ8, ICC8, ICT4, ICP6, IC20, IC10, IC15 C= MELE : 69 , 70 , 71 , 72 , 73 , 74 , 75 , 76 , 77 , 78 IF (MFR.EQ.31) THEN CALL BBCALC(MELE,NBNN,MFR,IDIM,XE, 1 NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU, 2 NSTRS,LRE,IFOUR,NHRM,A,BBX,SHPTOT,SHPWRK, 3 BGENE,XDPGE,YDPGE,PP,NOER) IF (NOER.NE.0) THEN CALL ERREUR(noer) RETURN ENDIF ENDIF C DO 594 IGAU=1,NBPTEL CALL BMATST(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU, 1 MELE,MFR,NBNN,LRE,IFOUR,NSTRS,NHRM,1.D0,XE, 2 SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) IF (DJAC.EQ.0.D0) THEN kerr=259 if (noer.eq.0) THEN INTERR(1)=IB CALL ERREUR(259) endif GOTO 9964 ENDIF IF (DJAC.LT.0.D0) ISDJC=ISDJC+1 C En cas d'elements incompressibles : BGENE selon la methode B-BARRE IF (MFR.EQ.31) THEN CALL BBAR(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU, & MELE,NBNN,LRE,IFOUR,NSTRS,XE,DJAC,A,BBX,BGENE) ENDIF C CALL BST(BGENE,XDDL,LRE,NSTRS,XSTRS) C dans xstrs on a les deformations II sur config mi pas C on va calculer grad u/2 puis decomposition polaire puis rtens C on retravaille sur config initiale r_z=XZero iipdpg=0 CALL BGRMAS(iGau,NOELE,NBNO,LRE,IFOUR,NGRA,NIFOUR,XE1, . r_z,SHPTOT,SHPWRK,BB,BGR,DJAC,IIPDPG) do iou=1,lre xddls2(iou)= 0.5D0 * xddl(iou) enddo CALL BGRDEP(BGR,NGRA,XDDLs2,LRE,GRADI) C on ajoute l'identite au gradient if(idim.eQ.2) then gradi(1)=gradi(1)+1.D0 gradi(4)=gradi(4)+1.D0 ELSE IF(IDIM.EQ.3) THEN gradi(1)=gradi(1)+1.D0 gradi(5)=gradi(5)+1.D0 gradi(9)=gradi(9)+1.D0 ENDIF CALL POLA2(gradi,R,U,IDIM) C fait le rtens Rt.A.R on utilise u pour mettre Rt C et on met le tenseur dans le tableau tens C attention vu le stockage R est en fait Rt if(idim.eq.2) then U(1)=r(1) u(2)=r(3) U(3)=R(2) u(4)=R(4) tens(1)=xstrs(1) tens(2)=xstrs(4)*0.5d0 tens(3)=xstrs(4)*0.5d0 tens(4)=xstrs(2) elseif(idim.eq.3) then U(1)=r(1) u(2)=r(4) U(3)=R(7) u(4)=R(2) u(5)=r(5) u(6)=r(8) u(7)=r(3) u(8)=r(6) u(9)=r(9) tens(1)=xstrs(1) tens(2)=xstrs(4)*0.5D0 tens(3)=xstrs(5)*0.5D0 tens(4)=tens(2) tens(5)=xstrs(2) tens(6)=xstrs(6)*0.5D0 tens(7)=tens(3) tens(8)=tens(6) tens(9)=xstrs(3) else write(6,*)' idim est ni 2 ni 3 stop' stop endif CALL MULMAT(tentra,tens,U,IDIM,IDIM,IDIM) CALL MULMAT(tens,R,Tentra,IDIM,IDIM,IDIM) C tens contient le nouveau tenseur on va remplir xstrs C en 2 D epzz ne change pas if(idim.eq.2) then xstrs(1)=tens(1) xstrs(2)=tens(4) xstrs(4)=tens(2)*2.D0 else xstrs(1)=tens(1) xstrs(2)=tens(5) xstrs(3)=tens(9) xstrs(4)=tens(2)*2.D0 xstrs(5)=tens(3)*2.D0 xstrs(6)=tens(6)*2.D0 endif C C remplissage du segment contenant les deformations C MPTVAL=IVAEPS DO ICOMP=1,NSTRS MELVAL=IVAL(ICOMP) VELCHE(IGAU,IB)=XSTRS(ICOMP) ENDDO C 594 CONTINUE C C fin de la boucle sur les points de gauss C IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN INTERR(1)=IB if (noer.eq.1) then kerr=195 else CALL ERREUR(195) GOTO 9964 endif ENDIF 394 CONTINUE C C fin de la boucle sur les elements C 9964 CONTINUE endif C GOTO 510 C Elements massifs en FORMULATIONs 'ELECTROSTATIQUE' et 'DIFFUSION' C ----------------------------------------------------------------- 97173 CONTINUE NBNO = NBNN NDDD = NDEP SEGINI,MWRK1 C------------------------- C Boucle sur les elements C------------------------- DO IEL = 1, NBELEM C - Recuperation des coordonnees des noeuds de l element IEL CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE) C - Recuperation des inconnues primales aux noeuds de l element IEL MPTVAL = IVADEP IE = 1 DO IGAU = 1, NBNN DO ICOMP = 1, NDDD MELVAL = IVAL(ICOMP) IGMN = MIN(IGAU,VELCHE(/1)) IEMN = MIN(IEL ,VELCHE(/2)) XDDL(IE) = VELCHE(IGMN,IEMN) IE = IE+1 ENDDO ENDDO C-- -- -- -- -- -- -- -- -- C - Boucle sur les points de Gauss C-- -- -- -- -- -- -- -- -- ISDJC=0 DO IGAU = 1, NBPTEL C -- Calcul de la matrice B et du jacobien au point de Gauss IGAU IF (MFR.EQ.71) THEN CALL BELEC(XE,SHPTOT(1,1,IGAU),NBNN,LHOOK,-1, & SHPWRK,BGENE,DJAC) ELSE IF (MFR.EQ.73) THEN CALL BDIFF(XE,SHPTOT(1,1,IGAU),NBNN,LHOOK,-1, & SHPWRK,BGENE,DJAC) ENDIF IF (DJAC.EQ.0.D0) THEN kerr=259 if (noer.eq.0) THEN INTERR(1)=IEL CALL ERREUR(259) endif GOTO 98173 ENDIF IF (DJAC.LT.0.D0) ISDJC = ISDJC+1 CALL BST(BGENE,XDDL,LRE,NSTRS,XSTRS) C -- Remplissage du segment contenant les "deformations" = -grad(.) MPTVAL = IVAEPS DO ICOMP = 1, NSTRS MELVAL = IVAL(ICOMP) VELCHE(IGAU,IEL) = XSTRS(ICOMP) ENDDO C-- -- -- -- -- -- -- -- -- ENDDO C-- -- -- -- -- -- -- -- -- IF (ISDJC.NE.0 .AND. ISDJC.NE.NBPGAU) THEN kerr=195 if (noer.eq.0) THEN INTERR(1)=IEL CALL ERREUR(195) GOTO 98173 endif ENDIF C------------------------- ENDDO C------------------------- 98173 CONTINUE GOTO 510 C_______________________________________________________________________ C C milieux poreux C_______________________________________________________________________ C 79 CONTINUE C C pour ces elements nbbb = nombre de noeuds C nbno = nombre de fonctions de forme C NBNO=IPORE NSTN=1 LRN=NBNO-NBBB LRB=LRE-LRN C SEGINI,MWRK1,MWRK5 C Initialisation de MTRACE necessaire mais inutilise pour ces elements IF (IREPS2.EQ.1) SEGINI MTRACE C DO 3079 IB=1,NBELEM C C on cherche les coordonnees des noeuds de l element ib C CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) C C on cherche les deplacements C MPTVAL=IVADEP IE=1 DO IGAU=1,NBNN DO ICOMP=1,NDEP-1 MELVAL=IVAL(ICOMP) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) XDDL(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO ENDDO C C puis les pressions C MELVAL=IVAL(NDEP) IBMN=MIN(IB,VELCHE(/2)) DO IGAU=1,LRN IGAUSO=IBSOM(NSPOS(IELE)+IGAU-1) IGMN=MIN(IGAUSO,VELCHE(/1)) XDDL(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO C C boucle sur les points de gauss C ISDJC=0 C DO 5079 IGAU=1,NBPTEL C CALL BNPORE(IGAU,NBNO,NBBB,LRE,IFOUR,LHOOK,NSTN,NHRM, & 1.D0,XE,SHPTOT,SHPWRK,BGENE,XGENE,DJAC,1) C IF (DJAC.EQ.0.D0) THEN INTERR(1)=IB if (noer.eq.0) CALL ERREUR(259) kerr=259 GOTO 9979 ENDIF IF (DJAC.LT.0.D0) ISDJC=ISDJC+1 C CALL BST(BGENE,XDDL,LRE,LHOOK,XSTRS) C C calcul des eps 2 C IF (IREPS2.EQ.1) & CALL BST2(SHPWRK,XDDL,XE,NBNN,IFOUR,XSTRS,TRACE, & IGAU,XDPGE,YDPGE,UDPGE,NHRM) C C calcul de la pression C XP=0.D0 DO ID=1,LRN XP=XP+XGENE(1,ID)*XDDL(LRB+ID) ENDDO XSTRS(NSTRS)=XP C C remplissage du segment contenant les deformations C MPTVAL=IVAEPS DO 7079 ICOMP=1,NSTRS MELVAL=IVAL(ICOMP) VELCHE(IGAU,IB)=XSTRS(ICOMP) 7079 CONTINUE C 5079 CONTINUE C C fin de la boucle sur les points de gauss C IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN INTERR(1)=IB if (noer.eq.1) then kerr=195 else CALL ERREUR(195) GOTO 9979 endif ENDIF C 3079 CONTINUE C C fin de la boucle sur les elements C 9979 CONTINUE C GOTO 510 C_______________________________________________________________________ C C milieux poreux - SUITE C_______________________________________________________________________ C 173 CONTINUE C C pour ces elements nbbb = nombre de noeuds C nbno = nombre de fonctions de forme C IF (MELE.GE.173.AND.MELE.LE.177) THEN IDECAP = 2 ELSE IF (MELE.GE.178.AND.MELE.LE.182) THEN IDECAP = 3 ENDIF C NBNO=IPORE NSTN=IDECAP NSTB=4 IF(IFOUR.EQ.1.OR.IFOUR.EQ.2) NSTB=6 C LPP=NBNO-NBBB LRN=IDECAP*LPP LRB=LRE-LRN C SEGINI,MWRK1,MWRK5 C Initialise de MTRACE necessaire mais inutilise pour cet element IF (IREPS2.EQ.1) SEGINI MTRACE C DO 3173 IB=1,NBELEM C C on cherche les coordonnees des noeuds de l element ib C CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) C C on cherche les deplacements C MPTVAL=IVADEP IE=1 DO IGAU=1,NBNN DO ICOMP=1,NDEP-IDECAP MELVAL=IVAL(ICOMP) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) XDDL(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO ENDDO C C puis les pressions C DO IPR = 1,IDECAP MELVAL=IVAL(NDEP-IDECAP+IPR) DO IGAU=1,LPP IGAUSO=IBSOM(NSPOS(IELE)+IGAU-1) IGMN=MIN(IGAUSO,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) XDDL(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO ENDDO C C boucle sur les points de gauss C ISDJC=0 C DO 5173 IGAU=1,NBPTEL C CALL BNQORE(IGAU,NBNO,NBBB,LRE,IFOUR,NSTB,NSTN,NHRM, & 1.D0,XE,SHPTOT,SHPWRK,BGENE,XGENE,DJAC,IDECAP,LHOOK,1) C IF (DJAC.EQ.0.D0) THEN INTERR(1)=IB if (noer.eq.0) CALL ERREUR(259) kerr=259 GOTO 9173 ENDIF IF (DJAC.LT.0.D0) ISDJC=ISDJC+1 C CALL BST(BGENE,XDDL,LRE,LHOOK,XSTRS) C C calcul des eps 2 C IF (IREPS2.EQ.1) & CALL BST2(SHPWRK,XDDL,XE,NBNN,IFOUR,XSTRS,TRACE, & IGAU,XDPGE,YDPGE,UDPGE,NHRM) C C calcul des pressions C IE=LRB DO IPR=1,IDECAP XP=0.D0 IPR1=(IPR-1)*LPP DO ID=1,LPP IE=IE+1 XP=XP+XGENE(IPR,ID+IPR1)*XDDL(IE) ENDDO XSTRS(NSTRS-IDECAP+IPR)=XP ENDDO C C remplissage du segment contenant les deformations C MPTVAL=IVAEPS DO 7173 ICOMP=1,NSTRS MELVAL=IVAL(ICOMP) VELCHE(IGAU,IB)=XSTRS(ICOMP) 7173 CONTINUE C 5173 CONTINUE C C fin de la boucle sur les points de gauss C IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN INTERR(1)=IB if (noer.eq.1) then kerr=195 else CALL ERREUR(195) GOTO 9173 endif ENDIF C 3173 CONTINUE C C fin de la boucle sur les elements C 9173 CONTINUE C GOTO 510 C_______________________________________________________________________ C C joints poreux C_______________________________________________________________________ C 80 CONTINUE C C pour ces elements nbbb = nombre de noeuds C nbno = nombre de fonctions de forme C NBNO=IPORE NSTN=1 LRN=(NBNO-NBBB)*3/2 LPP = LRN LRB=LRE-LRN NFAC=(3*NBBB-NBNO)/2 C SEGINI,MWRK1,MWRK3,MWRK5 C DO 3080 IB=1,NBELEM C C on cherche d'abord les deplacements C MPTVAL=IVADEP IE=1 DO 4180 IGAU=1,NFAC DO 4280 ICOMP=1,NDEP-1 MELVAL=IVAL(ICOMP) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) XDDL(IE)=VELCHE(IGMN,IBMN) IE=IE+1 4280 CONTINUE 4180 CONTINUE C C puis les pressions C MELVAL=IVAL(NDEP) DO 4080 IGAU=1,NBNN DO 4190 INSOM=1,NBSOM(IELE) IF (IGAU.EQ.IBSOM(NSPOS(IELE)+INSOM-1)) GO TO 4191 4190 CONTINUE IF (IGAU.GT.NFAC) GO TO 4191 GO TO 4080 4191 CONTINUE IBMN=MIN(IB ,VELCHE(/2)) IGMN=MIN(IGAU,VELCHE(/1)) XDDL(IE)=VELCHE(IGMN,IBMN) IE=IE+1 4080 CONTINUE C C on cherche les coordonnees des noeuds de l element ib C CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) C C calcul des exes locaux et des coordonnees locales C CALL JOPLOC(XE,SHPTOT,NBBB,NBNO,IFOUR,XEL,BPSS) C CALL INTDEL(XNTH,XNTB,XNTT,LRN,MELE) C C boucle sur les points de gauss C ISDJC=0 C DO 5080 IGAU=1,NBPTEL C CALL BNPORJ(IGAU,NBNO,NBBB,LRE,IFOUR,LHOOK,NSTN,XE,XEL, & SHPTOT,SHPWRK,BPSS,BGENE,XGENE,DJAC,1) C IF (DJAC.EQ.0.D0) THEN INTERR(1)=IB if (noer.eq.0) CALL ERREUR(259) kerr=259 GOTO 9980 ENDIF IF (DJAC.LT.0.D0) ISDJC=ISDJC+1 C CALL BST(BGENE,XDDL,LRB,LHOOK,XSTRS) C C calcul de la pression C XP=0.D0 DO 4480 ID=1,LRN XP=XP+XNTT(ID)*XGENE(1,ID)*XDDL(LRB+ID) 4480 CONTINUE XSTRS(NSTRS)=XP C C remplissage du segment contenant les deformations C MPTVAL=IVAEPS DO 7080 ICOMP=1,NSTRS MELVAL=IVAL(ICOMP) VELCHE(IGAU,IB)=XSTRS(ICOMP) 7080 CONTINUE C 5080 CONTINUE C IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN INTERR(1)=IB if (noer.eq.1) then kerr=195 else CALL ERREUR(195) GOTO 9980 endif ENDIF C 3080 CONTINUE C 9980 CONTINUE GOTO 510 C_______________________________________________________________________ C C joints poreux - SUITE C_______________________________________________________________________ C 185 CONTINUE C C pour ces elements nbbb = nombre de noeuds C nbno = nombre de fonctions de forme C IF (MELE.GE.185.AND.MELE.LE.187) THEN IDECAP = 2 ELSE IF (MELE.GE.188.AND.MELE.LE.190) THEN IDECAP = 3 ENDIF NBNO=IPORE NSTN=IDECAP NSTB=2 IF(IFOUR.EQ.1.OR.IFOUR.EQ.2) NSTB=3 LPP=(NBNO-NBBB)*3/2 LRN=IDECAP*LPP LRB=LRE-LRN NFAC=(3*NBBB-NBNO)/2 SEGINI,MWRK1,MWRK3,MWRK5 DO 3185 IB=1,NBELEM C C on cherche d'abord les deplacements C MPTVAL=IVADEP IE=1 DO 4185 IGAU=1,NFAC DO 4285 ICOMP=1,NDEP-IDECAP MELVAL=IVAL(ICOMP) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) XDDL(IE)=VELCHE(IGMN,IBMN) IE=IE+1 4285 CONTINUE 4185 CONTINUE C C puis les pressions C DO 4785 IPR=1,IDECAP MELVAL=IVAL(NDEP-IDECAP+IPR) DO 4085 IGAU=1,NBNN DO 4195 INSOM=1,NBSOM(IELE) IF (IGAU.EQ.IBSOM(NSPOS(IELE)+INSOM-1)) GO TO 4891 4195 CONTINUE IF (IGAU.GT.NFAC) GO TO 4891 GO TO 4085 4891 CONTINUE IBMN=MIN(IB ,VELCHE(/2)) IGMN=MIN(IGAU,VELCHE(/1)) XDDL(IE)=VELCHE(IGMN,IBMN) IE=IE+1 4085 CONTINUE 4785 CONTINUE C C on cherche les coordonnees des noeuds de l element ib C CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) C C calcul des exes locaux et des coordonnees locales C CALL JOPLOC(XE,SHPTOT,NBBB,NBNO,IFOUR,XEL,BPSS) C CALL INTDEL(XNTH,XNTB,XNTT,LPP,MELE) C C boucle sur les points de gauss C ISDJC=0 C DO 5185 IGAU=1,NBPTEL C CALL BNPQRJ(IGAU,NBNO,NBBB,LRE,IFOUR,LHOOK,NSTN,XE,XEL, & SHPTOT,SHPWRK,BPSS,BGENE,XGENE,DJAC,IDECAP,NSTB,1) C IF (DJAC.EQ.0.D0) THEN INTERR(1)=IB if (noer.eq.0) CALL ERREUR(259) kerr=259 GOTO 9985 ENDIF IF (DJAC.LT.0.D0) ISDJC=ISDJC+1 C CALL BST(BGENE,XDDL,LRB,LHOOK,XSTRS) C C calcul de la pression C IE=LRB DO 4985 IPR=1,IDECAP XP=0.D0 IPR1=(IPR-1)*LPP DO 4485 ID=1,LPP IE=IE+1 XP=XP+XNTT(ID)*XGENE(IPR,ID+IPR1)*XDDL(IE) 4485 CONTINUE XSTRS(NSTRS-IDECAP+IPR)=XP 4985 CONTINUE C C remplissage du segment contenant les deformations C MPTVAL=IVAEPS DO 7185 ICOMP=1,NSTRS MELVAL=IVAL(ICOMP) VELCHE(IGAU,IB)=XSTRS(ICOMP) 7185 CONTINUE C 5185 CONTINUE C IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN kerr=195 INTERR(1)=IB if (noer.eq.0) then CALL ERREUR(195) GOTO 9985 endif ENDIF C 3185 CONTINUE C 9985 CONTINUE C GOTO 510 C____________________________________________________________________ 99 CONTINUE MOTERR(1:4)=NOMTP(MELE) MOTERR(9:12)='EPSI' CALL ERREUR(86) 510 CONTINUE SEGSUP,MWRK1 IF (MWRK2.NE.0) SEGSUP,MWRK2 IF (MWRK3.NE.0) SEGSUP,MWRK3 IF (MWRK5.NE.0) SEGSUP,MWRK5 IF (MTRACE.NE.0) SEGSUP MTRACE c RETURN END