C NUMOP3 SOURCE GOUNAND 25/05/15 21:15:09 12268 SUBROUTINE NUMOP3(MELEME,ICPR,NODES) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER (I-N) C*********************************************************************** C NOM : NUMOP3 C DESCRIPTION : Numerotation des noeuds d'un maillage C par la methode de SLOAN C Reprise de NUMOP2 pour le graphe d'adjacence C et le placement des multiplicateurs C C C LANGAGE : ESOPE C AUTEUR : Stephane GOUNAND (CEA/DES/ISAS/DM2S/SEMT/LTA) C mel : gounand@semt2.smts.cea.fr C*********************************************************************** C APPELES : C APPELES (E/S) : C APPELES (BLAS) : C APPELES (CALCUL) : C APPELE PAR : C*********************************************************************** C SYNTAXE GIBIANE : C ENTREES : C ENTREES/SORTIES : C SORTIES : C*********************************************************************** C VERSION : v1, 06/05/2025, version initiale C HISTORIQUE : v1, 06/05/2025, creation C HISTORIQUE : C HISTORIQUE : C*********************************************************************** -INC PPARAM -INC CCOPTIO -INC CCREEL -INC SMELEME -INC SMCOORD LOGICAL CONNEC SEGMENT JMEM(NODES+1),JMEMN(NODES+1) C JMEM et JMEMN contiennent le nombre d'element auquel appartient un noeud SEGMENT JNT(NODES) C JNT contient la nouvelle numerotation SEGMENT IJNT(NODES) C IJNT contient l'inverse de la nouvelle numerotation SEGMENT JORDRE(NODES) C Ordre des noeuds SEGMENT IWORK(3*NODES+1) C Tableau de travail pour prsloa SEGMENT IJNT2(NODES) SEGMENT JOR2(NODES) C Tableau de travail pour trifus SEGMENT ICPR(nbpts) C ICPR au debut contient l'ancienne numerotation , C a la fin la nouvelle. SEGMENT IADJ(NODES+1) SEGMENT JADJC(0) C IADJ(i) pointe sur JADJC qui contient les voisins de i entre C IADJ(i) et IADJ(i+1)-1 SEGMENT LAGRAN(NB) C contient les noeud de lagrange et les noeuds les suivant directement C cf element de type 49 SEGMENT BOOLEEN LOGICAL BOOL(NODES) ENDSEGMENT C BOOL(i) = true si le noeud i a ete deja mentionne dans la liste C des voisins JADJC. SEGMENT IMEMOIR(NBV),LMEMOIR(NBV) C contient les elements appartenant a chaque noeud,sous forme de liste. INTEGER ELEM C nom d'un element LOGICAL LDBNUM,LVERIF * * Executable statements * LDBNUM=.FALSE. LVERIF=.FALSE. SEGACT ICPR*MOD NODES=ICPR(/1) SEGACT MELEME C icpr: numero des noeuds. C meleme: objet de maillage (cf assem2.eso) DO 10 I=1,ICPR(/1) ICPR(I)=0 10 CONTINUE IPT1=MELEME IKOU=0 NBV=0 NB1=0 NB2=0 DO 100 IO=1,MAX(1,LISOUS(/1)) IF (LISOUS(/1).GT.0) THEN IPT1=LISOUS(IO) SEGACT IPT1 ENDIF C on cree la numerotation des noeuds. C 'nb noeuds/element'=IPT1.NUM(/1) C 'nb element'=IPT1.NUM(/2) IF(abs(IPT1.ITYPEL).EQ.49) THEN NB1=NB1+IPT1.NUM(/2) NB2=MAX(NB2,IPT1.NUM(/1)) C NB1= nbre d'éléments de type 49. C NB2=nbre de noeuds/élément maximum parmi C les éléments de type 22. ENDIF DO J=1,IPT1.NUM(/2) DO I=1,IPT1.NUM(/1) IJ=IPT1.NUM(I,J) C IJ est le Ième noeud du Jème élément. IF (ICPR(IJ).EQ.0) THEN C s'il est déjà numéroté, on ne fait rien. IKOU=IKOU+1 ICPR(IJ)=IKOU ENDIF enddo enddo 100 CONTINUE NODES=IKOU NB=NB2*NB1 C***** initalisation des segments********* SEGINI,LAGRAN SEGINI IADJ,JADJC SEGINI,JMEM,JMEMN SEGINI BOOLEEN C****************************************** IPT1=MELEME IADJ(1)=1 INC=0 DO 200 IO=1,MAX(1,LISOUS(/1)) IF (LISOUS(/1).GT.0) IPT1=LISOUS(IO) DO 210 J=1,IPT1.NUM(/2) IF(abs(IPT1.ITYPEL).EQ.49) THEN is=sign(1,ipt1.itypel) DO 220 I=1,IPT1.NUM(/1) C chaque element de type 49 a au plus NB2 noeuds. LAGRAN(INC*NB2+I)=ICPR(IPT1.NUM(I,J))*is C les noeuds de l'elements de type 49 C sont ranges dans le vecteur LAGRAN. 220 CONTINUE DO 225 I=IPT1.NUM(/1)+1,NB2 LAGRAN(INC*NB2+I)=0 C comme on a alloue la place memoire maximale, C on remplit les cases restantes avec des 0. 225 CONTINUE INC=INC+1 C INC=le nbre d'elements de type 49. ENDIF DO 230 I=1,IPT1.NUM(/1) IJ=ICPR(IPT1.NUM(I,J))+1 JMEM(IJ)=JMEM(IJ)+1 C JMEM(I+1): nb elements auquel le noeud I appartient 230 CONTINUE 210 CONTINUE 200 CONTINUE JMEM(1)=1 DO 30 I=1,NODES JMEM(I+1)=JMEM(I)+JMEM(I+1) C JMEM(I+1)=indice de depart des elements C auxquels le noeud I appartient. 30 CONTINUE NBV=JMEM(NODES+1) C NBV= dimension de IMEMOIR. SEGINI IMEMOIR,LMEMOIR IPT1=MELEME DO 300 IO=1,MAX(1,LISOUS(/1)) IF (LISOUS(/1).GT.0) THEN IPT1=LISOUS(IO) ENDIF DO J=1,IPT1.NUM(/2) DO I=1,IPT1.NUM(/1) IJ=ICPR(IPT1.NUM(I,J)) JMEMN(IJ+1)=JMEMN(IJ+1)+1 IMEMOIR(JMEM(IJ)+JMEMN(IJ+1)-1)=J LMEMOIR(JMEM(IJ)+JMEMN(IJ+1)-1)=IO C on range dans IMEMOIR tous les elements des sous-objets C IO auxquels appartient le noeud ICPR(IPT1.NUM(I,J)). C On connait pour chaque element, le sous-objet auquel C il appartient grace a LMEMOIR enddo enddo 300 CONTINUE DO 410 J=1,NODES BOOL(J)=.FALSE. 410 CONTINUE DO 400 I=1,NODES IADJ(I+1)=IADJ(I) DO 420 J=JMEM(I),JMEM(I+1)-1 ELEM=IMEMOIR(J) C ELEM=element auquel appartient le noeud I. IPT1=MELEME IF (LISOUS(/1).GT.0) IPT1=LISOUS(LMEMOIR(J)) itype = abs(ipt1.itypel) * si element de type 49 ou 22, on ne connecte pas 2 noeuds non LX connec=.true. if (itype.eq.49) then do k=3,ipt1.num(/1) if (i.eq.icpr(ipt1.num(k,elem))) connec=.false. enddo endif DO 430 K=1,IPT1.NUM(/1) C k representatif du nb de noeuds par elements. IK=ICPR(IPT1.NUM(K,ELEM)) if (k.gt.2.and..not.connec) goto 430 IF ((I.NE.IK).AND. & (.NOT.(BOOL(IK)))) THEN C si i n'est pas egal a un des nouveaux numeros des noeuds C de l'element ELEM et si il n'appartient pas deja a l'ens des C voisins du noeud i(jadjc(i)),alors on le rajoute. C JADJC(IADJ(I+1))=IK JADJC(**)=IK IADJ(I+1)=IADJ(I+1)+1 BOOL(IK)=.TRUE. ENDIF 430 CONTINUE 420 CONTINUE ** pour ne pas avoir un tableau de longueur nulle, ce que esope n'aime pas * if(jadjc(/1).eq.0) jadjc(**)=0 * remise a faux de bool DO 412 J=IADJ(I),IADJ(I+1)-1 IK=JADJC(J) BOOL(IK)=.FALSE. 412 CONTINUE 400 CONTINUE SEGSUP JMEM,JMEMN SEGSUP IMEMOIR,LMEMOIR SEGSUP BOOLEEN * * A ce stade, le graphe d'adjacence est près * if (ldbnum) Write(ioimp,*) 'Graphe adjacence nodes=',nodes if (ldbnum) segprt,IADJ if (ldbnum) segprt,JADJC SEGINI JNT SEGINI IWORK E2=JADJC(/1) CALL LABEL(NODES,E2,JADJC(1),IADJ(1), $ JNT(1), $ IWORK(1), $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 SEGSUP IWORK SEGSUP IADJ,JADJC if (ldbnum) SEGPRT,JNT * * Maintenant, il faut mettre les Lagrange a la bonne place !!!! * On fait comme dans KRES24 * IF (NB1.EQ.0) GOTO 2000 SEGINI IJNT DO I=1,NODES IJNT(JNT(I))=I ENDDO * Il y a trois types différents de -1 à 1 : il faut pouvoir placer * les multiplicateurs de lagrange entre les autres noeuds NTYP=3 SEGINI JORDRE DO I=1,NODES JORDRE(I)=I*NTYP ENDDO JORMAX= (NODES+1)*NTYP if (ldbnum) write(ioimp,*) 'Avant mise a la bonne place' if (ldbnum) segprt,jordre do 700 J=0,NB1-1 ipaur=-igrand ipaus=igrand do 800 il=3,nb2 ip = abs(LAGRAN(J*NB2+il)) * if (ip.eq.0) write (6,*) ' prob numop3 ' if (ip.eq.0) goto 710 jddln=JNT(IP) jordre(jddln)=-abs(jordre(jddln)) ipaur=max(ipaur,jordre(jddln)) ipaus=min(ipaus,jordre(jddln)) 800 continue 710 continue * * On va laisser comme ça * Ce cas peut arriver après élimination * Cela devrait revenir à placer les multiplicateurs en fin ou début * de matrice if (ipaur.eq.-igrand.or.ipaus.eq.igrand) then * Write(ioimp,*) 'mulag sans relations pas ok' * goto 9999 ipaur=0 ipaus=-jormax endif if (ldbnum) write(ioimp,*) 'j,ipaur,ipaus=',j,ipaur $ ,ipaus * * le premier mult avant le premier noeud ip=abs(LAGRAN(J*NB2+1)) iddln=JNT(IP) jordre(iddln)=ipaur+1 * * le deuxieme mult apres le dernier noeud ip=abs(LAGRAN(J*NB2+2)) iddln=JNT(IP) jordre(iddln)=ipaus-1 * 700 continue * WRITE(IOIMP,*) 'Avant chgt signe' * segprt,jordre DO I=1,nodes JORDRE(I)=-JORDRE(I) ENDDO * Avant tri if (ldbnum) WRITE(IOIMP,*) 'Avant TRIFUS' if (ldbnum) SEGPRT,IJNT if (ldbnum) SEGPRT,JORDRE SEGINI JOR2 SEGINI IJNT2 * ok maintenant on trie CALL TRIFUS(nodes,JORDRE(1),IJNT(1),JOR2(1),IJNT2(1)) SEGSUP IJNT2 SEGSUP JOR2 * Apres tri if (ldbnum) WRITE(IOIMP,*) 'Apres TRIFUS' if (ldbnum) SEGPRT,IJNT if (ldbnum) SEGPRT,JORDRE * permutation inverse DO I=1,nodes JNT(IJNT(I))=I ENDDO * Verification que dans la nouvelle numerotation les multiplicateurs * sont correctement places... IF (LVERIF) THEN write(ioimp,*) 'VERIF NUMOP3' do 1700 J=0,NB1-1 ipaur=-igrand ipaus=igrand do 1800 il=3,nb2 ip = abs(LAGRAN(J*NB2+il)) * if (ip.eq.0) write (6,*) ' prob numop3 ' if (ip.eq.0) goto 1710 jddln=JNT(IP) ipaur=max(ipaur,jddln) ipaus=min(ipaus,jddln) 1800 continue 1710 continue if (ipaur.eq.-igrand.or.ipaus.eq.igrand) then goto 1700 * Write(ioimp,*) 'mulag sans relations pas ok' * goto 9999 endif ip=abs(LAGRAN(J*NB2+1)) iddln=JNT(IP) if (ipaus.le.iddln) then write(ioimp,*) 'Erreur numerotation' write(ioimp,*) 'ip,iddln,ipaus=',ip,iddln,ipaus goto 9999 endif ip=abs(LAGRAN(J*NB2+2)) iddln=JNT(IP) if (ipaur.ge.iddln) then write(ioimp,*) 'Erreur numerotation' write(ioimp,*) 'ip,iddln,ipaur=',ip,iddln,ipaur goto 9999 endif 1700 continue ENDIF C*************************************************************************** SEGSUP JORDRE SEGSUP IJNT 2000 CONTINUE DO 860 I=1,ICPR(/1) IF(ICPR(I).NE.0) ICPR(I)=JNT(ICPR(I)) C numerotation finale. 860 CONTINUE SEGSUP JNT SEGSUP LAGRAN * * Normal termination * RETURN * * Format handling * * * Error handling * 9999 CONTINUE WRITE(IOIMP,*) 'An error was detected in subroutine numop3' RETURN * * End of subroutine NUMOP3 * END