Télécharger la plaquette Cast3M


fr en

Télécharger numopt.eso

Retour à la liste

Numérotation des lignes : 

numopt

  1. C NUMOPT    SOURCE    GOUNAND   25/05/15    21:15:10     12268          
  2. C   RACINE DE LA NUMEROTATION      POUR LA SORTIE SUR FAC
  3. C
  4.       SUBROUTINE NUMOPT(MELEME,ICPR,NODES)
  5.       IMPLICIT INTEGER(I-N)
  6.  
  7. -INC PPARAM
  8. -INC CCOPTIO
  9. -INC CCPRECO
  10.  
  11. -INC SMELEME
  12. -INC SMCOORD
  13.  
  14.       CHARACTER*4 MVAL
  15.       DATA MVAL/'NOOP'/
  16.       SEGMENT MEMJT1(NKON)
  17.       SEGMENT MEMJT2(NKON)
  18.       SEGMENT IPOME(NODES+1)
  19.       SEGMENT JNT(NODES)
  20.       SEGMENT JMEM(NODES)
  21.       SEGMENT ICPR(nbpts)
  22.       SEGMENT ICPREN(nbpts)
  23.       SEGMENT IDCP(nbpts)
  24.       SEGMENT JOINT(NODES)
  25.       SEGMENT NEWJT(NODES)
  26.  
  27.       icpren=0
  28.       segact meleme
  29. *  cas trivial ou le maillage est fait de points isoles
  30.       CALL LIRMOT(MVAL,1,IRET,0)
  31.       itypto=0
  32.       ipt1=meleme
  33.       do k=1,max(1,lisous(/1))
  34.          if (lisous(/1).ne.0) ipt1=lisous(k)
  35.          segact ipt1
  36.          itypto = itypto+ipt1.itypel-1
  37.       enddo
  38.  
  39.       if (itypto.eq.0.or.iret.eq.1.or.nucrou.eq.-1) then
  40.          segact icpr*MOD
  41.          ia=0
  42.          ipt1=meleme
  43.          do k=1,max(1,lisous(/1))
  44.             if (lisous(/1).ne.0) ipt1=lisous(k)
  45.             segact ipt1
  46.             do j=  ipt1.num(/2),1,-1
  47.                do i=1,ipt1.num(/1)
  48.                   if (icpr(ipt1.num(i,j)).eq.0) then
  49.                      ia=ia+1
  50.                      icpr(ipt1.num(i,j))=ia
  51.                   endif
  52.                enddo
  53.             enddo
  54.          enddo
  55.          nodes=ia
  56.          return
  57.       endif
  58.  
  59.       nucro=0
  60. *  si plus de 4000 noeuds dans le meleme numérotation NS
  61.       segact icpr*mod
  62.       do i=1,icpr(/1)
  63.          icpr(i)=0
  64.       enddo
  65.       ia=0
  66.       ib=0
  67.       ipd=0
  68.       ipt1=meleme
  69.       ifrot=0
  70.       do 60 isous=1,max(1,lisous(/1))
  71.          if (lisous(/1).ne.0) ipt1=lisous(isous)
  72.          segact ipt1
  73.          if (ipt1.itypel.eq.-49) then
  74. *         write (6,*) ' frottement '
  75.             ifrot=1
  76.          endif
  77.          NBNN1=ipt1.num(/1)
  78.          NBEL1=ipt1.num(/2)
  79.          do  i=1,NBEL1
  80.             do  j=1,NBNN1
  81.                iel1=ipt1.num(j,i)
  82.                if (icpr(iel1).eq.0) then
  83.                   ia=ia+1
  84.                   icpr(iel1)=ia
  85.                endif
  86.             enddo
  87.          enddo
  88.  
  89.          if (NBNN1.gt.ib .AND. NBEL1.GT.0.and.abs(ipt1.itypel).ne.49)
  90.      $        then
  91.             ib=NBNN1
  92.             ipd=icpr(ipt1.num(1,1))
  93.          endif
  94.  60   continue
  95.       if (ia.lt.ib*2 .and.ia.gt.100) then
  96. **     write (6,*) ' superelement detecte ',ia,ib
  97.          ia=1
  98.       else
  99.          ipd=0
  100.       endif
  101.       do i=1,icpr(/1)
  102.          icpr(i)=0
  103.       enddo
  104. *  choix de la numerotation en fonction du nombre de noeuds
  105. **    write(6,*) ' numopt ia ',ia
  106.       if (ia.gt.4000) nucro=2
  107.       if (ia.le.4000) nucro=0
  108. *  Reverse CM en iteratif
  109.       IF(NUCROU.EQ.1) nucro=0
  110.       IF(NUCROU.EQ.2) nucro=2
  111.       IF(NUCROU.EQ.3) nucro=3
  112. *  en cas de frottement forcément nested dissection pour placer correctement les mult de frot
  113. *     if (ifrot.eq.1) nucro=2
  114. *  ajuster le nombre de zero consecutif
  115.       if (nucro.eq.0) izrosf=16
  116.       if (nucro.eq.2) izrosf=48
  117.       if(ipd.ne.0) izrosf=48
  118.       CALL LIRMOT(MVAL,1,IRET,0)
  119.       if (iret.eq.0) then
  120. *  a t'on deja fait l'operation ???
  121.          segini icpren,idcp
  122.          incren=0
  123.          do isous=2,lisous(/1)
  124.             ipt3=lisous(isous)
  125.             segact ipt3
  126. * SG 2011/10/07 on a échangé l'ordre des boucles pour le remettre dans
  127. *               l'ordre "usuel"
  128.             do j=1,ipt3.num(/2)
  129.                do i=1,ipt3.num(/1)
  130.                   ip=ipt3.num(i,j)
  131.                   if (icpren(ip).eq.0) then
  132.                      incren=incren+1
  133.                      icpren(ip)=incren
  134.                   endif
  135.                enddo
  136.             enddo
  137.          enddo
  138.          do i=1,icpren(/1)
  139.             if (icpren(i).ne.0) idcp(icpren(i))=i
  140.          enddo
  141.  
  142.          do 150 ip=1,npreco
  143.             if (prenum(ip).ne.0) then
  144.                ipt1=prenum(ip)
  145.                segact ipt1
  146. *      write (6,*) ' prenum ',ip,lisous(/1),ipt1.lisous(/1)
  147.                if (lisous(/1)+1.eq.ipt1.lisous(/1)) then
  148. *SG Attention ! On a osé stocker le type de renumérotation qui avait été
  149. *SG   calculée dans le itypel
  150.                   ipt2=ipt1.lisous(ipt1.lisous(/1))
  151.                   segact ipt2
  152.                   nucroa=ipt2.itypel
  153.                   segdes ipt2
  154.                   if (nucro.ne.nucroa) goto 11
  155. *  verif maillage identique
  156. *      write (6,*) ' maillage 1 ',lisous(/1)
  157. *      write (6,*) ' maillage 2 ',ipt1.lisous(/1)
  158.                   do isous=2,lisous(/1)
  159.                      ipt3=lisous(isous)
  160.                      ipt4=ipt1.lisous(isous)
  161.                      segact ipt3,ipt4
  162. *       write (6,*) ipt3.num(/1),ipt4.num(/1),ipt3.num(/2),ipt4.num(/2)
  163. *       write (6,*) ' itypel ',ipt3.itypel,ipt4.itypel
  164.                      if (ipt3.itypel.ne.ipt4.itypel) goto 10
  165.                      if (ipt3.num(/1).ne.ipt4.num(/1)) goto 10
  166.                      if (ipt3.num(/2).ne.ipt4.num(/2)) goto 10
  167. * SG 2011/10/07 on a échangé l'ordre des boucles pour le remettre dans
  168. *               l'ordre "usuel"
  169.                      do j=1,ipt3.num(/2)
  170.                         do i=1,ipt3.num(/1)
  171. *        write (6,*) i,j,ipt3.num(i,j),
  172. *    >     icpren(ipt3.num(i,j)),ipt4.num(i,j)
  173.                            if (icpren(ipt3.num(i,j)).ne.ipt4.num(i,j))
  174.      $                          goto 10
  175.                         enddo
  176.                      enddo
  177.                      segdes ipt3,ipt4
  178.                   enddo
  179.                   if (iimpi.eq.2022) write (6,*)
  180.      >                 ' preconditionnement de la numerotation ',ip,
  181.      >                 (prenum(i),i=1,30)
  182.                   ipt2=ipt1.lisous(ipt1.lisous(/1))
  183.                   segact ipt2
  184.                   segact icpr*MOD
  185.                   nodes=ipt2.num(/2)
  186. *       write (6,*) ' icpr(/1),nodes ',icpr(/1),nodes
  187.                   do i=1,nodes
  188.                      if (idcp(ipt2.num(1,i)).ne.0) icpr(idcp(ipt2.num(1
  189.      $                    ,i)))=i
  190.                   enddo
  191.                   ia=nodes
  192.                   segdes ipt1,ipt2
  193. *     write (6,*) ' apres precond dans numopt '
  194. *     write (6,*) (icpr(i),i=1,icpr(/1))
  195.                   return
  196.  
  197.  10               continue
  198.                   segdes ipt3,ipt4
  199.  11               continue
  200.                endif
  201.             endif
  202.  150     continue
  203.          if (nucro.eq.1) then
  204.             call numop1(MELEME,ICPR,NODES)
  205.             goto 300
  206.          elseif (nucro.eq.2) then
  207.             if (iimpi.ne.0) write(ioimp,*)
  208.      $           'Numerotation Nested Dissection'
  209.             call numop2(MELEME,ICPR,NODES)
  210.             goto 300
  211.          elseif (nucro.eq.3) then
  212.             if (iimpi.ne.0) write(ioimp,*) 'Numerotation Sloan'
  213.             call numop3(MELEME,ICPR,NODES)
  214.             goto 300
  215.          endif
  216.       endif
  217.       if (iret.eq.1) call refus
  218.       NODES=nbpts
  219.       SEGACT ICPR*MOD
  220.       SEGACT MELEME
  221.       DO  I=1,ICPR(/1)
  222.          ICPR(I)=0
  223.       ENDDO
  224.       IPT1=MELEME
  225.       IKOU=0
  226.       DO 202 IO=1,MAX(1,LISOUS(/1))
  227.          IF (LISOUS(/1).NE.0) THEN
  228.             IPT1=LISOUS(IO)
  229.             SEGACT IPT1
  230.          ENDIF
  231.          DO 203 J=1,IPT1.NUM(/2)
  232. *  inversion car on reinverse a la fin
  233.             DO 204 I=1,IPT1.NUM(/1)
  234.                IJ=IPT1.NUM(I,J)
  235.                IF (ICPR(IJ).NE.0) GOTO 204
  236.                IKOU=IKOU+1
  237.                ICPR(IJ)=IKOU
  238.  204        CONTINUE
  239.  203     CONTINUE
  240.  202  CONTINUE
  241.       NODES=IKOU
  242.       SEGINI JNT,JMEM,NEWJT,JOINT,IPOME
  243.       IPT1=MELEME
  244.       DO 3 IO=1,MAX(1,LISOUS(/1))
  245.          IF (LISOUS(/1).NE.0) IPT1=LISOUS(IO)
  246.          DO  J=1,IPT1.NUM(/2)
  247.             DO  I=1,IPT1.NUM(/1)
  248.                JMEM(ICPR(IPT1.NUM(I,J)))=JMEM(ICPR(IPT1.NUM(I,J)))+1
  249.             enddo
  250.          enddo
  251.  3    CONTINUE
  252.       IPOME(1)=0
  253.       DO 6 I=1,NODES
  254.          IPOME(I+1)=IPOME (I) + JMEM(I)
  255.    6  CONTINUE
  256.       NKON=IPOME(NODES+1)
  257.       SEGINI MEMJT1,memjt2
  258.       DO  J=1,NODES
  259.          JMEM(J)=0
  260.       enddo
  261.       IPT1=MELEME
  262.       DO 102 IPASS=1,1
  263.          DO 101 IO=1,MAX(1,LISOUS(/1))
  264.             IF (LISOUS(/1).NE.0) IPT1=LISOUS(IO)
  265. **    if (ipt1.itypel.eq.-49) then
  266. **     if (ipass.eq.1) goto 101
  267. **    else
  268. **     if (ipass.eq.2) goto 101
  269. **    endif
  270.             DO  I=1,IPT1.NUM(/2)
  271.                DO  J=1,IPT1.NUM(/1)
  272.                   IND=ICPR(IPT1.NUM(J,I))
  273.                   JMEM(IND)=JMEM(IND)+1
  274.                   MEMJT1(IPOME(IND)+JMEM(IND))=I
  275.                   MEMJT2(IPOME(IND)+JMEM(IND))=IO
  276.                enddo
  277.             enddo
  278.  101     CONTINUE
  279.  102  CONTINUE
  280.       IDIFF=1
  281.       if (iimpi.ne.0) write(ioimp,*)
  282.      $     'Numerotation Reverse Cuthill McKee'
  283.       CALL OPTNUM(MELEME,MEMJT1,memjt2,JMEM,JNT,NEWJT,JOINT,IDIFF,IPOME,
  284.      # ICPR,NODES,NOOPTI,ipd,nucro)
  285. *  ne pas inverse si superelement pour conserver sa numerotation
  286.       if (ipd.eq.0) then
  287.          DO 110 I=1,NODES
  288.             JNT(I)=NODES-JNT(I)+1
  289.  110     CONTINUE
  290.       endif
  291. C  PERMUTER LES COORDONNEES ET CORRIGER IPAD
  292.       DO 111 I=1,icpr(/1)
  293.          IF (ICPR(I).EQ.0) GOTO 111
  294.          ICPR(I)=JNT(ICPR(I))
  295.  111  CONTINUE
  296.       IF(NOOPTI .EQ.0)  GO TO 116
  297.       IA=0
  298.       DO 115 I=1,icpr(/1)
  299.          IF(ICPR(I).EQ.0) GO TO 115
  300.          IA=IA+1
  301.          ICPR(I)=IA
  302.  115  CONTINUE
  303.  116  CONTINUE
  304.       SEGSUP MEMJT1,MEMJT2,JMEM,NEWJT,JOINT,IPOME
  305.       SEGSUP JNT
  306. * Mise en commentaire du return par SG car comme cela on preconditionne
  307. * aussi le Cuthill-McKee
  308. *      RETURN
  309.  300  continue
  310. *  sauver le resultat au cas ou on veuille recommencer
  311.       segini,ipt1=meleme
  312.       nbsous=ipt1.lisous(/1)+1
  313.       nbref=ipt1.lisref(/1)
  314.       nbnn=ipt1.num(/1)
  315.       nbelem=ipt1.num(/2)
  316.       segadj ipt1
  317.       nbnn=1
  318.       nbsous=0
  319.       nbref=0
  320.       nbelem=nodes
  321.       segini ipt2
  322. *SG Attention ! On ose stocker le type de renumérotation qui a été
  323. *SG   calculée dans le itypel
  324.       ipt2.itypel=nucro
  325.       ia=0
  326.       do 120 i=1,icpr(/1)
  327.          if (icpr(i).gt.nodes.or.icpr(i).eq.0) goto 120
  328.          ia=ia+1
  329.          ipt2.num(1,icpr(i))=icpren(i)
  330.  120  continue
  331. *     call ecmail(ipt2)
  332.       ipt1.lisous(ipt1.lisous(/1))=ipt2
  333.       segdes ipt2
  334.       do 125 is=1,ipt1.lisous(/1)-1
  335.          ipt2=ipt1.lisous(is)
  336.          segini,ipt3=ipt2
  337.          segdes ipt2
  338. * SG 2011/10/07 on a échangé l'ordre des boucles pour le remettre dans
  339. *               l'ordre "usuel"
  340.          do il=1,ipt3.num(/2)
  341.             do ip=1,ipt3.num(/1)
  342.                ipt3.num(ip,il)=icpren(ipt3.num(ip,il))
  343.             enddo
  344.          enddo
  345.          segdes ipt3
  346.          ipt1.lisous(is)=ipt3
  347.  125  continue
  348.       if (ipt1.lisous(/1).eq.1) goto 131
  349.       segdes ipt1
  350.       do 130 ip=npreco-1,1,-1
  351.          prenum(ip+1)=prenum(ip)
  352.  130  continue
  353.       prenum(1)=ipt1
  354.  131  continue
  355.       if (icpren.ne.0) segsup icpren,idcp
  356.       segact icpr*mod
  357. *     write (6,*) ' fin de numopt '
  358. *     write (6,*) (icpr(i),i=1,icpr(/1))
  359. *
  360.       RETURN
  361.       END
  362.  
  363.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales