C EPSI2 SOURCE PV 13/11/08 21:15:24 7855 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 calcul des deformations c c massif, poreux, joints poreux, incompressibles *---------------------------------------------------------------------* * * * entrees : * * ________ * * * * ipmail pointeur sur un segment meleme * * ipmint pointeur sur un segment minte * * mele numero de l'element fini * * iele numero geometrique de l'element * nbpgau nombre de point d'integration pour la rigidite * * ivadep pointeur sur le chamelem de deplacements * * nbptel nombre de points par element * * lre nombre de ddl dans la matrice de rigidite * * nstrs nombre de composante de contraintes/deformations * * pour une matrice de hooke * * lhook dimension de la matrice de hooke * * mfr numero de la formulation de l'element fini * * ndep nombre de composantes de deplacements * * ipore nombre de fonctions de forme * * iresp2 flag pour indiquer si on veut les contraintes * * de piola-kirchhoff * * uzdpg = deformation au point nsdpge support de la * * rydpf = deformation plane generalisee * * rxdpg = * * * * sorties : * * ________ * * * * ivaeps pointeur sur un segment mptval contenant les * * les melvals de deformations * *---------------------------------------------------------------------* * Pour MEMOIRE : Si MELE element incompressible alors MFR = 31 *---------------------------------------------------------------------* IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC CCOPTIO -INC CCREEL -INC CCHAMP -INC CCGEOME -INC SMCHAML -INC SMCHPOI -INC SMELEME -INC SMCOORD -INC SMMODEL -INC SMINTE -INC SMLREEL SEGMENT WRK1 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 WRK2 REAL*8 BGR(NGRA,LRE),BB(2,NGRA),gradi(ngra),R(ngra),u(ngra) REAL*8 TENS(9),tentra(9),xddls2(lre) ENDSEGMENT c SEGMENT WRK3 REAL*8 BPSS(3,3),XEL(3,NBBB) REAL*8 XNTH(LRN,LRN),XNTB(LRN,LRN),XNTT(LRN) ENDSEGMENT c SEGMENT WRK5 REAL*8 XGENE(NSTN,LRN) ENDSEGMENT c SEGMENT NOTYPE CHARACTER*16 TYPE(NBTYPE) ENDSEGMENT c SEGMENT MPTVAL INTEGER IPOS(NS),NSOF(NS) INTEGER IVAL(NCOSOU) CHARACTER*16 TYVAL(NCOSOU) ENDSEGMENT c SEGMENT MTRACE REAL*8 TRACE(NBPTEL) ENDSEGMENT CHARACTER*(NCONCH) CONM PARAMETER (NINF=3) DIMENSION A(4,60),BBX(3,60),UDPGE(3) INTEGER INFOS(NINF) LOGICAL BDPGE CHARACTER*8 CMATE DIMENSION IN(6),JN(6),ITAB(3,3) C DATA IN/1,2,3,1,1,2/ DATA JN/1,2,3,2,3,3/ C 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/ c c Introduction du point autour duquel se fait le mouvement c de la section en defo plane generalisee C En 1D : pas de rotation kerr=0 BDPGE=.FALSE. NDPGE=0 XDPGE=XZero YDPGE=XZero IF (IFOUR.EQ.-3) THEN BDPGE=.TRUE. NDPGE=3 UDPGE(1)=UZDPG UDPGE(2)=RYDPG UDPGE(3)=RXDPG SEGACT,MCOORD IREF=(IIPDPG-1)*(IDIM+1) XDPGE=XCOOR(IREF+1) YDPGE=XCOOR(IREF+2) ELSE IF (IDIM.EQ.1) THEN IF (IFOUR.EQ. 7 .OR. IFOUR.EQ. 8 .OR. IFOUR.EQ. 9 .OR. . IFOUR.EQ.10 .OR. IFOUR.EQ.14) THEN BDPGE=.TRUE. NDPGE=1 UDPGE(1)=UZDPG ELSE IF (IFOUR.EQ.11) THEN BDPGE=.TRUE. NDPGE=2 UDPGE(1)=UZDPG UDPGE(2)=RXDPG ENDIF ENDIF c MELEME=IPMAIL NBNN=NUM(/1) NBELEM=NUM(/2) c NHRM=NIFOUR MINTE=IPMINT NBBB=NBNN c 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 GOTO (99,99,99, 4,99, 4,99, 4,99, 4,99,99,99, 4, 4, 4, 4,99,99,99, 1 99,99, 4, 4, 4, 4,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 2 99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 3 99,99,99,99,99,99,99,99, 4, 4, 4, 4, 4,99,99, 4,99,99,79,79, 4 79,79,79,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 5 99,99,99,99,99,99,99,80,80,80, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 6 4, 4),MELE c IF (MELE.EQ.183.OR.MELE.EQ.184.OR. . MELE.EQ.193.OR.MELE.EQ.194) GOTO 4 IF (MELE.GE.173.AND.MELE.LE.182) GOTO 173 c GOTO 99 c_______________________________________________________________________ c c elements massifs et elements incompressibles MECANIQUE c_______________________________________________________________________ c 4 CONTINUE IF (MFR.EQ.71 .OR. MFR.EQ.73) GOTO 97173 IF( ideri.le.2.or.ideri.eq.5 ) then * ideri le 2 est pour lineaire et quadratique et =5 pour utlisateur C Elements massifs en FORMULATION 'MECANIQUE' C ------------------------------------------- NBNO=NBNN c* NDDD=NDEP c* IF (IFOUR.EQ.-3) NDDD=NDEP-3 NDDD=NDEP-NDPGE c c Donnees liees a l'element de reference c SEGACT MINTE c SEGINI,WRK1 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 4004 IGAU=1,NBNN DO 4004 ICOMP=1,NDDD MELVAL=IVAL(ICOMP) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) XDDL(IE)=VELCHE(IGMN,IBMN) IE=IE+1 4004 CONTINUE IF (BDPGE) 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) if( ideri.eq.5) then * on se emet à mi-pas do iou=1,idim do iyu=1,nbnn XE(iou,iyu)= xe(iou,iyu) + xddl( iou+ (iyu-1)*nddd)/2.D0 enddo enddo endif c c Calcul des coeff de modification de b-barre (elts incompres) c IF (MFR.EQ.31) THEN CALL BBCALC(XE,MELE,NBNN,IDIM,NBPGAU, 1 POIGAU,QSIGAU,ETAGAU,DZEGAU,NSTRS, 2 LRE,IFOUR,A,BBx,NHRM,SHPTOT,SHPWRK,XDPGE,YDPGE) ENDIF c c boucle sur les points de gauss c ISDJC=0 C DO 5004 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 INTERR(1)=IB CALL ERREUR(259) 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 7004 ICOMP=1,NSTRS MELVAL=IVAL(ICOMP) IBMN=MIN(IB,VELCHE(/2)) VELCHE(IGAU,IBMN)=XSTRS(ICOMP) 7004 CONTINUE c 5004 CONTINUE c c fin de la boucle sur les points de gauss c IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN if (noer.eq.1) then kerr=195 else 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) L=2 IF (IDIM.EQ.3 .OR. IFOUR.EQ.0) L=3 DO 5005 ICOMP=1,L MELVAL=IVAL(ICOMP) IBMN=MIN(IB ,VELCHE(/2)) DO 5006 IGAU=1,NBPTEL VELCHE(IGAU,IBMN)=VELCHE(IGAU,IBMN)+TRACE(IGAU) 5006 CONTINUE 5005 CONTINUE ENDIF ENDIF 3004 CONTINUE c c fin de la boucle sur les elements c 9904 CONTINUE SEGDES MINTE SEGSUP WRK1 IF (IREPS2.EQ.1) SEGSUP MTRACE * cas de la dérivée de Truesdell elseif( ideri.eq.3) then NBNO=NBNN NDDD=NDEP-NDPGE SEGACT MINTE SEGINI,WRK1 * IF (IREPS2.EQ.1) SEGINI,MTRACE c c boucle sur les elements c DO 34 IB=1,NBELEM c c on cherche les deplacements c MPTVAL=IVADEP IE=1 DO 44 IGAU=1,NBNN DO 44 ICOMP=1,NDDD MELVAL=IVAL(ICOMP) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) XDDL(IE)=VELCHE(IGMN,IBMN) IE=IE+1 44 CONTINUE IF (BDPGE) THEN DO i=1,NDPGE XDDL(IE)=UDPGE(i) IE=IE+1 ENDDO ENDIF c on cherche le max des variations des champs pour savoir s'il faut c appeler hookis plusieurs fois c nbgmat=0 nelmat=0 mptval=ivamat segact mptval DO IM=1,NMATT IF (IVAL(IM).NE.0) THEN MELVAL=IVAL(IM) nelmat=Max(nelmat ,VELCHE(/2)) nbgmat=Max(nbgmat,VELCHE(/1)) ENDIF ENDDO segdes mptval 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 coordonnées la moitié du deplacements pour faire C la configuration à mi-pas do iou=1,idim do iyu=1,nbnn XE(iou,iyu)= xe(iou,iyu) + xddl( iou+ (iyu-1)*nddd)/2.D0 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) c IF (DJAC.EQ.0.D0) THEN INTERR(1)=IB CALL ERREUR(259) GOTO 994 ELSE IF (DJAC.LT.0.D0) THEN ISDJC=ISDJC+1 ENDIF c c on cherche les matrices de Hooke c mptval=ivamat segact mptval DO 94 IM=1,NMATT IF (IVAL(IM).NE.0) THEN MELVAL=IVAL(IM) IBMN=MIN(IB ,VELCHE(/2)) IGMN=MIN(IGAU,VELCHE(/1)) VALMAT(IM)=VELCHE(IGMN,IBMN) ELSE VALMAT(IM)=0.D0 ENDIF 94 CONTINUE segdes mptval kcas=2 if(nbgmat+nelmat.gt.2 . or . ib+igau.eq.2) then 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 iou=1,nstrs do iyu=1,nstrs ddhomu(iyu,iou)=ddhook(iyu,iou) enddo enddo CALL DBST(BGENE,DDHomu,XDDL,LRE,NSTRS,XSTRS) c xstrs contient la contrainte on va faire pica xstrs zdep05 DO 220 INO=1,NBNN DO 220 ID=1,IDIM XE1(ID,INO)=XE(ID,INO) XE2(ID,INO)=XE(ID,INO)-xddl( id+ (ino-1)*nddd)/2.D0 220 CONTINUE DO iou=1,3 DO IYU=1,3 XJAC(iou,iyu)=0.D0 enddo enddo CALL ZERO(XJAC,3,3) CALL HPRIME(XE1,NBNN,IDIM,SHPtot,IGAU,SHpwrk,DJAC) C C CALCUL DE LA MATRICE F C DO 140 ID=1,NBNN DO 140 IE=1,IDIM DO 140 IF=1,IDIM XJAC(IE,IF)=SHpwrk(IF+1,ID)*XE2(IE,ID)+XJAC(IE,IF) 140 CONTINUE 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.2) THEN DETF=XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1) DETF = DETF * XJAC (3,3) ENDIF 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)) ENDIF DETF=1.D0/(DETF+1.D-20) C C CALCUL DES CONTRAINTES DE CAUCHY C DO 160 ID=1,NSTRS IND=IN(ID) JND=JN(ID) xstrs2(ID)=0.D0 DO 170 IE=1,IDIM DO 170 IF=1,IDIM ICO=ITAB(IE,IF) XsTRS2(ID)=XsTRS(ICO)*XJAC(IND,IE)*XJAC(JND,IF)*DETF 1 +xstrs2(ID) 170 CONTINUE 160 CONTINUE 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 * fin du caldul de capi dans dans xstrs2 la contrainte de kirchoff * on continu en calculant epsi sur config initiale DO 221 INO=1,NBNN DO 221 ID=1,IDIM XE(ID,INO)=XE2(ID,INO)+xddl( id+ (ino-1)*nddd)/2.D0 221 CONTINUE * inversion loi de hook CALL INVALM(DDHOMU,LHOOK,LHOOK,KERRE,0.D0) DO 6000 I=1,LHOOK xstrs(I)=0.D0 DO 60001 J=1,LHOOK xstrs(I)=xstrs(I)+DDHOMU(I,J)*xstrs2(J) 60001 CONTINUE 6000 CONTINUE c c remplissage du segment contenant les deformations c MPTVAL=IVAEPS DO 74 ICOMP=1,NSTRS MELVAL=IVAL(ICOMP) IBMN=MIN(IB,VELCHE(/2)) VELCHE(IGAU,IBMN)=XSTRS(ICOMP) 74 CONTINUE 54 continue IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN if (noer.eq.1) then kerr=195 else INTERR(1)=IB CALL ERREUR(195) GOTO 994 endif ENDIF 34 CONTINUE 994 CONTINUE SEGSUP WRK1 segdes minte * fin du truesdell * debut du jaumann elseif(ideri.eq.4) then C========================== NBNO=NBNN c* NDDD=NDEP c* IF (IFOUR.EQ.-3) NDDD=NDEP-3 NDDD=NDEP-NDPGE c c Donnees liees a l'element de reference c SEGACT MINTE c SEGINI,WRK1 SEGINI,MTRACE segini wrk2 c boucle sur les elements c DO 394 IB=1,NBELEM c c on cherche les deplacements c MPTVAL=IVADEP IE=1 DO 494 IGAU=1,NBNN DO 494 ICOMP=1,NDDD MELVAL=IVAL(ICOMP) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) XDDL(IE)=VELCHE(IGMN,IBMN) IE=IE+1 494 CONTINUE IF (BDPGE) 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) * * on se met sur la config à mi pas (XE) xe1 est la config initiale * do iou=1,idim do iyu=1,nbnn XE1(iou,iyu)= xe(iou,iyu) XE(iou,iyu)= xe(iou,iyu) + xddl( iou+ (iyu-1)*nddd)/2.D0 enddo enddo c c Calcul des coeff de modification de b-barre (elts incompres) c IF (MFR.EQ.31) THEN CALL BBCALC(XE,MELE,NBNN,IDIM,NBPGAU, 1 POIGAU,QSIGAU,ETAGAU,DZEGAU,NSTRS, 2 LRE,IFOUR,A,BBX,NHRM,SHPTOT,SHPWRK,XDPGE,YDPGE) ENDIF c c boucle sur les points de gauss c ISDJC=0 C DO 594 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 INTERR(1)=IB CALL ERREUR(259) 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) * dans xstrs on a les deformations II sur config mi pas * on va calculer grad u/2 puis decomposition polaire puis rtens * on retravaille sur config initiale xxzero=0.d0 iipdpg=0 CALL BGRMAS(iGau,NOELE,NBNO,LRE,IFOUR,NGRA,NIFOUR,XE1, . xXZero,SHPTOT,SHPWRK,BB,BGR,DJAC,IIPDPG) do iou=1,lre xddls2(iou)= xddl(iou) / 2.D0 enddo CALL BGRDEP(BGR,NGRA,XDDLs2,LRE,GRADI) * on ajoute l'identité 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) * fait le rtens Rt.A.R on utilise u pour mettre Rt * et on met le tenseur dans le tableau tens * 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)/2.d0 tens(3)=xstrs(4)/2.D0 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)/2.D0 tens(3)=xstrs(5)/2.D0 tens(4)=tens(2) tens(5)=xstrs(2) tens(6)=xstrs(6)/2.D0 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) * tens contient le nouveau tenseur on va remplir xstrs * 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 794 ICOMP=1,NSTRS MELVAL=IVAL(ICOMP) IBMN=MIN(IB,VELCHE(/2)) VELCHE(IGAU,IBMN)=XSTRS(ICOMP) 794 CONTINUE 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 SEGDES MINTE SEGSUP WRK1,wrk2 SEGSUP MTRACE endif c GOTO 510 C Elements massifs en FORMULATIONs 'ELECTROSTATIQUE' et 'DIFFUSION' C ----------------------------------------------------------------- 97173 CONTINUE SEGACT,MINTE NBNO = NBNN NDDD = NDEP SEGINI,WRK1 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 INTERR(1) = IEL CALL ERREUR(259) 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) IEMN = MIN(IEL,VELCHE(/2)) VELCHE(IGAU,IEMN) = XSTRS(ICOMP) ENDDO C-- -- -- -- -- -- -- -- -- ENDDO C-- -- -- -- -- -- -- -- -- IF (ISDJC.NE.0 .AND. ISDJC.NE.NBPGAU) THEN INTERR(1) = IEL if (noer.eq.1) then kerr=195 else CALL ERREUR(195) GOTO 98173 endif ENDIF C------------------------- ENDDO C------------------------- 98173 CONTINUE SEGDES,MINTE SEGSUP,WRK1 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 WRK1,WRK5 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 4079 IGAU=1,NBNN DO 4079 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 4079 CONTINUE c c puis les pressions c MELVAL=IVAL(NDEP) DO 4179 IGAU=1,LRN IGAUSO=IBSOM(NSPOS(IELE)+IGAU-1) IGMN=MIN(IGAUSO,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) XDDL(IE)=VELCHE(IGMN,IBMN) IE=IE+1 4179 CONTINUE 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 CALL ERREUR(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 4279 ID=1,LRN XP=XP+XGENE(1,ID)*XDDL(LRB+ID) 4279 CONTINUE XSTRS(NSTRS)=XP c c remplissage du segment contenant les deformations c MPTVAL=IVAEPS DO 7079 ICOMP=1,NSTRS MELVAL=IVAL(ICOMP) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) VELCHE(IGMN,IBMN)=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 SEGSUP WRK1,WRK5 IF (IREPS2.EQ.1) SEGSUP MTRACE 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 c LPP=NBNO-NBBB LRN=IDECAP*LPP LRB=LRE-LRN c SEGINI WRK1,WRK5 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 4173 IGAU=1,NBNN DO 4173 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 4173 CONTINUE c c puis les pressions c DO 4473 IPR = 1,IDECAP MELVAL=IVAL(NDEP-IDECAP+IPR) DO 4273 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 4273 CONTINUE 4473 CONTINUE 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,LHOOK,NSTN,NHRM, & 1.D0,XE,SHPTOT,SHPWRK,BGENE,XGENE,DJAC,IDECAP,1) c IF (DJAC.EQ.0.D0) THEN INTERR(1)=IB CALL ERREUR(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 4673 IPR=1,IDECAP XP=0.D0 IPR1=(IPR-1)*LPP DO 4373 ID=1,LPP IE=IE+1 XP=XP+XGENE(IPR,ID+IPR1)*XDDL(IE) 4373 CONTINUE XSTRS(NSTRS-IDECAP+IPR)=XP 4673 CONTINUE c c remplissage du segment contenant les deformations c MPTVAL=IVAEPS DO 7173 ICOMP=1,NSTRS MELVAL=IVAL(ICOMP) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) VELCHE(IGMN,IBMN)=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 SEGSUP WRK1,WRK5 IF (IREPS2.EQ.1) SEGSUP MTRACE 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 LRB=LRE-LRN NFAC=(3*NBBB-NBNO)/2 c SEGINI WRK1,WRK3,WRK5 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 CALL ERREUR(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) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) VELCHE(IGMN,IBMN)=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 SEGSUP WRK1,WRK3,WRK5 c GOTO 510 c____________________________________________________________________ 99 CONTINUE MOTERR(1:4)=NOMTP(MELE) MOTERR(9:12)='EPSI' CALL ERREUR(86) 510 CONTINUE RETURN END