Télécharger la plaquette Cast3M


fr en

Télécharger numop3.eso

Retour à la liste

Numérotation des lignes : 

numop3

  1. C NUMOP3    SOURCE    GOUNAND   25/05/15    21:15:09     12268          
  2.       SUBROUTINE NUMOP3(MELEME,ICPR,NODES)
  3.       IMPLICIT REAL*8 (A-H,O-Z)
  4.       IMPLICIT INTEGER (I-N)
  5. C***********************************************************************
  6. C NOM         : NUMOP3
  7. C DESCRIPTION : Numerotation des noeuds d'un maillage
  8. C               par la methode de SLOAN
  9. C               Reprise de NUMOP2 pour le graphe d'adjacence
  10. C               et le placement des multiplicateurs
  11. C
  12. C
  13. C LANGAGE     : ESOPE
  14. C AUTEUR      : Stephane GOUNAND (CEA/DES/ISAS/DM2S/SEMT/LTA)
  15. C               mel : gounand@semt2.smts.cea.fr
  16. C***********************************************************************
  17. C APPELES          :
  18. C APPELES (E/S)    :
  19. C APPELES (BLAS)   :
  20. C APPELES (CALCUL) :
  21. C APPELE PAR       :
  22. C***********************************************************************
  23. C SYNTAXE GIBIANE    :
  24. C ENTREES            :
  25. C ENTREES/SORTIES    :
  26. C SORTIES            :
  27. C***********************************************************************
  28. C VERSION    : v1, 06/05/2025, version initiale
  29. C HISTORIQUE : v1, 06/05/2025, creation
  30. C HISTORIQUE :
  31. C HISTORIQUE :
  32. C***********************************************************************
  33. -INC PPARAM
  34. -INC CCOPTIO
  35. -INC CCREEL
  36. -INC SMELEME
  37. -INC SMCOORD
  38.  
  39.  
  40.       LOGICAL CONNEC
  41.       SEGMENT JMEM(NODES+1),JMEMN(NODES+1)
  42. C       JMEM et JMEMN contiennent le nombre d'element auquel appartient un noeud
  43.  
  44.       SEGMENT JNT(NODES)
  45. C       JNT contient la nouvelle numerotation
  46.       SEGMENT IJNT(NODES)
  47. C     IJNT contient l'inverse de la nouvelle numerotation
  48.       SEGMENT JORDRE(NODES)
  49. C     Ordre des noeuds
  50.       SEGMENT IWORK(3*NODES+1)
  51. C       Tableau de travail pour prsloa
  52.       SEGMENT IJNT2(NODES)
  53.       SEGMENT JOR2(NODES)
  54. C       Tableau de travail pour trifus
  55.  
  56.       SEGMENT ICPR(nbpts)
  57. C       ICPR au debut contient l'ancienne numerotation ,
  58. C                               a la fin la nouvelle.
  59.  
  60.       SEGMENT IADJ(NODES+1)
  61.       SEGMENT JADJC(0)
  62. C       IADJ(i) pointe sur JADJC qui contient les voisins de i entre
  63. C       IADJ(i) et IADJ(i+1)-1
  64.  
  65.       SEGMENT LAGRAN(NB)
  66. C       contient les noeud de lagrange et les noeuds les suivant directement
  67. C       cf element de type 49
  68.  
  69.       SEGMENT BOOLEEN
  70.          LOGICAL BOOL(NODES)
  71.       ENDSEGMENT
  72. C       BOOL(i) = true si le noeud i a ete deja mentionne dans la liste
  73. C       des voisins JADJC.
  74.  
  75.       SEGMENT IMEMOIR(NBV),LMEMOIR(NBV)
  76. C       contient les elements appartenant a chaque noeud,sous forme de liste.
  77.  
  78.       INTEGER ELEM
  79. C       nom d'un element
  80.       LOGICAL LDBNUM,LVERIF
  81. *
  82. * Executable statements
  83. *
  84.       LDBNUM=.FALSE.
  85.       LVERIF=.FALSE.
  86.       SEGACT ICPR*MOD
  87.       NODES=ICPR(/1)
  88.       SEGACT MELEME
  89. C       icpr: numero des noeuds.
  90. C       meleme: objet de maillage (cf assem2.eso)
  91.  
  92.       DO 10 I=1,ICPR(/1)
  93.          ICPR(I)=0
  94.  10   CONTINUE
  95.  
  96.       IPT1=MELEME
  97.       IKOU=0
  98.       NBV=0
  99.       NB1=0
  100.       NB2=0
  101.  
  102.       DO 100 IO=1,MAX(1,LISOUS(/1))
  103.          IF (LISOUS(/1).GT.0) THEN
  104.             IPT1=LISOUS(IO)
  105.             SEGACT IPT1
  106.          ENDIF
  107. C          on cree la numerotation des noeuds.
  108. C          'nb noeuds/element'=IPT1.NUM(/1)
  109. C          'nb element'=IPT1.NUM(/2)
  110.          IF(abs(IPT1.ITYPEL).EQ.49) THEN
  111.             NB1=NB1+IPT1.NUM(/2)
  112.             NB2=MAX(NB2,IPT1.NUM(/1))
  113. C            NB1= nbre d'éléments de type 49.
  114. C            NB2=nbre de noeuds/élément maximum parmi
  115. C                 les éléments de type 22.
  116.          ENDIF
  117.          DO J=1,IPT1.NUM(/2)
  118.             DO I=1,IPT1.NUM(/1)
  119.                IJ=IPT1.NUM(I,J)
  120. C              IJ est le Ième noeud du Jème élément.
  121.                IF (ICPR(IJ).EQ.0) THEN
  122. C              s'il est déjà numéroté, on ne fait rien.
  123.                   IKOU=IKOU+1
  124.                   ICPR(IJ)=IKOU
  125.                ENDIF
  126.             enddo
  127.          enddo
  128.  100  CONTINUE
  129.  
  130.       NODES=IKOU
  131.       NB=NB2*NB1
  132.  
  133. C*****  initalisation des segments*********
  134.  
  135.       SEGINI,LAGRAN
  136.       SEGINI IADJ,JADJC
  137.       SEGINI,JMEM,JMEMN
  138.       SEGINI BOOLEEN
  139.  
  140. C******************************************
  141.  
  142.       IPT1=MELEME
  143.       IADJ(1)=1
  144.       INC=0
  145.       DO 200 IO=1,MAX(1,LISOUS(/1))
  146.          IF (LISOUS(/1).GT.0) IPT1=LISOUS(IO)
  147.          DO 210 J=1,IPT1.NUM(/2)
  148.             IF(abs(IPT1.ITYPEL).EQ.49) THEN
  149.                is=sign(1,ipt1.itypel)
  150.                DO 220 I=1,IPT1.NUM(/1)
  151. C               chaque element de type 49 a  au plus NB2 noeuds.
  152.                   LAGRAN(INC*NB2+I)=ICPR(IPT1.NUM(I,J))*is
  153. C                  les noeuds de l'elements de type 49
  154. C                  sont ranges dans le vecteur LAGRAN.
  155.  220           CONTINUE
  156.                DO 225 I=IPT1.NUM(/1)+1,NB2
  157.                   LAGRAN(INC*NB2+I)=0
  158. C                 comme on a alloue la place memoire maximale,
  159. C                 on remplit les cases restantes avec des 0.
  160.  225           CONTINUE
  161.                INC=INC+1
  162. C               INC=le nbre d'elements de type 49.
  163.             ENDIF
  164.             DO 230 I=1,IPT1.NUM(/1)
  165.                IJ=ICPR(IPT1.NUM(I,J))+1
  166.                JMEM(IJ)=JMEM(IJ)+1
  167. C              JMEM(I+1): nb elements auquel le noeud I appartient
  168.  230        CONTINUE
  169.  210     CONTINUE
  170.  200  CONTINUE
  171.  
  172.  
  173.       JMEM(1)=1
  174.       DO 30 I=1,NODES
  175.          JMEM(I+1)=JMEM(I)+JMEM(I+1)
  176. C           JMEM(I+1)=indice de depart des elements
  177. C           auxquels le noeud I appartient.
  178.  30   CONTINUE
  179.       NBV=JMEM(NODES+1)
  180. C       NBV= dimension de IMEMOIR.
  181.       SEGINI IMEMOIR,LMEMOIR
  182.       IPT1=MELEME
  183.       DO 300 IO=1,MAX(1,LISOUS(/1))
  184.          IF (LISOUS(/1).GT.0) THEN
  185.             IPT1=LISOUS(IO)
  186.          ENDIF
  187.          DO J=1,IPT1.NUM(/2)
  188.             DO I=1,IPT1.NUM(/1)
  189.                IJ=ICPR(IPT1.NUM(I,J))
  190.                JMEMN(IJ+1)=JMEMN(IJ+1)+1
  191.                IMEMOIR(JMEM(IJ)+JMEMN(IJ+1)-1)=J
  192.                LMEMOIR(JMEM(IJ)+JMEMN(IJ+1)-1)=IO
  193. C             on range dans IMEMOIR tous les elements des sous-objets
  194. C             IO auxquels appartient le noeud ICPR(IPT1.NUM(I,J)).
  195. C             On connait pour chaque element, le sous-objet auquel
  196. C             il appartient grace a LMEMOIR
  197.             enddo
  198.          enddo
  199.  300  CONTINUE
  200.       DO 410 J=1,NODES
  201.          BOOL(J)=.FALSE.
  202.  410  CONTINUE
  203.       DO 400 I=1,NODES
  204.          IADJ(I+1)=IADJ(I)
  205.          DO 420 J=JMEM(I),JMEM(I+1)-1
  206.             ELEM=IMEMOIR(J)
  207. C              ELEM=element auquel appartient le noeud I.
  208.             IPT1=MELEME
  209.             IF (LISOUS(/1).GT.0) IPT1=LISOUS(LMEMOIR(J))
  210.             itype = abs(ipt1.itypel)
  211. *  si element de type 49 ou 22, on ne connecte pas 2 noeuds non LX
  212.             connec=.true.
  213.             if (itype.eq.49) then
  214.                do k=3,ipt1.num(/1)
  215.                   if (i.eq.icpr(ipt1.num(k,elem))) connec=.false.
  216.                enddo
  217.             endif
  218.             DO 430 K=1,IPT1.NUM(/1)
  219. C              k representatif du nb de noeuds par elements.
  220.                IK=ICPR(IPT1.NUM(K,ELEM))
  221.                if (k.gt.2.and..not.connec) goto 430
  222.                IF ((I.NE.IK).AND.
  223.      &              (.NOT.(BOOL(IK)))) THEN
  224. C                 si i n'est pas  egal a un des nouveaux numeros des noeuds
  225. C                 de l'element ELEM et si il n'appartient pas deja  a l'ens des
  226. C                 voisins du noeud i(jadjc(i)),alors on le rajoute.
  227. C                   JADJC(IADJ(I+1))=IK
  228.                   JADJC(**)=IK
  229.                   IADJ(I+1)=IADJ(I+1)+1
  230.                   BOOL(IK)=.TRUE.
  231.                ENDIF
  232.  430        CONTINUE
  233.  420     CONTINUE
  234. **  pour ne pas avoir un tableau de longueur nulle, ce que esope n'aime pas
  235. *         if(jadjc(/1).eq.0) jadjc(**)=0
  236. *  remise a faux de bool
  237.          DO 412 J=IADJ(I),IADJ(I+1)-1
  238.             IK=JADJC(J)
  239.             BOOL(IK)=.FALSE.
  240.  412     CONTINUE
  241.  400  CONTINUE
  242.       SEGSUP JMEM,JMEMN
  243.       SEGSUP IMEMOIR,LMEMOIR
  244.       SEGSUP BOOLEEN
  245. *
  246. * A ce stade, le graphe d'adjacence est près
  247. *
  248.       if (ldbnum) Write(ioimp,*) 'Graphe adjacence nodes=',nodes
  249.       if (ldbnum) segprt,IADJ
  250.       if (ldbnum) segprt,JADJC
  251.       SEGINI JNT
  252.       SEGINI IWORK
  253.       E2=JADJC(/1)
  254.       CALL LABEL(NODES,E2,JADJC(1),IADJ(1),
  255.      $     JNT(1),
  256.      $     IWORK(1),
  257.      $     IMPR,IRET)
  258.       IF (IRET.NE.0) GOTO 9999
  259.       SEGSUP IWORK
  260.       SEGSUP IADJ,JADJC
  261.       if (ldbnum) SEGPRT,JNT
  262. *
  263. * Maintenant, il faut mettre les Lagrange a la bonne place !!!!
  264. * On fait comme dans KRES24
  265. *
  266.       IF (NB1.EQ.0) GOTO 2000
  267.       SEGINI IJNT
  268.       DO I=1,NODES
  269.          IJNT(JNT(I))=I
  270.       ENDDO
  271. *     Il y a trois types différents de -1 à 1 : il faut pouvoir placer
  272. *     les multiplicateurs de lagrange entre les autres noeuds
  273.       NTYP=3
  274.       SEGINI JORDRE
  275.       DO I=1,NODES
  276.          JORDRE(I)=I*NTYP
  277.       ENDDO
  278.       JORMAX= (NODES+1)*NTYP
  279.  
  280.       if (ldbnum) write(ioimp,*) 'Avant mise a la bonne place'
  281.       if (ldbnum) segprt,jordre
  282.  
  283.       do 700 J=0,NB1-1
  284.          ipaur=-igrand
  285.          ipaus=igrand
  286.          do 800 il=3,nb2
  287.             ip = abs(LAGRAN(J*NB2+il))
  288. *            if (ip.eq.0) write (6,*) ' prob numop3 '
  289.             if (ip.eq.0) goto 710
  290.             jddln=JNT(IP)
  291.             jordre(jddln)=-abs(jordre(jddln))
  292.             ipaur=max(ipaur,jordre(jddln))
  293.             ipaus=min(ipaus,jordre(jddln))
  294.  800     continue
  295.  710     continue
  296. *
  297. *     On va laisser comme ça
  298. *     Ce cas peut arriver après élimination
  299. *     Cela devrait revenir à placer les multiplicateurs en fin ou début
  300. *     de matrice
  301.          if (ipaur.eq.-igrand.or.ipaus.eq.igrand) then
  302. *            Write(ioimp,*) 'mulag sans relations pas ok'
  303. *            goto 9999
  304.             ipaur=0
  305.             ipaus=-jormax
  306.          endif
  307.          if (ldbnum) write(ioimp,*) 'j,ipaur,ipaus=',j,ipaur
  308.      $        ,ipaus
  309.  
  310. *
  311. *  le premier mult avant le premier noeud
  312.          ip=abs(LAGRAN(J*NB2+1))
  313.          iddln=JNT(IP)
  314.          jordre(iddln)=ipaur+1
  315. *
  316. *  le deuxieme mult apres le dernier noeud
  317.          ip=abs(LAGRAN(J*NB2+2))
  318.          iddln=JNT(IP)
  319.          jordre(iddln)=ipaus-1
  320. *
  321.  700  continue
  322. *     WRITE(IOIMP,*) 'Avant chgt signe'
  323. *     segprt,jordre
  324.       DO I=1,nodes
  325.          JORDRE(I)=-JORDRE(I)
  326.       ENDDO
  327. *     Avant tri
  328.       if (ldbnum) WRITE(IOIMP,*) 'Avant TRIFUS'
  329.       if (ldbnum) SEGPRT,IJNT
  330.       if (ldbnum) SEGPRT,JORDRE
  331.       SEGINI JOR2
  332.       SEGINI IJNT2
  333. *     ok maintenant on trie
  334.       CALL TRIFUS(nodes,JORDRE(1),IJNT(1),JOR2(1),IJNT2(1))
  335.       SEGSUP IJNT2
  336.       SEGSUP JOR2
  337. *     Apres tri
  338.       if (ldbnum) WRITE(IOIMP,*) 'Apres TRIFUS'
  339.       if (ldbnum) SEGPRT,IJNT
  340.       if (ldbnum) SEGPRT,JORDRE
  341. *     permutation inverse
  342.       DO I=1,nodes
  343.          JNT(IJNT(I))=I
  344.       ENDDO
  345. *     Verification que dans la nouvelle numerotation les multiplicateurs
  346. *     sont correctement places...
  347.       IF (LVERIF) THEN
  348.          write(ioimp,*) 'VERIF NUMOP3'
  349.          do 1700 J=0,NB1-1
  350.             ipaur=-igrand
  351.             ipaus=igrand
  352.             do 1800  il=3,nb2
  353.                ip = abs(LAGRAN(J*NB2+il))
  354. *            if (ip.eq.0) write (6,*) ' prob numop3 '
  355.                if (ip.eq.0) goto 1710
  356.                jddln=JNT(IP)
  357.                ipaur=max(ipaur,jddln)
  358.                ipaus=min(ipaus,jddln)
  359.  1800       continue
  360.  1710       continue
  361.             if (ipaur.eq.-igrand.or.ipaus.eq.igrand) then
  362.                goto 1700
  363. *     Write(ioimp,*) 'mulag sans relations pas ok'
  364. *     goto 9999
  365.             endif
  366.             ip=abs(LAGRAN(J*NB2+1))
  367.             iddln=JNT(IP)
  368.             if (ipaus.le.iddln) then
  369.                write(ioimp,*) 'Erreur numerotation'
  370.                write(ioimp,*) 'ip,iddln,ipaus=',ip,iddln,ipaus
  371.                goto 9999
  372.             endif
  373.             ip=abs(LAGRAN(J*NB2+2))
  374.             iddln=JNT(IP)
  375.             if (ipaur.ge.iddln) then
  376.                write(ioimp,*) 'Erreur numerotation'
  377.                write(ioimp,*) 'ip,iddln,ipaur=',ip,iddln,ipaur
  378.                goto 9999
  379.             endif
  380.  1700    continue
  381.       ENDIF
  382.  
  383. C***************************************************************************
  384.       SEGSUP JORDRE
  385.       SEGSUP IJNT
  386.  2000 CONTINUE
  387.  
  388.       DO 860 I=1,ICPR(/1)
  389.          IF(ICPR(I).NE.0) ICPR(I)=JNT(ICPR(I))
  390. C          numerotation finale.
  391.  860  CONTINUE
  392.       SEGSUP JNT
  393.       SEGSUP LAGRAN
  394. *
  395. * Normal termination
  396. *
  397.       RETURN
  398. *
  399. * Format handling
  400. *
  401. *
  402. * Error handling
  403. *
  404.  9999 CONTINUE
  405.       WRITE(IOIMP,*) 'An error was detected in subroutine numop3'
  406.       RETURN
  407. *
  408. * End of subroutine NUMOP3
  409. *
  410.       END
  411.  
  412.  

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