C SIGMAX SOURCE MB234859 25/09/08 21:16:11 12358 subroutine SIGMAX (MATE,IMAT,NBGMAT,NELMAT,NMATT,CMATE, & IVAMAT,IMODEL,IREPS2,IVADEP, & IVASTR,UZDPG,RYDPG,RXDPG,IIPDPG,IRETER) c C PROCEDURE UTILISEE DANS LE CAS D'ELEMENTS XFEM c POUR LE CALCUL DE la contrainte (élastique) C C C********************************************************* C PARTIE DECLARATIVE C********************************************************* IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC SMCOORD -INC SMMODEL -INC SMCHAML -INC SMINTE -INC SMELEME -INC SMLREEL C POINTEUR MCHEX1.MCHELM,MINT1.MINTE,MINT2.MINTE -INC TMPTVAL c SEGMENT WRK1 REAL*8 XE(3,NBBB) REAL*8 DDHOOK(LHOOK,LHOOK) REAL*8 XDDL(LRE),XSTRS(LHOOK) ENDSEGMENT c SEGMENT WRK2 REAL*8 SHPWRK(6,NBNO),BGENE(LHOOK,LRE) c REAL*8 LV7WRK(NBENRMA2,2,6,NBNO) REAL*8 LV7WRK(NBENRMA2,2,6,NBBB) ENDSEGMENT C SEGMENT,MVELCH REAL*8 VALMAT(NV1) ENDSEGMENT C SEGMENT MRACC INTEGER TLREEL(NBENRMA2,NBI) ENDSEGMENT C SEGMENT MTRACE REAL*8 TRACE(3,NBPTEL) ENDSEGMENT C PARAMETER (NDDLMAX=30,NBNIMAX=10) CTY PARAMETER (NDDLMAX=20,NBNIMAX=10) PARAMETER (NBENRMAX=5) DIMENSION MLRE(NBENRMAX+1) C CHARACTER*8 CMATE DIMENSION UDPGE(3) LOGICAL BDPGE C c write (*,*) '#############################' c write (*,*) '##### DEBUT DE SIGMAX #####' c write (*,*) '#############################' C C********************************************************* c Introduction du point autour duquel se fait le mouvement c de la section en defo plane generalisee C********************************************************* C En 1D : pas de rotation BDPGE=.FALSE. NDPGE=0 XDPGE=0.D0 YDPGE=0.D0 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 NV1=NMATT SEGINI,MVELCH C C C********************************************************* c C RECUPERATION + ACTIVATIONS + VALEURS UTILES c C********************************************************* C SEGACT MMODEL,IMODEL deja activé dans epsi1 c C++++ Activation au cas ou ++++++++++++++++++++++++++++++ SEGACT MCOORD C++++ Recup + Activation de la geometrie ++++++++++++++++ MELEME= IMAMOD c SEGACT MELEME deja activé dans epsi1 C++++ RECUP DES INFOS EF ++++++++++++++++++++++++++++++++ c + OBTENUES PAR ELQUOI DANS RIGI1 PENDANT PHASE 1 C segment INFO deja actif dans RIGI1 c + rigi1 n appelle pas elquoi, c'est modeli qui l'a fait c mais du coup, on na pas de segment minte c (car depend de si pt de g pour rigi, pour sigma....) c c'est + simple de rappeler elquoi ici MELE = NEFMOD NGAU1 = INFELE(6) LRE = INFELE(9) LHOOK = INFELE(10) MINT1 = INFMOD(6) segact,MINT1 MINT2 = INFMOD(10) if(MINT2.ne.0) segact,MINT2 MFR = INFELE(13) IELE = INFELE(14) NDDL = INFELE(15) NSTRS = INFELE(16) c write(6,*)'-> EPSIX INFELE',(INFELE(iou),iou=1,16) c + AUTRES INFOS C nbre de noeuds par element NBNN1 = NUM(/1) C nbre d elements NBEL1 = NUM(/2) c REM: pour se passer du dimensionnement du nbre d'enrichissement dans c elquoi et le realiser localement , on pourrait ecrire: c LRE = NDDLMAX*NBNN1 c NDDL= NDDLMAX C sous decoupage et points de Gauss de l element geometrique de base CTY if(MELE.eq.263) then if(MELE.EQ.263.OR.MELE.EQ.264) then NGAU2 = MINT2.POIGAU(/1) endif c write(*,*) 'dim de MINT2=6,',(MINT2.SHPTOT(/2)),(MINT2.SHPTOT(/3)) c write(*,*) 'MINT2',(MINT2.QSIGAU(iou),iou=1,NGAU) c nbre maxi de fonction de forme par noeud (fonction std comprise) c NBNI = NDDL/IDIM inutile! C++++ Recup des infos d enrichissement +++++++++++++++++++ c recup du MCHEX1 d'enrichissement NOBMO1 = IVAMOD(/1) if(NOBMO1.ne.0) then do iobmo1=1,NOBMO1 if((TYMODE(iobmo1)).eq.'MCHAML') then MCHEX1 = IVAMOD(iobmo1) segact,MCHEX1 if((MCHEX1.TITCHE).eq.'ENRICHIS') then MCHAM1 = MCHEX1.ICHAML(1) segact,MCHAM1 goto 1000 endif endif enddo write(*,*) 'Le modele est vide (absence d enrichissement)' * return else write(*,*) 'Il n y a pas de MCHEML enrichissement dans le Modele' * return endif 1000 continue c niveau d enrichissement(s) du modele (ddls std u exclus) c NBENR1= 0 si std, 1 si H seul, 2 si H+F1, 3 si H+F1+F2, etc... if(NOBMO1.ne.0) then NBENR1= MCHAM1.IELVAL(/1) else NBENR1 = 0 endif c write(*,*) 'niveau d enrichissement(s) du modele',NBENR1 c C********************************************************* C INITIALISATIONS... C********************************************************* IRETER = 0 c c preparation des tables avec: if(NBENR1.ne.0) then do ienr=1,NBENR1 c -le nombre d'inconnues de chaque sous-zone c determinee depuis le nombre de fonction de forme c ienr= 1: U+H(1+1=2), 2: U+H+F1(2+4=6), 3: U+H+F1+F2(6+4=10) nbniJ = 2 + ((ienr-1)*4) MLRE(1+ienr) = IDIM*NBNN1*nbniJ enddo endif C Tables + longues car 1er indice correspond au fontion de forme std MLRE(1) = IDIM*NBNN1*1 if(NBENR1.lt.(NBENRMAX+1)) then do ienr=(NBENR1+1),(NBENRMAX) MLRE(1+ienr) = 0 enddo endif c c ...DU SEGMENT WRK1 NBENRMA2 = NBENRMAX NBBB = NBNN1 SEGINI,WRK1 c ...DU SEGMENT WRK2 c NBNO = NBNI NBNO = LRE/IDIM SEGINI,WRK2 C c ...DU SEGMENT MRACC NBENRMA2 = NBENRMAX NBI = NBNN1 segini,MRACC C C du nombre d erreur sur le noms de composantes NBERR1=0 C********************************************************* C C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> BOUCLE SUR LES ELEMENTS C DO 2000 J=1,NBEL1 c write(*,*) '===================================' c write(*,*) '========= EF',J,' /NBEL1 =========' C C C********************************************************* C POUR CHAQUE ELEMENT, ... C********************************************************* C C++++ ON RECUPERE LES COORDONNEES DES NOEUDS DE L ELEMENT CALL DOXE(XCOOR,IDIM,NBNN1,NUM,J,XE) C c C++++ NBENRJ = niveau d enrichissement de l element ++++ C =0 si EF std =1 si U+H =2 si U+H+F1 =3 si U+H+F1+F2 NBENRJ=0 if(NBENR1.ne.0) then do IENR=1,NBENR1 MELVA1 = MCHAM1.IELVAL(IENR) segact,MELVA1 do I=1,NBNN1 mlree1 = MELVA1.IELCHE(I,J) C on en profite pour remplir MRACC table de raccourcis pour cet element TLREEL(IENR,I) = mlree1 if(mlree1.ne.0) then NBENRJ = max(NBENRJ,IENR) C et on active la listreel segact,mlree1 endif enddo enddo endif if(NBENRMA2.gt.NBENR1) then do IENR2=(NBENR1+1),NBENRMA2 do I=1,NBNN1 TLREEL(IENR2,I) = 0 enddo enddo endif C c c++++ quelques indices pour dimensionner au plus juste c nbre total de ddl de l'élément considéré NLIGRD = MLRE(1+NBENRJ) NLIGRP = MLRE(1+NBENRJ) c nbre fonction de forme=((Ni_std+Ni_enrichi)*nbnoeud)=Ninconnues/idim NBNI = (MLRE(1+NBENRJ)) / IDIM c write(*,*) 'EF',J,' NBENRJ=',NBENRJ, c &'donc',NLIGRD,' ddls et ',NBNI,' fonctions de forme' c C++++ ON CALCULE XDDL ++++ MPTVAL = IVADEP NCOSOU = IVAL(/1) INITOT = 0 C-->> BOUCLE SUR LES niveaux d'ENRICHISSEMENTS (0:U 1:A 2:BCDE1 3:BCDE2) DO IENR=0,NBENRJ *nbre de fonction(s) de ce niveau d'enrichissement (=1 si std ou H, =4 pour F1,2,...) if(IENR .le. 1) then NBNIENR = 1 else NBNIENR = 4 endif C---->> BOUCLE SUR LES fonctions de forme de ce type d'enrichissement DO INI=1,NBNIENR INITOT = INITOT + 1 C------>> BOUCLE SUR LA DIMENSION DO 2220 KDIM=1,IDIM ICOMP = (INITOT-1)*IDIM + KDIM c --cas ou on n'a pas trouvé assez de composantes-- if(ICOMP.GT.NCOSOU) GOTO 2221 MELVAL = IVAL(ICOMP) c --cas ou on a pas trouvé cette composante dans cette zone du c chpoint solution devenu mchaml -- if(MELVAL.eq.0)then c Avait on besoin de cette composante? c oui, si c'est une composante obligatoire if(IENR.eq.0) goto 2991 c oui, si l'un des noeuds est enrichi do I=1,NBNN1 if(TLREEL(IENR,I).ne.0) goto 2991 enddo c non, si c'est facultatif et qu'on n'est pas enrichi -> on saute goto 2220 c ->AVERTISSEMENT puis on saute 2991 NBERR1=NBERR1+1 if(IIMPI.lt.1) goto 2220 c write(IOIMP,*) 'PB OPERATEUR SIGM :' write(IOIMP,991) ICOMP,IENR,INI,KDIM 991 format(2X,'ABSENCE DANS LE CHPOINT DEPLACEMENT DE LA ',I3, $ 'ieme COMPOSANTE (enrichissement',I3, $ ', fonction',I3,', direction ',I3, $ ') NECESSAIRE POUR L UN DES NOEUDS SUIVANTS :') write(IOIMP,*)' noeuds :',(NUM(iou,J),iou=1,NBNN1) goto 2220 endif C---------->> BOUCLE SUR LES NOEUDS DO I=1,NBNN1 IQ = ((INITOT-1)*NBNN1*IDIM) + ((I-1)*IDIM) + KDIM XDDL(IQ) = VELCHE(I,J) ENDDO 2220 CONTINUE ENDDO ENDDO c --cas normal (toutes les composantes souhaitees etaient presentes)-- GOTO 2223 c --cas ou on n'a pas trouvé assez de composantes-- 2221 CONTINUE if (IIMPI.ge.1) then WRITE(IOIMP,2222) J,NCOSOU,ICOMP 2222 FORMAT(2X,'ATTENTION : ELEMENT ',I6, & ' LE CHAMP DE DEPLACEMENT CONTIENT ',I3,' COMPOSANTES', & ' PAR NOEUD AU LIEU DE ',I3,' ATTENDUES') endif NDDL=NCOSOU*NBNN1 NBENRJ=IENR 2223 CONTINUE c --cas ou on a une ou des erreurs-- IF (NBERR1.gt.0.and.J.eq.NBEL1) THEN write(IOIMP,*) 'OPERATEUR GRAD : ABSENCE DANS LE CHPOINT ', & 'DEPLACEMENT DE CERTAINES INCONNUES ATTENDUES PAR LE MODELE' ENDIF c c rem: il serait probablement interessant au niveau du tems cpu c d'utiliser moins de pts de Gauss lorsque l element est élastique c On pourrait par ex. utiliser MINT2 = infele(12) qui contient c le segment d'integration de l'EF std (QUA4 par ex.) * if((NBENRJ.eq.0).and.(MINT2.ne.0)) then * MINTE = MINT2 * NBPGAU= NGAU2 * else MINTE = MINT1 NBPGAU= NGAU1 * endif c c pour les def quadratiques ISDJC=0 NBPTEL=NBPGAU IF (IREPS2.EQ.1) SEGINI MTRACE c C C********************************************************* C C>>>>>>>>>> BOUCLE SUR LES POINTS DE GAUSS C DO 2500 KGAU=1,NBPGAU C C********************************************************* C Initialisation à 0 C********************************************************* c ZERO ne serait-il pas facultatif? CALL ZERO(SHPWRK,6,NBNO) C i6zz = 3 IF (IDIM.EQ.3) i6zz = 4 c do ienr7=1,NBENRMAX do ienr7=1,NBENRJ do inod7=1,NBNN1 c do i6=1,6 do i6=1,i6zz LV7WRK(ienr7,1,i6,inod7) = 0 LV7WRK(ienr7,2,i6,inod7) = 0 enddo enddo enddo c********************************************************* c Calcul des fonction de forme std dans repere local c********************************************************* ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccccc c (et donc sur les Ni std) DO 2510 I=1,NBNN1 C++++ Calcul des Ni std c (rappel: 1:Ni 2:Ni,qsi 3:Ni,eta avec i=1,4) SHPWRK(1,I) = SHPTOT(1,I,KGAU) SHPWRK(2,I) = SHPTOT(2,I,KGAU) SHPWRK(3,I) = SHPTOT(3,I,KGAU) IF (IDIM.EQ.3) SHPWRK(4,I) = SHPTOT(4,I,KGAU) 2510 CONTINUE ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc c********************************************************* c Passage des fonctions de forme std dans repere global c********************************************************* C++++ CALCUL DES Ni,x Ni,y (i=1,4) + CALCUL DE det(J) CALL JACOBI(XE,SHPWRK,IDIM,NBNN1,DJAC) c if(J.eq.1.and.KGAU.eq.1) c &write(*,*) 'Ni(i=1,4)=',(SHPWRK(1,iou),iou=1,NBNN1) c********************************************************* c Si on est pas du tout enrichi on peut sauter une grosse partie c********************************************************* if(NBENRJ.eq.0) goto 2999 c********************************************************* c Calcul des level set + leurs derivees dans repere global c********************************************************* ccccc BOUCLE SUR LES enrichissements ccccccccccccccccccc do 2520 ienr=1,NBENRJ c MELVA1=MCHAM1.IELVAL(IENR) c segact,MELVA1 ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccc DO 2521 I=1,NBNN1 C++++ Le I eme noeud est-il ienr-enrichi? mlree1= TLREEL(IENR,I) if(mlree1.eq.0) goto 2521 c Calcul du repere local de fissure(=PSI,PHI) c (rappel: 1,1:psi 1,2:phi 2,1 psi,x ...etc...) do 2522 inode=1,NBNN1 c pour le H-enrichissement, on n a pas gardé PSI (inutile) if(ienr.ne.1) then c valeur de PSI au inode^ieme noeud xpsi1 = mlree1.PROG(inode) c qu on multiplie par la valeur de Ni^std au pt de G considéré LV7WRK(ienr,1,1,I)= LV7WRK(ienr,1,1,I) & + (SHPWRK(1,inode)*xpsi1) LV7WRK(ienr,1,2,I)= LV7WRK(ienr,1,2,I) & + (SHPWRK(2,inode)*xpsi1) LV7WRK(ienr,1,3,I)= LV7WRK(ienr,1,3,I) & + (SHPWRK(3,inode)*xpsi1) IF (IDIM.EQ.3) LV7WRK(ienr,1,4,I)= LV7WRK(ienr,1,4,I) & + (SHPWRK(4,inode)*xpsi1) c valeur de PHI au inode^ieme noeud xphi1 = mlree1.PROG(NBNN1+inode) else xphi1 = mlree1.PROG(inode) endif LV7WRK(ienr,2,1,I)= LV7WRK(ienr,2,1,I) & + (SHPWRK(1,inode)*xphi1) LV7WRK(ienr,2,2,I)= LV7WRK(ienr,2,2,I) & + (SHPWRK(2,inode)*xphi1) LV7WRK(ienr,2,3,I)= LV7WRK(ienr,2,3,I) & + (SHPWRK(3,inode)*xphi1) IF (IDIM.EQ.3) LV7WRK(ienr,2,4,I)= LV7WRK(ienr,2,4,I) & + (SHPWRK(4,inode)*xphi1) 2522 continue 2521 continue ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc 2520 CONTINUE ccccc fin de BOUCLE SUR LES enrichissements cccccccccccccccc c on a construit C LV7WRK(ienr, PSI/PHI, valeur/deriveeparqsi/pareta, iNOEUD) c********************************************************* c Ajout des fonctions de forme d enrichissement c + leurs derivees dans repere global c********************************************************* CALL SHAPX(MELE,LV7WRK,NBENRMAX,NBENRJ,TLREEL,SHPWRK,IRET) c retour a la partie commune aux EF enrichi et non enrichi 2999 continue C********************************************************* C CALCUL DE B' C********************************************************* c ZERO ne serait-il pas facultatif? c call ZERO(BGENE,LHOOK,NLIGRP) KB=1 C boucle sur tous les Ni DO 3001 II=1,NBNI BGENE(1,KB) = SHPWRK(2,II) BGENE(2,KB+1) = SHPWRK(3,II) BGENE(4,KB) = SHPWRK(3,II) BGENE(4,KB+1) = SHPWRK(2,II) IF(IDIM.EQ.3) THEN BGENE(3,KB+2)=SHPWRK(4,II) BGENE(5,KB)=SHPWRK(4,II) BGENE(5,KB+2)=SHPWRK(2,II) BGENE(6,KB+1)=SHPWRK(4,II) BGENE(6,KB+2)=SHPWRK(3,II) ENDIF KB = KB + IDIM 3001 CONTINUE C c if(J.eq.5.and.KGAU.eq.1) then c if(KGAU.eq.1) then c write(*,*) 'BGENE(1,..)=',(BGENE(1,iou),iou=1,2*NBNI) c write(*,*) 'BGENE(2,..)=',(BGENE(2,iou),iou=1,2*NBNI) c write(*,*) 'BGENE(4,..)=',(BGENE(3,iou),iou=1,2*NBNI) c endif C********************************************************* C CALCUL DE D C********************************************************* c on cherche les matrices de Hooke MPTVAL=IVAMAT IF (IMAT.EQ.2) THEN MELVAL=IVAL(1) IBMN=MIN(J ,IELCHE(/2)) IGMN=MIN(KGAU,IELCHE(/1)) MLREEL=IELCHE(IGMN,IBMN) SEGACT MLREEL IF (KGAU.LE.NBGMAT.AND.(J.LE.NELMAT.OR.NBGMAT.GT.1)) 1 CALL DOHOOO(PROG,LHOOK,DDHOOK) ELSE IF (IMAT.EQ.1) THEN DO 9004 IM=1,NMATT IF (IVAL(IM).NE.0) THEN MELVAL=IVAL(IM) IBMN=MIN(J ,VELCHE(/2)) IGMN=MIN(KGAU,VELCHE(/1)) VALMAT(IM)=VELCHE(IGMN,IBMN) ELSE VALMAT(IM)=0.D0 ENDIF 9004 CONTINUE c IF (MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN * IF (IGAU.LE.NBGMAT) * 1 CALL DOHMAO(VALMAT,CMATE,IFOUR,IDIM,TXR,XLOC,XGLOB,D1HO, * 2 ROTH,DDHOOK,LHOOK,1,IRTD) write(*,*) 'cas non prevu a ce jour pour les EF xfem' return ELSE IF (KGAU.LE.NBGMAT.AND.(J.LE.NELMAT.OR.NBGMAT.GT.1)) 1 CALL DOHMAS(VALMAT,CMATE,IFOUR,LHOOK,1,DDHOOK,IRTD) ENDIF ENDIF c c C********************************************************* C CALCUL DE sigma = D*B*u C********************************************************* c CALL DBST(BGENE,DDHOOK,XDDL,(NBNI*IDIM),LHOOK,XSTRS) c c calcul des eps 2 (termes quadratiques de la deformation) c IF (IREPS2.EQ.1) 1 CALL DBST2(SHPWRK,DDHOOK,XDDL,XE,NBNI,IFOUR,LHOOK,XSTRS, 2 TRACE,KGAU,XDPGE,YDPGE,UDPGE,NIFOUR) c c remplissage du segment contenant les contraintes c MPTVAL=IVASTR DO 7004 ICOMP=1,LHOOK MELVAL=IVAL(ICOMP) IBMN=MIN(J,VELCHE(/2)) VELCHE(KGAU,IBMN)=XSTRS(ICOMP) 7004 CONTINUE c * if(J.eq.5 .and. KGAU.eq.1) then c if(KGAU.eq.1) then c write(*,*) J,KGAU,'sigma(..)=',(XSTRS(iou),iou=1,LHOOK) c endif c 2500 CONTINUE C FIN DE BOUCLE SUR LES POINTS DE GAUSS <<<<<<<<<<<<<< C c quelques suppressions c (ici element non-incompressible=> pas besoin de MTRACE (cf epsi2) IF (IREPS2.EQ.1) THEN SEGSUP MTRACE ENDIF c 2000 CONTINUE C FIN DE BOUCLE SUR LES ELEMENTS <<<<<<<<<<<<<<<<<<<<< c C********************************************************* C SUPPRESSION ET DESACTIVATION DE SEGMENTS C********************************************************* C SEGSUP,WRK1,WRK2,MVELCH SEGSUP,MRACC RETURN END