C RELXFE SOURCE PV090527 26/04/30 21:16:12 12529 C RELXFE SOURCE BP208322 17/04/18 21:15:13 9395 SUBROUTINE RELXFE(IMODEL, MRIGID) C*********************************************************************** C cet operateur créé d'une matrice élémentaire de rigidité c pour imposer à zéro les ddl xfem crée mais non actif c (dans les éléments de transition) C C ENTREES : C ________ C C IMODEL pointeur sur le modele élémentaire C C ENTREES/SORTIE : C ________ C c MRIGID rigidité chapeu dans laquelle on va écrire c (dans la dernière sous matrice) la rigidité voulue C--------------------------------------------------------------------- C*********************************************************************** IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*8 CMATE PARAMETER (NDDLMAX=30,NBNIMAX=10) INTEGER LOCIRI(10,(1+NDDLMAX)) DIMENSION MLRE(NDDLMAX+1) CHARACTER*4 MOTINC(NDDLMAX),MOTDUA(NDDLMAX) DATA MOTINC/'UX ','UY ','UZ ','AX ','AY ','AZ ', >'B1X ','B1Y ','B1Z ','C1X ','C1Y ','C1Z ','D1X ','D1Y ', >'D1Z ','E1X ','E1Y ','E1Z ','B2X ','B2Y ','B2Z ','C2X ', >'C2Y ','C2Z ','D2X ','D2Y ','D2Z ','E2X ','E2Y ','E2Z '/ DATA MOTDUA/'FX ','FY ', 'FZ ','FAX ','FAY ','FAZ ', >'FB1X','FB1Y','FB1Z','FC1X','FC1Y','FC1Z','FD1X','FD1Y', >'FD1Z','FE1X','FE1Y','FE1Z','FB2X','FB2Y','FB2Z','FC2X', >'FC2Y','FC2Z','FD2X','FD2Y','FD2Z','FE2X','FE2Y','FE2Z'/ -INC PPARAM -INC CCOPTIO -INC CCGEOME -INC SMELEME -INC SMCOORD -INC SMRIGID -INC SMMODEL -INC CCHAMP -INC SMCHAML -INC SMLREEL -INC SMINTE POINTEUR MCHEX1.MCHELM c sement raccourcis par éléments SEGMENT MRACC INTEGER TLREEL(NBENRMA2,NBI) INTEGER MELRIG(NBELEM) ENDSEGMENT c Segment contenant l'info DDL a mettre à 0 SEGMENT TBLOQ INTEGER MBLOQ(3,NBPTB) INTEGER NBLOQ(3) ENDSEGMENT C++++ Recup + Activation de la geometrie ++++++++++++++++ IPT1= IMODEL.IMAMOD SEGACT, IPT1 C nbre de noeuds par element NBNN1 = IPT1.NUM(/1) C nbre d elements NBEL1 = IPT1.NUM(/2) C++++ RECUP DES INFOS EF ++++++++++++++++++++++++++++++++ MELE = IMODEL.NEFMOD MFR = INFELE(13) IELE = INFELE(14) NDDL = INFELE(15) NSTRS = INFELE(16) C+++++nombre de sous matrice de mrigid (va être ammené a changé) NRIGEL = MRIGID.IRIGEL(/2) C++++ Recup des infos d enrichissement +++++++++++++++++++ c recup du MCHEX1 d'enrichissement NBENR1=0 MCHAM1=0 NOBMO1=IMODEL.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(ioimp,*) 'Le modele est vide (absence d enrichissement)' * return else write(ioimp,*) 'Aucun MCHAML 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 if (MCHAM1.ne.0) NBENR1=MCHAM1.IELVAL(/1) C++++ INITIALISATIONS... C c ... des tables LOCIRI et MLRE c MLRE contient le nombre d'inconnues de chaque sous-zone c determinee depuis le nombre de fonctions 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) if (NBENR1.ne.0) then do ienr=1,NBENR1 nbniJ = 2 + ((ienr-1)*4) MLRE(1+ienr) = IDIM*NBNN1*nbniJ c write(*,*) 'ienr', ienr, 'mlre', MLRE(1+ienr) c -LOCIRI LOCIRI(5,1+ienr)= NIFOUR enddo endif C Tables + longues car 1er indice correspond au fontion de forme std MLRE(1) = IDIM*NBNN1*1 LOCIRI(5,1)= NIFOUR c on complete avec des 0 if (NBENR1.lt.(NDDLMAX+1)) then do ienr=(NBENR1+1),(NDDLMAX) MLRE(1+ienr) = 0 enddo endif c c ...DU SEGMENT MRACC NBENRMA2 = NDDLMAX NBI = NBNN1 NBELEM = NBEL1 segini , MRACC C initialisation du tableau des ddl a mettre à zéro SEGACT,MCOORD*MOD NBPTB= nbpts SEGINI,TBLOQ C++++ TBLOQ.MBLOQ(NBENRJ, INUM) = faut il mettre a 0 les ddl ++++ C =0 si pas encore passé dans les ddl C =1 si déja passé dans les ddl pas de mise à 0 c =2 si déja passé dans les ddl mise à 0 nécéssaire C ++++ TBLOQ.NBLOQ(NBENRJ) Compteur de ddl de type NBENRJ à mettre à 0 C********************************************************* C C>>>>>>>>>>>>>>>>>>>>>>>>>>> 1ere BOUCLE SUR LES ELEMENTS C NBENR = 0 DO 2000 J=1,NBEL1 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) NBENRJ=max(NBENRJ,IENR) enddo enddo endif NBENR=max(NBENRJ,NBENR) NDDLE = MLRE(NBENRJ+1) c if (NBENRJ.ne.0) then c write(*,*) '***********************************************' c write(*,*) 'ELEMENT', J c write(*,* ) 'Niveau d enrichssement', NBENRJ c write(*,* ) 'Nb de ddl', NDDLE c endif C********************************************************* C C>>>>>>>>>>>>>>>>>>>>>>>>>>> BOUCLE SUR LES DDL De L' ELEMENT C DO 3000 II=1,NDDLE C********************************************************************** C On cherche les noeuds qui ne sont pas effectivement enrichis C pour forcer à 0 les DDL correspondants C********************************************************************** c on trouve le type de fonction de ce ddl: IENR=0 si U, =1 si A, c =2 si B1, =3 si C1, = 4 si D1, =5 si E1, =6 si B2, ... =9 si E2 IENR = ((II-1)/IDIM) / NBNN1 c on trouve le niveau d enrichissement KENR de ce ddl si nonstd if (IENR.eq.0) then go to 3001 elseif(IENR.ge.2) then C write(*,*) 'DDL' , II, 'Niveu d enr', IENR KENR = ((IENR - 2) / 4) + 2 c ci dessus: 4 represente le nombre de fonction de la base d'enrichissement c et 2 est le decalage du a U et H else KENR = IENR endif c on trouve le noeud correspondant au ddl CYT INODE = ((II+1)/IDIM) - ((IENR)*NBNN1) INODE = 1 + ((II-1)/IDIM) - ((IENR)*NBNN1) c numero global du noeud INUM = IPT1.NUM(INODE , J) c est ont déja passé dans ce ddl ? c write(*,*) 'INUM', INUM, 'Kenr', KENR c write(*,*) 'mlree1', mlree1,'Mbloq',Tbloq.Mbloq(KENR, INUM) if (Tbloq.Mbloq(KENR, INUM).gt.0) GOTO 3001 c est-ce un noeud vraiment enrichi? mlree1 = TLREEL(KENR,INODE) Tbloq.Mbloq(KENR, INUM)=Tbloq.Mbloq(KENR, INUM)+1 if (mlree1.eq.0) then Tbloq.Mbloq(KENR, INUM)=Tbloq.Mbloq(KENR, INUM)+1 Tbloq.Nbloq(KENR) = Tbloq.Nbloq(KENR)+1 endif c write(*,*) 'Nouveau Mbloq', Tbloq.Mbloq(KENR, INUM) 3001 CONTINUE 3000 CONTINUE c c 2000 CONTINUE C FIN DE BOUCLE SUR LES ELEMENTS C********************************************************* C Creation des matrices de bloquage C pour les DDL non effectivement enrichis C********************************************************* NLIGRD = 2 NLIGRP = 2 NBNN = 2 NBPTS1=NBPTB C********************************************************* C C>>>>>>>>>>>>>>>>>>>>>>>>>>> Boucle sur les types d'enrichissement C C********************************************************* IDDL = 3 DO 4000 IENR =1, NBENR1 C Maillage des noeuds à bloquer C nombre de noeuds a bloquer -> nombre d'éléments du maillage de blocage if (IENR.eq.1) then NBELEM = TBLOQ.NBLOQ(1) elseif (IENR.EQ.2) then NBELEM = TBLOQ.NBLOQ(2) else NBELEM = TBLOQ.NBLOQ(3) endif NFON=1 IF (IENR.gt.1) NFON=4 C********************************************************* C C>>>>>>>>>>>>>>>>>>>>>>>>>>> Boucle sur les fonctions de formes d'enrichissement C C********************************************************* DO 4001 IFON = 1 , NFON C********************************************************* C C>>>>>>>>>>>>>>>>>>>>>>>>>> Boucle sur les composantes C C********************************************************* DO 4002 ICOMP = 1, IDIM IDDL=IDDL+1 C si aucun noeud a bloquer on saute ce type d'enrichissement if (NBELEM.EQ.0) goto 4000 C********************************************************* C Maillage du blocage NBREF = 0 NBSOUS = 0 SEGINI IPT2 IPT2.ITYPEL = 22 Nelem = 0 C ajustement de XMCOORD pour ajouter les noeud des multiplicateurs NBPTS = NBPTS1 + NBELEM SEGADJ MCOORD C>>>>>>>>>>>>>>>>>>>>>>>>>> Boucle sur les noeuds DO 4010 INUM = 1, NBPTB IF (TBLOQ.MBLOQ(IENR, INUM).LT.2) goto 4010 Nelem = Nelem + 1 IPT2.NUM(2,Nelem)=INUM NBPTS1 = NBPTS1 + 1 IPT2.NUM(1,Nelem)=NBPTS1 IPT2.icolor(Nelem)=1 4010 CONTINUE C Fin de boucle sur les noeuds C coordonées des nouveaux noeud C>>>>>>>>>>>>>>>>>>>>>>>>>> Boucle sur les éléments de ipt2 DO 4011 IELE = 1, NBELEM INUM2 = IPT2.NUM(2,IELE) NBPTS2 = NBPTS1 - NBELEM + IELE INEW = (IDIM+1)*(NBPTS2-1) IOLD = (IDIM+1)*(INUM2-1) DO ID = 1, IDIM XCOOR(INEW +ID ) = XCOOR(IOLD +ID ) ENDDO 4011 CONTINUE C********************************************************* C Matrice de blocage NELRIG = NELEM rigrel=0 SEGINI XMATR1 DO i=1,NELRIG XMATR1.RE(1,1,i)=0.D0 XMATR1.RE(2,1,i)=1.D0 XMATR1.RE(2,2,i)=0.D0 XMATR1.RE(1,2,i)=1.D0 ENDDO C********************************************************* C Segment Descripteur SEGINI DES1 DES1.LISINC(1)='LX' DES1.LISDUA(1)='FLX' DES1.LISINC(2)=MOTINC(IDDL) DES1.LISDUA(2)=MOTDUA(IDDL) DES1. NOELEP(1)=1 DES1. NOELEP(2)=2 DES1. NOELED(1)=1 DES1. NOELED(2)=2 C********************************************************* C stockage de la rigidité construite dans MRIGID NRIGEL = NRIGEL+1 SEGADJ MRIGID SEGACT, MRIGID*MOD MRIGID.IRIGEL(1, NRIGEL)= IPT2 MRIGID.IRIGEL(3, NRIGEL)=DES1 MRIGID.IRIGEL(4, NRIGEL)=XMATR1 MRIGID.COERIG(NRIGEL)=1.D0 c desactivations SEGDES XMATR1, DES1 c write(*,*) 'SOUS MATRICE', NRIGEL, MOTINC(IDDL) c WRITE(*,*) 'NB de blocages ', NBELEM c WRITE(*,*) 'maillage ', IPT2 4002 CONTINUE IF (IDIM.EQ.2) IDDL=IDDL+1 4001 CONTINUE 4000 CONTINUE C write (*,*) '***************NRIGEL*******************', NRIGEL C Fin des boucles sur les niveau d'enrichissement et composantes c C********************************************************* C SUPPRESSION ET DESACTIVATION DE SEGMENTS C********************************************************* SEGSUP,MRACC,TBLOQ segdes, MCOORD END