C NUMOP2 SOURCE PV090527 26/03/06 21:15:06 12394 C RACINE DE LA NUMEROTATION POUR LA SORTIE SUR FAC C C methode utilisee: NESTED DISSECTION. C ce programme refait l'indicage de la matrice afin de minimiser C le profil. SUBROUTINE NUMOP2(MELEME,ICPR,NODES) IMPLICIT INTEGER(I-N) -INC SMELEME -INC SMCOORD -INC PPARAM -INC CCOPTIO -INC CCASSIS -INC CCREEL 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 ICPR(nbpts) segment idcp(nodes) 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 INVLAG(NB1) 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 INTEGER N SEGMENT MASQUE LOGICAL MASQ(NODES) ENDSEGMENT C MASQ(X)=.TRUE. si le noeud X n'a pas ete renumerote; C .FALSE. si il l'a ete. INTEGER DIM,DIMSEP C DIM= nombre de noeuds renumerotes. INTEGER PIVOT C PIVOT est le noeud utile a la division du domaine. SEGMENT IPOS(NODES*3+5) C est le vecteur contenant la numerotation dans les deux sens,de 1 a NODES C puis de NODES+1 a 2* NODES, cf la subroutine SEPAR C C segments utilisés dans sepa2 C SEGMENT NRELONG(NODES*nbthr) C NRELONG contient pour chaque noeud sa profondeur. SEGMENT NOELON(NODES*nbthr) SEGMENT NOEL2(NODES) SEGMENT LONDIM(NODES*nbthr) C NOELON contient les noeuds de profondeur LONG. C DIMLON= dimension de NOELON. C********************************** C debut du program C********************************** * pour que izero soit de la bonne longueur izero=0 C initialisation C******************************* C norme d'erreur 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 * reordonner suivant la numerotation initiale * pour ne pas etre perturbe par l'elimination des inconnues * qui change l'ordre des elements segini idcp ib=0 do i=1,icpr(/1) if(icpr(i).ne.0) then if(idcp(icpr(i)).ne.0) then else ib=ib+1 idcp(icpr(i))=i icpr(i)=ib endif endif enddo segsup idcp C***** initalisation des segments********* SEGINI IADJ,JADJC,JMEM,JMEMN,LAGRAN SEGINI BOOLEEN,JNT DO 20 I=1,NODES+1 IADJ(I)=0 JMEM(I)=0 JMEMN(I)=0 20 CONTINUE 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 BOOL(IK)=.TRUE. ENDIF 430 CONTINUE 420 CONTINUE IADJ(I+1)=JADJC(/1)+1 * remise a faux de bool DO 412 J=IADJ(I),IADJ(I+1)-1 IK=JADJC(J) BOOL(IK)=.FALSE. 412 CONTINUE * tri suivant numero croissant? A priori pas utile ** call trient(jadjc(iadj(i)),jnt(1),iadj(i+1)-iadj(i)) 400 CONTINUE * pour ne pas avoir un tableau de longueur nulle, ce que esope n'aime pas if(jadjc(/1).eq.0) jadjc(**)=0 SEGSUP JMEMN,IMEMOIR,LMEMOIR,BOOLEEN C************************************************************************** C affectation C************************ if (nbthrs.gt.1) call threadii SEGINI IPOS,MASQUE IPOSMAX=1 ** write (6,*) ' nodes ',nodes DO 50 I=1,NODES MASQ(I)=.TRUE. IPOS(I)=0 IPOS(NODES+I)=0 IPOS(2*NODES+I)=0 50 CONTINUE C initialement, les noeuds ne sont pas masques,ont donc C une position nulle. DIM=0 C le nombre de noeuds renumerotes DIM est initialement egal a zero. C on initialise un premier separateur avec lisous(1) et lisous(2) C qui contiennent les noeuds maitres du super element (si appele par assem4) mdomn=0 iposv=ipos(mdomn+1) iposmax=iposmax+3 ipos(iposmax+1-2)=iposv+1 ipos(iposmax+1-1)=iposv+2 ipos(iposmax+1-0)=iposv+3 dimsep=0 do io=1,2 ipt1=lisous(io) if (ipt1.itypel.eq.0 ) then do j=1,ipt1.num(/2) do i=1,ipt1.num(/1) ip=icpr(ipt1.num(i,j)) if (masq(ip)) then masq(ip)=.false. ipos(ip+nodes)=iposmax-2 dimsep=dimsep+1 endif enddo enddo endif enddo do i=1,nodes if (masq(i)) then ipos(i+nodes)=iposmax ipos(i+2*nodes)=ipos(i+2*nodes)+1 endif enddo ** write(6,*) ' dimsep dans numop2 ',dimsep dim=dim+dimsep C **************************************** C boucle principale NS=NODES ns=ns-dimsep nbthr=min(127,nbthrs) 2020 continue nrelong=0 noelon=0 noel2=0 londim=0 SEGINI/err=2000/ NRELONG SEGINI/err=2000/ NOELON SEGINI/err=2000/ noel2 SEGINI/err=2000/ londim goto 2010 2000 continue if(nrelong.ne.0) segsup nrelong if(noelon.ne.0) segsup noelon if(noel2.ne.0) segsup noel2 if(londim.ne.0) segsup londim nbthr=(nbthr-1)/2+1 if (nbthr.ne.1) goto 2020 call erreur(48) return 2010 continue ** write (6,*) ' avant appel sepa2 ' DO 500 I=1,NODES 550 IF(.NOT.MASQ(I)) GOTO 500 C si le noeud est masque alors ne rien faire: il est deja C renumerote. On passe au noeud suivant. PIVOT=I CALL SEPA2(IADJ,JADJC,PIVOT,MASQUE,DIMSEP,NS, > IPOS,NODES,IPOSMAX,nrelong,noelon,noel2, > londim,nbthr,izero) if (ierr.ne.0) then if (nbthrs.gt.1) call threadis return endif C separe le domaine d'etude en 2 parties. C on decrit le domaine d'etude a partir du pivot et on cherche la C longueur maximale en decrivant les voisins de pivot, et leurs C voisins... jusqu'a rencontrer un voisin masque. On cree alors C une nouvelle separation. C les noeuds masques delimitent la separation du domaine. DIM=DIM+DIMSEP NS=NS-DIMSEP C la dimension de noeuds renumerotes est augmente de DIMSEP. C Celle de noeuds a renumeroter est diminue de DIMSEP. * IF (DIM.GE.NODES) GOTO 600 C si tous les noeuds ont ete renumerotes, on arrete. GOTO 550 500 CONTINUE ** write (6,*) ' apres appel sepa2 ' SEGSUP NRELONG,NOELON,noel2,londim,jmem 600 CONTINUE if (nbthrs.gt.1) call threadis * * tri dans chaque zone * je ne sais pas trop pourquoi ca marche iposmx=0 do 610 lpoint=1,nodes mdomn=ipos(lpoint+nodes) iposv=ipos(mdomn+1) iposi=0 do 620 kk=iadj(lpoint),iadj(lpoint+1)-1 k=jadjc(kk) if(k.eq.lpoint) goto 620 iposk=ipos(ipos(k+nodes)+1) if (iposk.ne.iposv) then iposi=max(iposi,iposk) endif 620 continue ipos(lpoint+2*nodes)=5*(lpoint+iposi*nodes) iposmx=max(iposmx,ipos(lpoint+2*nodes)) 610 continue iposmx=iposmx+5 * * mise a la bonne place des multiplicateurs de Lagrange do 700 J=0,NB1-1 jp1=j+1 iposvs=igrand iposvr=-igrand mdomnr=0 mdomns=0 ipaur=0 ipaus=0 * write (6,*) 'numop2 ',(lagran(J*NB2+il),il=1,nb2) do 800 il=3,nb2 ip = abs(LAGRAN(J*NB2+il)) * if (ip.eq.0) write (6,*) ' prob numop2 ' if (ip.eq.0) goto 710 mdomn=ipos(ip+nodes) iposv=ipos(mdomn+1) * deplacer les noeuds en relation en fin de zone ipos(ip+2*nodes)=-abs(ipos(ip+2*nodes)) ipos(ip+2*nodes)=mod(ipos(ip+2*nodes),iposmx)-iposmx if (iposv.gt.iposvr) then iposvr=iposv mdomnr=mdomn ipaur=ipos(ip+2*nodes) elseif (iposv.eq.iposvr) then ipaur=max(ipaur,ipos(ip+2*nodes)) endif if (iposv.lt.iposvs) then iposvs=iposv mdomns=mdomn ipaus=ipos(ip+2*nodes) elseif (iposv.eq.iposvs) then ipaus=min(ipaus,ipos(ip+2*nodes)) endif 800 continue 710 continue * * le premier mult avant le premier noeud ip=abs(LAGRAN(J*NB2+1)) IPOS(IP+2*NODES)= ipaur+1 * numeroter le frottement apres le contact if (LAGRAN(J*NB2+1).lt.0) > IPOS(IP+2*NODES)=ipaur+2 IPOS(IP+nodes)=mdomnr ** write (6,*) 'premier mult ',ip,ipos(ip+nodes),ipos(ip+2*nodes) * * le deuxieme mult apres le dernier noeud ip=abs(LAGRAN(J*NB2+2)) IPOS(IP+2*NODES)= ipaus-1 * numeroter le frottement avant le contact if (LAGRAN(J*NB2+2).lt.0) > IPOS(IP+2*NODES)=ipaus-2 ** IPOS(IP+nodes)=mdomns * pv deuxieme mult en bout de maillage IPOS(IP+nodes)=0 ** write (6,*) 'deuxieme mult ',ip,ipos(ip+nodes),ipos(ip+2*nodes) * 700 continue SEGSUP IADJ,JADJC,MASQUE * ok maintenant on trie CALL SORTI2(IPOS,JNT,NODES) C*************************************************************************** DO 860 I=1,ICPR(/1) IF(ICPR(I).NE.0) ICPR(I)=JNT(ICPR(I)) C numerotation finale. 860 CONTINUE SEGSUP JNT,IPOS,LAGRAN RETURN END