Télécharger la plaquette Cast3M


fr en

Télécharger numop2.eso

Retour à la liste

Numérotation des lignes : 

numop2

  1. C NUMOP2    SOURCE    PV090527  26/09/15    06:38:59     12641          
  2. C   RACINE DE LA NUMEROTATION      POUR LA SORTIE SUR FAC
  3. C
  4. C       methode utilisee: NESTED DISSECTION.
  5.  
  6. C       ce programme refait l'indicage de la matrice afin de minimiser
  7. C       le profil.
  8.  
  9.  
  10.       SUBROUTINE NUMOP2(MELEME,ICPR,NODES)
  11.       IMPLICIT INTEGER(I-N)
  12. -INC SMELEME
  13. -INC SMCOORD
  14.  
  15. -INC PPARAM
  16. -INC CCOPTIO
  17. -INC CCASSIS
  18. -INC CCREEL
  19.  
  20.         LOGICAL CONNEC
  21.         SEGMENT JMEM(NODES+1),JMEMN(NODES+1)
  22. C       JMEM et JMEMN contiennent le nombre d'element auquel appartient un noeud
  23.  
  24.         SEGMENT JNT(NODES)
  25. C       JNT contient la nouvelle numerotation
  26.  
  27.         SEGMENT ICPR(nbpts)
  28.         segment idcp(nodes)
  29. C       ICPR au debut contient l'ancienne numerotation ,
  30. C                               a la fin la nouvelle.
  31.  
  32.         SEGMENT IADJ(NODES+1)
  33.         SEGMENT JADJC(0)
  34. C       IADJ(i) pointe sur JADJC qui contient les voisins de i entre
  35. C       IADJ(i) et IADJ(i+1)-1
  36.  
  37.         SEGMENT LAGRAN(NB)
  38. C       contient les noeud de lagrange et les noeuds les suivant directement
  39. C       cf element de type 49
  40.  
  41.         SEGMENT BOOLEEN
  42.           LOGICAL BOOL(NODES)
  43.         ENDSEGMENT
  44. C       BOOL(i) = true si le noeud i a ete deja mentionne dans la liste
  45. C       des voisins JADJC.
  46.  
  47.         SEGMENT IMEMOIR(NBV),LMEMOIR(NBV)
  48. C       contient les elements appartenant a chaque noeud,sous forme de liste.
  49.  
  50.         INTEGER ELEM
  51. C       nom d'un element
  52.  
  53.         INTEGER N
  54.  
  55.  
  56.         SEGMENT MASQUE
  57.           LOGICAL MASQ(NODES)
  58.         ENDSEGMENT
  59. C       MASQ(X)=.TRUE. si le noeud X n'a pas ete renumerote;
  60. C       .FALSE. si il l'a ete.
  61.  
  62.         INTEGER DIM,DIMSEP
  63. C       DIM= nombre de noeuds renumerotes.
  64.  
  65.         INTEGER PIVOT
  66. C       PIVOT est le noeud utile a la division du domaine.
  67.  
  68.         SEGMENT IPOS(NODES*3+5)
  69. C       est le vecteur contenant la numerotation dans les deux sens,de 1 a NODES
  70. C       puis de NODES+1 a 2* NODES, cf la subroutine SEPAR
  71. C
  72. C  segments utilisés dans sepa2
  73. C
  74.         SEGMENT NRELONG(NODES*nbthr)
  75. C       NRELONG contient pour chaque noeud sa profondeur.
  76.  
  77.         SEGMENT NOELON(NODES*nbthr)
  78.         SEGMENT NOEL2(NODES)
  79.         SEGMENT LONDIM(NODES*nbthr)
  80. C       NOELON contient les noeuds de profondeur LONG.
  81. C       DIMLON= dimension de NOELON.
  82.  
  83.  
  84. C**********************************
  85.  
  86. C debut du program
  87.  
  88. C**********************************
  89.  
  90.  
  91. *  pour que izero soit de la bonne longueur
  92.         izero=0
  93. C       initialisation
  94. C*******************************
  95.  
  96.  
  97. C       norme d'erreur
  98.         SEGACT ICPR*MOD
  99.         NODES=ICPR(/1)
  100.         SEGACT MELEME
  101. C       icpr: numero des noeuds.
  102. C       meleme: objet de maillage (cf assem2.eso)
  103.  
  104.         DO 10 I=1,ICPR(/1)
  105.           ICPR(I)=0
  106. 10      CONTINUE
  107.  
  108.         IPT1=MELEME
  109.         IKOU=0
  110.         NBV=0
  111.         NB1=0
  112.         NB2=0
  113.  
  114.         DO 100 IO=1,MAX(1,LISOUS(/1))
  115.            IF (LISOUS(/1).GT.0) THEN
  116.                 IPT1=LISOUS(IO)
  117.                 SEGACT IPT1
  118.            ENDIF
  119. C          on cree la numerotation des noeuds.
  120. C          'nb noeuds/element'=IPT1.NUM(/1)
  121. C          'nb element'=IPT1.NUM(/2)
  122.            itypl=abs(ipt1.itypel)
  123.            IF(itypl.EQ.49.or.itypl.eq.22) then
  124.              NB1=NB1+IPT1.NUM(/2)
  125.              NB2=MAX(NB2,IPT1.NUM(/1))
  126. C            NB1= nbre d'éléments de type 49.
  127. C            NB2=nbre de noeuds/élément maximum parmi
  128. C                 les éléments de type 22.
  129.            ENDIF
  130.            DO  J=1,IPT1.NUM(/2)
  131.            DO  I=1,IPT1.NUM(/1)
  132.                IJ=IPT1.NUM(I,J)
  133. C              IJ est le Ième noeud du Jème élément.
  134.                IF (ICPR(IJ).EQ.0) THEN
  135. C              s'il est déjà numéroté, on ne fait rien.
  136.                  IKOU=IKOU+1
  137.                  ICPR(IJ)=IKOU
  138.                ENDIF
  139.            enddo
  140.            enddo
  141. 100     CONTINUE
  142.  
  143.         NODES=IKOU
  144.         NB=NB2*NB1
  145.  
  146. *  reordonner suivant la numerotation initiale
  147. *  pour ne pas etre perturbe par l'elimination des inconnues
  148. *  qui change l'ordre des elements
  149.        segini idcp
  150.        ib=0
  151.        do i=1,icpr(/1)
  152.          if(icpr(i).ne.0) then
  153.            if(idcp(icpr(i)).ne.0) then
  154.            else
  155.             ib=ib+1
  156.             idcp(icpr(i))=i
  157.             icpr(i)=ib
  158.            endif
  159.          endif
  160.        enddo
  161.        segsup idcp
  162. C*****  initalisation des segments*********
  163.  
  164.         SEGINI IADJ,JADJC,JMEM,JMEMN,LAGRAN
  165.         SEGINI BOOLEEN,JNT
  166.  
  167.         DO 20 I=1,NODES+1
  168.                 IADJ(I)=0
  169.                 JMEM(I)=0
  170.                 JMEMN(I)=0
  171. 20      CONTINUE
  172.  
  173. C******************************************
  174.  
  175.         IPT1=MELEME
  176.         IADJ(1)=1
  177.         INC=0
  178.         DO 200 IO=1,MAX(1,LISOUS(/1))
  179.             IF (LISOUS(/1).GT.0) IPT1=LISOUS(IO)
  180.  
  181.             itypl=abs(ipt1.itypel)
  182.             DO 210 J=1,IPT1.NUM(/2)
  183.               IF(ITYPL.EQ.49.or.itypl.eq.22) THEN
  184.                is=sign(1,ipt1.itypel)
  185.                 DO 220 I=1,IPT1.NUM(/1)
  186. C               chaque element de type 49 a  au plus NB2 noeuds.
  187.                    LAGRAN(INC*NB2+I)=ICPR(IPT1.NUM(I,J))*is
  188. C                  les noeuds de l'elements de type 49
  189. C                  sont ranges dans le vecteur LAGRAN.
  190. 220             CONTINUE
  191.                 DO 225 I=IPT1.NUM(/1)+1,NB2
  192.                   LAGRAN(INC*NB2+I)=0
  193. C                 comme on a alloue la place memoire maximale,
  194. C                 on remplit les cases restantes avec des 0.
  195. 225             CONTINUE
  196.                 INC=INC+1
  197. C               INC=le nbre d'elements de type 49.
  198.               ENDIF
  199.  
  200.             DO  230 I=1,IPT1.NUM(/1)
  201.                IJ=ICPR(IPT1.NUM(I,J))+1
  202.                JMEM(IJ)=JMEM(IJ)+1
  203. C              JMEM(I+1): nb elements auquel le noeud I appartient
  204. 230         CONTINUE
  205. 210         CONTINUE
  206.  
  207. 200     CONTINUE
  208.  
  209.  
  210.         JMEM(1)=1
  211.         DO 30 I=1,NODES
  212.             JMEM(I+1)=JMEM(I)+JMEM(I+1)
  213. C           JMEM(I+1)=indice de depart des elements
  214. C           auxquels le noeud I appartient.
  215. 30      CONTINUE
  216.         NBV=JMEM(NODES+1)
  217. C       NBV= dimension de IMEMOIR.
  218.         SEGINI IMEMOIR,LMEMOIR
  219.  
  220.  
  221.  
  222.         IPT1=MELEME
  223.  
  224.         DO 300 IO=1,MAX(1,LISOUS(/1))
  225.            IF (LISOUS(/1).GT.0) THEN
  226.                 IPT1=LISOUS(IO)
  227.            ENDIF
  228.            DO  J=1,IPT1.NUM(/2)
  229.            DO  I=1,IPT1.NUM(/1)
  230.                IJ=ICPR(IPT1.NUM(I,J))
  231.               JMEMN(IJ+1)=JMEMN(IJ+1)+1
  232.               IMEMOIR(JMEM(IJ)+JMEMN(IJ+1)-1)=J
  233.               LMEMOIR(JMEM(IJ)+JMEMN(IJ+1)-1)=IO
  234. C             on range dans IMEMOIR tous les elements des sous-objets
  235. C             IO auxquels appartient le noeud ICPR(IPT1.NUM(I,J)).
  236. C             On connait pour chaque element, le sous-objet auquel
  237. C             il appartient grace a LMEMOIR
  238.            enddo
  239.            enddo
  240. 300     CONTINUE
  241.  
  242.            DO 410 J=1,NODES
  243.                 BOOL(J)=.FALSE.
  244. 410        CONTINUE
  245.         DO 400 I=1,NODES
  246.            IADJ(I+1)=IADJ(I)
  247.            DO 420 J=JMEM(I),JMEM(I+1)-1
  248.                ELEM=IMEMOIR(J)
  249. C              ELEM=element auquel appartient le noeud I.
  250.  
  251.                IPT1=MELEME
  252.                IF (LISOUS(/1).GT.0) IPT1=LISOUS(LMEMOIR(J))
  253.                itype = abs(ipt1.itypel)
  254. *  si element de type 49 ou 22, on ne connecte pas 2 noeuds non LX
  255.                connec=.true.
  256.                if (itype.eq.49.or.itype.eq.22) then
  257.                 kd=3
  258.                 if(itype.eq.22) kd=2
  259.                 do k=kd,ipt1.num(/1)
  260.                   if (i.eq.icpr(ipt1.num(k,elem))) connec=.false.
  261.                 enddo
  262.                endif
  263.                DO 430 K=1,IPT1.NUM(/1)
  264. C              k representatif du nb de noeuds par elements.
  265.                   IK=ICPR(IPT1.NUM(K,ELEM))
  266.                   if (k.ge.kd.and..not.connec) goto 430
  267.                   IF ((I.NE.IK).AND.
  268.      &                 (.NOT.(BOOL(IK)))) THEN
  269. C                 si i n'est pas  egal a un des nouveaux numeros des noeuds
  270. C                 de l'element ELEM et si il n'appartient pas deja  a l'ens des
  271. C                 voisins du noeud i(jadjc(i)),alors on le rajoute.
  272. C                   JADJC(IADJ(I+1))=IK
  273.                     JADJC(**)=IK
  274.                     BOOL(IK)=.TRUE.
  275.                   ENDIF
  276. 430            CONTINUE
  277. 420        CONTINUE
  278.            IADJ(I+1)=JADJC(/1)+1
  279. *  remise a faux de bool
  280.            DO 412 J=IADJ(I),IADJ(I+1)-1
  281.              IK=JADJC(J)
  282.              BOOL(IK)=.FALSE.
  283. 412        CONTINUE
  284. *  tri suivant numero croissant? A priori pas utile
  285. **         call trient(jadjc(iadj(i)),jnt(1),iadj(i+1)-iadj(i))
  286.  
  287. 400     CONTINUE
  288. *  pour ne pas avoir un tableau de longueur nulle, ce que esope n'aime pas
  289.            if(jadjc(/1).eq.0) jadjc(**)=0
  290.  
  291.         SEGSUP      JMEMN,IMEMOIR,LMEMOIR,BOOLEEN
  292.  
  293.  
  294.  
  295. C**************************************************************************
  296.  
  297.  
  298. C       affectation
  299. C************************
  300.  
  301.  
  302.         if (nbthrs.gt.1) call threadii
  303.         SEGINI IPOS,MASQUE
  304.         IPOSMAX=1
  305.  
  306. **      write (6,*) ' nodes ',nodes
  307.         DO 50 I=1,NODES
  308.            MASQ(I)=.TRUE.
  309.            IPOS(I)=0
  310.            IPOS(NODES+I)=0
  311.            IPOS(2*NODES+I)=0
  312. 50      CONTINUE
  313. C       initialement, les noeuds ne sont pas masques,ont donc
  314. C       une position nulle.
  315.  
  316.         DIM=0
  317. C       le nombre de noeuds renumerotes  DIM est initialement egal a zero.
  318. C on initialise un premier separateur avec lisous(1) et lisous(2)
  319. C qui contiennent les noeuds maitres du super element (si appele par assem4)
  320.         mdomn=0
  321.         iposv=ipos(mdomn+1)
  322.         iposmax=iposmax+3
  323.         ipos(iposmax+1-2)=iposv+1
  324.         ipos(iposmax+1-1)=iposv+2
  325.         ipos(iposmax+1-0)=iposv+3
  326.         dimsep=0
  327.         do io=1,2
  328.          ipt1=lisous(io)
  329.           if (ipt1.itypel.eq.0 ) then
  330.           do j=1,ipt1.num(/2)
  331.           do i=1,ipt1.num(/1)
  332.            ip=icpr(ipt1.num(i,j))
  333.            if (masq(ip)) then
  334.              masq(ip)=.false.
  335.              ipos(ip+nodes)=iposmax-2
  336.              dimsep=dimsep+1
  337.            endif
  338.           enddo
  339.           enddo
  340.           endif
  341.         enddo
  342.         do i=1,nodes
  343.           if (masq(i)) then
  344.             ipos(i+nodes)=iposmax
  345.             ipos(i+2*nodes)=ipos(i+2*nodes)+1
  346.           endif
  347.         enddo
  348. **      write(6,*) ' dimsep dans numop2 ',dimsep
  349.         dim=dim+dimsep
  350. C               ****************************************
  351. C                       boucle principale
  352.         NS=NODES
  353.         ns=ns-dimsep
  354.  
  355.         nbthr=min(127,nbthrs)
  356.  2020   continue
  357.         nrelong=0
  358.         noelon=0
  359.         noel2=0
  360.         londim=0
  361.         SEGINI/err=2000/ NRELONG
  362.         SEGINI/err=2000/ NOELON
  363.         SEGINI/err=2000/ noel2
  364.         SEGINI/err=2000/ londim
  365.         goto 2010
  366.  2000   continue
  367.          if(nrelong.ne.0) segsup nrelong
  368.          if(noelon.ne.0) segsup noelon
  369.          if(noel2.ne.0) segsup noel2
  370.          if(londim.ne.0) segsup londim
  371.          nbthr=(nbthr-1)/2+1
  372.          if (nbthr.ne.1) goto 2020
  373.          call erreur(48)
  374.          return
  375.  2010   continue
  376. **      write (6,*) ' avant appel sepa2 '
  377.         DO 500 I=1,NODES
  378. 550       IF(.NOT.MASQ(I)) GOTO 500
  379. C         si le noeud est masque alors ne rien faire: il est deja
  380. C         renumerote. On passe au noeud suivant.
  381.  
  382.           PIVOT=I
  383.  
  384.           CALL SEPA2(IADJ,JADJC,PIVOT,MASQUE,DIMSEP,NS,
  385.      >      IPOS,NODES,IPOSMAX,nrelong,noelon,noel2,
  386.      >      londim,nbthr,izero)
  387.           if (ierr.ne.0) then
  388.             if (nbthrs.gt.1) call threadis
  389.             return
  390.           endif
  391. C         separe le domaine d'etude en 2 parties.
  392. C         on decrit le domaine d'etude a partir du pivot et on cherche la
  393. C         longueur maximale en decrivant les voisins de pivot, et leurs
  394. C         voisins... jusqu'a rencontrer un voisin masque. On cree alors
  395. C         une nouvelle separation.
  396. C         les noeuds masques delimitent la separation du domaine.
  397.  
  398.  
  399.           DIM=DIM+DIMSEP
  400.           NS=NS-DIMSEP
  401. C         la dimension de noeuds renumerotes est augmente de DIMSEP.
  402. C         Celle de noeuds a renumeroter est diminue de DIMSEP.
  403.  
  404. *         IF (DIM.GE.NODES) GOTO 600
  405. C         si tous les noeuds ont ete renumerotes, on arrete.
  406.  
  407.           GOTO 550
  408.  
  409. 500     CONTINUE
  410. **      write (6,*) ' apres appel sepa2 '
  411.  
  412.         SEGSUP NRELONG,NOELON,noel2,londim,jmem
  413.  
  414. 600     CONTINUE
  415.         if (nbthrs.gt.1) call threadis
  416. *
  417. *  tri dans chaque zone
  418. *  je ne sais pas trop pourquoi ca marche
  419.          iposmx=0
  420.          do 610 lpoint=1,nodes
  421.           mdomn=ipos(lpoint+nodes)
  422.           iposv=ipos(mdomn+1)
  423.           iposi=0
  424.          do 620 kk=iadj(lpoint),iadj(lpoint+1)-1
  425.           k=jadjc(kk)
  426.           if(k.eq.lpoint) goto 620
  427.           iposk=ipos(ipos(k+nodes)+1)
  428.            if (iposk.ne.iposv) then
  429.              iposi=max(iposi,iposk)
  430.            endif
  431.  620     continue
  432.          ipos(lpoint+2*nodes)=5*(lpoint+iposi*nodes)
  433.          iposmx=max(iposmx,ipos(lpoint+2*nodes))
  434.  610     continue
  435.          iposmx=iposmx+5
  436. *
  437. *  mise a la bonne place des multiplicateurs de Lagrange
  438.         do 700 J=0,NB1-1
  439.           jp1=j+1
  440.           iposvs=igrand
  441.           iposvr=-igrand
  442.           mdomnr=0
  443.           mdomns=0
  444.           ipaur=0
  445.           ipaus=0
  446. *         write (6,*) 'numop2 ',(lagran(J*NB2+il),il=1,nb2)
  447.         do 800 il=3,nb2
  448.           ip = abs(LAGRAN(J*NB2+il))
  449. *         if (ip.eq.0) write (6,*) ' prob numop2 '
  450.           if (ip.eq.0) goto 800
  451.           mdomn=ipos(ip+nodes)
  452.           iposv=ipos(mdomn+1)
  453. *  deplacer les noeuds en relation en fin de zone
  454.           ipos(ip+2*nodes)=-abs(ipos(ip+2*nodes))
  455.           ipos(ip+2*nodes)=mod(ipos(ip+2*nodes),iposmx)-iposmx
  456.           if (iposv.gt.iposvr) then
  457.               iposvr=iposv
  458.               mdomnr=mdomn
  459.               ipaur=ipos(ip+2*nodes)
  460.           elseif (iposv.eq.iposvr) then
  461.               ipaur=max(ipaur,ipos(ip+2*nodes))
  462.           endif
  463.           if (iposv.lt.iposvs) then
  464.               iposvs=iposv
  465.               mdomns=mdomn
  466.               ipaus=ipos(ip+2*nodes)
  467.           elseif (iposv.eq.iposvs) then
  468.                 ipaus=min(ipaus,ipos(ip+2*nodes))
  469.           endif
  470.  800    continue
  471.  710    continue
  472. *
  473. *  le premier mult avant le premier noeud
  474.         ip=abs(LAGRAN(J*NB2+1))
  475.         IPOS(IP+2*NODES)= ipaur+1
  476. * numeroter le frottement apres le contact
  477.          if (LAGRAN(J*NB2+1).lt.0) 
  478.      >       IPOS(IP+2*NODES)=ipaur+2
  479.         IPOS(IP+nodes)=mdomnr
  480. **      write (6,*) 'premier  mult ',ip,ipos(ip+nodes),ipos(ip+2*nodes)
  481. *
  482. *  le deuxieme mult apres le dernier noeud
  483.         ip=abs(LAGRAN(J*NB2+2))
  484.         if (ip.eq.0) goto 700
  485.         IPOS(IP+2*NODES)= ipaus-1
  486. * numeroter le frottement avant le contact
  487.          if (LAGRAN(J*NB2+2).lt.0) 
  488.      >       IPOS(IP+2*NODES)=ipaus-2
  489. **      IPOS(IP+nodes)=mdomns
  490. * pv deuxieme mult en bout de maillage
  491.         IPOS(IP+nodes)=0                
  492. **      write (6,*) 'deuxieme mult ',ip,ipos(ip+nodes),ipos(ip+2*nodes)
  493. *
  494.  700    continue
  495.         SEGSUP IADJ,JADJC,MASQUE
  496. *  ok maintenant on trie
  497.         CALL SORTI2(IPOS,JNT,NODES)
  498.  
  499. C***************************************************************************
  500.  
  501.         DO 860 I=1,ICPR(/1)
  502.            IF(ICPR(I).NE.0) then
  503.              ICPR(I)=JNT(ICPR(I))
  504. *            write(6,*) 'i icpr ',i,icpr(i)
  505.            endif
  506. C          numerotation finale.
  507. 860     CONTINUE
  508.  
  509.  
  510.         SEGSUP JNT,IPOS,LAGRAN
  511.  
  512.  
  513.         RETURN
  514.         END
  515.  
  516.  
  517.  
  518.  
  519.  
  520.  
  521.  
  522.  
  523.  
  524.  
  525.  
  526.  
  527.  
  528.  
  529.  
  530.  
  531.  
  532.  
  533.  
  534.  
  535.  
  536.  
  537.  
  538.  
  539.  
  540.  
  541.  
  542.  
  543.  
  544.  
  545.  
  546.  
  547.  
  548.  
  549.  
  550.  
  551.  
  552.  
  553.  
  554.  
  555.  
  556.  
  557.  
  558.  
  559.  
  560.  
  561.  
  562.  
  563.  
  564.  

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