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

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