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

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