C TRIELE SOURCE PV090527 26/04/30 21:16:43 12529 C======================================================================= C CREATION : bp, le 18.10.2005 C MODIFS : -modif de la projection : bp, le 07/02/2006 c -refonte de l'operateur pour integration dans Cast3M (aout2008) c -mise à 0. de Phi et Psi moins aleatoire (septembre 2009) c -syntaxe 2 (bp, juin 2010) c -ajout option 'DESENRICHISSEMENT' (bp, avril 2012) c -ajout SURE, GG 2017-2019 c C======================================================================= C FONCTION : C TRIE LES ELEMENTS SELON LES ISOZEROS DES CHPOINTS LEVEL SET C C======================================================================= c SYNTAXE 1 : c chelX_n+1 | | = TRIE mod_n psi_n+1 phi_n+1 | | ; c | | | 'SAUT' | c | relaX | | 'DESENR' | c c (mod_n est actualisé --> mod_n+1) C C E: phi_n+1 = MCHPOI distance au plan de fissure a linstant n+1 C E: psi_n+1 = MCHPOI distance au front de fissure a linstant n+1 C E: mod_n = MMODEL a linstant n C (contient eventuellement une zone avec des element X-FEM C et IVAMOD pointe alors vers chelX_n) C E: 'SAUT' = mot cle optionnel permettant de limiter le niveau C d'enrichissement a H. La pointe doit etre sur un bord C d'element C E: 'DESENR' = mot cle optionnel permettant de 'desenrichir' le modele C S: mod_n+1 = MMODEL a linstant n+1 (mise a jour de IVAMOD) C S: chelX_n+1 = MCHELM d'enrichissement a linstant n+1 C (construit d'apres les iso-zero de phi et psi) C S: relaX : relation associée au desenrichissement c (mise a zero des 'pointes devenues obsolètes') C C======================================================================= c SYNTAXE 2 : c TRIE mod_? chelX_n+1 c mod_n+1 c C E: mod_? = MMODEL contenant ou pas des enrichissements C E: chelX_n+1 = MCHELM d'enrichissement a l instant n+1 C S: mod_n+1 = MMODEL pointant vers chelX_n+1 (inscrit dans IVAMOD) C C======================================================================= C SUBROUTINE TRIELE c C======================================================================= c IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C SEGMENTS INCLUDE -INC PPARAM -INC CCOPTIO -INC SMCOORD -INC SMELEME -INC SMCHPOI -INC SMMODEL -INC SMCHAML -INC SMLREEL -INC SMINTE -INC SMRIGID LOGICAL LOG1 CHARACTER*(LOCOMP) nomc1,nomc2 C CHARACTER*(LOCOMP) char1 PARAMETER (NMOT0=2) CHARACTER*4 MOTCLE(NMOT0) DATA MOTCLE/'SAUT','DESE'/ PARAMETER (NDDL1=12) CHARACTER*4 MOTPR1(NDDL1),MOTDU1(NDDL1) DATA MOTPR1/'B1X ','B1Y ','B1Z ','C1X ','C1Y ','C1Z ', >'D1X ','D1Y ','D1Z ','E1X ','E1Y ','E1Z '/ DATA MOTDU1/'FB1X','FB1Y','FB1Z','FC1X','FC1Y','FC1Z', >'FD1X','FD1Y','FD1Z','FE1X','FE1Y','FE1Z'/ CHARACTER*4 MOTPR2(NDDL1),MOTDU2(NDDL1) DATA MOTPR2/'B2X ','B2Y ','B2Z ','C2X ','C2Y ','C2Z ', >'D2X ','D2Y ','D2Z ','E2X ','E2Y ','E2Z '/ DATA MOTDU2/'FB2X','FB2Y','FB2Z','FC2X','FC2Y','FC2Z', >'FD2X','FD2Y','FD2Z','FE2X','FE2Y','FE2Z'/ C SEGMENT ITYPND(NI) SEGMENT NDPT1(NI) SEGMENT NDPT2(NI) SEGMENT NUMND(NI) SEGMENT KKGAU(NI) SEGMENT IDEJVU(NBPTI) c segment pour tableau de travail de reel SEGMENT XPSI(NI) SEGMENT XPHI(NI) c autrespointeurs POINTEUR MCHEX1.MCHELM,MCHEX2.MCHELM POINTEUR IPTX1.MELEME POINTEUR MINT1.MINTE,MINT2.MINTE,MINT0.MINTE C======================================================================= c NDPT1=0 NDPT2=0 IDIM1 = IDIM + 1 segact,MCOORD*mod C======================================================================= c c LECTURE DES ARGUMENTS c C======================================================================= C++++ le modele a l'instant n CALL LIROBJ('MMODEL ',MMODE1,1,IRETOU) IF (IERR .NE. 0) RETURN CALL ACTOBJ('MMODEL ',MMODE1,1) c++++ L'option pour le niveau d'enrichissement max souhaité *as 2010_01_20 : selon l'option le niveau d'enrichissement est limité à 1 ou non: lsaut=100 c CALL LIRCHA(char1,0,iretou) c if (iretou.ne.0) then c if (char1.ne.'SAUT') then c MOTERR=char1 c CALL ERREUR(7) c RETURN c else c lsaut=1000 c endif c endif *fin as 2010_01_20 CALL LIRMOT(MOTCLE,NMOT0,IMOT0,0) if(IMOT0.eq.1) lsaut=1000 c++++ lecture d un MCHEX1 d'enrichissement a brancher CALL LIROBJ('MCHAML',IPIN,0,ISYNTAX) IF (IERR .NE. 0) RETURN ISYNTAX=ISYNTAX+1 C (syntaxe 2) if (ISYNTAX .eq. 2) then CALL ACTOBJ('MCHAML',IPIN,1) CALL REDUAF(IPIN,MMODE1,MCHEX1,0,IR,KER) IF(IR .NE. 1) CALL ERREUR(KER) IF(IERR .NE. 0) RETURN segact,MCHEX1 c verif : il doit etre de type enrichissement if ((MCHEX1.TITCHE).ne.'ENRICHIS') then MOTERR='MCHAML' MOTERR(9:16)='ENRICHIS' c "On attendait un %M1:8 de %m9:24" CALL ERREUR(1025) return endif c verif : il doit comporter au moins 1 sous-zone if(MCHEX1.ICHAML(/1) .EQ. 0) THEN if (impi.ge.1) write(ioimp,*) & 'MCHAML doit avoir au moins 1 zone !' c "On n'a pas trouve, dans un MCHAML, de zone geometrique ou de c constituant correspondant a l'objet MODELE" CALL ERREUR(472) RETURN ENDIF C++++ OU lecture des Level Set (syntaxe 1) else CALL LIROBJ('CHPOINT ',MCHPO1,1,IRETOU) CALL ACTOBJ('CHPOINT ',MCHPO1,1) CALL LIROBJ('CHPOINT ',MCHPO2,1,IRETOU) CALL ACTOBJ('CHPOINT ',MCHPO2,1) if (IERR.ne.0) then CALL ERREUR(181) c write(*,*) 'Besoin des chpoints Level Sets' c return endif C qq verif de base sur les chpoints segact MCHPO1,MCHPO2 nsoup1 = MCHPO1.IPCHP(/1) nsoup2 = MCHPO2.IPCHP(/1) if ((nsoup1.ne.1).or.(nsoup2.ne.1)) then INTERR = MCHPO1 c CALL ERREUR(377) c write(*,*) 'Chpoint avec plus d 1 zone: Non traite test1' c return endif MSOUP1 = MCHPO1.IPCHP(1) MSOUP2 = MCHPO2.IPCHP(1) C composante des chpoints: 1=PSI 2=PHI segact MSOUP1,MSOUP2 nc1 = MSOUP1.NOCOMP(/2) nc2 = MSOUP2.NOCOMP(/2) if ((nc1.ne.1).or.(nc2.ne.1)) then CALL ERREUR(181) c write(*,*) 'Chpoint avec plus d 1 composante: Non traite' c return endif nomc1 = MSOUP1.NOCOMP(1) nomc2 = MSOUP2.NOCOMP(1) if (nomc1.ne.'PSI ') then MOTERR = 'PSI' CALL ERREUR(243) c write(*,*) 'Chpoint 1 doit avoir pour composante PSI ' c return endif if (nomc2.ne.'PHI ') then MOTERR = 'PHI' CALL ERREUR(243) c write(*,*) 'Chpoint 2 doit avoir pour composante PHI ' c return endif C recuperation des maillages et valeurs IPT1 = MSOUP1.IGEOC IPT2 = MSOUP2.IGEOC c if (IPT1.ne.IPT2) then c write(*,*) 'Chpoints doivent avoir le meme support maillage' c return c bp (07/07/2011) : trop contraignant on introduit ndpt2 et on c va faire des test + tard c endif MPOVA1 = MSOUP1.IPOVAL MPOVA2 = MSOUP2.IPOVAL C on crée MPOVA3 et MPOVA4 = versions "simplifiées" de MPOVA1 et MPOVA2 SEGINI,MPOVA3=MPOVA1 SEGINI,MPOVA4=MPOVA2 endif C++++ fin distinction syntaxe 2/1 C Initialisation des tableaux de sortie info erichissement a C chaques noeuds NI = nbpts SEGINI,ITYPND C======================================================================= c C PREMIERE BOUCLE SUR LES SOUS MODELES c c cas syntaxe 1 : on cree un nouveau MCHEX2+MCHAM2 a partir de : c - copie de l'ancien MCHEX1+MCHAM1 existant c - ou nouveau MCHEX1+MCHAM1 vide c si desenrichissement, on cree les relations IPRIG2 c c cas syntaxe 2 : on branche le MCHEX1 fourni c C======================================================================= SEGACT,MMODE1 NBSMOD = MMODE1.KMODEL(/1) DO 100 ISMOD1 = 1,NBSMOD c CcccccC ouverture du sous modele IMODEL = MMODE1.KMODEL(ISMOD1) segact IMODEL*mod c on travaille seulement si EF XFEM if(NEFMOD.ne.263.and.NEFMOD.ne.264) goto 100 NOBMO1 = IVAMOD(/1) ccccccC branchement d 1 eventuel MCHEX1 d'enrichissement fourni (syntaxe 2) if (ISYNTAX.EQ.2) then NOBMO1 = IVAMOD(/1) c NOBMOD = NOBMO1 + 1 c if(NOBMO1.ne.0) write(6,*)'Attention: Les objets contenus', c & ' dans le IVAMOD du modele seront oubliés !!!' NOBMOD=1 cbp: attention: on ecrase tout pour l instant MN3 = INFMOD(/1) NFOR = FORMOD(/2) NMAT = MATMOD(/2) segadj,IMODEL TYMODE(NOBMOD)='MCHAML ' IVAMOD(NOBMOD)= MCHEX1 goto 100 endif ccccccC (syntaxe 1) : ccccccC recup du MCHEX1 d'enrichissement pre-existant 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 isouX=1 MCHAM1 = MCHEX1.ICHAML(isouX) cbp: on suppose donc l'exitence d'autant de MCHEX1 (avec 1 seul MCHAM1) c que de zones elementaires, au lieu de 1 MCHEX1 avec autant de MCHAM1 c que de zones elementaires... il faudrait corriger ca tres vite... segact,MCHAM1 goto 101 endif endif enddo endif ccccccC OU creation d un MCHEX1 d'enrichissement vierge c write(*,*) 'Le modele ne contient pas le MCHELM ENRICHIS' c write(*,*) '=> On initialise a 0 le MCHELM ENRICHIS' L1 = 8 N1 = 1 N3 = 6 SEGINI,MCHEX1 MCHEX1.TITCHE = 'ENRICHIS' MCHEX1.IFOCHE = IFOUR MCHEX1.IMACHE(1)= IMAMOD MCHEX1.CONCHE(1)= IMODEL.CONMOD MCHEX1.INFCHE(1,2) = 0 MCHEX1.INFCHE(1,3) = NIFOUR MCHEX1.INFCHE(1,6) = 1 *as 2010_01_20: * N2 = 2 N2 = 1 SEGINI,MCHAM1 MCHEX1.ICHAML(1) = MCHAM1 MELEME=IMAMOD SEGACT,MELEME N1PTEL = 0 N1EL = 0 N2PTEL = NUM(/1) N2EL = NUM(/2) MCHAM1.NOMCHE(1)= 'H ' *as 2010_01_20: * MCHAM1.NOMCHE(2)= 'F1 ' DO I2=1,N2 MCHAM1.TYPCHE(I2)= 'POINTEURLISTREEL' SEGINI,MELVA1 MCHAM1.IELVAL(I2) = MELVA1 ENDDO c on ajoute une petite place dans IVAMOD NOBMOD=NOBMO1+1 MN3=INFMOD(/1) NFOR=FORMOD(/2) NMAT=MATMOD(/2) segadj,IMODEL iobmo1=NOBMOD ccccccC point de rassemblement une fois obtenu MCHEX1 et MCHAM1 101 continue c qq verif de base NCHAM1 = MCHEX1.ICHAML(/1) if (NCHAM1.ne.1) then INTERR(1) = MCHEX1.ICHAML c CALL ERREUR(377) c write(*,*) 'MCHELM doit avoir 1 seule zone test2' endif c recup du maillage et du MCHAML associé a cette zone MELEME = MCHEX1.IMACHE(1) if (MELEME.ne.IMAMOD) then if (impi.ge.1) write(ioimp,*) & 'MCHAML et MODELE doivent avoir le meme support !' c "Geometrie incompatible avec le MCHAML" call ERREUR(706) return endif segact,MELEME MCHAM1 = MCHEX1.ICHAML(1) c recup de l enrichissement précédent (ITIP1) et calcul suivant (ITIP2) ITIP1=MCHEX1.INFCHE(1,2) if(ITIP1.lt.0.or.ITIP1.gt.2) then if (impi.ge.1) write(ioimp,*) 'Pb : MCHEX1.INFCHE(1,2)=',ITIP1 call erreur(5) return endif c cas option 'SAUT' (H-only) if(IMOT0.eq.1) then ITIP2=0 c autre cas (H, F1 et F2 enrichissement) : on permute else if(ITIP1.eq.1) then ITIP2=2 else c cas ou ITIP1 = 0 (initial) ou 2 ITIP2=1 endif endif C======================================================= c write(*,*) 'CREATION D UN JOLI MCHEX2 identique au MCHEX1' c et eventuellement d une matrice de relation pour désenrichir C======================================================= SEGINI,MCHEX2=MCHEX1 IVAMOD(iobmo1)=MCHEX2 TYMODE(iobmo1)='MCHAML ' MELE = NEFMOD SEGINI,MCHAM2=MCHAM1 MCHEX2.ICHAML(1)=MCHAM2 MCHEX2.INFCHE(1,2)=ITIP2 N2ENR = MCHAM2.IELVAL(/1) do i2enr=1,N2ENR MELVA1 = MCHAM1.IELVAL(i2enr) SEGINI,MELVA2=MELVA1 MCHAM2.IELVAL(i2enr) = MELVA2 enddo c --- DESENRICHISSEMENT --- IPRIG2=0 if(IMOT0.eq.2) then c -On "oublie" les très vieux enrichissements devenus obsolètes c et qu on va réutiliser pour la nouvelle pointe de fissure if(MCHAM2.IELVAL(/1).ge.ITIP2+1) then MELVA2=MCHAM2.IELVAL(ITIP2+1) segact,MELVA2*mod do j=1,MELVA2.IELCHE(/2) do i=1,MELVA2.IELCHE(/1) MELVA2.IELCHE(i,j)=0 enddo enddo endif c -On met a 0 (avec une relation) les vieux enrichissements if(MCHAM2.IELVAL(/1).ge.ITIP1+1) then c On commence par creer un maillage prototype de la rigidite (IPT6) c en reperant les noeuds a traiter via idejvu NBPTI=NBPTS segini,idejvu NBSOUS=0 NBREF=0 NBNN=2 NBELEM=1000 segini,IPT6 IPT6.ITYPEL=22 c on noublie pas de créer des noeuds pour les LX NBPTS0=NBPTS iilx=NBPTS0 NBPTS=NBPTS+NBELEM segadj,MCOORD c on travaille sur le melva2 de la pointe précédente MELVA2=MCHAM2.IELVAL(ITIP1+1) segact,MELVA2*mod ii2=0 do 120 j=1,MELVA2.IELCHE(/2) do 121 i=1,MELVA2.IELCHE(/1) i120 = NUM(i,j) c rem: commme l enrichissement est nodal et que l on travaille sur 1 seul c enrichissement a annuler on se permet les lignes ci dessous if(idejvu(i120).gt.0) goto 121 idejvu(i120)=1 if (MELVA2.IELCHE(i,j).ne.0) then ii2=ii2+1 if(ii2.gt.NBELEM) then NBELEM=NBELEM + 1000 NBPTS=NBPTS + 1000 segadj,IPT6,MCOORD endif c on crée un noeud pour le multiplicateur de lagrange LX c iilx = NBPTS0+ii2 iilx = iilx+1 XCOOR((iilx-1)*(IDIM1)+1) = XCOOR((i120-1)*(IDIM1)+1) XCOOR((iilx-1)*(IDIM1)+2) = XCOOR((i120-1)*(IDIM1)+2) if(idim.eq.3) & XCOOR((iilx-1)*(IDIM1)+3) = XCOOR((i120-1)*(IDIM1)+3) XCOOR(iilx*(IDIM1)) = XCOOR(i120*(IDIM1)) IPT6.NUM(1,ii2)=iilx IPT6.NUM(2,ii2)=i120 endif 121 continue 120 continue c on ajuste IPT6 et on supprime NBELEM=ii2 segadj,IPT6 segsup,idejvu c on finit d ajuster MCOORD NRIGEL=idim*4 NBPTS = NBPTS0+(ii2*NRIGEL) segadj,MCOORD c On crée la relation qui permet de désenrichir segini,MRIGID IPRIG2=MRIGID IFORIG = IFOUR MTYMAT ='RIGIDITE' NLIGRD=2 NLIGRP=2 NELRIG=ii2 jb=0 jdim=0 do iriri=1,NRIGEL COERIG(iriri)=1.d0 rigrel=0 segini,DESCR,XMATRI c IRIGEL(1,iriri)=IPT6 c On duplique le maillage prototype de la rigidite if(iriri.eq.1) then IPT7=IPT6 else segini,IPT7=IPT6 c on crée un noeud pour le multiplicateur de lagrange LX c pour chaque element de ipt7 do ii2=1,NBELEM iilx = iilx+1 c iilx6=IPT6.NUM(1,ii2) on prend ipt7 actif et pas ipt6 iilx6=IPT7.NUM(1,ii2) XCOOR((iilx-1)*(IDIM1)+1) = XCOOR((iilx6-1)*(IDIM1)+1) XCOOR((iilx-1)*(IDIM1)+2) = XCOOR((iilx6-1)*(IDIM1)+2) if(idim.eq.3) & XCOOR((iilx-1)*(IDIM1)+3) = XCOOR((iilx6-1)*(IDIM1)+3) XCOOR(iilx*(IDIM1)) = XCOOR(iilx6*(IDIM1)) IPT7.NUM(1,ii2)=iilx enddo endif IRIGEL(1,iriri)=IPT7 IRIGEL(2,iriri)=0 IRIGEL(3,iriri)=DESCR IRIGEL(4,iriri)=XMATRI IRIGEL(5,iriri)=NIFOUR IRIGEL(6,iriri)=0 IRIGEL(7,iriri)=0 IRIGEL(8,iriri)=0 c remplissage du DESCR LISINC(1)='LX' LISDUA(1)='FLX' jb=jb+1 jdim=jdim+1 if(ITIP1.eq.1) then LISINC(2)=MOTPR1(jb) LISDUA(2)=MOTDU1(jb) elseif(ITIP1.eq.2) then LISINC(2)=MOTPR2(jb) LISDUA(2)=MOTDU2(jb) endif c petit ajustement pour "sauter" les composantes en Z if(jdim.eq.idim) then jdim=0 if(idim.eq.2) jb=jb+1 endif do jinc=1,NLIGRP NOELEP(jinc)=jinc NOELED(jinc)=jinc enddo c remplissage de XMATRI do jel=1,NELRIG RE(1,1,jel)=0.d0 RE(1,2,jel)=1.d0 RE(2,1,jel)=1.d0 RE(2,2,jel)=0.d0 enddo segdes,DESCR,XMATRI enddo segdes,MRIGID 139 continue endif endif c --- fin du cas DESENRICHISSEMENT --- C======================================================= c BOUCLE SUR LES POINTS DES CHPOINTS c pour creer une numerotation locale c --> NDPT1/2(#global) = #chpo1/2 C======================================================= c -chpoint psi segact,IPT1 NBNO1 = IPT1.NUM(/2) C on crée la table qui depuis num de noeud donne point dans chpoint NI = nbpts SEGINI,NDPT1 cccccc Vrai Debut de boucle DO I1=1,NBNO1 C on recupere le numero de noeud INUM1 = IPT1.NUM(1,I1) C on remplit la table NDPT1(INUM1) = I1 ENDDO c -chpoint phi if (IPT1.eq.IPT2) then NDPT2=NDPT1 else segact,IPT2 NBNO2 = IPT2.NUM(/2) if (NBNO1.ne.NBNO2) then MOTERR='CHPOINT POINTS !' c "Les deux objets de type %m1:8 n'ont pas le meme nombre de %m9:16" call ERREUR(403) return endif SEGINI,NDPT2 DO I2=1,NBNO2 C on recupere le numero de noeud INUM2 = IPT2.NUM(1,I2) if (NDPT1(INUM2).eq.0) then IF(IMPI.GE.1) write(IOIMP,*) 'Les chpoints PSI et PHI ', & 'doivent avoir les memes points suuports !' c "Les champs par point ont des supports geometriques incompatibles" call ERREUR(348) return endif C on remplit la table NDPT2(INUM2) = I2 ENDDO endif C=======FIN DE BOUCLE SUR LES POINTS DES CHPOINTS======= C======================================================= c 1ere BOUCLE SUR LES ELEMENTS c POUR LA "SIMPLIFICATION" des level sets : c MPOVA1/2 --> MPOVA3/4 C======================================================= C Ouverture du maillage et des chpoints level set c SEGACT MELEME,MPOVA1,MPOVA2 SEGACT MPOVA1,MPOVA2 NBNN1 = NUM(/1) NBEL1 = NUM(/2) c on recupere les coord dans l element parents des noeuds MINT0 = INFMOD(3) segact,MINT0 c on recupere les coord dans l element parents des points de gauss MINT1 = INFMOD(5) segact,MINT1 c segment integration ef std correspondant -> inutile pour l instant c MINT2 = INFELE(12) c if(MINT2.ne.0) segact,MINT2 c recherche du num du pt de gauss KKGAU le plus proche du noeud inode1 NGAU1 = INFELE(6) NI = NBNN1 segini,KKGAU do inode1=1,NBNN1 kgau1 = 1 DXi1 = ((MINT1.QSIGAU(kgau1)) - (MINT0.QSIGAU(inode1)))**2 & + ((MINT1.ETAGAU(kgau1)) - (MINT0.ETAGAU(inode1)))**2 & + ((MINT1.DZEGAU(kgau1)) - (MINT0.DZEGAU(inode1)))**2 do kgau2=2,NGAU1 DXi2 = ((MINT1.QSIGAU(kgau2)) - (MINT0.QSIGAU(inode1)))**2 & + ((MINT1.ETAGAU(kgau2)) - (MINT0.ETAGAU(inode1)))**2 & + ((MINT1.DZEGAU(kgau2)) - (MINT0.DZEGAU(inode1)))**2 if (DXi2.lt.DXi1) then DXi1 = DXi2 kgau1 = kgau2 endif enddo KKGAU(inode1) = kgau1 enddo NI = NBEL1 C Initialisation des tableaux de travail NI = NBNN1 segini,NUMND,XPSI,XPHI c++++++ Debut de boucle sur les elements DO 200 J=1,NBEL1 C++++++++ Mini-Boucle(S) sur les noeuds de l element DO K=1,NBNN1 C numero INUM du Kieme noeud du Jieme element INUM = NUM(K,J) C numero du point correspondant des chpoints I1 = NDPT1(INUM) I2 = NDPT2(INUM) c on verifie que le chpoint y est bien defini sinon on saute tout! if(I1.eq.0) goto 200 C valeurs de psi et phi au Kieme noeud XPSI(K) = MPOVA1.VPOCHA(I1,1) XPHI(K) = MPOVA2.VPOCHA(I2,1) ENDDO C++++++++ Fin de la Mini-Boucle sur les noeuds de l element C++++++++ Mini-Boucle(S) sur les noeuds de l element DO K=1,NBNN1 C on calcule phi au pts de gauss le + proche du noeud considéré kgau1 = KKGAU(K) XPSIGAU = 0.d0 XPHIGAU = 0.d0 do ii1=1,NBNN1 XPSIGAU = XPSIGAU & + ( (MINT1.SHPTOT(1,ii1,kgau1)) * (XPSI(ii1)) ) XPHIGAU = XPHIGAU & + ( (MINT1.SHPTOT(1,ii1,kgau1)) * (XPHI(ii1)) ) enddo C numero INUM du Kieme noeud du Jieme element INUM = NUM(K,J) C numero du point correspondant des chpoints I1 = NDPT1(INUM) I2 = NDPT2(INUM) c C simplification c IF (((XPSI(K))*XPSIGAU).LE.0.d0) THEN c c XPSI(K) = 0.d0 c MPOVA3.VPOCHA(I1,1) = 0.d0 c ELSEIF (ABS(XPSI(K)).LE.XTOL) THEN c XPSI(K) = 0.d0 c MPOVA3.VPOCHA(I1,1) = 0.d0 c ENDIF c IF (((XPHI(K))*XPHIGAU).LE.0.d0) THEN c c XPHI(K) = 0.d0 c MPOVA4.VPOCHA(I2,1) = 0.d0 c ELSEIF (ABS(XPHI(K)).LE.XTOL) THEN c c XPHI(K) = 0.d0 c MPOVA4.VPOCHA(I2,1) = 0.d0 c ENDIF C simplification (bp: 13./07/2011) IF(((XPSI(K))*XPSIGAU).LE.0.d0) & MPOVA3.VPOCHA(I1,1) = 0.d0 IF(((XPHI(K))*XPHIGAU).LE.0.d0) & MPOVA4.VPOCHA(I2,1) = 0.d0 ENDDO C++++++++ Fin de la Mini-Boucle sur les noeuds de l element 200 CONTINUE C=======FIN DE la 1ere BOUCLE SUR LES ELEMENTS========== C======================================================= c 2EME BOUCLE SUR LES ELEMENTS C POUR LE TRI DES NOEUDS ("TEST EN PAGAILLE") c --> ITYPND(#global) = 0 : noeud pas vu c 1 : standard (non enrichi) c 10 : H-enrichi c 100 : F-enrichi C======================================================= c++++++ Debut de boucle sur les elements DO 222 J=1,NBEL1 C++++++++ Mini-Boucle(S) sur les noeuds de l element DO K=1,NBNN1 C numero INUM du Kieme noeud du Jieme element INUM = NUM(K,J) C numero du point correspondant des chpoints I1 = NDPT1(INUM) I2 = NDPT2(INUM) c on verifie que le chpoint y est bien defini sinon on saute tout! if(I1.eq.0) goto 222 C numero INUM du Kieme noeud de l'element en cours (le Jieme) NUMND(K) = INUM C valeurs de psi et phi au Kieme noeud XPSIK = MPOVA3.VPOCHA(I1,1) XPHIK = MPOVA4.VPOCHA(I1,1) XPSI(K) = XPSIK XPHI(K) = XPHIK C Calcul de valeurs utiles POUR LE TRI if (K.eq.1) then PSIMAX = XPSIK PSIMIN = XPSIK PSIMIA = ABS(XPSIK) PHIMAX = XPHIK PHIMIN = XPHIK PHIMIA = ABS(XPHIK) else PSIMAX = MAX(PSIMAX,XPSIK) PSIMIN = MIN(PSIMIN,XPSIK) PSIMIA = MIN(PSIMIA,ABS(XPSIK)) PHIMAX = MAX(PHIMAX,XPHIK) PHIMIN = MIN(PHIMIN,XPHIK) PHIMIA = MIN(PHIMIA,ABS(XPHIK)) endif ENDDO C++++++++ Fin de la Mini-Boucle sur les noeuds de l element PSIPRD = PSIMAX * PSIMIN PHIPRD = PHIMAX * PHIMIN c C C++++++++ Tests en pagaille pour trier c LOG1 = Tous les noeuds de lelement sont situés en arriere du front? LOG1 = (PSIMAX.LE.0.d0) c ...noeuds en arriere du front + element coupé par la fissure IF (LOG1.AND.(PHIPRD.LT.0.d0)) THEN C => tous les noeud doivent etre H-enrichi (=10) DO K=1,NBNN1 INUM = NUMND(K) ITYPND(INUM) = MAX(10,ITYPND(INUM)) ENDDO C c ...noeuds en arriere du front + element tangent avec la fissure ELSEIF (LOG1.AND.(PHIMIA.EQ.0.d0)) THEN C => seuls les noeud tangents doivent etre H-enrichi (=10) DO K=1,NBNN1 IF (XPHI(K).EQ.0.d0) THEN INUM = NUMND(K) ITYPND(INUM) = MAX(10,ITYPND(INUM)) ENDIF ENDDO C c Element contenant la pointe de fissure ELSEIF ((PSIPRD.LT.0.d0).AND.(PHIPRD.LT.0.d0)) THEN C => tous les noeud doivent etre F-enrichi (=100) DO K=1,NBNN1 INUM = NUMND(K) ITYPND(INUM) = 100 ENDDO C c Element dont un bord contient la pointe de fissure ELSEIF ((PSIPRD.LT.0.d0).AND.(PHIMIA.EQ.0.d0)) THEN C => tous les noeud doivent etre F-enrichi (=100) DO K=1,NBNN1 INUM = NUMND(K) ITYPND(INUM) = 100 ENDDO C c La pointe arrive sur un bord de l element ELSEIF ((PSIMIA.EQ.0.d0).AND.(PHIPRD.LT.0.d0)) THEN C => choix: seuls les noeud en avant du front sont F-enrichi (=100) DO K=1,NBNN1 INUM = NUMND(K) if (XPSI(K).GE.0.d0) ITYPND(INUM) = 100 ENDDO C c Element dont un coin contient la pointe de fissure ELSEIF ((PSIMIA.EQ.0.d0).AND.(PHIMIA.EQ.0.d0)) THEN C => choix: seuls les noeud en avant du front sont F-enrichi (=100) DO K=1,NBNN1 INUM = NUMND(K) if (XPSI(K).GE.0.d0) ITYPND(INUM) = 100 ENDDO C c Les autres elements ne sont pas directement coupés ni tangents avec la fissure ELSE C => on met a 1 pour dire qu on les a vu DO K=1,NBNN1 INUM = NUMND(K) ITYPND(INUM) = MAX(1,ITYPND(INUM)) ENDDO C ENDIF C++++++++ Fin des Tests en pagaille pour trier 222 CONTINUE C=======FIN DE la 2eme BOUCLE SUR LES ELEMENTS========== segsup,NUMND,XPSI,XPHI,KKGAU 100 CONTINUE C=====FIN DE LA PREMIERE BOUCLE SUR LES SOUS MODELES================ C A ce stade, pour la syntaxe 1, c on a cree le champ par point ITYPND qui vaut : c 1 -> noeud non enrichi c 10 -> noeud H enrichi c 100 -> noeud F1 ou F2 enrichi cbp, TODO : surement a revoir ... IF (ISYNTAX.EQ.2) GOTO 999 C======================================================================= c C 2EME BOUCLE SUR LES SOUS MODELES c c --> remplissage de MCHAM2.IELVAL(1/2/3) = MELVA2 c avec MELVA2.IELCHE(K,J) = listreel de phi ou psi associe c a l'enrichissement du Kieme noeud du Jieme element c c --> ITYPND =101/102 selon F1/F2-enrichi c c --> creation de MCHAM3.IELVAL(1/2/3) = MELVA3 c associee aux elements SURE c avec MELVA3.VELCHE(1,J) = 1. si H/F1/F2-enrichi , 0. sinon c c --> ITYPND = ITYPND + 200/2000/20000 si le hanging node appartient c a une relation ou l'un des noeuds est H/F1/F2-enrichi c C======================================================================= nexsur = 0 DO 400 ISMOD1 = 1,NBSMOD c write (*,*) 'boucle 400 : sous modele ', ISMOD1,'/',NBSMOD c CcccccC ouverture du sous modele IMODEL = MMODE1.KMODEL(ISMOD1) segact IMODEL*mod IPT3 = IMODEL.IMAMOD SEGACT IPT3 c travail specifique pour les elements SURE if (NEFMOD.EQ.259) GOTO 410 c on travaille seulement si EF XFEM if(NEFMOD.ne.263.and.NEFMOD.ne.264) goto 400 C======================================================= c 3EME BOUCLE SUR LES ELEMENTS C POUR ATTRIBUER L'ENRICHISSEMENT C======================================================= c write(*,*) 'MCHAM2=',MCHAM2 cbp : il n'y en a qu'1 --> deja en segini segact,MCHAM2 c on evite d'ouvrir les MELVA2 pour chaque element DO I2ENR=1,N2ENR MELVA2 = MCHAM2.IELVAL(I2ENR) if(MELVA2.ne.0) segact,MELVA2*mod ENDDO segact,MPOVA1,MPOVA2 c++++++ Debut de boucle sur les elements DO 300 J=1,NBEL1 c write(*,*) '***** ELEM', J,'/', NBEL1 c pour eviter de creer bcp de listreels on utilise l astuce : mphi =0 mpsiph=0 c++++++++ Boucle sur les noeuds de l element j DO 310 K=1,NBNN1 C Kieme noeud du Jieme element a pour numero INUM INUM = NUM(K,J) c write(*,*) 'INUM', INUM, 'ITYPND(INUM)', ITYPND(INUM) c on saute les noeuds non enrichis if((ITYPND(INUM)).le.1) goto 310 Ccccccccccc noeud H-enrichi if ((ITYPND(INUM)).eq.10) then MELVA2 = MCHAM2.IELVAL(1) c segact,MELVA2*mod MLREE1 = MELVA2.IELCHE(K,J) c on H-enrichit ssi ce noeud n etait pas deja H-enrichi if (MLREE1.eq.0) then if (mphi.ne.0) then MELVA2.IELCHE(K,J) = mphi else JG = NBNN1 segini,MLREE2 MELVA2.IELCHE(K,J) = MLREE2 c on remplit la listreel des valeurs de PHI aux noeuds de l element do inode=1,NBNN1 INUM2 = NUM(inode,J) I2 = NDPT2(INUM2) MLREE2.PROG(inode) = MPOVA2.VPOCHA(I2,1) enddo mphi=MLREE2 endif endif endif Ccccccccccc fin des noeud H-enrichi Ccccccccccc noeud F-enrichi * if((ITYPND(INUM)).eq.100) then *as 2010_01_20: lsaut =1000 si option 'SAUT' , =100 sinon c if (((ITYPND(INUM)).eq.lsaut).OR.((ITYPND(INUM)).eq.101)) then c write(*,*) 'ITYPND(INUM)', ITYPND(INUM), 'lsaut', lsaut c write(*,*) 'Nelem', j , 'INUM', inum if (((ITYPND(INUM)).ge.lsaut)) then c si F1 pas prevu on le rajoute if (N2ENR.lt.2) then N2ENR = 2 N2 = N2ENR segadj,MCHAM2 MCHAM2.NOMCHE(N2ENR) = 'F1' MCHAM2.TYPCHE(N2ENR) = 'POINTEURLISTREEL' N1PTEL=0 N1EL =0 N2PTEL=NBNN1 N2EL =NBEL1 segini,MELVA2 MCHAM2.IELVAL(N2ENR) = MELVA2 else MELVA2 = MCHAM2.IELVAL(2) c segact,MELVA2*mod endif c - on F1-enrichit si ce noeud n etait pas deja F1-enrichi c OU si on est dans une stratégie de désenrichissement c et que c'est le tour de F1 (on alterne F1, F2, F1, etc.) MLREE1 = MELVA2.IELCHE(K,J) c if (MLREE1.eq.0) then if ((MLREE1.eq.0.and.IMOT0.ne.2).OR & .(ITIP2.eq.1.and.IMOT0.eq.2)) then ITYPND(INUM)=101; if (mpsiph.ne.0) then MELVA2.IELCHE(K,J) = mpsiph else NBLV7=2 JG = NBLV7*NBNN1 segini,MLREE2 MELVA2.IELCHE(K,J) = MLREE2 c on remplit la listreel des valeurs de PSI et PHI do inode=1,NBNN1 INUM2 = NUM(inode,J) I1 = NDPT1(INUM2) I2 = NDPT2(INUM2) MLREE2.PROG(inode) = MPOVA1.VPOCHA(I1,1) MLREE2.PROG(NBNN1+inode) = MPOVA2.VPOCHA(I2,1) enddo mpsiph=MLREE2 endif c - si noeud deja F1-enrichi, alors on essaie de le F2-enrichir c OU si on est dans une stratégie de désenrichissement c et que c'est le tour de F2 else ITYPND(INUM)=102; c si F2 pas prevu on le rajoute if (N2ENR.lt.3) then N2ENR = 3 N2 = N2ENR segadj,MCHAM2 MCHAM2.NOMCHE(N2ENR) = 'F2' MCHAM2.TYPCHE(N2ENR) = 'POINTEURLISTREEL' N1PTEL=0 N1EL =0 N2PTEL=NBNN1 N2EL =NBEL1 segini,MELVA2 MCHAM2.IELVAL(N2ENR) = MELVA2 else MELVA2 = MCHAM2.IELVAL(3) segact,MELVA2*mod endif c +-- on F2-enrichit si ce noeud n etait pas deja F1-enrichi MLREE1 = MELVA2.IELCHE(K,J) if (MLREE1.eq.0) then if (mpsiph.ne.0) then MELVA2.IELCHE(K,J) = mpsiph else NBLV7=2 JG = NBLV7*NBNN1 segini,MLREE2 MELVA2.IELCHE(K,J) = MLREE2 c on remplit la listreel des valeurs de PSI et PHI do inode=1,NBNN1 INUM2 = NUM(inode,J) I1 = NDPT1(INUM2) I2 = NDPT2(INUM2) MLREE2.PROG(inode) = MPOVA1.VPOCHA(I1,1) MLREE2.PROG(NBNN1+inode) = MPOVA2.VPOCHA(I2,1) enddo mpsiph=MLREE2 endif c +-- si on est la, c est que le noeud est deja F1 et F2 enrichi c => pour l instant on genere une erreur, c mais pour le futur ajouter la possiblite d avoir N2ENR >3 else if(impi.ge.1) write(ioimp,*) & 'Trop d enrichissement (> 2 pointes) pour le ', & 'noeud',INUM,' = ',K,'ieme noeud du',J,'ieme element' c "Nombre maximum de fonctions de forme atteint au noeud %i1." INTERR(1)=INUM CALL ERREUR(1112) RETURN endif c +-- fin de la distinction noeud deja F2-enrichi ou pas endif c - fin de la distinction noeud deja F1-enrichi ou pas endif Ccccccccccc fin du cas noeud F-enrichi 310 CONTINUE c++++++++ Fin de Boucle sur les noeuds de l element j 300 CONTINUE C=======FIN DE la 3eme BOUCLE SUR LES ELEMENTS=========== C A ce stade, c ITYPND(#global) = 1 -> noeud non enrichi c 10 -> noeud H enrichi c 101 -> noeud F1 enrichi c 102 -> noeud F2 enrichi c C On a construit le champ par element MCHAM2 à 3 composantes : C 1. MCHAM2.IELVAL(1) : C champ par elements qui vaut phi sur les elements h-enrichis C 2. MCHAM2.IELVAL(2) C champ par elements contenant les valeurs de psi et phi C aux noeuds F1-enrichis C 3. MCHAM2.IELVAL(3) C champ par elements contenant les valeurs de psi et phi C aux noeuds F2-enrichis c on n' plus besoin des level-sets simplifiees SEGSUP,MPOVA3,MPOVA4 C if(NDPT1.ne.NDPT2) SEGSUP,NDPT2 C SEGSUP,NDPT1 C======================================================= c 4eme BOUCLE SUR LES ELEMENTS C construction d'un champ d'enrichissement associée aux SURE C======================================================= 410 CONTINUE IF(NEFMOD.NE.259 ) GOTO 400 IF(infele(13).NE.63) GOTO 400 C gg 2017_07_05 : traitement des elements sure c (pour les relations de compatibilite issues de RAFF) nexsur = 1 c recuperation du maillage support du sous-modele de sure NBNN3 = IPT3.NUM(/1) NBEL3 = IPT3.NUM(/2) MN3=IMODEL.INFMOD(/1) NFOR=IMODEL.FORMOD(/2) NMAT=IMODEL.MATMOD(/2) * on ajoute au sous-modele . IVAMOD(1) = MCHAM3 NOBMOD = 1 SEGADJ, IMODEL C write(*,*) 'Creation du champ d enrichissement de sure' C++++ INITIALISATION d'un nouveau MCHAML vierge MCHAM3 c L1 = 8 C nombre de sous champs c N1 = 1 c N3 = 6 c SEGINI, MCHEX3 c MCHEX3.TITCHE = 'ENRICHIS' c MCHEX3.IFOCHE = IFOUR c MCHEX3.IMACHE(1)= IMAMOD c MCHEX3.CONCHE(1)= IMODEL.CONMOD c MCHEX3.INFCHE(1,2) = 0 c MCHEX3.INFCHE(1,3) = NIFOUR c MCHEX3.INFCHE(1,6) = 1 C nombre de composantes N2 = 3 SEGINI, MCHAM3 c MCHEX3.ICHAML(1)=MCHAM3 MCHAM3.NOMCHE(1)= 'H ' MCHAM3.NOMCHE(2)= 'F1 ' MCHAM3.NOMCHE(3)= 'F2 ' N1PTEL=NBNN3 N1EL =NBEL3 N2PTEL=0 N2EL =0 SEGINI, MELVA3 DO 409 NCOMP2=1,N2 MCHAM3.TYPCHE(NCOMP2)= 'REAL*8' SEGINI, MELVA3 MCHAM3.IELVAL(NCOMP2)=MELVA3 409 CONTINUE c on branche dans le sous-modele le nouveau MCHAM3 cbp, TODO : ligne ci-dessous me semble etre une bidouille --> a verifier TYMODE(1)='MCHAMLL ' IVAMOD(1)=MCHAM3 c ajout du nouveau sous-champ d'enrichissement de sure MCHAM3 c dans le champ d'enrichissement de sortie MCHEX2 SEGACT MCHEX2 L1 = 8 N3 = 6 N1 = MCHEX2.IMACHE(/1) + 1 SEGADJ,MCHEX2 SEGACT MCHEX2*mod MCHEX2.IMACHE(N1)= IPT3 MCHEX2.CONCHE(N1)= IMODEL.CONMOD MCHEX2.ICHAML(N1)= MCHAM3 MCHEX2.INFCHE(N1,2) = 0 MCHEX2.INFCHE(N1,3) = NIFOUR MCHEX2.INFCHE(N1,6) = 1 c++++++ Debut de boucle sur les elements DO 411 J=1,NBEL3 c compteurs d'enrichissements dans l'element Nenrh = 0 Nenrf1 = 0 Nenrf2 = 0 c++++++++ Boucle sur les noeuds de l element j DO 412 K=1,NBNN3 C Kieme noeud du Jieme element a pour numero INUM INUM = IPT3.NUM(K,J) c on saute les noeuds non enrichis if((ITYPND(INUM)).le.1) goto 412 Ccccccccccc noeud H-enrichi if ((ITYPND(INUM)).eq.10) then MELVA3 = MCHAM3.IELVAL(1) segact,MELVA3*mod MELVA3.VELCHE(K,J) = 1.0 Nenrh = Nenrh +1 Ccccccccccc noeud F1-enrichi elseif ((ITYPND(INUM)).eq.101) then MELVA3 = MCHAM3.IELVAL(2) segact,MELVA3*mod MELVA3.VELCHE(K,J) = 1.0 Nenrf1 = Nenrf1 + 1 Ccccccccccc noeud F2-enrichi elseif ((ITYPND(INUM)).eq.102) then MELVA3 = MCHAM3.IELVAL(3) segact,MELVA3*mod MELVA3.VELCHE(K,J) = 1.0 Nenrf2 = Nenrf2 + 1 endif 412 CONTINUE Ccccccccccc reperage des hanging nodes a enrichir if (ITYPND(IPT3.NUM(1,J)).le.1) then if (nenrh.gt.0) then ITYPND(IPT3.NUM(1,J))= ITYPND(IPT3.NUM(1,J)) + 200 endif endif if (ITYPND(IPT3.NUM(1,J)).lt.11) then if (nenrf1.gt.0) then ITYPND(IPT3.NUM(1,J))= ITYPND(IPT3.NUM(1,J)) + 2000 endif if (nenrf2.gt.0) then ITYPND(IPT3.NUM(1,J))= ITYPND(IPT3.NUM(1,J)) + 20000 endif endif 411 CONTINUE c++++++ fin de boucle sur les elements 400 CONTINUE C=====FIN DE LA DEUXIEME BOUCLE SUR LES SOUS MODELES================ C A ce stade, c ITYPND(#global) = c 1 -> noeud non enrichi c 10 -> noeud H enrichi c 101 -> noeud F1 enrichi c 102 -> noeud F2 enrichi c 201 -> noeud non enrichi à H enrichir c 301 -> noeud F1 enrichi à H enrichir c 302 -> noeud F2 enrichi à H enrichir c 2001 -> noeud non enrichi à F1 enrichir c 2010 -> noeud H enrichi à F1 enrichier c 2102 -> noeud F2 enrichi à F1 enrichir c 2201 -> noeud non enrichi à F1 enrichir et H enrichir c 2302 -> noeud F2 enrichi à F1 enrichir et H enrichir c 20001 -> noeud non enrichi à F2 enrichir c 20010 -> noeud H enrichi à F2 enrichir c 20101 -> noeud F1 enrichi à F2 enrichir c 20201 -> noeud non enrichi à F2 enrichir et H enrichir c 20301 -> noeud F1 enrichi à F2 enrichir et H enrichir C 22001 -> noeud non enrichi à F2 enrichir et F1 enrichir c 22010 -> noeud H enrichi à F2 enrichir et F1 enrichir c 22201 -> noeud non enrichi à F2 enrichir, F1 enrichir et H enrichir c C On a construit le champ par element à 3 composantes C Portant sur le maillage de XC8R ou XQ4R C MCHAM2.IELVAL(1) C-> champ par elements qui vaut phi sur les elements h-enrichis C MCHAM2.IELVAL(2) C-> champ par elements contenant les valeurs de psi et phi C aux noeds F1-enrichis C MCHAM2.IELVAL(3) C-> champ par elements contenant les valeurs de psi et phi C aux noeds F2-enrichis C On a construit le champ par element à 3 composantes C Portant sur le maillage de SURE C MCHAM3.IELVAL(1) C-> champ par elements qui vaut 1 sur les noeud h-enrichis C MCHAM3.IELVAL(2) C--> champ par elements qui vaut 1 sur les noeud F1-enrichis C MCHAM3.IELVAL(3) C--> champ par elements qui vaut 1 sur les noeud F2-enrichis if (nexsur.eq.0) goto 999 C======================================================================= c C 3EME BOUCLE SUR LES SOUS MODELES c C======================================================================= DO 500 ISMOD1 = 1,NBSMOD c CcccccC ouverture du sous modele IMODEL = MMODE1.KMODEL(ISMOD1) segact IMODEL*mod IPT4 = IMODEL.IMAMOD SEGACT IPT4 NBNN4 = IPT4.NUM(/1) NBEL4 = IPT4.NUM(/2) c on travaille seulement si EF XFEM if(NEFMOD.ne.263.and.NEFMOD.ne.264) goto 500 C======================================================= c 5EME BOUCLE SUR LES ELEMENTS C pour attribuer un enrichissement aux hanging nodes non enrichis c dont les sure vont etre enrichis C======================================================= segact,MCHAM2 c on evite d'ouvrir les MELVA2 pour chaque element c DO I2ENR=1,N2ENR c MELVA2 = MCHAM2.IELVAL(I2ENR) c if(MELVA2.ne.0) segact,MELVA2*mod c ENDDO segact,MPOVA1,MPOVA2 c++++++ Debut de boucle sur les elements DO 600 J=1,NBEL4 c pour eviter de creer bcp de listreels on utilise l astuce : mphi =0 mpsiph=0 c++++++++ Boucle sur les noeuds de l element j DO 610 K=1,NBNN4 C Kieme noeud du Jieme element a pour numero INUM INUM = IPT4.NUM(K,J) if (INUM.EQ.132) then endif c on saute les noeuds qui ne nous interessent pas if((ITYPND(INUM)).le.200) goto 610 Ccccccccccc noeud H-enrichi if ((ITYPND(INUM)).eq.201) GOTO 611 if ((ITYPND(INUM)).eq.301) GOTO 611 if ((ITYPND(INUM)).eq.302) GOTO 611 if ((ITYPND(INUM)).eq.2201) GOTO 611 if ((ITYPND(INUM)).eq.2302) GOTO 611 if ((ITYPND(INUM)).eq.20201) GOTO 611 if ((ITYPND(INUM)).eq.20301) GOTO 611 if ((ITYPND(INUM)).eq.22201) GOTO 611 GOTO 612 611 CONTINUE MELVA2 = MCHAM2.IELVAL(1) segact,MELVA2*mod MLREE1 = MELVA2.IELCHE(K,J) c on H-enrichit ssi ce noeud n etait pas deja H-enrichi if (MLREE1.eq.0) then if (mphi.ne.0) then MELVA2.IELCHE(K,J) = mphi else JG = NBNN1 segini,MLREE2 MELVA2.IELCHE(K,J) = MLREE2 c on remplit la listreel des valeurs de PHI aux noeuds de l element do inode=1,NBNN1 INUM2 = IPT4.NUM(inode,J) I2 = NDPT2(INUM2) MLREE2.PROG(inode) = MPOVA2.VPOCHA(I2,1) enddo mphi=MLREE2 endif endif Ccccccccccc fin des noeud H-enrichi 612 CONTINUE Ccccccccccc noeud F1-enrichi if ((ITYPND(INUM)).eq.2001) GOTO 613 if ((ITYPND(INUM)).eq.2010) GOTO 613 if ((ITYPND(INUM)).eq.2102) GOTO 613 if ((ITYPND(INUM)).eq.2201) GOTO 613 if ((ITYPND(INUM)).eq.2302) GOTO 613 if ((ITYPND(INUM)).eq.22001) GOTO 613 if ((ITYPND(INUM)).eq.22010) GOTO 613 if ((ITYPND(INUM)).eq.22201) GOTO 613 goto 614 613 CONTINUE MELVA2 = MCHAM2.IELVAL(2) segact,MELVA2*mod if (mpsiph.ne.0) then MELVA2.IELCHE(K,J) = mpsiph else NBLV7=2 JG = NBLV7*NBNN1 segini,MLREE2 MELVA2.IELCHE(K,J) = MLREE2 c on remplit la listreel des valeurs de PSI et PHI do inode=1,NBNN4 INUM2 = IPT4.NUM(inode,J) I1 = NDPT1(INUM2) I2 = NDPT2(INUM2) MLREE2.PROG(inode) = MPOVA1.VPOCHA(I1,1) MLREE2.PROG(NBNN1+inode) = MPOVA2.VPOCHA(I2,1) enddo mpsiph=MLREE2 endif Ccccccccccc fin des noeud F1-enrichi 614 CONTINUE Ccccccccccc noeud F2-enrichi if ((ITYPND(INUM)).eq.20001) GOTO 615 if ((ITYPND(INUM)).eq.20010) GOTO 615 if ((ITYPND(INUM)).eq.20101) GOTO 615 if ((ITYPND(INUM)).eq.20201) GOTO 615 if ((ITYPND(INUM)).eq.20301) GOTO 615 if ((ITYPND(INUM)).eq.22001) GOTO 615 if ((ITYPND(INUM)).eq.22010) GOTO 615 if ((ITYPND(INUM)).eq.22201) GOTO 615 goto 610 615 CONTINUE MELVA2 = MCHAM2.IELVAL(3) segact,MELVA2*mod if (mpsiph.ne.0) then MELVA2.IELCHE(K,J) = mpsiph else NBLV7=2 JG = NBLV7*NBNN1 segini,MLREE2 MELVA2.IELCHE(K,J) = MLREE2 c on remplit la listreel des valeurs de PSI et PHI do inode=1,NBNN1 INUM2 = IPT4.NUM(inode,J) I1 = NDPT1(INUM2) I2 = NDPT2(INUM2) MLREE2.PROG(inode) = MPOVA1.VPOCHA(I1,1) MLREE2.PROG(NBNN1+inode) = MPOVA2.VPOCHA(I2,1) enddo mpsiph=MLREE2 endif Ccccccccccc fin du cas noeud F2-enrichi 610 CONTINUE 600 CONTINUE c++++++++ Fin de Boucle sur les noeuds de l element j 500 CONTINUE c fin de la boucle sur les sous modeles 999 CONTINUE C======================================================================= c MENAGE C======================================================================= if(NDPT1.ne.0) then if(NDPT1.ne.NDPT2.and.NDPT2.ne.0) SEGSUP,NDPT2 SEGSUP,NDPT1 endif SEGSUP,ITYPND C======================================================================= c ECRITURE DES RESULTATS DE L'OPERATEUR C======================================================================= IF(ISYNTAX.EQ.1) THEN if(IPRIG2.ne.0) CALL ECROBJ('RIGIDITE',IPRIG2) CALL ACTOBJ('MCHAML ',MCHEX2,1) CALL ECROBJ('MCHAML ',MCHEX2) ENDIF c rem : pour la SYNTAXE 2, c'est une directive --> rien a faire END