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  23/03/14    21:15:10     11629          
  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. C       ICPR au debut contient l'ancienne numerotation ,
  29. C                               a la fin la nouvelle.
  30.  
  31.         SEGMENT IADJ(NODES+1)
  32.         SEGMENT JADJC(0)
  33. C       IADJ(i) pointe sur JADJC qui contient les voisins de i entre
  34. C       IADJ(i) et IADJ(i+1)-1
  35.  
  36.         SEGMENT LAGRAN(NB)
  37. C       contient les noeud de lagrange et les noeuds les suivant directement
  38. C       cf element de type 49
  39.         SEGMENT INVLAG(NB1)
  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)
  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.            IF(abs(IPT1.ITYPEL).EQ.49) THEN
  123.              NB1=NB1+IPT1.NUM(/2)
  124.              NB2=MAX(NB2,IPT1.NUM(/1))
  125. C            NB1= nbre d'éléments de type 49.
  126. C            NB2=nbre de noeuds/élément maximum parmi
  127. C                 les éléments de type 22.
  128.            ENDIF
  129.            DO  J=1,IPT1.NUM(/2)
  130.            DO  I=1,IPT1.NUM(/1)
  131.                IJ=IPT1.NUM(I,J)
  132. C              IJ est le Ième noeud du Jème élément.
  133.                IF (ICPR(IJ).EQ.0) THEN
  134. C              s'il est déjà numéroté, on ne fait rien.
  135.                  IKOU=IKOU+1
  136.                  ICPR(IJ)=IKOU
  137.                ENDIF
  138.            enddo
  139.            enddo
  140. 100     CONTINUE
  141.  
  142.         NODES=IKOU
  143.         NB=NB2*NB1
  144.  
  145. C*****  initalisation des segments*********
  146.  
  147.         SEGINI IADJ,JADJC,JMEM,JMEMN,LAGRAN
  148.         SEGINI BOOLEEN,JNT
  149.  
  150.         DO 20 I=1,NODES+1
  151.                 IADJ(I)=0
  152.                 JMEM(I)=0
  153.                 JMEMN(I)=0
  154. 20      CONTINUE
  155.  
  156. C******************************************
  157.  
  158.         IPT1=MELEME
  159.         IADJ(1)=1
  160.         INC=0
  161.         DO 200 IO=1,MAX(1,LISOUS(/1))
  162.             IF (LISOUS(/1).GT.0) IPT1=LISOUS(IO)
  163.  
  164.             DO 210 J=1,IPT1.NUM(/2)
  165.               IF(abs(IPT1.ITYPEL).EQ.49) THEN
  166.                is=sign(1,ipt1.itypel)
  167.                 DO 220 I=1,IPT1.NUM(/1)
  168. C               chaque element de type 49 a  au plus NB2 noeuds.
  169.                    LAGRAN(INC*NB2+I)=ICPR(IPT1.NUM(I,J))*is
  170. C                  les noeuds de l'elements de type 49
  171. C                  sont ranges dans le vecteur LAGRAN.
  172. 220             CONTINUE
  173.                 DO 225 I=IPT1.NUM(/1)+1,NB2
  174.                   LAGRAN(INC*NB2+I)=0
  175. C                 comme on a alloue la place memoire maximale,
  176. C                 on remplit les cases restantes avec des 0.
  177. 225             CONTINUE
  178.                 INC=INC+1
  179. C               INC=le nbre d'elements de type 49.
  180.               ENDIF
  181.  
  182.             DO  230 I=1,IPT1.NUM(/1)
  183.                IJ=ICPR(IPT1.NUM(I,J))+1
  184.                JMEM(IJ)=JMEM(IJ)+1
  185. C              JMEM(I+1): nb elements auquel le noeud I appartient
  186. 230         CONTINUE
  187. 210         CONTINUE
  188.  
  189. 200     CONTINUE
  190.  
  191.  
  192.         JMEM(1)=1
  193.         DO 30 I=1,NODES
  194.             JMEM(I+1)=JMEM(I)+JMEM(I+1)
  195. C           JMEM(I+1)=indice de depart des elements
  196. C           auxquels le noeud I appartient.
  197. 30      CONTINUE
  198.         NBV=JMEM(NODES+1)
  199. C       NBV= dimension de IMEMOIR.
  200.         SEGINI IMEMOIR,LMEMOIR
  201.  
  202.  
  203.  
  204.         IPT1=MELEME
  205.  
  206.         DO 300 IO=1,MAX(1,LISOUS(/1))
  207.            IF (LISOUS(/1).GT.0) THEN
  208.                 IPT1=LISOUS(IO)
  209.            ENDIF
  210.            DO  J=1,IPT1.NUM(/2)
  211.            DO  I=1,IPT1.NUM(/1)
  212.                IJ=ICPR(IPT1.NUM(I,J))
  213.               JMEMN(IJ+1)=JMEMN(IJ+1)+1
  214.               IMEMOIR(JMEM(IJ)+JMEMN(IJ+1)-1)=J
  215.               LMEMOIR(JMEM(IJ)+JMEMN(IJ+1)-1)=IO
  216. C             on range dans IMEMOIR tous les elements des sous-objets
  217. C             IO auxquels appartient le noeud ICPR(IPT1.NUM(I,J)).
  218. C             On connait pour chaque element, le sous-objet auquel
  219. C             il appartient grace a LMEMOIR
  220.            enddo
  221.            enddo
  222. 300     CONTINUE
  223.  
  224.            DO 410 J=1,NODES
  225.                 BOOL(J)=.FALSE.
  226. 410        CONTINUE
  227.         DO 400 I=1,NODES
  228.            IADJ(I+1)=IADJ(I)
  229.            DO 420 J=JMEM(I),JMEM(I+1)-1
  230.                ELEM=IMEMOIR(J)
  231. C              ELEM=element auquel appartient le noeud I.
  232.  
  233.                IPT1=MELEME
  234.                IF (LISOUS(/1).GT.0) IPT1=LISOUS(LMEMOIR(J))
  235.                itype = abs(ipt1.itypel)
  236. *  si element de type 49 ou 22, on ne connecte pas 2 noeuds non LX
  237.                connec=.true.
  238.                if (itype.eq.49) then
  239.                 do k=3,ipt1.num(/1)
  240.                   if (i.eq.icpr(ipt1.num(k,elem))) connec=.false.
  241.                 enddo
  242.                endif
  243.                DO 430 K=1,IPT1.NUM(/1)
  244. C              k representatif du nb de noeuds par elements.
  245.                   IK=ICPR(IPT1.NUM(K,ELEM))
  246.                   if (k.gt.2.and..not.connec) goto 430
  247.                   IF ((I.NE.IK).AND.
  248.      &                 (.NOT.(BOOL(IK)))) THEN
  249. C                 si i n'est pas  egal a un des nouveaux numeros des noeuds
  250. C                 de l'element ELEM et si il n'appartient pas deja  a l'ens des
  251. C                 voisins du noeud i(jadjc(i)),alors on le rajoute.
  252. C                   JADJC(IADJ(I+1))=IK
  253.                     JADJC(**)=IK
  254.                     IADJ(I+1)=IADJ(I+1)+1
  255.                     BOOL(IK)=.TRUE.
  256.                   ENDIF
  257. 430            CONTINUE
  258. 420        CONTINUE
  259. *  remise a faux de bool
  260.            DO 412 J=IADJ(I),IADJ(I+1)-1
  261.              IK=JADJC(J)
  262.              BOOL(IK)=.FALSE.
  263. 412        CONTINUE
  264. 400     CONTINUE
  265.  
  266.         SEGSUP      JMEMN,IMEMOIR,LMEMOIR,BOOLEEN
  267.  
  268.  
  269.  
  270. C**************************************************************************
  271.  
  272.  
  273. C       affectation
  274. C************************
  275.  
  276.  
  277.         if (nbthrs.gt.1) call threadii
  278.         SEGINI IPOS,MASQUE
  279.         IPOSMAX=1
  280.  
  281. **      write (6,*) ' nodes ',nodes
  282.         DO 50 I=1,NODES
  283.            MASQ(I)=.TRUE.
  284.            IPOS(I)=0
  285.            IPOS(NODES+I)=0
  286.            IPOS(2*NODES+I)=0
  287. 50      CONTINUE
  288. C       initialement, les noeuds ne sont pas masques,ont donc
  289. C       une position nulle.
  290.  
  291.         DIM=0
  292. C       le nombre de noeuds renumerotes  DIM est initialement egal a zero.
  293. C on initialise un premier separateur avec lisous(1) et lisous(2)
  294. C qui contiennent les noeuds maitres du super element (si appele par assem4)
  295.         mdomn=0
  296.         iposv=ipos(mdomn+1)
  297.         iposmax=iposmax+3
  298.         ipos(iposmax+1-2)=iposv+1
  299.         ipos(iposmax+1-1)=iposv+2
  300.         ipos(iposmax+1-0)=iposv+3
  301.         dimsep=0
  302.         do io=1,2
  303.          ipt1=lisous(io)
  304.           if (ipt1.itypel.eq.0 ) then
  305.           do j=1,ipt1.num(/2)
  306.           do i=1,ipt1.num(/1)
  307.            ip=icpr(ipt1.num(i,j))
  308.            if (masq(ip)) then
  309.              masq(ip)=.false.
  310.              ipos(ip+nodes)=iposmax-2
  311.              dimsep=dimsep+1
  312.            endif
  313.           enddo
  314.           enddo
  315.           endif
  316.         enddo
  317.         do i=1,nodes
  318.           if (masq(i)) then
  319.             ipos(i+nodes)=iposmax
  320.             ipos(i+2*nodes)=ipos(i+2*nodes)+1
  321.           endif
  322.         enddo
  323. **      write(6,*) ' dimsep dans numop2 ',dimsep
  324.         dim=dim+dimsep
  325. C               ****************************************
  326. C                       boucle principale
  327.         NS=NODES
  328.         ns=ns-dimsep
  329.  
  330.         nbthr=min(64,nbthrs)
  331.         SEGINI NRELONG,NOELON,noel2,londim
  332. **      write (6,*) ' avant appel sepa2 '
  333.         DO 500 I=1,NODES
  334. 550       IF(.NOT.MASQ(I)) GOTO 500
  335. C         si le noeud est masque alors ne rien faire: il est deja
  336. C         renumerote. On passe au noeud suivant.
  337.  
  338.           PIVOT=I
  339.  
  340.           CALL SEPA2(IADJ,JADJC,PIVOT,MASQUE,DIMSEP,NS,
  341.      >      IPOS,NODES,IPOSMAX,nrelong,noelon,noel2,
  342.      >      londim,nbthr,izero)
  343.           if (ierr.ne.0) then
  344.             if (nbthrs.gt.1) call threadis
  345.             return
  346.           endif
  347. C         separe le domaine d'etude en 2 parties.
  348. C         on decrit le domaine d'etude a partir du pivot et on cherche la
  349. C         longueur maximale en decrivant les voisins de pivot, et leurs
  350. C         voisins... jusqu'a rencontrer un voisin masque. On cree alors
  351. C         une nouvelle separation.
  352. C         les noeuds masques delimitent la separation du domaine.
  353.  
  354.  
  355.           DIM=DIM+DIMSEP
  356.           NS=NS-DIMSEP
  357. C         la dimension de noeuds renumerotes est augmente de DIMSEP.
  358. C         Celle de noeuds a renumeroter est diminue de DIMSEP.
  359.  
  360. *         IF (DIM.GE.NODES) GOTO 600
  361. C         si tous les noeuds ont ete renumerotes, on arrete.
  362.  
  363.           GOTO 550
  364.  
  365. 500     CONTINUE
  366. **      write (6,*) ' apres appel sepa2 '
  367.  
  368.         SEGSUP NRELONG,NOELON,noel2,londim,jmem
  369.  
  370. 600     CONTINUE
  371.         if (nbthrs.gt.1) call threadis
  372. *
  373. *  tri dans chaque zone
  374. *  je ne sais pas trop pourquoi ca marche
  375.          iposmx=0
  376.          do 610 lpoint=1,nodes
  377.           mdomn=ipos(lpoint+nodes)
  378.           iposv=ipos(mdomn+1)
  379.           iposi=0
  380.          do 620 kk=iadj(lpoint),iadj(lpoint+1)-1
  381.           k=jadjc(kk)
  382.           if(k.eq.lpoint) goto 620
  383.           iposk=ipos(ipos(k+nodes)+1)
  384.            if (iposk.ne.iposv) then
  385.              iposi=max(iposi,iposk)
  386.            endif
  387.  620     continue
  388.          ipos(lpoint+2*nodes)=5*(lpoint+iposi*nodes)
  389.          iposmx=max(iposmx,ipos(lpoint+2*nodes))
  390.  610     continue
  391.          iposmx=iposmx+5
  392. *
  393. *  mise a la bonne place des multiplicateurs de Lagrange
  394.         do 700 J=0,NB1-1
  395.           jp1=j+1
  396.           iposvs=igrand
  397.           iposvr=-igrand
  398.           mdomnr=0
  399.           mdomns=0
  400.           ipaur=0
  401.           ipaus=0
  402. *         write (6,*) 'numop2 ',(lagran(J*NB2+il),il=1,nb2)
  403.         do 800 il=3,nb2
  404.           ip = abs(LAGRAN(J*NB2+il))
  405. *         if (ip.eq.0) write (6,*) ' prob numop2 '
  406.           if (ip.eq.0) goto 710
  407.           mdomn=ipos(ip+nodes)
  408.           iposv=ipos(mdomn+1)
  409. *  deplacer les noeuds en relation en fin de zone
  410.           ipos(ip+2*nodes)=-abs(ipos(ip+2*nodes))
  411.           ipos(ip+2*nodes)=mod(ipos(ip+2*nodes),iposmx)-iposmx
  412.           if (iposv.gt.iposvr) then
  413.               iposvr=iposv
  414.               mdomnr=mdomn
  415.               ipaur=ipos(ip+2*nodes)
  416.           elseif (iposv.eq.iposvr) then
  417.               ipaur=max(ipaur,ipos(ip+2*nodes))
  418.           endif
  419.           if (iposv.lt.iposvs) then
  420.               iposvs=iposv
  421.               mdomns=mdomn
  422.               ipaus=ipos(ip+2*nodes)
  423.           elseif (iposv.eq.iposvs) then
  424.                 ipaus=min(ipaus,ipos(ip+2*nodes))
  425.           endif
  426.  800    continue
  427.  710    continue
  428. *
  429. *  le premier mult avant le premier noeud
  430.         ip=abs(LAGRAN(J*NB2+1))
  431.         IPOS(IP+2*NODES)= ipaur+1
  432. * numeroter le frottement apres le contact
  433.          if (LAGRAN(J*NB2+1).lt.0) 
  434.      >       IPOS(IP+2*NODES)=ipaur+2
  435.         IPOS(IP+nodes)=mdomnr
  436. **      write (6,*) 'premier  mult ',ip,ipos(ip+nodes),ipos(ip+2*nodes)
  437. *
  438. *  le deuxieme mult apres le dernier noeud
  439.         ip=abs(LAGRAN(J*NB2+2))
  440.         IPOS(IP+2*NODES)= ipaus-1
  441. * numeroter le frottement avant le contact
  442.          if (LAGRAN(J*NB2+2).lt.0) 
  443.      >       IPOS(IP+2*NODES)=ipaus-2
  444. **      IPOS(IP+nodes)=mdomns
  445. * pv deuxieme mult en bout de maillage
  446.         IPOS(IP+nodes)=0                
  447. **      write (6,*) 'deuxieme mult ',ip,ipos(ip+nodes),ipos(ip+2*nodes)
  448. *
  449.  700    continue
  450.         SEGSUP IADJ,JADJC,MASQUE
  451. *  ok maintenant on trie
  452.         CALL SORTI2(IPOS,JNT,NODES)
  453.  
  454. C***************************************************************************
  455.  
  456.         DO 860 I=1,ICPR(/1)
  457.            IF(ICPR(I).NE.0) ICPR(I)=JNT(ICPR(I))
  458. C          numerotation finale.
  459. 860     CONTINUE
  460.  
  461.  
  462.         SEGSUP JNT,IPOS,LAGRAN
  463.  
  464.  
  465.         RETURN
  466.         END
  467.  
  468.  
  469.  
  470.  
  471.  
  472.  
  473.  
  474.  
  475.  
  476.  
  477.  
  478.  
  479.  
  480.  
  481.  
  482.  
  483.  
  484.  
  485.  
  486.  
  487.  
  488.  
  489.  
  490.  
  491.  
  492.  
  493.  
  494.  
  495.  
  496.  
  497.  
  498.  
  499.  
  500.  
  501.  
  502.  
  503.  
  504.  
  505.  
  506.  
  507.  
  508.  
  509.  

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