Télécharger la plaquette Cast3M


fr en

Télécharger noepr2.eso

Retour à la liste

Numérotation des lignes : 

noepr2

  1. C NOEPR2    SOURCE    PV090527  25/09/14    21:15:02     12360          
  2. C****************************************************************************
  3. C****************************************************************************
  4. C*************NOEPERI ..NOEuds PERIpheriques*********************************
  5. C****************************************************************************
  6. C****************************************************************************
  7.  
  8. C  NOEPERI part de PIVOT,lui associe l=1, associe l=l+1 a ses voisins directs,
  9. C  repart des voisins directs pour associer un l a leur voisins....
  10. C  LONG=max(l).
  11. C  NRELONG(I)=l ,NOELON contient les noeuds tels que l=LONG.
  12. *  noepr2 rend en prime la distance a la frontiere, sa taille et le
  13. *  desequilibre des domaines
  14. *  si nbtot n'est pas nul, noepr2 s'arrete des qu'il a trouve la
  15. *  frontiere (il a remplis nbtot/2) pts. nrelong est alors remis a
  16. *  zero
  17. *  Octobre 2014: croissance simultanee de deux zones basees sur les pivots impairs ou pairs
  18. *
  19.         SUBROUTINE NOEPR2(IADJ,JADJC,PIVOT,NRELONG,NOELON,isens,dimlon,
  20.      >  NODES,ipon,NBTOT,lfront,lfron,londim,fcout,lpiv,iccon,icouch)
  21.       IMPLICIT INTEGER(I-N)
  22.       IMPLICIT REAL*8 (A-H,O-Z)
  23. -INC CCREEL
  24.         integer iadj(*),jadjc(*)
  25.         integer pivot(lpiv)
  26.  
  27.         INTEGER LONG,DIMLON,Y,X,dimini,diminj
  28.         integer nrelong(*),noelon(*)
  29.  
  30.         integer ipon(*)
  31.         integer londim(*)
  32.         INTEGER NODES
  33.  
  34.         INTEGER diml1,DIML2,diml3,diml1i,diml2f,diml1f,diml2i
  35.         LOGICAL bool1,bool2,pretest,pretest2
  36.  
  37.         pretest=.false.
  38.         MDOMN=-1
  39.         do i=1,lpiv
  40.           if (pivot(i).ne.0) then
  41.            MDOMN=ipon(PIVOT(i))
  42.            goto 1
  43.           endif
  44.         enddo
  45.    1    continue
  46.         mode=1
  47.         if (nbtot.eq.0) mode=2
  48.         if (nbtot.eq.-1) mode=3
  49.         if (nbtot .gt. nodes)  call erreur(5)
  50.         nodz=nodes
  51.         if (nbtot.gt.0) nodz=nbtot
  52.         DIMLON=0
  53.         LONG=3
  54. *       write (6,*) ' noepr2 lpiv pivot ',lpiv,(pivot(i),i=1,lpiv)
  55.         nbz1=0
  56.         nbz2=0
  57. *  insertion pivot impair zone 1
  58.         DIML1=0
  59.         do 10 i=1,lpiv,2
  60. **       write (6,*) ' i pivot(i) ',i,pivot(i)
  61.          if (pivot(i).le.0) goto 10
  62.          if (ipon(pivot(i)).ne.mdomn) goto 10
  63. *        if (i.eq.1) write (6,*) 'noepr2 mode i pivot nrelong',
  64. *    >    i,pivot(i),nrelong(pivot(i))
  65.          if (nrelong(pivot(i)).ne.0) goto 10
  66.          diml1=diml1+1
  67.          noelon(diml1)=pivot(i)
  68.          nrelong(pivot(i))=1
  69.          nbz1=nbz1+1
  70.   10   continue
  71. *  insertion pivot pair zone 2
  72.         DIML2=nodz+1
  73.         do 11 i=2,lpiv,2
  74.          if (pivot(i).le.0) goto 11
  75. *        write (6,*) 'noepr2 mode i pivot nrelong',
  76. *    >    i,pivot(i),nrelong(pivot(i))
  77.          if (ipon(pivot(i)).ne.mdomn) goto 11
  78.          if (nrelong(pivot(i)).ne.0) goto 11
  79.          diml2=diml2-1
  80.          noelon(diml2)=pivot(i)
  81.          nrelong(pivot(i))=3
  82.          nbz2=nbz2+1
  83.   11   continue
  84.        if (diml1.eq.0.and.diml2.eq.nodz+1) then
  85.          fcout=xgrand
  86.          lfront=1
  87.          return
  88.        endif
  89. *
  90. *  croissance des deux zones
  91. *
  92.        icr1=2
  93.        icr2=3
  94.        bool1=.true.
  95.        bool2=.true.
  96.        diml1i=1
  97.        diml2f=nodz
  98.        icouch=0
  99.   20   continue
  100. ** evaluation frontiere
  101.        if (bool1.and.(icr1*nbz1.le.icr2*nbz2.or..not.bool2)) then
  102.        diml1f=diml1
  103.        if (pretest.and..false.) then
  104.        do 22 i=diml1i,diml1f
  105.             nz1=0
  106.             nz2=0
  107.             x=noelon(i)
  108.             if (nrelong(x).eq.2) goto 22
  109.             DO 29 J=IADJ(X),IADJ(X+1)-1
  110.               Y=JADJC(J)
  111. **            if (mdomn.ne.ipon(y))goto 29
  112.               nrly=nrelong(y)
  113.               if (nrly.eq.1) nz1=nz1+1
  114.               if (nrly.eq.3) nz2=nz2+1
  115.   29        continue
  116.               if (nz2.gt.nz1) then
  117.                 nbz1=nbz1-1
  118.                 nrelong(x)=2
  119.               endif
  120. *             write (6,*) 'frontiere 1 => 2',nz1,nz2
  121.   22   continue
  122.        endif
  123.         bool1=.false.
  124. *       write (6,*) ' diml1i diml1f ',diml1i,diml1f
  125. **      diml1f=diml1
  126.         pretest =.false.
  127.         pretest2=.false.
  128.         do 30 i=diml1i,diml1f
  129.             nz1=0
  130.             nz2=0
  131.             x=noelon(i)
  132.             if(nrelong(x).eq.2) goto 30
  133.             DO 40 J=IADJ(X),IADJ(X+1)-1
  134.                Y=JADJC(J)
  135.                if (mdomn.ne.ipon(y))goto 40
  136.                nrly=nrelong(y)
  137.                if (pretest) then
  138.                  if (nrly.eq.1) nz1=nz1+1
  139.                  if (nrly.eq.3) nz2=nz2+1
  140.                endif
  141.                  IF(NRLY.EQ.0) then
  142. *  insertion zone 1
  143.                    diml1=diml1+1
  144.                    if (.not.bool1) icouch=icouch+1
  145.                    noelon(diml1)=y
  146.                    nrelong(y)=1
  147.                    nbz1=nbz1+1
  148.                    BOOL1=.true.
  149. *     if (mode.ne.2) write (6,*) ' insertion dans 1 de ',diml1,y
  150.                    goto 40
  151.                  elseif (NRLY.eq.1) then
  152. *  zone 1
  153. *     if (mode.ne.2) write (6,*) ' deja en zone 1      ',y
  154.                    goto 40
  155.                  elseif (nry.eq.2) then
  156. *  frontiere
  157. *     if (mode.ne.2) write (6,*) ' deja en frontiere 1 ',y
  158.                    goto 40
  159.                  elseif (nrly.eq.3) then
  160. *  zone 2 <<> frontiere
  161. *     if (mode.ne.2) write (6,*) ' passage 1 sur frontiere de ',y
  162.                    nrelong(y)=2
  163.                    nbz2=nbz2-1
  164.                    pretest2=.true.
  165.                    goto 40
  166.                  endif
  167.   40        continue
  168.               if (pretest.and.nz2.gt.nz1.and.nrelong(x).eq.1) then
  169.                 nbz1=nbz1-1
  170.                 nrelong(x)=2
  171.               endif
  172.   30    continue
  173.         pretest=pretest2
  174.         diml1i=diml1f+1
  175.        endif
  176.        if (bool2.and.(icr1*nbz2.le.icr2*nbz1.or..not.bool1)) then
  177.        diml2i=diml2
  178.        if (pretest.and..false.) then
  179.        do 23 i=diml2f,diml2i,-1
  180.             nz1=0
  181.             nz2=0
  182.             x=noelon(i)
  183.             if (nrelong(x).eq.2) goto 23
  184.             DO 28 J=IADJ(X+1)-1,IADJ(X),-1
  185.               Y=JADJC(J)
  186.               nrly=nrelong(y)
  187. **            if (mdomn.ne.ipon(y))goto 28
  188.               if (nrly.eq.1) nz1=nz1+1
  189.               if (nrly.eq.3) nz2=nz2+1
  190.   28        continue
  191.               if (nz1.gt.nz2) then
  192.                 nbz2=nbz2-1
  193.                 nrelong(x)=2
  194.               endif
  195. *             write (6,*) 'frontiere 3 => 2',nz1,nz2
  196.   23   continue
  197.        endif
  198.         bool2=.false.
  199. **      diml2i=diml2
  200. *       write (6,*) ' diml2f diml2i ',diml2f,diml2i
  201.         do 50 i=diml2f,diml2i,-1
  202.             nz1=0
  203.             nz2=0
  204.             x=noelon(i)
  205.             if (nrelong(x).eq.2) goto 50
  206.             DO 60 J=IADJ(X+1)-1,IADJ(X),-1
  207.                Y=JADJC(J)
  208.                if (mdomn.ne.ipon(y))goto 60
  209.                  nrly=nrelong(y)
  210.                if (pretest) then
  211.                  if (nrly.eq.1) nz1=nz1+1
  212.                  if (nrly.eq.3) nz2=nz2+1
  213.                endif
  214.                  IF(NRLY.EQ.0) then
  215. *  insertion zone 3
  216.                    diml2=diml2-1
  217.                    if (.not.bool2) icouch=icouch+1
  218.                    noelon(diml2)=y
  219.                    nrelong(y)=3
  220.                    nbz2=nbz2+1
  221.                    BOOL2=.true.
  222. *     if (mode.ne.2) write (6,*) ' insertion dans 2 de ',diml2,y
  223.                    goto 60
  224.                  elseif (NRLY.eq.3) then
  225. *  zone 3
  226. *     if (mode.ne.2) write (6,*) ' deja en zone 2      ',y
  227.                    goto 60
  228.                  elseif (nrly.eq.2) then
  229. *  frontiere
  230. *     if (mode.ne.2) write (6,*) ' deja en frontiere 2 ',y
  231.                    goto 60
  232.                  elseif (nrly.eq.1) then
  233. *  zone 1 ==> frontiere
  234. *     if (mode.ne.2) write (6,*) ' passage 2 sur frontiere de ',y
  235.                    nrelong(y)=2
  236.                    nbz1=nbz1-1
  237.                    pretest2=.true.
  238.                    goto 60
  239.                  endif
  240.   60        continue
  241.               if (pretest.and.nz1.gt.nz2.and.nrelong(x).eq.3) then
  242.                 nbz2=nbz2-1
  243.                 nrelong(x)=2
  244.               endif
  245.   50    continue
  246.         pretest=pretest2
  247.         diml2f=diml2i-1
  248.        endif
  249.        if (bool1.or.bool2) goto 20
  250.        nbtotn=diml1+nodz+1-diml2
  251. *  dans le cas ou la zone n'est pas connexe on va completer par
  252. *  la partie non connexe
  253.         if (mode.ne.2.and.nbtot.gt.nbtotn) then
  254. *       write (6,*) ' ajout autres composantes connexes ',diml1,diml2,
  255. *    >  nbtotn,nbtot,nodz
  256. *  si pas trop de noeuds
  257. *  on commence par examiner le voisinage de la frontiere car c'est moins cher
  258. ****    if (mode.ne.3) goto 21
  259.         diml1i=diml1
  260.         do 200 i=1,nbtot
  261.          x=noelon(i)
  262.          if (x.eq.0) goto 200
  263.          if (nrelong(x).ne.2) goto 200
  264.             DO 210 J=IADJ(X),IADJ(X+1)-1
  265.                Y=JADJC(J)
  266.                if (mdomn.ne.ipon(y))goto 210
  267.                  IF(NRELONG(Y).EQ.0) then
  268.                     DIML1=DIML1+1
  269.                     NOELON(DIML1)=Y
  270.                     NRELONG(Y)=1
  271.                     nbtotn=nbtotn+1
  272.                     if (nbtot.le.nbtotn) goto 21
  273.                endif
  274.  210        continue
  275.  200    continue
  276.           if (diml1.eq.diml1i) goto 220
  277. *  Si on a rajoute des noeuds, on fait aussi leurs voisinage
  278.  230    continue
  279.           nodi=diml1i+1
  280.           diml1i=diml1
  281.           nodf=diml1
  282.         do 240 i=nodi,nodf
  283.          x=noelon(i)
  284.             DO 250 J=IADJ(X),IADJ(X+1)-1
  285.                Y=JADJC(J)
  286.                if (mdomn.ne.ipon(y))goto 250
  287.                  IF(NRELONG(Y).EQ.0) then
  288.                     DIML1=DIML1+1
  289.                     NOELON(DIML1)=Y
  290.                     NRELONG(Y)=1
  291.                     nbtotn=nbtotn+1
  292.                     if (nbtot.le.nbtotn) goto 21
  293.                endif
  294.  250        continue
  295.  240    continue
  296.         if (diml1.ne.diml1i) goto 230
  297. 220     continue
  298.         if (mode.ne.3) goto 21
  299.  
  300. *       write (6,*) ' ajout 2autres composantes connexes ',diml1,diml2,
  301. *    >  nbtotn,nbtot,nodz
  302.   21    continue
  303.         endif
  304. * Si on n'a toujours pas notre compte, on balaye tout
  305.         if ((iccon.eq.1.and.mode.eq.2).or. 
  306.      >       (mode.eq.1.and.nbtot.gt.nbtotn)) then
  307.         do Y=1,NODes
  308.           IF((NRELONG(Y).EQ.0).AND.(MDOMN.EQ.ipon(Y))) THEN
  309.                     DIML1=DIML1+1
  310.                     NOELON(DIML1)=Y
  311.                     NRELONG(Y)=1
  312.                     nbtotn=nbtotn+1
  313.                     if (mode.eq.1.and.nbtot.le.nbtotn) goto 215
  314.           ENDIF
  315.         enddo
  316.   215   continue
  317.         endif
  318.        nbtot=nbtotn
  319. *
  320. *  remise au carre des deux zones et de la frontiere
  321. *
  322.        nbtot=diml1+nodz+1-diml2
  323. *      write (6,*) ' nodes diml1 diml2 ',nodes,diml1,diml2,nbtot
  324. *  classer a la fin par permutation dans zone 1 la frontiere
  325.        diml1f=diml1
  326.        do 100 i=1,diml1-1
  327.          x=noelon(i)
  328.            if (nrelong(x).eq.2) then
  329.              do 105 j=diml1f,i+1,-1
  330.                y=noelon(j)
  331.                if (nrelong(y).eq.1) goto 106
  332.  105         continue
  333.              goto 107
  334.  106        continue
  335.             diml1f=j-1
  336.             noelon(i)=y
  337.             noelon(j)=x
  338.            endif
  339.  100   continue
  340.  107   continue
  341. *  classer au debut par permutation dans zone 2 la frontiere
  342.        diml2i=diml2
  343.        do 110 i=nodz,diml2+1,-1
  344.          x=noelon(i)
  345.            if (nrelong(x).eq.2) then
  346.              do 115 j=diml2i,i-1
  347.                y=noelon(j)
  348.                if (nrelong(y).eq.3) goto 116
  349.  115         continue
  350.              goto 117
  351.  116       continue
  352.            diml2i=j+1
  353.            noelon(i)=y
  354.            noelon(j)=x
  355.            endif
  356.  110   continue
  357.  117   continue
  358. * suprimer le trou entre les deux
  359.            if (nodz.ne.nbtot) then
  360.             do 120 i=diml2,nodz
  361.              noelon(i-nodz+nbtot)=noelon(i)
  362.  120        continue
  363.            endif
  364.  
  365. * decompte des longueurs
  366.            n1bz1=0
  367.            n1bz2=0
  368.            nbfr=0
  369.            do 130 i=1,nbtot
  370.             x=noelon(i)
  371.             nrlx=nrelong(x)
  372.             if (nrlx.eq.1) n1bz1=n1bz1+1
  373.             if (nrlx.eq.2) nbfr=nbfr+1
  374.             if (nrlx.eq.3) n1bz2=n1bz2+1
  375.  130       continue
  376. **         if (nbfr.ne.nbfra) write (6,*) ' nbfr nbfra ',nbfr,nbfra
  377. *           write (6,*) 'nbtot n1bz1 nbfr n1bz2',nbtot,n1bz1,nbfr,n1bz2
  378. *
  379. *  finalisation des infos
  380. *
  381.        londim(1)=n1bz1
  382.        londim(2)=n1bz1+nbfr
  383.        londim(3)=n1bz1+nbfr+n1bz2
  384.        long=3
  385.        dimlon=londim(3)
  386.         ideseq=abs(n1bz1-n1bz2)
  387.         xbtot=nbtot
  388.         xbfr=nbfr
  389.         xdeseq=ideseq
  390. **      fcout=(xbfr        )**4  +  (xdeseq )**3
  391.  
  392.         fcout= 2*xbtot*xbfr+xbfr*xbfr*xbfr/3.0 + xdeseq**2 
  393.      >     -icouch**3.0
  394.         fcout=max(0.d0,fcout)
  395. **      write (6,*) 'xbtot',nbtot,'xbfr',nbfr,'xdeseq',ideseq,'fcout',
  396. **   >   fcout,'icouch',icouch
  397.  
  398.         if (nbtot/(xdeseq+1).le.2)
  399.      >     fcout=fcout+xbtot*xbtot*xbtot*xbtot
  400.         if (nbfr.eq.0.and.lpiv.ne.1.and.pivot(1).ne.0.and.
  401.      >        pivot(2).ne.0) then
  402. **       write (6,*) ' frontiere vide ',n1bz1,n1bz2,nbtot,mdomn
  403.          if (mode.ne.2.and.n1bz1*n1bz2.ne.0) fcout=1
  404.         endif
  405.  
  406. *       write (6,*) ' noepr2 ',n1bz1,n1bz2,nbfr,nbtot,
  407. *    >        noelon(1),noelon(nbtot),fcout
  408.         isens=2
  409.         if (n1bz1.lt.n1bz2) isens=1
  410. *  mise a zero eventuelle de nrelong
  411.         if (mode.ne.2) then
  412.          do i=1,dimlon
  413.           x=noelon(i)
  414.           if (x.ne.0) nrelong(x)=0
  415.          enddo
  416.         endif
  417.         lfron = long
  418.         lfront=2
  419.         if (n1bz2.eq.0.or.n1bz1.eq.0) lfront=1
  420.         if (mode.eq.1) nbtot=nodz
  421.         return
  422.         end
  423.  
  424.  
  425.  
  426.  
  427.  
  428.  
  429.  
  430.  
  431.  
  432.  
  433.  
  434.  
  435.  
  436.  
  437.  
  438.  
  439.  
  440.  
  441.  
  442.  
  443.  
  444.  
  445.  
  446.  
  447.  
  448.  
  449.  
  450.  
  451.  
  452.  
  453.  
  454.  
  455.  
  456.  
  457.  
  458.  
  459.  
  460.  
  461.  
  462.  
  463.  
  464.  

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