C NUMOPT SOURCE GOUNAND 25/05/15 21:15:10 12268 C RACINE DE LA NUMEROTATION POUR LA SORTIE SUR FAC C SUBROUTINE NUMOPT(MELEME,ICPR,NODES) IMPLICIT INTEGER(I-N) -INC PPARAM -INC CCOPTIO -INC CCPRECO -INC SMELEME -INC SMCOORD CHARACTER*4 MVAL DATA MVAL/'NOOP'/ SEGMENT MEMJT1(NKON) SEGMENT MEMJT2(NKON) SEGMENT IPOME(NODES+1) SEGMENT JNT(NODES) SEGMENT JMEM(NODES) SEGMENT ICPR(nbpts) SEGMENT ICPREN(nbpts) SEGMENT IDCP(nbpts) SEGMENT JOINT(NODES) SEGMENT NEWJT(NODES) icpren=0 segact meleme * cas trivial ou le maillage est fait de points isoles CALL LIRMOT(MVAL,1,IRET,0) itypto=0 ipt1=meleme do k=1,max(1,lisous(/1)) if (lisous(/1).ne.0) ipt1=lisous(k) segact ipt1 itypto = itypto+ipt1.itypel-1 enddo if (itypto.eq.0.or.iret.eq.1.or.nucrou.eq.-1) then segact icpr*MOD ia=0 ipt1=meleme do k=1,max(1,lisous(/1)) if (lisous(/1).ne.0) ipt1=lisous(k) segact ipt1 do j= ipt1.num(/2),1,-1 do i=1,ipt1.num(/1) if (icpr(ipt1.num(i,j)).eq.0) then ia=ia+1 icpr(ipt1.num(i,j))=ia endif enddo enddo enddo nodes=ia return endif nucro=0 * si plus de 4000 noeuds dans le meleme numérotation NS segact icpr*mod do i=1,icpr(/1) icpr(i)=0 enddo ia=0 ib=0 ipd=0 ipt1=meleme ifrot=0 do 60 isous=1,max(1,lisous(/1)) if (lisous(/1).ne.0) ipt1=lisous(isous) segact ipt1 if (ipt1.itypel.eq.-49) then * write (6,*) ' frottement ' ifrot=1 endif NBNN1=ipt1.num(/1) NBEL1=ipt1.num(/2) do i=1,NBEL1 do j=1,NBNN1 iel1=ipt1.num(j,i) if (icpr(iel1).eq.0) then ia=ia+1 icpr(iel1)=ia endif enddo enddo if (NBNN1.gt.ib .AND. NBEL1.GT.0.and.abs(ipt1.itypel).ne.49) $ then ib=NBNN1 ipd=icpr(ipt1.num(1,1)) endif 60 continue if (ia.lt.ib*2 .and.ia.gt.100) then ** write (6,*) ' superelement detecte ',ia,ib ia=1 else ipd=0 endif do i=1,icpr(/1) icpr(i)=0 enddo * choix de la numerotation en fonction du nombre de noeuds ** write(6,*) ' numopt ia ',ia if (ia.gt.4000) nucro=2 if (ia.le.4000) nucro=0 * Reverse CM en iteratif IF(NUCROU.EQ.1) nucro=0 IF(NUCROU.EQ.2) nucro=2 IF(NUCROU.EQ.3) nucro=3 * en cas de frottement forcément nested dissection pour placer correctement les mult de frot * if (ifrot.eq.1) nucro=2 * ajuster le nombre de zero consecutif if (nucro.eq.0) izrosf=16 if (nucro.eq.2) izrosf=48 if(ipd.ne.0) izrosf=48 CALL LIRMOT(MVAL,1,IRET,0) if (iret.eq.0) then * a t'on deja fait l'operation ??? segini icpren,idcp incren=0 do isous=2,lisous(/1) ipt3=lisous(isous) segact ipt3 * SG 2011/10/07 on a échangé l'ordre des boucles pour le remettre dans * l'ordre "usuel" do j=1,ipt3.num(/2) do i=1,ipt3.num(/1) ip=ipt3.num(i,j) if (icpren(ip).eq.0) then incren=incren+1 icpren(ip)=incren endif enddo enddo enddo do i=1,icpren(/1) if (icpren(i).ne.0) idcp(icpren(i))=i enddo do 150 ip=1,npreco if (prenum(ip).ne.0) then ipt1=prenum(ip) segact ipt1 * write (6,*) ' prenum ',ip,lisous(/1),ipt1.lisous(/1) if (lisous(/1)+1.eq.ipt1.lisous(/1)) then *SG Attention ! On a osé stocker le type de renumérotation qui avait été *SG calculée dans le itypel ipt2=ipt1.lisous(ipt1.lisous(/1)) segact ipt2 nucroa=ipt2.itypel segdes ipt2 if (nucro.ne.nucroa) goto 11 * verif maillage identique * write (6,*) ' maillage 1 ',lisous(/1) * write (6,*) ' maillage 2 ',ipt1.lisous(/1) do isous=2,lisous(/1) ipt3=lisous(isous) ipt4=ipt1.lisous(isous) segact ipt3,ipt4 * write (6,*) ipt3.num(/1),ipt4.num(/1),ipt3.num(/2),ipt4.num(/2) * write (6,*) ' itypel ',ipt3.itypel,ipt4.itypel if (ipt3.itypel.ne.ipt4.itypel) goto 10 if (ipt3.num(/1).ne.ipt4.num(/1)) goto 10 if (ipt3.num(/2).ne.ipt4.num(/2)) goto 10 * SG 2011/10/07 on a échangé l'ordre des boucles pour le remettre dans * l'ordre "usuel" do j=1,ipt3.num(/2) do i=1,ipt3.num(/1) * write (6,*) i,j,ipt3.num(i,j), * > icpren(ipt3.num(i,j)),ipt4.num(i,j) if (icpren(ipt3.num(i,j)).ne.ipt4.num(i,j)) $ goto 10 enddo enddo segdes ipt3,ipt4 enddo if (iimpi.eq.2022) write (6,*) > ' preconditionnement de la numerotation ',ip, > (prenum(i),i=1,30) ipt2=ipt1.lisous(ipt1.lisous(/1)) segact ipt2 segact icpr*MOD nodes=ipt2.num(/2) * write (6,*) ' icpr(/1),nodes ',icpr(/1),nodes do i=1,nodes if (idcp(ipt2.num(1,i)).ne.0) icpr(idcp(ipt2.num(1 $ ,i)))=i enddo ia=nodes segdes ipt1,ipt2 * write (6,*) ' apres precond dans numopt ' * write (6,*) (icpr(i),i=1,icpr(/1)) return 10 continue segdes ipt3,ipt4 11 continue endif endif 150 continue if (nucro.eq.1) then call numop1(MELEME,ICPR,NODES) goto 300 elseif (nucro.eq.2) then if (iimpi.ne.0) write(ioimp,*) $ 'Numerotation Nested Dissection' call numop2(MELEME,ICPR,NODES) goto 300 elseif (nucro.eq.3) then if (iimpi.ne.0) write(ioimp,*) 'Numerotation Sloan' call numop3(MELEME,ICPR,NODES) goto 300 endif endif if (iret.eq.1) call refus NODES=nbpts SEGACT ICPR*MOD SEGACT MELEME DO I=1,ICPR(/1) ICPR(I)=0 ENDDO IPT1=MELEME IKOU=0 DO 202 IO=1,MAX(1,LISOUS(/1)) IF (LISOUS(/1).NE.0) THEN IPT1=LISOUS(IO) SEGACT IPT1 ENDIF DO 203 J=1,IPT1.NUM(/2) * inversion car on reinverse a la fin DO 204 I=1,IPT1.NUM(/1) IJ=IPT1.NUM(I,J) IF (ICPR(IJ).NE.0) GOTO 204 IKOU=IKOU+1 ICPR(IJ)=IKOU 204 CONTINUE 203 CONTINUE 202 CONTINUE NODES=IKOU SEGINI JNT,JMEM,NEWJT,JOINT,IPOME IPT1=MELEME DO 3 IO=1,MAX(1,LISOUS(/1)) IF (LISOUS(/1).NE.0) IPT1=LISOUS(IO) DO J=1,IPT1.NUM(/2) DO I=1,IPT1.NUM(/1) JMEM(ICPR(IPT1.NUM(I,J)))=JMEM(ICPR(IPT1.NUM(I,J)))+1 enddo enddo 3 CONTINUE IPOME(1)=0 DO 6 I=1,NODES IPOME(I+1)=IPOME (I) + JMEM(I) 6 CONTINUE NKON=IPOME(NODES+1) SEGINI MEMJT1,memjt2 DO J=1,NODES JMEM(J)=0 enddo IPT1=MELEME DO 102 IPASS=1,1 DO 101 IO=1,MAX(1,LISOUS(/1)) IF (LISOUS(/1).NE.0) IPT1=LISOUS(IO) ** if (ipt1.itypel.eq.-49) then ** if (ipass.eq.1) goto 101 ** else ** if (ipass.eq.2) goto 101 ** endif DO I=1,IPT1.NUM(/2) DO J=1,IPT1.NUM(/1) IND=ICPR(IPT1.NUM(J,I)) JMEM(IND)=JMEM(IND)+1 MEMJT1(IPOME(IND)+JMEM(IND))=I MEMJT2(IPOME(IND)+JMEM(IND))=IO enddo enddo 101 CONTINUE 102 CONTINUE IDIFF=1 if (iimpi.ne.0) write(ioimp,*) $ 'Numerotation Reverse Cuthill McKee' CALL OPTNUM(MELEME,MEMJT1,memjt2,JMEM,JNT,NEWJT,JOINT,IDIFF,IPOME, # ICPR,NODES,NOOPTI,ipd,nucro) * ne pas inverse si superelement pour conserver sa numerotation if (ipd.eq.0) then DO 110 I=1,NODES JNT(I)=NODES-JNT(I)+1 110 CONTINUE endif C PERMUTER LES COORDONNEES ET CORRIGER IPAD DO 111 I=1,icpr(/1) IF (ICPR(I).EQ.0) GOTO 111 ICPR(I)=JNT(ICPR(I)) 111 CONTINUE IF(NOOPTI .EQ.0) GO TO 116 IA=0 DO 115 I=1,icpr(/1) IF(ICPR(I).EQ.0) GO TO 115 IA=IA+1 ICPR(I)=IA 115 CONTINUE 116 CONTINUE SEGSUP MEMJT1,MEMJT2,JMEM,NEWJT,JOINT,IPOME SEGSUP JNT * Mise en commentaire du return par SG car comme cela on preconditionne * aussi le Cuthill-McKee * RETURN 300 continue * sauver le resultat au cas ou on veuille recommencer segini,ipt1=meleme nbsous=ipt1.lisous(/1)+1 nbref=ipt1.lisref(/1) nbnn=ipt1.num(/1) nbelem=ipt1.num(/2) segadj ipt1 nbnn=1 nbsous=0 nbref=0 nbelem=nodes segini ipt2 *SG Attention ! On ose stocker le type de renumérotation qui a été *SG calculée dans le itypel ipt2.itypel=nucro ia=0 do 120 i=1,icpr(/1) if (icpr(i).gt.nodes.or.icpr(i).eq.0) goto 120 ia=ia+1 ipt2.num(1,icpr(i))=icpren(i) 120 continue * call ecmail(ipt2) ipt1.lisous(ipt1.lisous(/1))=ipt2 segdes ipt2 do 125 is=1,ipt1.lisous(/1)-1 ipt2=ipt1.lisous(is) segini,ipt3=ipt2 segdes ipt2 * SG 2011/10/07 on a échangé l'ordre des boucles pour le remettre dans * l'ordre "usuel" do il=1,ipt3.num(/2) do ip=1,ipt3.num(/1) ipt3.num(ip,il)=icpren(ipt3.num(ip,il)) enddo enddo segdes ipt3 ipt1.lisous(is)=ipt3 125 continue if (ipt1.lisous(/1).eq.1) goto 131 segdes ipt1 do 130 ip=npreco-1,1,-1 prenum(ip+1)=prenum(ip) 130 continue prenum(1)=ipt1 131 continue if (icpren.ne.0) segsup icpren,idcp segact icpr*mod * write (6,*) ' fin de numopt ' * write (6,*) (icpr(i),i=1,icpr(/1)) * RETURN END