Télécharger la plaquette Cast3M


fr en

Télécharger numop1.eso

Retour à la liste

Numérotation des lignes : 

numop1

  1. C NUMOP1    SOURCE    PV        20/03/30    21:21:40     10567          
  2. C   RACINE DE LA NUMEROTATION      POUR LA SORTIE SUR FAC
  3. C
  4. C       methode utilisee: NESTED DISSECTION & MULTIFRONTALE
  5.  
  6. C       ce programme decompose un domaine d'etude par elements finis
  7. C       et utilise ensuite la methode multifrontale.
  8.  
  9.  
  10.       SUBROUTINE NUMOP1(MELEME,ICPR,NODES)
  11.       IMPLICIT INTEGER(I-N)
  12. -INC SMELEME
  13. -INC SMCOORD
  14.  
  15. -INC PPARAM
  16. -INC CCOPTIO
  17.  
  18.         SEGMENT JMEM(NODES+1),JMEMN(NODES+1)
  19. C       JMEM et JMEMN contiennent le nombre d'elements auquel appartient
  20. C       un noeud
  21.  
  22.         SEGMENT JNT(NODES),NUMERO(NODES)
  23. C       JNT contient la nouvelle numerotation
  24.  
  25.         SEGMENT ICPR(nbpts)
  26. C       ICPR au debut contient l'ancienne numerotation,
  27. C       a la fin la nouvelle.
  28.  
  29.         SEGMENT IADJ(NODES+1),IVOIS(0)
  30.         SEGMENT NPVOIS(2*NBENS),NVOIS(0)
  31. C       IADJ(i) pointe sur IVOIS qui contient les voisins de i entre
  32. C       IADJ(i) et IADJ(i+1)-1
  33.  
  34.         SEGMENT LAGRAN(NB)
  35. C       contient les noeud de lagrange et les noeuds les suivant directement
  36. C       cf element de type 22
  37.         SEGMENT INVLAG(NB1)
  38.  
  39.         SEGMENT BOOLEEN
  40.           LOGICAL BOOL(NODES)
  41.         ENDSEGMENT
  42. C       BOOL(i) = true si le noeud i a etait deja memtioner dans la liste
  43. C       des voisin IVOIS
  44.  
  45.         SEGMENT IMEMOIR(NBV)
  46. C       contient les element appartenant a chaque noeud,sous forme de liste,
  47. C       au premier, puis au second, etc...
  48.  
  49.         INTEGER ELEM
  50. C       nom d'un element
  51.  
  52.         INTEGER N,NBENS
  53.  
  54.         SEGMENT MASQUE
  55.           LOGICAL MASQ(NODES)
  56.         ENDSEGMENT
  57. C       MASQ(X)=.TRUE. si le noeud X n'a pas ete renumerote;
  58. C               .FALSE. si il l'a ete.
  59.  
  60.         INTEGER DIM,DIMSEP
  61. C       DIM= nombre de noeuds renumerotes.
  62.  
  63.         INTEGER PIVOT
  64. C       PIVOT est le noeud utile a la division du domaine.
  65.  
  66.         SEGMENT IPOS(NODES*2)
  67. C       est le vecteur contenant la numerotation dans les deux sens,de 1 a NODES
  68. C       puis de NODES+1 a 2* NODES, cf la subroutine SEPAR
  69.  
  70.         SEGMENT ICOMPT(NODES),IORDRE(2*NBENS)
  71.  
  72.  
  73.  
  74. C       initialisation
  75. C*******************************
  76.  
  77.  
  78.  
  79.         IENORM=2000000000
  80. C       norme d'erreur.
  81.  
  82.         NODES=nbpts
  83.         SEGACT ICPR*MOD
  84.         SEGACT MELEME
  85. C       icpr: numero des noeuds.
  86. C       meleme: objet de maillage (cf assem2.eso)
  87.  
  88.         DO 10 I=1,ICPR(/1)
  89. 10        ICPR(I)=0
  90.  
  91.         IPT1=MELEME
  92.         IKOU=0
  93.         NBV=0
  94.         NB1=0
  95.         NB2=0
  96.  
  97.         DO 100 IO=1,MAX(1,LISOUS(/1))
  98.            IF (LISOUS(/1).NE.0) THEN
  99.               IPT1=LISOUS(IO)
  100.               SEGACT IPT1
  101.            ENDIF
  102. C          on cree la numerotation des noeuds.
  103. C          'nb noeuds/element'=IPT1.NUM(/1)
  104. C          'nb element'=IPT1.NUM(/2)
  105.            IF(IPT1.ITYPEL.EQ.22)  THEN
  106.              NB1=NB1+IPT1.NUM(/2)
  107.              NB2=MAX(NB2,IPT1.NUM(/1))
  108. C            NB1= nbre d'elements de type 22.
  109. C            NB2=nbre de noeuds/element maximum parmi
  110. C                 les elements de type 22.
  111.            ENDIF
  112.            DO 150 I=1,IPT1.NUM(/1)
  113.            DO 150 J=1,IPT1.NUM(/2)
  114.                IJ=IPT1.NUM(I,J)
  115. C              IJ est le Ieme noeud du Jeme element.
  116.                IF (ICPR(IJ).NE.0) GOTO 150
  117. C              s'il est deja numerote, on ne fait rien.
  118.                IKOU=IKOU+1
  119.                ICPR(IJ)=IKOU
  120. 150        CONTINUE
  121.         NODES=IKOU
  122. C       NODES= nbre de noeuds.
  123. 100     CONTINUE
  124.  
  125.         NB=NB2*NB1
  126. C       NB= dimension(LAGRAN).
  127.  
  128. C*****  initalisation des segments*********
  129.  
  130.         SEGINI IADJ,IVOIS,JMEM,JMEMN,LAGRAN
  131.         SEGINI BOOLEEN,JNT
  132.  
  133.         DO 20 I=1,NODES+1
  134.                 IADJ(I)=0
  135.                 JMEM(I)=0
  136.                 JMEMN(I)=0
  137. 20      CONTINUE
  138. C******************************************
  139.  
  140.  
  141.  
  142.         IPT1=MELEME
  143.         NGRAND=0
  144.         IADJ(1)=1
  145.         INC=0
  146.  
  147.  
  148.         DO 200 IO=1,MAX(1,LISOUS(/1))
  149.  
  150.             IF (LISOUS(/1).NE.0) THEN
  151.                IPT1=LISOUS(IO)
  152.             ENDIF
  153.  
  154.             DO 210 J=1,IPT1.NUM(/2)
  155.  
  156.              IF(IPT1.ITYPEL.EQ.22) THEN
  157.                DO 220 I=1,IPT1.NUM(/1)
  158. C               chaque element de type 22 a  au plus NB2 noeuds.
  159.                 LAGRAN(INC*NB2+I)=ICPR(IPT1.NUM(I,J))
  160. C               les noeuds de l'elements de type 22
  161. C               sont ranges dans le vecteur LAGRAN.
  162. 220            CONTINUE
  163.                DO 225 I=IPT1.NUM(/1)+1,NB2
  164.                  LAGRAN(INC*NB2+I)=0
  165. C                comme on a alloue la place memeoire maximale,
  166. C                on remplit les cases restantes avec des 0.
  167. 225            CONTINUE
  168.                INC=INC+1
  169. C              INC=le nbre d'elements de type 22.
  170.              ENDIF
  171.  
  172.              DO  230 I=1,IPT1.NUM(/1)
  173.               IJ=ICPR(IPT1.NUM(I,J))+1
  174.               JMEM(IJ)=JMEM(IJ)+1
  175. C             JMEM(I+1): nb elements auquel le noeud I appartient
  176. 230          CONTINUE
  177.  
  178. 210         CONTINUE
  179.  
  180.             NGRAND=MAX(NGRAND,IPT1.NUM(/2))
  181.  
  182. 200     CONTINUE
  183.  
  184.  
  185.  
  186.         NGRAND=NGRAND+1
  187.         JMEM(1)=1
  188.         DO 30 I=1,NODES
  189.             JMEM(I+1)=JMEM(I)+JMEM(I+1)
  190. C           JMEM(I+1)=indice de depart des elements
  191. C           auxquels le noeud I appartient.
  192. 30      CONTINUE
  193.         NBV=JMEM(NODES+1)
  194. C       NBV= dimension de IMEMOIR.
  195.         SEGINI IMEMOIR
  196.  
  197.  
  198.         IPT1=MELEME
  199.  
  200.         IMATOT=LISOUS(/1)*(NGRAND+1)
  201.         IF(IMATOT.GT.IENORM) THEN
  202.            CALL ERREUR(382)
  203.            RETURN
  204.         ENDIF
  205.  
  206.  
  207.  
  208.         DO 300 IO=1,MAX(1,LISOUS(/1))
  209.            IF (LISOUS(/1).NE.0) THEN
  210.                 IPT1=LISOUS(IO)
  211.            ENDIF
  212.            DO 350 I=1,IPT1.NUM(/1)
  213.            DO 350 J=1,IPT1.NUM(/2)
  214.              JMEMN(ICPR(IPT1.NUM(I,J))+1)=JMEMN(ICPR(IPT1.NUM(I,J))+1)+1
  215.              IMEMOIR(JMEM(ICPR(IPT1.NUM(I,J)))+
  216.      &       JMEMN(ICPR(IPT1.NUM(I,J))+1)-1)=J+IO*NGRAND
  217. C             on range dans IMEMOIR tous les elements des sous-objets
  218. C             IO auxquels appartient le noeud ICPR(IPT1.NUM(I,J)).
  219. C             On connait pour chaque element, le sous-objet auquel
  220. C             il appartient.
  221. 350        CONTINUE
  222. 300     CONTINUE
  223.  
  224.  
  225.         DO 400 I=1,NODES
  226.            IADJ(I+1)=IADJ(I)
  227.            DO 410 J=1,NODES
  228.                 BOOL(J)=.FALSE.
  229. 410        CONTINUE
  230.            DO 420 J=JMEM(I),JMEM(I+1)-1
  231.                ELEM=IMEMOIR(J)
  232. C              ELEM=element auquel appartient le noeud I.
  233.  
  234.                IPT1=MELEME
  235.                IF (LISOUS(/1).NE.0) IPT1=LISOUS(ELEM/NGRAND)
  236.                ELEM=MOD(ELEM,NGRAND)
  237.                DO 430 K=1,IPT1.NUM(/1)
  238. C              k representatif du nb de noeuds par elements.
  239.                   IK=ICPR(IPT1.NUM(K,ELEM))
  240.                   IF ((I.NE.IK).AND.
  241.      &                 (.NOT.(BOOL(IK)))) THEN
  242. C                 si i n'est pas  egal a un des nouveaux numeros des noeuds
  243. C                 de l'element ELEM et si il n'appartient pas deja  a l'ens des
  244. C                 voisins du noeud i(jadjc(i)),alors on le rajoute.
  245. C                   IVOIS(IADJ(I+1))=IK
  246.                     IVOIS(**)=IK
  247.                     IADJ(I+1)=IADJ(I+1)+1
  248.                     BOOL(IK)=.TRUE.
  249.                   ENDIF
  250. 430            CONTINUE
  251.  
  252.  
  253. 420        CONTINUE
  254. 400     CONTINUE
  255.  
  256.         SEGSUP JMEM,JMEMN,BOOLEEN
  257.         SEGSUP IMEMOIR
  258.  
  259.  
  260. C       affectation
  261. C************************
  262.  
  263.  
  264.         N=NODES
  265.         SEGINI IPOS,MASQUE
  266.  
  267.         DO 50 I=1,N
  268.            MASQ(I)=.TRUE.
  269.            IPOS(I)=0
  270.            IPOS(NODES+I)=0
  271. 50      CONTINUE
  272. C       initialement, les noeuds ne sont pas masques,ont donc
  273. C       une position nulle.
  274.  
  275.         DIM=0
  276. C       le nombre de noeuds renumerotes  DIM est initialement egal a zero.
  277.  
  278.  
  279. C               ****************************************
  280. C                       boucle principale
  281.         MPOS=0
  282.         DO 500 I=1,N
  283.  
  284. 550        IF(.NOT.MASQ(I)) GOTO 500
  285. C          si le noeud a ete renumerote, on passe au suivant.
  286.  
  287.            PIVOT=I
  288. C          sinon, il devient PIVOT.
  289.  
  290.           CALL SEPAR(IADJ,IVOIS,PIVOT,MASQUE,DIMSEP,N,IPOS,NODES)
  291. C         separe le domaine d'etude en 2 parties.
  292. C         on decrit le domaine d'etude a partir du pivot et on cherche la
  293. C         longueur maximale en decrivant les voisins de pivot, et leurs
  294. C         voisins... jusqu'a rencontrer un voisin masque. On cree alors
  295. C         une nouvelle separation.
  296. C         les noeuds masques delimitent la separation du domaine.
  297.  
  298.           DIM=DIM+DIMSEP
  299.  
  300.           N=N-DIMSEP
  301.  
  302.           IF (DIM.GE.NODES) GOTO 600
  303. C         si tous les noeuds ont une position non nulle,on arrete.
  304.  
  305.           GOTO 550
  306.  
  307. 500     CONTINUE
  308.  
  309. 600     SEGSUP MASQUE
  310.         SEGINI NUMERO
  311.  
  312.         CALL PREPA(IPOS,NPVOIS,NVOIS,NBENS,NUMERO,NODES,IADJ,IVOIS)
  313. C       procedure  donnant une nouvelle gestion de donnees.
  314. C       NBENS= nbre d'ensembles.
  315. C       un ensemble est compose de noeuds ayant meme NUMERO.
  316. C       NPVOIS(I,2)-NPVOIS(I,1)-1=nbre d'ensemble voisins a l'ensemble I.
  317. C       NVOIS(J), pour J compris entre NPVOIS(I,1) et NPVOIS(I,2)-1,
  318. C       contient les ensembles voisins de I.
  319.  
  320.         CALL MINDEG(NPVOIS,NVOIS,NBENS,IORDRE)
  321. C       on recherche l'ensemble de minimum degre.
  322. C       on le numerote et on recommence.
  323.  
  324.         CALL SORTID(IADJ,IVOIS,IORDRE,NUMERO,NODES,NBENS,JNT)
  325. C       JNT(I) est le nouveau numero du NOEUD I.
  326.  
  327.         SEGSUP IORDRE
  328.  
  329. C******************************************************
  330. C       insertion des lagrangiens.
  331. C******************************************************
  332. C       principe: on cherche le premier noeud auquel s'applique les 2
  333. C       coefficients de lagrange, on les insere de part et d'autre de
  334. C       ce noeud dans la numerotation.
  335.  
  336.         DO 1 I=1,NODES
  337.                 IPOS(JNT(I))=I
  338. 1       CONTINUE
  339.  
  340.         IF(NB2.NE.0) THEN
  341.            N=NB/NB2
  342.         ELSE
  343.            N=0
  344.         ENDIF
  345.  
  347.         DO 700 I=1,NODES
  348.                 NOEUD=IPOS(I)
  349.                 IF (NOEUD.EQ.0) GOTO 700
  350.                 DO 750 J=0,N-1
  351.                    IF(NOEUD.EQ.LAGRAN(J*NB2+1)) THEN
  353.                         IPOS(NODES+NOEUD)=0
  354.                         GOTO 700
  355.                    ENDIF
  356.  
  357.                    IF(NOEUD.EQ.LAGRAN(J*NB2+2))  THEN
  359.                         IPOS(NODES+NOEUD)=0
  360.                         GOTO 700
  361.                    ENDIF
  362. 750             CONTINUE
  363.                 IPOS(NODES+NOEUD)=I-NORME
  364. C               IPOS(NODES+NOEUD)= numero d'elimination obtenu en tenant
  365. C                                       compte des lagrangiens.
  366.  
  367. 700     CONTINUE
  368.  
  369.         NCOMPT=0
  370.         NORM=0
  371.         SEGINI ICOMPT,INVLAG
  372.  
  373.         DO 800 I=1,NODES
  374.  
  375.            NOEUD=IPOS(I)
  376. C          NOEUD =Ieme noeud elimine.
  377.  
  378.            IF(IPOS(NODES+NOEUD).EQ.0) GOTO 800
  379.  
  380.            ICOMPT(NOEUD)=0
  381.  
  382.            DO 850 J=0,NB1-1
  383.             IF(LAGRAN(J*NB2+1).NE.0) THEN
  384.              K=3
  385. 855          IF((K.LE.NB2).AND.(LAGRAN(J*NB2+K).NE.0)) THEN
  386.  
  387.                 IF(NOEUD.EQ.LAGRAN(J*NB2+K)) THEN
  388.  
  389.                    INVLAG(ICOMPT(NOEUD)+1)=J
  390. C                  INVLAG mémorise l'élément de type 22.
  391.  
  392. C                  on change le numero d'élimination de noeud,
  393. C                  en insérant de part et d'autre les coef de lagrange.
  394.  
  395.                    JNT(LAGRAN(J*NB2+1))=
  396.      &                IPOS(NOEUD+NODES)+NORM-ICOMPT(NOEUD)
  397.  
  398.                    LAGRAN(J*NB2+1)=0
  399. C                  permet de savoir si le lagrangien a ete insere.
  400.  
  401.                    NORM=NORM+1
  402.                    JNT(NOEUD)=
  403.      &               IPOS(NODES+NOEUD)+NORM-ICOMPT(NOEUD)
  404. C                  renumerotation de noeud.
  405.  
  406.                    DO 852 IJ=1,ICOMPT(NOEUD)
  407.                      IL=INVLAG(IJ)
  408.                      JNT(LAGRAN(IL*NB2+2))=JNT(LAGRAN(IL*NB2+2))+1
  409. 852                CONTINUE
  410.  
  411.                    NORM=NORM+1
  412.                    JNT(LAGRAN(J*NB2+2))=IPOS(NODES+NOEUD)+NORM
  413.                    ICOMPT(NOEUD)=ICOMPT(NOEUD)+1
  414.                    GOTO 850
  415.                 ELSE
  416.                    K=K+1
  417.                    GOTO 855
  418.                 ENDIF
  419.  
  420.              ENDIF
  421.  
  422.             ENDIF
  423.  
  424. 850        CONTINUE
  425.  
  426.            IF(ICOMPT(NOEUD).EQ.0) THEN
  427.               JNT(NOEUD)=IPOS(NODES+NOEUD)+NORM
  428.            ENDIF
  429.  
  430. 800     CONTINUE
  431.  
  432.         SEGSUP INVLAG
  433. C******************************************************************
  434.         DO 860 I=1,nbpts
  435.            IF(ICPR(I).NE.0) THEN
  436.                  ICPR(I)=JNT(ICPR(I))
  437. C                numerotation finale.
  438.            ENDIF
  439. 860     CONTINUE
  440.  
  441.         SEGSUP IADJ,IVOIS,JNT,IPOS,ICOMPT,NUMERO
  442.         SEGSUP LAGRAN,NPVOIS,NVOIS
  443.  
  444.         RETURN
  445.         END
  446.  
  447.  
  448.  
  449.  

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