Télécharger la plaquette Cast3M


fr en

Télécharger voro1.eso

Retour à la liste

Numérotation des lignes : 

voro1

  1. C VORO1     SOURCE    PV        20/04/01    21:16:44     10569          
  2. C
  3.       SUBROUTINE VORO1(IPT1,IPT2,LNBOIT,ILEA1,ITRC1,DDMAX,
  4.      &                 MTAB1)
  5. C----------------------------------------------------------------------C
  6. C     Partition de Voronoi (non coupee) d'un nuage de points (IPT1) a  C
  7. C     partir de la triangulation de Delaunay (IPT2)                    C
  8. C                                                                      C
  9. C     IPT1   = pointeur sur le maillage de points, centres des         C
  10. C              polyedres (elements POI1)                               C
  11. C     IPT2   = pointeur sur le maillage de la triangulation de         C
  12. C              Delaunay de IPT1 (avec les coins de la boite de         C
  13. C              triangulation) (elements TRI3 ou TET4)                  C
  14. C     LNBOIT = tableau contenant les numeros des noeuds formant les    C
  15. C              coins de la boite de triangulation                      C
  16. C     ILEA1  = pointeur sur le segment MADJACEN contenant la table des C
  17. C              connectivite des triangles de IPT1                      C
  18. C     ITRC1  = pointeur sur le segment MCIRCONS contenant la table des C
  19. C              coordonnees des cercles circonscrits des triangles de   C
  20. C              IPT1                                                    C
  21. C     XZERO   = critere de distance pour considerer deux centres de    C
  22. C              cercles circonscrits confondus                          C
  23. C                                                                      C
  24. C     MTAB1  = table de maillage des polyedres non coupes              C
  25. C----------------------------------------------------------------------C
  26.       IMPLICIT INTEGER(I-N)
  27.       IMPLICIT REAL*8 (A-H,O-Z)
  28. c
  29.  
  30. -INC PPARAM
  31. -INC CCOPTIO
  32. -INC SMELEME
  33. -INC SMCOORD
  34. -INC CCGEOME
  35. -INC SMLENTI
  36. -INC SMTABLE
  37.       DIMENSION LNBOIT(8)
  38.       DIMENSION XPI1(3)
  39.       LOGICAL BOOL0,BOOL1,BOOL2,BOOL3
  40. C
  41.       POINTEUR LCCC1.MLENTI,LAR1.MLENTI,LAR2.MLENTI
  42. C     Tables, maillages et listes pour la table de sortie
  43.       POINTEUR IPTG.MELEME,MTABC.MTABLE,MTABI1.MTABLE
  44.       POINTEUR IPTI1.MELEME,LFI1.MLENTI,IPTVI1.MELEME,MTABF.MTABLE
  45.       POINTEUR IPTF1.MELEME,LA1.MLENTI,MTABA.MTABLE,IPTA1.MELEME
  46. C
  47.       SEGMENT,MADJACEN
  48.          INTEGER LEAC(NBL1,IDIM+1,2)
  49.       ENDSEGMENT
  50.       POINTEUR ILEA1.MADJACEN
  51. C
  52.       SEGMENT,MCIRCONS
  53.          REAL*8 TRC(NBE1,4)
  54.       ENDSEGMENT
  55.       POINTEUR ITRC1.MCIRCONS
  56. C
  57.       SEGMENT,MVORO
  58.          INTEGER LACT(N1,N2)
  59.       ENDSEGMENT
  60.       POINTEUR ITL1.MVORO,ITCEL1.MVORO,ITCEL2.MVORO,ITARE1.MVORO
  61.       POINTEUR ITFAC1.MVORO,ITFAC2.MVORO
  62.  
  63. C
  64.       XZERO = DDMAX*1.D-10
  65.       segact mcoord*mod
  66. C
  67. C     1) Pre-traitement : creation des points representant les centres
  68. C        cercles circonscrits de chaque triangle (le meme numero est
  69. C        attribue aux centres de cercles circonscrits superposes)
  70. C
  71. C     Ajustement initial de MCOORD au nombre de cercles circonscrits
  72.       NBTRI=ITRC1.TRC(/1)
  73.       NBPTS0=nbpts
  74.       NBPTS=NBPTS0+NBTRI
  75.       SEGADJ,MCOORD
  76. C     LCCC1.LECT(i)=j le triangle i a pour Centre de Cerclce Circonscrit
  77. C                     le point numero j
  78.       JG=NBTRI
  79.       SEGINI,LCCC1
  80.       I1=0
  81. C     Boucle sur les centres des cercles circonscrits des triangles
  82.       DO I=1,NBTRI
  83.          DO K=1,IDIM
  84.             XPI1(K)=ITRC1.TRC(I,K+1)
  85.          ENDDO
  86.          DO J=1,NBTRI
  87.             NCJ=0
  88.             IF (J.GE.I) GOTO 101
  89. C           Test si les coordonnees de I et J sont egales
  90.             DIJ1=0.
  91.             DO K=1,IDIM
  92.                XPJ=ITRC1.TRC(J,K+1)
  93.                DIJ1=DIJ1+(XPI1(K)-XPJ)*(XPI1(K)-XPJ)
  94. c               XDEN=XPI1(K)
  95. c               IF (XDEN.EQ.0.) XDEN=1.
  96. c               IF (ABS((XPI1(K)-XPJ)/XDEN).LT.XZERO) NCJ=NCJ+1
  97.             ENDDO
  98.             IF(DIJ1.LT.XZERO) NCJ=IDIM
  99. C           Si les centres I et J sont superposes
  100.             IF (NCJ.EQ.IDIM) THEN
  101.                LCCC1.LECT(I)=LCCC1.LECT(J)
  102.                GOTO 102
  103.             ENDIF
  104. 101         CONTINUE
  105.          ENDDO
  106. C        Creation d'un nouveau point dans MCOORD
  107.          I1=I1+1
  108.          NREFI=NBPTS0+I1
  109.          LCCC1.LECT(I)=NREFI
  110. C        Ecriture des ses coordonnees
  111.          DO J=1,IDIM
  112.             XCOOR(((NREFI-1)*(IDIM+1))+J)=XPI1(J)
  113.          ENDDO
  114. C        et de la densite associee, choisie nulle
  115.          XCOOR(((NREFI-1)*(IDIM+1))+(IDIM+1))=0.
  116. 102      CONTINUE
  117.       ENDDO
  118. C     Ajustement final de MCOORD (on a ajoute I1 points)
  119.       IF (I1.NE.NBTRI) THEN
  120.          NBPTS=NBPTS0+I1
  121.          SEGADJ,MCOORD
  122.       ENDIF
  123. C
  124. C
  125. C     2.0) Initialisations avant creation du maillage des polyedres
  126. C          IPTG : Maillage des aretes des polyedres (elements SEG2)
  127.       NBSOUS=0
  128.       NBREF=0
  129.       NBNN=2
  130.       NBELEM=I1*4
  131.       SEGINI,IPTG
  132.       IPTG.ITYPEL=2
  133.       NBARETES=0
  134. C     Table MTAB1 de sortie contenant les maillages des polyedres
  135.       CALL CRTABL(MTAB1)
  136.       CALL ECCTAB(MTAB1,'MOT     ',0,0,'VISU',0,0,
  137.      &                  'MAILLAGE',0,0,0,0,IPTG)
  138.       CALL CRTABL(MTABC)
  139.       CALL ECCTAB(MTAB1,'MOT     ',0,0,'CELL',0,0,
  140.      &                  'TABLE   ',0,0,0,0,MTABC)
  141.       IF (IDIM.EQ.3) THEN
  142.          CALL CRTABL(MTABF)
  143.          CALL ECCTAB(MTAB1,'MOT     ',0,0,'FACS',0,0,
  144.      &                     'TABLE   ',0,0,0,0,MTABF)
  145.       ENDIF
  146.       CALL CRTABL(MTABA)
  147.       CALL ECCTAB(MTAB1,'MOT     ',0,0,'ARTS',0,0,
  148.      &                  'TABLE   ',0,0,0,0,MTABA)
  149. C     ITL1 : Table indiquant si la ligne entre deux centres de cercle
  150. C     circonscrits existe deja ou non
  151.       N1=I1
  152.       N2=I1
  153.       SEGINI,ITL1
  154. C     ITCEL1 : Table des polyedres voisins d'un polyedre
  155. C     ITCEL2 : Table des aretes d'un polyedre
  156.       N1=10
  157.       N2=IPT1.NUM(/2)
  158.       SEGINI,ITCEL1,ITCEL2
  159. C     ITARE1 : Table des cellules appuyees sur une arete, on la stocke
  160. C     dans la table de sortie pour pouvoir la retrouver lors du
  161. C     decoupage dans VORO2
  162. C     ITARE1.LACT(ICE,NLA)=I
  163. C     NLA : numero (dans IPTG) de l'arete
  164. C     ICE : index des cellules appuyees sur NLA (de 1 a n, puis 0)
  165. C     I : numero (dans IPT1) de la ICE-ieme cellule appuyee sur NLA
  166.       N1=10
  167.       N2=IPTG.NUM(/2)
  168.       SEGINI,ITARE1
  169.       CALL ECCTAB(MTAB1,'MOT     ',0,0,'TAR',0,0,
  170.      &                  'MVORO   ',0,0,0,0,ITARE1)
  171. C
  172. C
  173. C
  174. C     2.1) Maillage des cellules de Voronoi en dimension 2
  175. C     Boucle sur les centres des polyedres de Voronoi
  176.       IF (IDIM.EQ.2) THEN
  177.          DO I=1,IPT1.NUM(/2)
  178.             NREF_I=IPT1.NUM(1,I)
  179. C           Initialisation des indices pour la cellule NREF_I
  180.             CALL CRTABL(MTABI1)
  181.             CALL ECCTAB(MTABC,'POINT   ',0,0,' ',0,NREF_I,
  182.      &                        'TABLE   ',0,0,0,0,MTABI1)
  183. C           --> maillage de la cellule
  184.             NBNN=2
  185.             NBELEM=I1*4
  186.             SEGINI,IPTI1
  187.             IPTI1.ITYPEL=2
  188.             NBELEI1=0
  189.             CALL ECCTAB(MTABI1,'MOT     ',0,0,'VISU',0,0,
  190.      &                         'MAILLAGE',0,0,0,0,IPTI1)
  191. C           --> liste des aretes
  192.             JG=I1*4
  193.             SEGINI,LA1
  194.             CALL ECCTAB(MTABI1,'MOT     ',0,0,'ARTS',0,0,
  195.      &                         'LISTENTI',0,0,0,0,LA1)
  196. C           --> maillage des voisns
  197.             NBNN=1
  198.             NBELEM=I1*4
  199.             SEGINI,IPTVI1
  200.             IPTVI1.ITYPEL=1
  201.             NVI1=0
  202.             CALL ECCTAB(MTABI1,'MOT     ',0,0,'VOIS',0,0,
  203.      &                         'MAILLAGE',0,0,0,0,IPTVI1)
  204. C           Recherche du premier triangle appuye sur NREFA
  205.             NREFA=IPT1.NUM(1,I)
  206.             DO J=1,IPT2.NUM(/2)
  207.                DO K=1,IPT2.NUM(/1)
  208.                   IF (IPT2.NUM(K,J).EQ.NREFA) GOTO 211
  209.                ENDDO
  210.             ENDDO
  211. 211         NT0=J
  212.             NA=K
  213.             NB=NA+1
  214.             IF (NB.EQ.4) NB=1
  215.             NC=6-NA-NB
  216.             NREFB=IPT2.NUM(NB,NT0)
  217.             NREFC=IPT2.NUM(NC,NT0)
  218. C           Balayage des triangles appuyes sur NREFA en tournant autour
  219. C           de ce point, a patrir du triangle NT0
  220. C           La boucle s'arrete quand on est revenu au triangle NT0
  221.             NT1=NT0
  222.             NRCC1=LCCC1.LECT(NT1)
  223.             NCVI=0
  224.             NAI=0
  225.             DO J=1,I1*4
  226.                IJOK=0
  227. C              Determination du point NRCC2 a relier a NRCC1
  228.                NT2=ILEA1.LEAC(NT1,NB,1)
  229.                NRCC2=LCCC1.LECT(NT2)
  230. C              Si la liaison entre CC1 et CC2 existe, on la recupere et
  231. C              on passe au triangle suivant
  232.                IF (ITL1.LACT(NRCC1-NBPTS0,NRCC2-NBPTS0).EQ.1) THEN
  233. C                 On determine le numero de l'arete en question
  234.                   DO K=1,ITCEL2.LACT(/1)
  235.                      NVI=0
  236.                      DO L=1,IPT1.NUM(/2)
  237.                         IF (IPT1.NUM(1,L).EQ.NREFC) NVI=L
  238.                      ENDDO
  239.                      IF (NVI.EQ.0) THEN
  240.                         DO L=1,4
  241.                            IF (LNBOIT(L).EQ.NREFC) NVI=-1*L
  242.                         ENDDO
  243.                      ENDIF
  244.                      NAR1=ITCEL2.LACT(K,NVI)
  245.                      NPA=IPTG.NUM(1,NAR1)
  246.                      NPB=IPTG.NUM(2,NAR1)
  247.                      IF (((NPA.EQ.NRCC1).AND.(NPB.EQ.NRCC2)).OR.
  248.      &                   ((NPA.EQ.NRCC2).AND.(NPB.EQ.NRCC1))) THEN
  249.                         GOTO 212
  250.                      ENDIF
  251.                   ENDDO
  252.                ENDIF
  253. C              Si les centres CC1 et CC2 sont confondus, il n'y a pas de
  254. C              liaison a creer, on passe au triangle suivant
  255.                IF (NRCC1.EQ.NRCC2) THEN
  256.                   GOTO 214
  257.                ENDIF
  258. C              Dans les autres cas, il y a bien creation d'une nouvelle
  259. C              arete dans le maillage IPTG de la cellule I
  260.                IJOK=1
  261.                NBARETES=NBARETES+1
  262.                IPTG.NUM(1,NBARETES)=NRCC1
  263.                IPTG.NUM(2,NBARETES)=NRCC2
  264.                NAR1=NBARETES
  265. C              Remplissage de la table ILT1 des liaisons entre centres
  266. C              de cercles circonscrits
  267.                ITL1.LACT(NRCC1-NBPTS0,NRCC2-NBPTS0)=1
  268.                ITL1.LACT(NRCC2-NBPTS0,NRCC1-NBPTS0)=1
  269. C              Remplissage de la table ITCEL1 des numeros des cellules
  270. C              voisines de la cellule I
  271.                NVI=0
  272.                DO K=1,IPT1.NUM(/2)
  273.                   IF (IPT1.NUM(1,K).EQ.NREFC) NVI=K
  274.                ENDDO
  275.                IF (NVI.EQ.0) THEN
  276.                   DO K=1,4
  277.                      IF (LNBOIT(K).EQ.NREFC) NVI=-1*K
  278.                   ENDDO
  279.                ENDIF
  280. 212            IF (NVI.GT.0) THEN
  281.                   NCVI=NCVI+1
  282.                   IF (NCVI.GT.ITCEL1.LACT(/1)) THEN
  283.                      N1=ITCEL1.LACT(/1)+10
  284.                      N2=ITCEL1.LACT(/2)
  285.                      SEGADJ,ITCEL1
  286.                   ENDIF
  287.                   ITCEL1.LACT(NCVI,I)=NVI
  288.                ENDIF
  289. C              Remplissage de la table ITCEL2 des numeros d'aretes
  290. C              appartenant a la cellule I
  291.                NAI=NAI+1
  292.                IF (NAI.GT.ITCEL2.LACT(/1)) THEN
  293.                   N1=ITCEL2.LACT(/1)+10
  294.                   N2=ITCEL2.LACT(/2)
  295.                   SEGADJ,ITCEL2
  296.                ENDIF
  297.                ITCEL2.LACT(NAI,I)=NAR1
  298. C              Remplissage de la table ITARE1 des numeros de cellules
  299. C              contenant l'arrete creee/recuperee
  300.                ICE=0
  301.                DO K=1,ITARE1.LACT(/1)
  302.                   IF (ITARE1.LACT(K,NAR1).EQ.0) THEN
  303.                      ICE=K
  304.                      GOTO 213
  305.                   ENDIF
  306.                ENDDO
  307.                IF (ICE.EQ.0) THEN
  308.                   ICE=ITARE1.LACT(/1)+1
  309.                   N1=ITARE1.LACT(/1)+10
  310.                   N2=ITARE1.LACT(/2)
  311.                   SEGADJ,ITARE1
  312.                ENDIF
  313. 213            ITARE1.LACT(ICE,NAR1)=I
  314. C              Ajout d'un element au maillage IPTI1 de la cellule I
  315.                NBELEI1=NBELEI1+1
  316.                IPTI1.NUM(1,NBELEI1)=NRCC1
  317.                IPTI1.NUM(2,NBELEI1)=NRCC2
  318. C              Ajout du numero d'arete dans la liste de la cellule I
  319.                LA1.LECT(NBELEI1)=NAR1
  320. C              Ajout d'un voisin a la cellule I
  321.                IF (NVI.GT.0) THEN
  322.                   NVI1=NVI1+1
  323.                   NREFVI=IPT1.NUM(1,NVI)
  324.                   IPTVI1.NUM(1,NVI1)=NREFVI
  325.                ENDIF
  326. C              Ajout de l'arete dans MTABA
  327.                IF (IJOK.EQ.1) THEN
  328.                   NBNN=2
  329.                   NBELEM=1
  330.                   SEGINI,IPTA1
  331.                   IPTA1.ITYPEL=2
  332.                   IPTA1.NUM(1,1)=NRCC1
  333.                   IPTA1.NUM(2,1)=NRCC2
  334.                   CALL ECCTAB(MTABA,'ENTIER  ',NBARETES,0,' ',0,0,
  335.      &                              'MAILLAGE',0,0,0,0,IPTA1)
  336.                ENDIF
  337. C              Test d'arret de la boucle (est on revenu sur le triangle
  338. C              de depart ?)
  339. 214            IF (NT2.EQ.NT0) THEN
  340.                   GOTO 215
  341.                ENDIF
  342. C              Nouveau points A, B et C pour le prochain triangle
  343.                DO K=1,IPT2.NUM(/1)
  344.                   IF (IPT2.NUM(K,NT2).EQ.NREFA) NA=K
  345.                   IF (IPT2.NUM(K,NT2).EQ.NREFC) NB=K
  346.                ENDDO
  347.                NC=6-NA-NB
  348.                NREFB=IPT2.NUM(NB,NT2)
  349.                NREFC=IPT2.NUM(NC,NT2)
  350.                NT1=NT2
  351.                NRCC1=LCCC1.LECT(NT1)
  352.             ENDDO
  353. 215         CONTINUE
  354. C           Ajustement des maillages de la cellule I, de ses voisins et
  355. C           de la liste des aretes
  356.             NBNN=IPTI1.NUM(/1)
  357.             NBELEM=NBELEI1
  358.             SEGADJ,IPTI1
  359.             JG=NBELEI1
  360.             SEGADJ,LA1
  361.             NBNN=IPTVI1.NUM(/1)
  362.             NBELEM=NVI1
  363.             SEGADJ,IPTVI1
  364.          ENDDO
  365. C
  366. C
  367. C
  368. C     2.2) Maillage des cellules de Voronoi en dimension 3
  369.       ELSEIF (IDIM.EQ.3) THEN
  370. C        ITFAC1 : tableau indiquant si la face (I,J) existe deja
  371. C        ou I et J sont des numeros de cellule
  372.          N1=IPT1.NUM(/2)
  373.          N2=N1
  374.          SEGINI,ITFAC1
  375. C        ITFAC2 : tableau indiquant si la face (I,J) existe deja
  376. C        ou I est un numero de cellule exterieures infinie et J un
  377. C        numero des coins de la boite de triangulation
  378.          N1=IPT1.NUM(/2)
  379.          N2=8
  380.          SEGINI,ITFAC2
  381. C        Boucle sur les centres des polyedres de Voronoi
  382.          NBFACES=0
  383.          DO I=1,IPT1.NUM(/2)
  384. C           On va chercher la table de la cellule I
  385.             NREF_I=IPT1.NUM(1,I)
  386. C           Initialisation des indices pour la cellule NREF_I
  387.             CALL CRTABL(MTABI1)
  388.             CALL ECCTAB(MTABC,'POINT   ',0,0,' ',0,NREF_I,
  389.      &                        'TABLE   ',0,0,0,0,MTABI1)
  390. C           --> maillage de la cellule
  391.             NBNN=2
  392.             NBELEM=I1*4
  393.             SEGINI,IPTI1
  394.             IPTI1.ITYPEL=2
  395.             NBELEI1=0
  396.             CALL ECCTAB(MTABI1,'MOT     ',0,0,'VISU',0,0,
  397.      &                         'MAILLAGE',0,0,0,0,IPTI1)
  398. C           --> liste des faces
  399.             JG=I1*4
  400.             SEGINI,LFI1
  401.             NFI1=0
  402.             CALL ECCTAB(MTABI1,'MOT     ',0,0,'FACS',0,0,
  403.      &                         'LISTENTI',0,0,0,0,LFI1)
  404. C           --> maillage des voisns
  405.             NBNN=1
  406.             NBELEM=I1*4
  407.             SEGINI,IPTVI1
  408.             IPTVI1.ITYPEL=1
  409.             IV1=0
  410.             CALL ECCTAB(MTABI1,'MOT     ',0,0,'VOIS',0,0,
  411.      &                         'MAILLAGE',0,0,0,0,IPTVI1)
  412. C           Recherche du premier tetraedre appuye sur NREFA
  413.             NCVI=0
  414.             NREFA=IPT1.NUM(1,I)
  415.             DO J=1,IPT2.NUM(/2)
  416.                DO K=1,IPT2.NUM(/1)
  417.                   IF (IPT2.NUM(K,J).EQ.NREFA) GOTO 221
  418.                ENDDO
  419.             ENDDO
  420. 221         NT0=J
  421.             NA=K
  422.             NB=NA+1
  423.             NC=NA-1
  424.             IF (NB.EQ.5) NB=1
  425.             IF (NC.EQ.0) NC=3
  426.             ND=10-NA-NB-NC
  427. C           Liste des aretes appuyees sur NREFA autour desquelles on
  428. C           va tourner. On incremente ces liste au fur et a mesure que
  429. C           l'on rencontre de nouveau tetraedres appuyes sur NT0
  430.             JG=IPTG.NUM(/2)
  431.             SEGINI,LAR1,LAR2
  432. C           On connait les trois premieres aretes appuyees sur NREFA
  433.             LAR1.LECT(1)=IPT2.NUM(NB,NT0)
  434.             LAR1.LECT(2)=IPT2.NUM(NC,NT0)
  435.             LAR1.LECT(3)=IPT2.NUM(ND,NT0)
  436. C           et les numeros des premiers tetraedres contenant ces aretes
  437.             LAR2.LECT(1)=NT0
  438.             LAR2.LECT(2)=NT0
  439.             LAR2.LECT(3)=NT0
  440. C           Balayage des aretes des tetraedres appuyes sur NREFA,
  441. C           a patrir du tetraedre NT0 (on parcoure LAR1)
  442. C           Chaque nouvelle arete definit potentielement une nouvelle
  443. C           face de la partition de Voronoi
  444. C           La liste LAR1 s'incremente, la boucle s'arrete lorsque l'on
  445. C           a parcouru toutes les aretes dans LAR1
  446.             DO J=1,LAR1.LECT(/1)
  447.                IJOK=0
  448. C              On tourne autour de l'arete J de la liste LAR1
  449.                NREFB=LAR1.LECT(J)
  450. C              Test d'arret de la boucle, si on a fini de parcourir la
  451. C              liste LAR1, on a fini !
  452.                IF (NREFB.EQ.0) THEN
  453.                   GOTO 229
  454.                ENDIF
  455. C              Calcul de NVI, numero de cette cellule voisine
  456.                NVI=0
  457.                DO K=1,IPT1.NUM(/2)
  458.                   IF (IPT1.NUM(1,K).EQ.NREFB) NVI=K
  459.                ENDDO
  460. C              ou bien s'il s'agit d'un sommet de la boite de
  461. C              triangulation
  462.                IF (NVI.EQ.0) THEN
  463.                   DO K=1,8
  464.                      IF (LNBOIT(K).EQ.NREFB) NVI=-1*K
  465.                   ENDDO
  466.                ENDIF
  467. C              NFIJ est le numero de la face entre I et J
  468. C              (si elle existe deja)
  469.                IF (NVI.GT.0) THEN
  470.                   NFIJ=ITFAC1.LACT(NVI,I)
  471.                ELSE
  472.                   NVII=-1*NVI
  473.                   NFIJ=ITFAC2.LACT(I,NVII)
  474.                ENDIF
  475. C              IJOK est mis a 1 quand on identifie une face non deja
  476. C              existante (mais elle peut etre vide !)
  477.                IF (NFIJ.EQ.0) IJOK=1
  478. C              Dans le cas d'une potentielle nouvelle face, on prepare
  479. C              son maillage et sa liste des aretes
  480.                IF (IJOK.EQ.1) THEN
  481.                   NBNN=2
  482.                   NBELEM=LAR1.LECT(/1)
  483.                   SEGINI,IPTF1
  484.                   IPTF1.ITYPEL=2
  485.                   JG=LAR1.LECT(/1)
  486.                   SEGINI,LA1
  487.                ENDIF
  488.                NBELEF1=0
  489. C              Le premier tetraedre balaye est NT0
  490.                NT0=LAR2.LECT(J)
  491. C              Calcul de NA,NB,NC et ND dans NT0
  492.                DO K=1,IPT2.NUM(/1)
  493.                   IF (IPT2.NUM(K,NT0).EQ.NREFA) NA=K
  494.                   IF (IPT2.NUM(K,NT0).EQ.NREFB) NB=K
  495.                ENDDO
  496.                IF ((NA.EQ.1).AND.(NB.EQ.2)) NC=3
  497.                IF ((NA.EQ.1).AND.(NB.EQ.3)) NC=2
  498.                IF ((NA.EQ.1).AND.(NB.EQ.4)) NC=2
  499.                IF ((NA.EQ.2).AND.(NB.EQ.1)) NC=3
  500.                IF ((NA.EQ.2).AND.(NB.EQ.3)) NC=1
  501.                IF ((NA.EQ.2).AND.(NB.EQ.4)) NC=1
  502.                IF ((NA.EQ.3).AND.(NB.EQ.1)) NC=2
  503.                IF ((NA.EQ.3).AND.(NB.EQ.2)) NC=1
  504.                IF ((NA.EQ.3).AND.(NB.EQ.4)) NC=1
  505.                IF ((NA.EQ.4).AND.(NB.EQ.1)) NC=2
  506.                IF ((NA.EQ.4).AND.(NB.EQ.2)) NC=1
  507.                IF ((NA.EQ.4).AND.(NB.EQ.3)) NC=1
  508.                ND=10-NA-NB-NC
  509.                NREFC=IPT2.NUM(NC,NT0)
  510.                NREFD=IPT2.NUM(ND,NT0)
  511. C              Boucle autour de l'arete NREFA-NREFB, en partant du
  512. C              tetraedre NT0, on construit ainsi la face (I,J) en
  513. C              tournant autour de ses aretes
  514. C              La boucle s'arrete quand on est revenu au tetraedre NT0
  515.                NT1=NT0
  516.                NRCC1=LCCC1.LECT(NT1)
  517.                DO K=1,I1*4
  518.                   NT2=ILEA1.LEAC(NT1,NC,1)
  519. C                 On ajoute NREFC et NREFD a la liste des aretes a
  520. C                 balayer plus tard (si non deja presents)
  521.                   DO L=1,LAR1.LECT(/1)
  522.                      IF (LAR1.LECT(L).EQ.NREFC) GOTO 222
  523.                      IF (LAR1.LECT(L).EQ.0) THEN
  524.                         LAR1.LECT(L)=NREFC
  525.                         LAR2.LECT(L)=NT2
  526.                         GOTO 222
  527.                      ENDIF
  528.                   ENDDO
  529. 222               CONTINUE
  530.                   DO L=1,LAR1.LECT(/1)
  531.                      IF (LAR1.LECT(L).EQ.NREFD) GOTO 223
  532.                      IF (LAR1.LECT(L).EQ.0) THEN
  533.                         LAR1.LECT(L)=NREFD
  534.                         LAR2.LECT(L)=NT2
  535.                         GOTO 223
  536.                      ENDIF
  537.                   ENDDO
  538. 223               CONTINUE
  539.                   NRCC2=LCCC1.LECT(NT2)
  540. C
  541. C                 Plusieurs cas sont possibles pour l'arete reliant
  542. C                 CC1 et CC2 :
  543. C
  544. C                 CAS 1 : la liaison entre CC1 et CC2 existe
  545. C                 => on la recupere et on passe au tetraedre suivant
  546.                   IF (ITL1.LACT(NRCC1-NBPTS0,NRCC2-NBPTS0).EQ.1) THEN
  547. C                    On cherche le numero de l'arete existante reliant
  548. C                    deja CC1 et CC2
  549.                      DO L=1,IPTG.NUM(/2)
  550.                         NPA=IPTG.NUM(1,L)
  551.                         NPB=IPTG.NUM(2,L)
  552.                         IF (((NPA.EQ.NRCC1).AND.(NPB.EQ.NRCC2)).OR.
  553.      &                      ((NPA.EQ.NRCC2).AND.(NPB.EQ.NRCC1))) THEN
  554.                            NLA=L
  555.                            GOTO 224
  556.                         ENDIF
  557.                      ENDDO
  558.                   ENDIF
  559. C                 CAS 2 : les centres CC1 et CC2 sont confondus
  560. C                 => il n'y a pas d'arete a faire ! on passe au
  561. C                    tetraedre suivant
  562.                   IF (NRCC1.EQ.NRCC2) THEN
  563.                      GOTO 227
  564.                   ENDIF
  565. C                 CAS 3 : creation d'une nouvelle arete dans IPTG
  566.                   NBARETES=NBARETES+1
  567.                   IPTG.NUM(1,NBARETES)=NRCC1
  568.                   IPTG.NUM(2,NBARETES)=NRCC2
  569.                   NLA=NBARETES
  570. C                 et ajout de cette arete dans MTABA
  571.                   NBNN=2
  572.                   NBELEM=1
  573.                   SEGINI,IPTA1
  574.                   IPTA1.ITYPEL=2
  575.                   IPTA1.NUM(1,1)=NRCC1
  576.                   IPTA1.NUM(2,1)=NRCC2
  577.                   CALL ECCTAB(MTABA,'ENTIER  ',NLA,0,' ',0,0,
  578.      &                              'MAILLAGE',0,0,0,0,IPTA1)
  579. C                 Remplissage de la table ILT1 des liaisons en centres
  580. C                 de cercles circonscrits
  581.                   ITL1.LACT(NRCC1-NBPTS0,NRCC2-NBPTS0)=1
  582.                   ITL1.LACT(NRCC2-NBPTS0,NRCC1-NBPTS0)=1
  583. C
  584. 224               CONTINUE
  585. C                 On incremente le nombre d'arete dans cette face
  586.                   NBELEF1=NBELEF1+1
  587. C                 Dans le cas d'une potentielle nouvelle face, on
  588. C                 incremente le maillage de la face et la liste de ses
  589. C                 aretes (si l'arete en question n'est pas presente)
  590.                   IF (IJOK.EQ.1) THEN
  591.                      DO L=1,LA1.LECT(/1)
  592.                         IF (LA1.LECT(L).EQ.NLA) GOTO 230
  593.                         IF (LA1.LECT(L).EQ.0) GOTO 231
  594.                      ENDDO
  595. 231                  IF (NBELEF1.GT.LA1.LECT(/1)) THEN
  596.                         JG=LA1.LECT(/1)+10
  597.                         SEGADJ,LA1
  598.                         NBNN=IPTF1.NUM(/1)
  599.                         NBELEM=IPTF1.NUM(/2)+10
  600.                         SEGADJ,IPTF1
  601.                      ENDIF
  602. C                    ajout dans la liste des aretes pour cette face
  603.                      LA1.LECT(NBELEF1)=NLA
  604. C                    et au maillage de la face
  605.                      IPTF1.NUM(1,NBELEF1)=NRCC1
  606.                      IPTF1.NUM(2,NBELEF1)=NRCC2
  607. 230                  CONTINUE
  608.                   ENDIF
  609. C                 Dans le cas d'une potentielle nouvelle face, si le
  610. C                 nombre d'aretes est egal a 3, il s'agit bel et bien
  611. C                 d'une face
  612.                   IF ((NBELEF1.EQ.3).AND.(IJOK.EQ.1)) THEN
  613. C                    on l'ajoute a MTABF
  614.                      NBFACES=NBFACES+1
  615.                      NFIJ=NBFACES
  616.                      CALL CRTABL(MTABF1)
  617.                      CALL ECCTAB(MTABF,'ENTIER  ',NFIJ,0,' ',0,0,
  618.      &                                 'TABLE   ',0,0,0,0,MTABF1)
  619.                      CALL ECCTAB(MTABF1,'MOT     ',0,0,'VISU',0,0,
  620.      &                                  'MAILLAGE',0,0,0,0,IPTF1)
  621.                      CALL ECCTAB(MTABF1,'MOT     ',0,0,'ARTS',0,0,
  622.      &                                  'LISTENTI',0,0,0,0,LA1)
  623.                   ENDIF
  624. C                 Si le nombre d'arete est egal a 3, le point J est bel
  625. C                 et bien un voisin de I
  626.                   IF (NBELEF1.EQ.3) THEN
  627. C                    S'il s'agit d'une face exterieure (J est sommet de
  628. C                    la boite de triangulation) on l'indique juste dans
  629. C                    ITFAC2
  630.                      IF (NVI.LT.0) THEN
  631.                         NVII=-1*NVI
  632.                         ITFAC2.LACT(I,NVII)=NFIJ
  633. C                    Sinon, on ajoute un voisin et on marque la face
  634. C                    dans ITFAC1
  635.                      ELSE
  636.                         NREFVI=IPT1.NUM(1,NVI)
  637.                         IV1=IV1+1
  638.                         IF (IV1.GT.IPTVI1.NUM(/2)) THEN
  639.                            NBNN=IPTVI1.NUM(/1)
  640.                            NBELEM=IPTVI1.NUM(/2)+10
  641.                            SEGADJ,IPTVI1
  642.                         ENDIF
  643.                         IPTVI1.NUM(1,IV1)=NREFVI
  644.                         ITFAC1.LACT(I,NVI)=NFIJ
  645.                         ITFAC1.LACT(NVI,I)=NFIJ
  646.                      ENDIF
  647. C                    Ajout du numero global de la face dans la liste des
  648. C                    faces de la cellule NREF_I
  649.                      NFI1=NFI1+1
  650.                      IF (NFI1.GT.LFI1.LECT(/1)) THEN
  651.                         JG=LFI1.LECT(/1)+10
  652.                         SEGADJ,LFI1
  653.                      ENDIF
  654.                      LFI1.LECT(NFI1)=NFIJ
  655.                   ENDIF
  656. C                 Remplissage de la table ITCEL2 des numeros d'aretes
  657. C                 appartenant a la cellule I
  658.                   DO L=1,ITCEL2.LACT(/1)
  659.                      IF (ITCEL2.LACT(L,I).EQ.NLA) THEN
  660.                         GOTO 227
  661.                      ENDIF
  662.                      IF (ITCEL2.LACT(L,I).EQ.0) GOTO 225
  663.                   ENDDO
  664. 225               NBELEI1=NBELEI1+1
  665.                   IF (NBELEI1.GT.ITCEL2.LACT(/1)) THEN
  666.                      N1=ITCEL2.LACT(/1)+10
  667.                      N2=ITCEL2.LACT(/2)
  668.                      SEGADJ,ITCEL2
  669.                   ENDIF
  670.                   ITCEL2.LACT(NBELEI1,I)=NLA
  671. C                 Ajout de l'arete au maillage de la cellule NREF_I
  672.                   IPTI1.NUM(1,NBELEI1)=NRCC1
  673.                   IPTI1.NUM(2,NBELEI1)=NRCC2
  674. C                 Remplissage de la table ITARE1 des numeros de cellules
  675. C                 contenant l'arrete creee/recuperee
  676.                   ICE=0
  677.                   DO L=1,ITARE1.LACT(/1)
  678.                      IF (ITARE1.LACT(L,NLA).EQ.0) THEN
  679.                         ICE=L
  680.                         GOTO 226
  681.                      ENDIF
  682.                   ENDDO
  683.                   IF (ICE.EQ.0) THEN
  684.                      ICE=ITARE1.LACT(/1)+1
  685.                      N1=ITARE1.LACT(/1)+10
  686.                      N2=ITARE1.LACT(/2)
  687.                      SEGADJ,ITARE1
  688.                   ENDIF
  689. 226               ITARE1.LACT(ICE,NLA)=I
  690. C                 Test d'arret de la boucle (est on revenu sur le
  691. C                 tetraedre de depart ?)
  692. 227               IF (NT2.EQ.NT0) THEN
  693. C                    Remplissage de la table ITCEL1 des numeros des
  694. C                    cellules voisines de la cellule I
  695.                      IF (NBELEF1.GE.3) THEN
  696.                         NCVI=NCVI+1
  697.                         IF (NCVI.GT.ITCEL1.LACT(/1)) THEN
  698.                            N1=ITCEL1.LACT(/1)+10
  699.                            N2=ITCEL1.LACT(/2)
  700.                            SEGADJ,ITCEL1
  701.                         ENDIF
  702.                         ITCEL1.LACT(NCVI,I)=NVI
  703.                         IF (IJOK.EQ.1) THEN
  704. C                          ajustement du maillage de la face (I,J)
  705.                            NBNN=IPTF1.NUM(/1)
  706.                            NBELEM=NBELEF1
  707.                            SEGADJ,IPTF1
  708. C                          et de la liste de ses aretes
  709.                            JG=NBELEF1
  710.                            SEGADJ,LA1
  711.                         ENDIF
  712.                      ENDIF
  713.                      GOTO 228
  714.                   ENDIF
  715. C                 Nouveau points A, B, C et D pour le prochain tetraedre
  716.                   DO L=1,IPT2.NUM(/1)
  717.                      IF (IPT2.NUM(L,NT2).EQ.NREFA) NA=L
  718.                      IF (IPT2.NUM(L,NT2).EQ.NREFB) NB=L
  719.                      IF (IPT2.NUM(L,NT2).EQ.NREFD) NC=L
  720.                   ENDDO
  721.                   ND=10-NA-NB-NC
  722.                   NREFC=IPT2.NUM(NC,NT2)
  723.                   NREFD=IPT2.NUM(ND,NT2)
  724.                   NT1=NT2
  725.                   NRCC1=LCCC1.LECT(NT1)
  726.                ENDDO
  727. 228            CONTINUE
  728.             ENDDO
  729. 229         CONTINUE
  730. C           Ajustement du maillage de la cellule I, de la liste de ses
  731. C           faces et de ses voisins
  732.             NBNN=IPTI1.NUM(/1)
  733.             NBELEM=NBELEI1
  734.             SEGADJ,IPTI1
  735.             JG=NFI1
  736.             SEGADJ,LFI1
  737.             NBNN=IPTVI1.NUM(/1)
  738.             NBELEM=IV1
  739.             SEGADJ,IPTVI1
  740.          ENDDO
  741.          SEGSUP,LAR1,LAR2,ITFAC1,ITFAC2
  742.       ENDIF
  743. C
  744. C     2.3) Ajustement final de IPTG et ITARE1
  745.       NBNN=IPTG.NUM(/1)
  746.       NBELEM=NBARETES
  747.       SEGADJ,IPTG
  748.       N1=ITARE1.LACT(/1)
  749.       N2=NBARETES
  750.       SEGADJ,ITARE1
  751. C
  752. C
  753. C     3) Menage avant de quitter
  754.       SEGSUP,LCCC1,ITL1,ITCEL1,ITCEL2,ILEA1,ITRC1
  755. 999   RETURN
  756.       END
  757.  
  758.  
  759.  
  760.  
  761.  
  762.  
  763.  
  764.  

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