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

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