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/02/17    21:15:09     12152          
  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+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.            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. *  pour ne pas avoir un tableau de longueur nulle, ce que esope n'aime pas
  260.            if(jadjc(/1).eq.0) jadjc(**)=0
  261. *  remise a faux de bool
  262.            DO 412 J=IADJ(I),IADJ(I+1)-1
  263.              IK=JADJC(J)
  264.              BOOL(IK)=.FALSE.
  265. 412        CONTINUE
  266. 400     CONTINUE
  267.  
  268.         SEGSUP      JMEMN,IMEMOIR,LMEMOIR,BOOLEEN
  269.  
  270.  
  271.  
  272. C**************************************************************************
  273.  
  274.  
  275. C       affectation
  276. C************************
  277.  
  278.  
  279.         if (nbthrs.gt.1) call threadii
  280.         SEGINI IPOS,MASQUE
  281.         IPOSMAX=1
  282.  
  283. **      write (6,*) ' nodes ',nodes
  284.         DO 50 I=1,NODES
  285.            MASQ(I)=.TRUE.
  286.            IPOS(I)=0
  287.            IPOS(NODES+I)=0
  288.            IPOS(2*NODES+I)=0
  289. 50      CONTINUE
  290. C       initialement, les noeuds ne sont pas masques,ont donc
  291. C       une position nulle.
  292.  
  293.         DIM=0
  294. C       le nombre de noeuds renumerotes  DIM est initialement egal a zero.
  295. C on initialise un premier separateur avec lisous(1) et lisous(2)
  296. C qui contiennent les noeuds maitres du super element (si appele par assem4)
  297.         mdomn=0
  298.         iposv=ipos(mdomn+1)
  299.         iposmax=iposmax+3
  300.         ipos(iposmax+1-2)=iposv+1
  301.         ipos(iposmax+1-1)=iposv+2
  302.         ipos(iposmax+1-0)=iposv+3
  303.         dimsep=0
  304.         do io=1,2
  305.          ipt1=lisous(io)
  306.           if (ipt1.itypel.eq.0 ) then
  307.           do j=1,ipt1.num(/2)
  308.           do i=1,ipt1.num(/1)
  309.            ip=icpr(ipt1.num(i,j))
  310.            if (masq(ip)) then
  311.              masq(ip)=.false.
  312.              ipos(ip+nodes)=iposmax-2
  313.              dimsep=dimsep+1
  314.            endif
  315.           enddo
  316.           enddo
  317.           endif
  318.         enddo
  319.         do i=1,nodes
  320.           if (masq(i)) then
  321.             ipos(i+nodes)=iposmax
  322.             ipos(i+2*nodes)=ipos(i+2*nodes)+1
  323.           endif
  324.         enddo
  325. **      write(6,*) ' dimsep dans numop2 ',dimsep
  326.         dim=dim+dimsep
  327. C               ****************************************
  328. C                       boucle principale
  329.         NS=NODES
  330.         ns=ns-dimsep
  331.  
  332.         nbthr=min(64,nbthrs)
  333.         SEGINI NRELONG,NOELON,noel2,londim
  334. **      write (6,*) ' avant appel sepa2 '
  335.         DO 500 I=1,NODES
  336. 550       IF(.NOT.MASQ(I)) GOTO 500
  337. C         si le noeud est masque alors ne rien faire: il est deja
  338. C         renumerote. On passe au noeud suivant.
  339.  
  340.           PIVOT=I
  341.  
  342.           CALL SEPA2(IADJ,JADJC,PIVOT,MASQUE,DIMSEP,NS,
  343.      >      IPOS,NODES,IPOSMAX,nrelong,noelon,noel2,
  344.      >      londim,nbthr,izero)
  345.           if (ierr.ne.0) then
  346.             if (nbthrs.gt.1) call threadis
  347.             return
  348.           endif
  349. C         separe le domaine d'etude en 2 parties.
  350. C         on decrit le domaine d'etude a partir du pivot et on cherche la
  351. C         longueur maximale en decrivant les voisins de pivot, et leurs
  352. C         voisins... jusqu'a rencontrer un voisin masque. On cree alors
  353. C         une nouvelle separation.
  354. C         les noeuds masques delimitent la separation du domaine.
  355.  
  356.  
  357.           DIM=DIM+DIMSEP
  358.           NS=NS-DIMSEP
  359. C         la dimension de noeuds renumerotes est augmente de DIMSEP.
  360. C         Celle de noeuds a renumeroter est diminue de DIMSEP.
  361.  
  362. *         IF (DIM.GE.NODES) GOTO 600
  363. C         si tous les noeuds ont ete renumerotes, on arrete.
  364.  
  365.           GOTO 550
  366.  
  367. 500     CONTINUE
  368. **      write (6,*) ' apres appel sepa2 '
  369.  
  370.         SEGSUP NRELONG,NOELON,noel2,londim,jmem
  371.  
  372. 600     CONTINUE
  373.         if (nbthrs.gt.1) call threadis
  374. *
  375. *  tri dans chaque zone
  376. *  je ne sais pas trop pourquoi ca marche
  377.          iposmx=0
  378.          do 610 lpoint=1,nodes
  379.           mdomn=ipos(lpoint+nodes)
  380.           iposv=ipos(mdomn+1)
  381.           iposi=0
  382.          do 620 kk=iadj(lpoint),iadj(lpoint+1)-1
  383.           k=jadjc(kk)
  384.           if(k.eq.lpoint) goto 620
  385.           iposk=ipos(ipos(k+nodes)+1)
  386.            if (iposk.ne.iposv) then
  387.              iposi=max(iposi,iposk)
  388.            endif
  389.  620     continue
  390.          ipos(lpoint+2*nodes)=5*(lpoint+iposi*nodes)
  391.          iposmx=max(iposmx,ipos(lpoint+2*nodes))
  392.  610     continue
  393.          iposmx=iposmx+5
  394. *
  395. *  mise a la bonne place des multiplicateurs de Lagrange
  396.         do 700 J=0,NB1-1
  397.           jp1=j+1
  398.           iposvs=igrand
  399.           iposvr=-igrand
  400.           mdomnr=0
  401.           mdomns=0
  402.           ipaur=0
  403.           ipaus=0
  404. *         write (6,*) 'numop2 ',(lagran(J*NB2+il),il=1,nb2)
  405.         do 800 il=3,nb2
  406.           ip = abs(LAGRAN(J*NB2+il))
  407. *         if (ip.eq.0) write (6,*) ' prob numop2 '
  408.           if (ip.eq.0) goto 710
  409.           mdomn=ipos(ip+nodes)
  410.           iposv=ipos(mdomn+1)
  411. *  deplacer les noeuds en relation en fin de zone
  412.           ipos(ip+2*nodes)=-abs(ipos(ip+2*nodes))
  413.           ipos(ip+2*nodes)=mod(ipos(ip+2*nodes),iposmx)-iposmx
  414.           if (iposv.gt.iposvr) then
  415.               iposvr=iposv
  416.               mdomnr=mdomn
  417.               ipaur=ipos(ip+2*nodes)
  418.           elseif (iposv.eq.iposvr) then
  419.               ipaur=max(ipaur,ipos(ip+2*nodes))
  420.           endif
  421.           if (iposv.lt.iposvs) then
  422.               iposvs=iposv
  423.               mdomns=mdomn
  424.               ipaus=ipos(ip+2*nodes)
  425.           elseif (iposv.eq.iposvs) then
  426.                 ipaus=min(ipaus,ipos(ip+2*nodes))
  427.           endif
  428.  800    continue
  429.  710    continue
  430. *
  431. *  le premier mult avant le premier noeud
  432.         ip=abs(LAGRAN(J*NB2+1))
  433.         IPOS(IP+2*NODES)= ipaur+1
  434. * numeroter le frottement apres le contact
  435.          if (LAGRAN(J*NB2+1).lt.0) 
  436.      >       IPOS(IP+2*NODES)=ipaur+2
  437.         IPOS(IP+nodes)=mdomnr
  438. **      write (6,*) 'premier  mult ',ip,ipos(ip+nodes),ipos(ip+2*nodes)
  439. *
  440. *  le deuxieme mult apres le dernier noeud
  441.         ip=abs(LAGRAN(J*NB2+2))
  442.         IPOS(IP+2*NODES)= ipaus-1
  443. * numeroter le frottement avant le contact
  444.          if (LAGRAN(J*NB2+2).lt.0) 
  445.      >       IPOS(IP+2*NODES)=ipaus-2
  446. **      IPOS(IP+nodes)=mdomns
  447. * pv deuxieme mult en bout de maillage
  448.         IPOS(IP+nodes)=0                
  449. **      write (6,*) 'deuxieme mult ',ip,ipos(ip+nodes),ipos(ip+2*nodes)
  450. *
  451.  700    continue
  452.         SEGSUP IADJ,JADJC,MASQUE
  453. *  ok maintenant on trie
  454.         CALL SORTI2(IPOS,JNT,NODES)
  455.  
  456. C***************************************************************************
  457.  
  458.         DO 860 I=1,ICPR(/1)
  459.            IF(ICPR(I).NE.0) ICPR(I)=JNT(ICPR(I))
  460. C          numerotation finale.
  461. 860     CONTINUE
  462.  
  463.  
  464.         SEGSUP JNT,IPOS,LAGRAN
  465.  
  466.  
  467.         RETURN
  468.         END
  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.  
  510.  
  511.  
  512.  

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