Télécharger la plaquette Cast3M


fr en

Télécharger delaun.eso

Retour à la liste

Numérotation des lignes : 

delaun

  1. C DELAUN    SOURCE    CB215821  20/11/25    13:24:28     10792          
  2.       SUBROUTINE DELAUN(MPOVAL,IPT1,XBOI,IBOI,IPT2,LNBOIT,ITRC1,ILEA1)
  3. C----------------------------------------------------------------------C
  4. C     Triangulation d'un maillage de points (POI1), 1D, 2D ou 3D,      C
  5. C     suivant l'agorithme de DELAUNAY donne dans P.-L. George,         C
  6. C     Delaunay triangulation and meshing, Hermes, 1998, p. 41.         C
  7. C                                                                      C
  8. C     IPT1 = pointeur sur le maillage de POI1 en entree ;              C
  9. C            suppose ACTIF en ENTREE et en SORTIE                      C
  10. C     XBOI = rapport entre le "diametre" du nuage de points (IPT1) et  C
  11. C            la taille de la boite de triangulation ;                  C
  12. C     IBOI = indicateur permettant d'indiquer si on veut "garder" la   C
  13. C            boite de triangulation :                                  C
  14. C            - IBOI = 0 : on ne la garde pas (defaut) ;                C
  15. C            - IBOI = 1 : on la garde ;                                C
  16. C                                                                      C
  17. C     IPT2   = pointeur sur le maillage triangule en sortie            C
  18. C     LNBOIT = liste des noeuds formant la boite de triangulation      C
  19. C     ITRC1  = table contenant les centres des cercles circonscrits    C
  20. C     ILEA1  = table des d'adjacence de la triangulation               C
  21. C                                                                      C
  22. C----------------------------------------------------------------------C
  23. C Version V0 : centre et rayon cercles circons. calcules chaque fois   C
  24. C Version V1 : centre et rayon cercles circons. stockes dans SEGMENT   C
  25. C              MCIRCONS + elements adjacents stockes dans SEGMENT      C
  26. C              ILEA1                                                   C
  27. C----------------------------------------------------------------------C
  28.       IMPLICIT INTEGER(I-N)
  29.       IMPLICIT REAL*8 (A-H,O-Z)
  30. c
  31. -INC SMCHPOI
  32.  
  33. -INC PPARAM
  34. -INC CCOPTIO
  35. -INC SMELEME
  36. -INC SMCOORD
  37. -INC CCGEOME
  38. -INC SMLENTI
  39. -INC SMTABLE
  40.       PARAMETER (ICMIN=1000)
  41.       PARAMETER (ZERO=1D-13)
  42.       DIMENSION LNBOIT(8)
  43.       DIMENSION XCMIMA(2,3),X0(4,3),XCMIL(3),XPI1(3),XCBJ(3),XG0(3)
  44.       LOGICAL BMOD,BHORS,BFF
  45. C
  46.       POINTEUR LCAVI.MLENTI,LBOUL.MLENTI,LBASE.MLENTI,LREF1.MLENTI
  47.       POINTEUR LEADJ.MLENTI,LNVADJ.MLENTI,LSUP.MLENTI,LSCP.MLENTI
  48. C
  49.       SEGMENT,MCIRCONS
  50.          REAL*8 TRC(NBE1,4)
  51.       ENDSEGMENT
  52.       POINTEUR ITRC1.MCIRCONS
  53. C
  54.       SEGMENT,MADJACEN
  55.          INTEGER LEAC(NBL1,IDIM+1,2)
  56.       ENDSEGMENT
  57.       POINTEUR ILEA1.MADJACEN
  58. C
  59.       SEGMENT,MLIAIS
  60.          INTEGER LIAS(NB1,NB1)
  61.       ENDSEGMENT
  62.       POINTEUR ITL1.MLIAIS
  63. C
  64.       SEGMENT RAYONOE
  65.          REAL*8 RAYN(NBNER)
  66.       ENDSEGMENT
  67.       POINTEUR IRN1.RAYONOE
  68. c
  69. C----------------------------------------------------------------------C
  70. C                               PARTIE 1                               C
  71. C          Construction de la boite englobant tous les points          C
  72. C----------------------------------------------------------------------C
  73. C
  74. C---- 1.1 Je determine la taille du maillage
  75.       SEGACT,MCOORD*mod
  76.       IDIMP1=IDIM+1
  77.       IPT2=IPT1
  78.       IF ((IBOI.EQ.1).AND.(IPT1.LISREF(/1).NE.0)) THEN
  79.          IPT2=IPT1.LISREF(1)
  80.          SEGACT,IPT2
  81.       ENDIF
  82.       IREF1=(IPT1.NUM(1,1)-1)*IDIMP1
  83.       DDMAX=0.D0
  84.       DO I=1,IDIM
  85.          XCMIMA(1,I)=XCOOR(IREF1+I)
  86.          XCMIMA(2,I)=XCOOR(IREF1+I)
  87.          DO J=2,IPT1.NUM(/2)
  88.             IREFJ=(IPT1.NUM(1,J)-1)*IDIMP1
  89.             XCNJ1=XCOOR(IREFJ+I)
  90.             IF (XCNJ1.LT.XCMIMA(1,I)) XCMIMA(1,I)=XCNJ1
  91.             IF (XCNJ1.GT.XCMIMA(2,I)) XCMIMA(2,I)=XCNJ1
  92.          ENDDO
  93.          IF ((IBOI.EQ.1).AND.(IPT1.LISREF(/1).NE.0)) THEN
  94.             DO J=1,IPT2.NUM(/2)
  95.                IREFJ=(IPT2.NUM(1,J)-1)*IDIMP1
  96.                XCNJ1=XCOOR(IREFJ+I)
  97.                IF (XCNJ1.LT.XCMIMA(1,I)) XCMIMA(1,I)=XCNJ1
  98.                IF (XCNJ1.GT.XCMIMA(2,I)) XCMIMA(2,I)=XCNJ1
  99.             ENDDO
  100.          ENDIF
  101.          DDMAXI=0.5D0*(XCMIMA(2,I)-XCMIMA(1,I))
  102.          IF (DDMAXI.GT.DDMAX) DDMAX=DDMAXI
  103.       ENDDO
  104. C
  105. C---- 1.2 On construit la boite
  106. C     1.2.1 Points coins de la boite
  107.       DBOIT=XBOI*DDMAX
  108.       NBPTS0=nbpts
  109.       NBPTS=NBPTS0+(2**IDIM)
  110.       SEGADJ,MCOORD
  111.       IF (IDIM.EQ.1) THEN
  112.          XCMIL(1)=0.5D0*(XCMIMA(1,1)+XCMIMA(2,1))
  113.          IREF1=(NBPTS-2)*IDIMP1
  114.          IREF2=(NBPTS-1)*IDIMP1
  115.          XCOOR(IREF1+1)=XCMIL(1)-DBOIT
  116.          XCOOR(IREF2+1)=XCMIL(1)+DBOIT
  117.          XCOOR(IREF1+2)=0.D0
  118.          XCOOR(IREF2+2)=0.D0
  119.       ELSEIF (IDIM.EQ.2) THEN
  120.          DO I=1,IDIM
  121.             XCMIL(I)=0.5D0*(XCMIMA(1,I)+XCMIMA(2,I))
  122.          ENDDO
  123.          IREF1=(NBPTS-4)*IDIMP1
  124.          IREF2=(NBPTS-3)*IDIMP1
  125.          IREF3=(NBPTS-2)*IDIMP1
  126.          IREF4=(NBPTS-1)*IDIMP1
  127.          XCOOR(IREF1+1)=XCMIL(1)-DBOIT
  128.          XCOOR(IREF2+1)=XCMIL(1)+DBOIT
  129.          XCOOR(IREF3+1)=XCMIL(1)+DBOIT
  130.          XCOOR(IREF4+1)=XCMIL(1)-DBOIT
  131.          XCOOR(IREF1+2)=XCMIL(2)-DBOIT
  132.          XCOOR(IREF2+2)=XCMIL(2)-DBOIT
  133.          XCOOR(IREF3+2)=XCMIL(2)+DBOIT
  134.          XCOOR(IREF4+2)=XCMIL(2)+DBOIT
  135.          XCOOR(IREF1+3)=0.D0
  136.          XCOOR(IREF2+3)=0.D0
  137.          XCOOR(IREF3+3)=0.D0
  138.          XCOOR(IREF4+3)=0.D0
  139.       ELSEIF (IDIM.EQ.3) THEN
  140.          DO I=1,IDIM
  141.             XCMIL(I)=0.5*(XCMIMA(1,I)+XCMIMA(2,I))
  142.          ENDDO
  143.          IREF1=(NBPTS-8)*IDIMP1
  144.          IREF2=(NBPTS-7)*IDIMP1
  145.          IREF3=(NBPTS-6)*IDIMP1
  146.          IREF4=(NBPTS-5)*IDIMP1
  147.          IREF5=(NBPTS-4)*IDIMP1
  148.          IREF6=(NBPTS-3)*IDIMP1
  149.          IREF7=(NBPTS-2)*IDIMP1
  150.          IREF8=(NBPTS-1)*IDIMP1
  151.          XCOOR(IREF1+1)=XCMIL(1)-DBOIT
  152.          XCOOR(IREF2+1)=XCMIL(1)+DBOIT
  153.          XCOOR(IREF3+1)=XCMIL(1)+DBOIT
  154.          XCOOR(IREF4+1)=XCMIL(1)-DBOIT
  155.          XCOOR(IREF5+1)=XCMIL(1)-DBOIT
  156.          XCOOR(IREF6+1)=XCMIL(1)+DBOIT
  157.          XCOOR(IREF7+1)=XCMIL(1)+DBOIT
  158.          XCOOR(IREF8+1)=XCMIL(1)-DBOIT
  159.          XCOOR(IREF1+2)=XCMIL(2)-DBOIT
  160.          XCOOR(IREF2+2)=XCMIL(2)-DBOIT
  161.          XCOOR(IREF3+2)=XCMIL(2)+DBOIT
  162.          XCOOR(IREF4+2)=XCMIL(2)+DBOIT
  163.          XCOOR(IREF5+2)=XCMIL(2)-DBOIT
  164.          XCOOR(IREF6+2)=XCMIL(2)-DBOIT
  165.          XCOOR(IREF7+2)=XCMIL(2)+DBOIT
  166.          XCOOR(IREF8+2)=XCMIL(2)+DBOIT
  167.          XCOOR(IREF1+3)=XCMIL(3)-DBOIT
  168.          XCOOR(IREF2+3)=XCMIL(3)-DBOIT
  169.          XCOOR(IREF3+3)=XCMIL(3)-DBOIT
  170.          XCOOR(IREF4+3)=XCMIL(3)-DBOIT
  171.          XCOOR(IREF5+3)=XCMIL(3)+DBOIT
  172.          XCOOR(IREF6+3)=XCMIL(3)+DBOIT
  173.          XCOOR(IREF7+3)=XCMIL(3)+DBOIT
  174.          XCOOR(IREF8+3)=XCMIL(3)+DBOIT
  175.          XCOOR(IREF1+4)=0.D0
  176.          XCOOR(IREF2+4)=0.D0
  177.          XCOOR(IREF3+4)=0.D0
  178.          XCOOR(IREF4+4)=0.D0
  179.          XCOOR(IREF5+4)=0.D0
  180.          XCOOR(IREF6+4)=0.D0
  181.          XCOOR(IREF7+4)=0.D0
  182.          XCOOR(IREF8+4)=0.D0
  183.       ENDIF
  184. C
  185. C     1.2.2 Maillage de la boite
  186.       NBSOUS=0
  187.       NBREF=0
  188.       NNOE1=NBPTS-(2**IDIM)+1
  189.       IF (IDIM.EQ.1) THEN
  190. C        Creation des elements SEG2 en definissant les noeuds
  191.          NBNN=2
  192.          NBELEM=IPT1.NUM(/2)+1
  193.          SEGINI,IPT2
  194.          IPT2.ITYPEL=2
  195.          NBELE2=1
  196.          IPT2.NUM(1,1)=NNOE1
  197.          IPT2.NUM(2,1)=NNOE1+1
  198. C        Liste des numeros des noeuds de la boite
  199.          NNBOIT=2
  200.          LNBOIT(1)=NNOE1
  201.          LNBOIT(2)=NNOE1+1
  202. C        Initialisation du segment ILEA1 contenant la table de
  203. C        connectivite des elements de IPT2
  204.          NBL1=NBELEM
  205.          SEGINI,ILEA1
  206.          ILEA1.LEAC(1,1,1)= 0
  207.          ILEA1.LEAC(1,1,2)= 0
  208.          ILEA1.LEAC(1,2,1)= 0
  209.          ILEA1.LEAC(1,2,2)= 0
  210.       ELSEIF (IDIM.EQ.2) THEN
  211. C        Creation des elements TRI3 en definissant les noeuds
  212.          NBNN=3
  213.          NBELEM=IPT1.NUM(/2)*3
  214.          SEGINI,IPT2
  215.          IPT2.ITYPEL=4
  216.          NBELE2=2
  217.          IPT2.NUM(1,1)=NNOE1
  218.          IPT2.NUM(2,1)=NNOE1+1
  219.          IPT2.NUM(3,1)=NNOE1+3
  220.          IPT2.NUM(1,2)=NNOE1+1
  221.          IPT2.NUM(2,2)=NNOE1+2
  222.          IPT2.NUM(3,2)=NNOE1+3
  223. C        Liste des numeros des noeuds de la boite
  224.          NNBOIT=4
  225.          LNBOIT(1)=NNOE1
  226.          LNBOIT(2)=NNOE1+1
  227.          LNBOIT(3)=NNOE1+2
  228.          LNBOIT(4)=NNOE1+3
  229. C        Initialisation du segment ILEA1 contenant la table de
  230. C        connectivite des elements de IPT2
  231.          NBL1=NBELEM
  232.          SEGINI,ILEA1
  233.          ILEA1.LEAC(1,1,1)= 2
  234.          ILEA1.LEAC(1,1,2)= 2
  235.          ILEA1.LEAC(1,2,1)= 0
  236.          ILEA1.LEAC(1,2,2)= 0
  237.          ILEA1.LEAC(1,3,1)= 0
  238.          ILEA1.LEAC(1,3,2)= 0
  239.          ILEA1.LEAC(2,1,1)= 0
  240.          ILEA1.LEAC(2,1,2)= 0
  241.          ILEA1.LEAC(2,2,1)= 1
  242.          ILEA1.LEAC(2,2,2)= 1
  243.          ILEA1.LEAC(2,3,1)= 0
  244.          ILEA1.LEAC(2,3,2)= 0
  245.       ELSEIF (IDIM.EQ.3) THEN
  246. C        Creation des elements TET4 en definissant les noeuds
  247.          NBNN=4
  248.          NBELEM=IPT1.NUM(/2)*7
  249.          SEGINI,IPT2
  250.          IPT2.ITYPEL=23
  251.          NBELE2=5
  252.          IPT2.NUM(1,1)=NNOE1
  253.          IPT2.NUM(2,1)=NNOE1+1
  254.          IPT2.NUM(3,1)=NNOE1+3
  255.          IPT2.NUM(4,1)=NNOE1+4
  256.          IPT2.NUM(1,2)=NNOE1+1
  257.          IPT2.NUM(2,2)=NNOE1+2
  258.          IPT2.NUM(3,2)=NNOE1+3
  259.          IPT2.NUM(4,2)=NNOE1+6
  260.          IPT2.NUM(1,3)=NNOE1+1
  261.          IPT2.NUM(2,3)=NNOE1+4
  262.          IPT2.NUM(3,3)=NNOE1+5
  263.          IPT2.NUM(4,3)=NNOE1+6
  264.          IPT2.NUM(1,4)=NNOE1+3
  265.          IPT2.NUM(2,4)=NNOE1+4
  266.          IPT2.NUM(3,4)=NNOE1+6
  267.          IPT2.NUM(4,4)=NNOE1+7
  268.          IPT2.NUM(1,5)=NNOE1+1
  269.          IPT2.NUM(2,5)=NNOE1+3
  270.          IPT2.NUM(3,5)=NNOE1+4
  271.          IPT2.NUM(4,5)=NNOE1+6
  272. C        Liste des numeros des noeuds de la boite
  273.          NNBOIT=8
  274.          LNBOIT(1)=NNOE1
  275.          LNBOIT(2)=NNOE1+1
  276.          LNBOIT(3)=NNOE1+2
  277.          LNBOIT(4)=NNOE1+3
  278.          LNBOIT(5)=NNOE1+4
  279.          LNBOIT(6)=NNOE1+5
  280.          LNBOIT(7)=NNOE1+6
  281.          LNBOIT(8)=NNOE1+7
  282. C        Initialisation du segment ILEA1 contenant la table de
  283. C        connectivite des elements de IPT2
  284.          NBL1=NBELEM
  285.          SEGINI,ILEA1
  286.          ILEA1.LEAC(1,1,1)= 5
  287.          ILEA1.LEAC(1,1,2)= 4
  288.          ILEA1.LEAC(1,2,1)= 0
  289.          ILEA1.LEAC(1,2,2)= 0
  290.          ILEA1.LEAC(1,3,1)= 0
  291.          ILEA1.LEAC(1,3,2)= 0
  292.          ILEA1.LEAC(1,4,1)= 0
  293.          ILEA1.LEAC(1,4,2)= 0
  294.          ILEA1.LEAC(2,1,1)= 0
  295.          ILEA1.LEAC(2,1,2)= 0
  296.          ILEA1.LEAC(2,2,1)= 5
  297.          ILEA1.LEAC(2,2,2)= 3
  298.          ILEA1.LEAC(2,3,1)= 0
  299.          ILEA1.LEAC(2,3,2)= 0
  300.          ILEA1.LEAC(2,4,1)= 0
  301.          ILEA1.LEAC(2,4,2)= 0
  302.          ILEA1.LEAC(3,1,1)= 0
  303.          ILEA1.LEAC(3,1,2)= 0
  304.          ILEA1.LEAC(3,2,1)= 0
  305.          ILEA1.LEAC(3,2,2)= 0
  306.          ILEA1.LEAC(3,3,1)= 5
  307.          ILEA1.LEAC(3,3,2)= 2
  308.          ILEA1.LEAC(3,4,1)= 0
  309.          ILEA1.LEAC(3,4,2)= 0
  310.          ILEA1.LEAC(4,1,1)= 0
  311.          ILEA1.LEAC(4,1,2)= 0
  312.          ILEA1.LEAC(4,2,1)= 0
  313.          ILEA1.LEAC(4,2,2)= 0
  314.          ILEA1.LEAC(4,3,1)= 0
  315.          ILEA1.LEAC(4,3,2)= 0
  316.          ILEA1.LEAC(4,4,1)= 5
  317.          ILEA1.LEAC(4,4,2)= 1
  318.          ILEA1.LEAC(5,1,1)= 4
  319.          ILEA1.LEAC(5,1,2)= 4
  320.          ILEA1.LEAC(5,2,1)= 3
  321.          ILEA1.LEAC(5,2,2)= 3
  322.          ILEA1.LEAC(5,3,1)= 2
  323.          ILEA1.LEAC(5,3,2)= 2
  324.          ILEA1.LEAC(5,4,1)= 1
  325.          ILEA1.LEAC(5,4,2)= 1
  326.       ENDIF
  327. C
  328. C---- 1.3 Initialisation du segment contenant les rayons et centres
  329. C         de la triangulation :
  330.       IF (IDIM.NE.1) THEN
  331.          NBE1=IPT2.NUM(/2)
  332.          SEGINI,ITRC1
  333.          NBNER  = LNBOIT(NNBOIT)
  334.          SEGINI,IRN1
  335. C
  336.          if (MPOVAL .ne. 0) then
  337.             DO I=1,IPT1.NUM(/2)
  338.                IRN1.RAYN(IPT1.NUM(1,I))= MPOVAL.VPOCHA(I,1)
  339.             END DO
  340.          endif
  341. C
  342.          DO I0=1,NBELE2
  343. C           BCIRCO renvoie le centre et le rayon de l element considere
  344. C
  345.             CALL BCIRCO(IRN1,IPT2,I0,XCBJ,XRBJ)
  346. C
  347.             ITRC1.TRC(I0,1)=XRBJ
  348.             DO J0=1,IDIM
  349.                ITRC1.TRC(I0,J0+1)=XCBJ(J0)
  350.             ENDDO
  351.          ENDDO
  352.       ENDIF
  353. C
  354. C
  355. C
  356. C
  357. C----------------------------------------------------------------------C
  358. C                               PARTIE 2                               C
  359. C                      Triangulation incrementale                      C
  360. C----------------------------------------------------------------------C
  361. C
  362.       JG=1
  363.       SEGINI,LSUP
  364. C     On boucle sur tous les points a trianguler
  365.       NBSOUS=0
  366.       NBREF=0
  367.       NBNN=IPT2.NUM(/1)
  368.       DO I=1,IPT1.NUM(/2)
  369. C         WRITE(6,*)
  370. C         WRITE(6,*) 'Iteration I=',I,' /',IPT1.NUM(/2)
  371. C
  372. C------- 2.1 Construction de la base
  373. C
  374. C        Je tire un point a ajouter a la triangulation
  375.          NPI1=IPT1.NUM(1,I)
  376. C        XPI1 : vecteur contenant les coordonnees du point
  377.          DO J=1,IDIM
  378.             XPI1(J)=XCOOR((NPI1-1)*IDIMP1+J)
  379.          ENDDO
  380. C
  381. C        2.1.1 Recherche d'un element de la base (NBASE1)
  382. C        Je demarre avec le dernier elem. de IPT2 (NEL0) et cherche si
  383. C        NPI1 est dedans. Si oui, c'est gagne, sinon :
  384. C        je determine le voisin de NEL0 oriente vers NPI1
  385. C        => nouveau NEL0, puis je recommence jusqu'a y arriver
  386.          NEL0=NBELE2
  387.          NBASE=1
  388. C        Cas particulier de la dimension 1
  389.          IF (IDIM.EQ.1) THEN
  390.             DO J=1,NBELE2
  391. C              Coordonnees des noeuds A et B de l'element NEL0
  392. C              A est a gauche et B est a droite
  393.                NREFA=IPT2.NUM(1,NEL0)
  394.                NREFB=IPT2.NUM(2,NEL0)
  395.                XA=XCOOR((NREFA-1)*IDIMP1+1)
  396.                XB=XCOOR((NREFB-1)*IDIMP1+1)
  397. C              Cas ou NPI1 est confondu avec A ou B
  398.                X1=ABS(XPI1(1))
  399.                IF (X1.LT.ZERO) X1=1.
  400.                DTEST1=ABS(XPI1(1)-XA)/X1
  401.                DTEST2=ABS(XPI1(1)-XB)/X1
  402.                IF ((DTEST1.LT.ZERO).OR.(DTEST2.LT.ZERO)) THEN
  403.                   NBNN=IPT2.NUM(/1)
  404.                   NBELEM=NBELE2
  405.                   NBSOUS=0
  406.                   NBREF=0
  407.                   SEGADJ,IPT2
  408. C                 Erreur de triangulation : vérifiez qu'il n'y a pas de
  409. C                 points confondus
  410.                   CALL ERREUR(846)
  411.                   GOTO 2501
  412. C              Cas ou le NPT1 est entre A et B, on a trouve la base
  413.                ELSEIF ((XPI1(1).GT.XA).AND.(XPI1(1).LT.XB)) THEN
  414.                   NBASE1=NEL0
  415.                   GOTO 2381
  416. C              Si NPI1 est a gauche de A, on va chercher le voisin
  417.                ELSEIF (XPI1(1).LT.XA) THEN
  418.                   NEL0=ILEA1.LEAC(NEL0,2,1)
  419. C              Si NPI1 est a droite de B, on va chercher le voisin
  420.                ELSEIF (XPI1(1).GT.XB) THEN
  421.                   NEL0=ILEA1.LEAC(NEL0,1,1)
  422.                ENDIF
  423.             ENDDO
  424. C           Si l'on passe ici, c'est que l'on a pas trouve NBASE1 !
  425.             NBNN=IPT2.NUM(/1)
  426.             NBELEM=NBELE2
  427.             NBSOUS=0
  428.             NBREF=0
  429.             SEGADJ,IPT2
  430. C           Erreur de triangulation : vérifiez qu'il n'y a pas de
  431. C           points confondus
  432.             CALL ERREUR(846)
  433.             GOTO 2501
  434.          ENDIF
  435. C        Cas generaux en dimension 2 et 3
  436.          JG=1
  437.          SEGINI,LSCP
  438.          NEL00=NEL0
  439.          NEL000=NEL0
  440.          DO J=1,NBELE2
  441.             NCP=0
  442. C           Coordonnees des noeuds de l'element NEL0
  443.             DO K=1,IDIMP1
  444.                NREF0=IPT2.NUM(K,NEL0)
  445.                DO KK=1,IDIM
  446.                   X0(K,KK)=XCOOR((NREF0-1)*IDIMP1+KK)
  447.                ENDDO
  448.             ENDDO
  449. C           Coordonnees du barycentre G de NEL0
  450.             DO KK=1,IDIM
  451.                XSOM=0.D0
  452.                DO K=1,IDIMP1
  453.                   XSOM=XSOM+X0(K,KK)
  454.                ENDDO
  455.                XG0(KK)=XSOM/(IDIMP1)
  456.             ENDDO
  457.             IK1=0
  458.             IK2=0
  459.             IK3=0
  460. C           En dimension 2
  461.             IF (IDIM.EQ.2) THEN
  462. C              On effectue le meme test sur chaque cote
  463.                DO K=1,3
  464.                   IF (K.EQ.1) THEN
  465.                      NA=2
  466.                      NB=3
  467.                   ELSEIF (K.EQ.2) THEN
  468.                      NA=1
  469.                      NB=3
  470.                   ELSEIF (K.EQ.3) THEN
  471.                      NA=1
  472.                      NB=2
  473.                   ENDIF
  474. C                 Normale au cote AB du triangle
  475.                   XN=X0(NA,2)-X0(NB,2)
  476.                   YN=X0(NB,1)-X0(NA,1)
  477.                   XXN=SQRT((XN**2)+(YN**2))
  478.                   XN=XN/XXN
  479.                   YN=YN/XXN
  480. C                 On teste si NPI1 est du meme cote que G par rapport a
  481. C                 la droite (AB)
  482.                   PS1=(XN*(XPI1(1)-X0(NA,1)))+(YN*(XPI1(2)-X0(NA,2)))
  483.                   XX0=SQRT(((XPI1(1)-X0(NA,1))**2)+
  484.      &                 ((XPI1(2)-X0(NA,2))**2))
  485.                   PS11=PS1/XX0
  486.                   IF (ABS(PS11).LT.ZERO) THEN
  487.                      NSCP=K
  488.                      NCP=1
  489.                      NBASE=2
  490.                      GOTO 2111
  491.                   ENDIF
  492.                   PS2=(XN*(XG0(1)-X0(NA,1)))+(YN*(XG0(2)-X0(NA,2)))
  493.                   IF ((SIGN(1D0,PS1)).EQ.(SIGN(1D0,PS2))) THEN
  494.                      GOTO 2111
  495.                   ELSE
  496.                      NADJ0=K
  497.                      GOTO 2112
  498.                   ENDIF
  499. 2111              CONTINUE
  500.                   IF (K.EQ.1) IK1=1
  501.                   IF (K.EQ.2) IK2=1
  502.                   IF (K.EQ.3) IK3=1
  503.                ENDDO
  504.                IF ((IK1.EQ.1).AND.(IK2.EQ.1).AND.(IK3.EQ.1)) THEN
  505. C                 NPI1 est dans NEL0, on a trouve un element de la base
  506.                   NBASE1=NEL0
  507.                   GOTO 2115
  508.                ENDIF
  509. 2112           CONTINUE
  510. C              NPI1 et G ne sont pas du meme cote par rapport au cote
  511. C              vu par le noeud NADJ0 de NEL0 => on passe a l'element
  512. C              voisin NEL1
  513.                IF (J.GE.3) NEL000=NEL00
  514.                IF (J.GE.2) NEL00=NEL0
  515.                NEL0=ILEA1.LEAC(NEL0,NADJ0,1)
  516.                IF (NEL0.EQ.NEL00) THEN
  517. C                 PRINT*,'Erreur dans la table des elements adjacents'
  518. C                 PRINT*,'L''element',NEL0,'est voisin de lui meme !'
  519.                   CALL ERREUR(846)
  520.                   GOTO 2501
  521.                ENDIF
  522.                IF (NEL0.EQ.NEL000) THEN
  523. C                 PRINT*,'Erreur dans la construction du maillage'
  524. C                 PRINT*,'Les elements',NEL0,NEL00,'se chevauchent !'
  525.                   CALL ERREUR(846)
  526.                   GOTO 2501
  527.                ENDIF
  528.             ENDIF
  529. C           En dimension 3
  530.             IF (IDIM.EQ.3) THEN
  531.                IK4=0
  532. C              On effectue le test sur chaque face
  533.                DO K=1,4
  534.                   IF (K.EQ.1) THEN
  535.                      NA=2
  536.                      NB=3
  537.                      NC=4
  538.                   ELSEIF (K.EQ.2) THEN
  539.                      NA=1
  540.                      NB=3
  541.                      NC=4
  542.                   ELSEIF (K.EQ.3) THEN
  543.                      NA=1
  544.                      NB=2
  545.                      NC=4
  546.                   ELSEIF (K.EQ.4) THEN
  547.                      NA=1
  548.                      NB=2
  549.                      NC=3
  550.                   ENDIF
  551. C                 Normale au plan ABC du triangle (n = AB^AC)
  552.                   XN=((X0(NB,2)-X0(NA,2))*(X0(NC,3)-X0(NA,3))) -
  553.      &               ((X0(NB,3)-X0(NA,3))*(X0(NC,2)-X0(NA,2)))
  554.                   YN=((X0(NB,3)-X0(NA,3))*(X0(NC,1)-X0(NA,1))) -
  555.      &               ((X0(NB,1)-X0(NA,1))*(X0(NC,3)-X0(NA,3)))
  556.                   ZN=((X0(NB,1)-X0(NA,1))*(X0(NC,2)-X0(NA,2))) -
  557.      &               ((X0(NB,2)-X0(NA,2))*(X0(NC,1)-X0(NA,1)))
  558.                   XXN=SQRT((XN**2)+(YN**2)+(ZN**2))
  559.                   XN=XN/XXN
  560.                   YN=YN/XXN
  561.                   ZN=ZN/XXN
  562. C                 On teste si NPI1 est du meme cote que G par rapport au
  563. C                 plan ABC
  564.                   PS1=(XN*(XPI1(1)-X0(NA,1))) +
  565.      &                (YN*(XPI1(2)-X0(NA,2))) +
  566.      &                (ZN*(XPI1(3)-X0(NA,3)))
  567.                   XX0=SQRT(((XPI1(1)-X0(NA,1))**2) +
  568.      &                     ((XPI1(2)-X0(NA,2))**2)+
  569.      &                     ((XPI1(3)-X0(NA,3))**2))
  570.                   PS11=PS1/XX0
  571.                   IF (ABS(PS11).LT.ZERO) THEN
  572.                      NCP=NCP+1
  573.                      JG=NCP
  574.                      SEGADJ,LSCP
  575.                      LSCP.LECT(NCP)=K
  576.                      GOTO 2113
  577.                   ENDIF
  578.                   PS2=(XN*(XG0(1)-X0(NA,1))) +
  579.      &                (YN*(XG0(2)-X0(NA,2))) +
  580.      &                (ZN*(XG0(3)-X0(NA,3)))
  581.                   IF ((SIGN(1D0,PS1)).EQ.(SIGN(1D0,PS2))) THEN
  582.                      GOTO 2113
  583.                   ELSE
  584.                      GOTO 2114
  585.                   ENDIF
  586. 2113              CONTINUE
  587.                   IF (K.EQ.1) IK1=1
  588.                   IF (K.EQ.2) IK2=1
  589.                   IF (K.EQ.3) IK3=1
  590.                   IF (K.EQ.4) IK4=1
  591.                ENDDO
  592.                IF ((IK1.EQ.1).AND.(IK2.EQ.1).AND.
  593.      &             (IK3.EQ.1).AND.(IK4.EQ.1)) THEN
  594. C                 NPI1 est bien dans NEL0, on a donc trouve la base
  595.                   NBASE1=NEL0
  596.                   GOTO 2115
  597.                ENDIF
  598. 2114           CONTINUE
  599. C              NPI1 et G ne sont pas du meme cote par rapport au plan
  600. C              vu par le noeud NADJ0 de NEL0 => on passe a l'element
  601. C              adjacent a NEL0 vu par le noeud K
  602.                IF (J.GE.3) NEL000=NEL00
  603.                IF (J.GE.2) NEL00=NEL0
  604.                NEL0=ILEA1.LEAC(NEL0,K,1)
  605.                IF (NEL0.EQ.NEL00) THEN
  606. C                 PRINT*,'Erreur dans la table des elements adjacents'
  607. C                 PRINT*,'L''element',NEL0,'est voisin de lui meme !'
  608.                   CALL ERREUR(846)
  609.                   GOTO 2501
  610.                ENDIF
  611.                IF (NEL0.EQ.NEL000) THEN
  612. C                 PRINT*,'Erreur dans la construction du maillage'
  613. C                 PRINT*,'Les elements',NEL0,NEL00,'se chevauchent !'
  614.                   CALL ERREUR(846)
  615.                   GOTO 2501
  616.                ENDIF
  617.             ENDIF
  618.          ENDDO
  619. C        Si l'on passe ici, c'est que l'on a pas trouve NBASE1 !
  620.          NBNN=IPT2.NUM(/1)
  621.          NBELEM=NBELE2
  622.          NBSOUS=0
  623.          NBREF=0
  624.          SEGADJ,IPT2
  625. C        Erreur de triangulation : vérifiez qu'il n'y a pas de points
  626. C        confondus
  627.          CALL ERREUR(846)
  628.          GOTO 2501
  629. 2115     CONTINUE
  630. C        Test si NPI1 n'est pas confondu avec un des noeuds de NBASE1
  631.          DO J=1,IDIMP1
  632.             NREFJ=IPT2.NUM(J,NBASE1)
  633.             NC=0
  634.             DO JJ=1,IDIM
  635.                X1=ABS(XPI1(JJ))
  636.                IF (X1.LT.ZERO) X1=1.
  637.                DTEST=ABS(XPI1(JJ)-XCOOR((NREFJ-1)*IDIMP1+JJ))/X1
  638.                IF (DTEST.LT.ZERO) NC=NC+1
  639.             ENDDO
  640.             IF (NC.EQ.IDIM) THEN
  641.                NBNN=IPT2.NUM(/1)
  642.                NBELEM=NBELE2
  643.                NBSOUS=0
  644.                NBREF=0
  645.                SEGADJ,IPT2
  646. C              Erreur de triangulation : vérifiez qu'il n'y a pas de
  647. C              points confondus
  648.                CALL ERREUR(846)
  649.                GOTO 2501
  650.             ENDIF
  651.          ENDDO
  652. C
  653. C        2.1.2 Ajout eventuel d'autres elements a la base
  654. C        LBASE : liste des numeros des elements de IPT2 appartenant a
  655. C                la base, taille initialisee a NBASE
  656.          JG=NBASE
  657.          SEGINI,LBASE
  658.          LBASE.LECT(1)=NBASE1
  659. C        Ajout des elments supplementaires a LBASE si necessaire
  660.          IF (NCP.NE.0) THEN
  661. C           En dimension 2, un seul adjacent est ajoute a LBASE, celui
  662. C           partageant le cote
  663.             IF (IDIM.EQ.2) THEN
  664.                LBASE.LECT(2)=ILEA1.LEAC(NBASE1,NSCP,1)
  665. C           En dimension 3 plusieurs cas sont possibles
  666.             ELSEIF (IDIM.EQ.3) THEN
  667. C              Cas ou NPI1 est sur un des 4 plans
  668.                IF (NCP.EQ.1) THEN
  669. C                 On recupere l'adjacent partageant le plan
  670.                   NBASE=2
  671.                   JG=NBASE
  672.                   SEGADJ,LBASE
  673.                   LBASE.LECT(2)=ILEA1.LEAC(NBASE1,LSCP.LECT(1),1)
  674. C              Cas ou NPI1 est sur une des aretes
  675.                ELSEIF (NCP.EQ.2) THEN
  676. C                 On recupere tous les elements partageant l'arete
  677.                   NA=1
  678.                   NB=2
  679.                   NC=LSCP.LECT(1)
  680.                   ND=LSCP.LECT(2)
  681.                   NREFC=IPT2.NUM(NC,NBASE1)
  682.                   NREFD=IPT2.NUM(ND,NBASE1)
  683.                   DO L=1,4
  684.                      IF ((NA.EQ.NC).OR.(NA.EQ.ND).OR.(NA.EQ.NB)) NA=NA+1
  685.                      IF ((NB.EQ.NC).OR.(NB.EQ.ND).OR.(NB.EQ.NA)) NB=NB+1
  686.                   ENDDO
  687.                   NREFA=IPT2.NUM(NA,NBASE1)
  688.                   NREFB=IPT2.NUM(NB,NBASE1)
  689.                   N1=NBASE1
  690.                   DO L=1,100
  691.                      N2=ILEA1.LEAC(N1,NC,1)
  692.                      ND=ILEA1.LEAC(N1,NC,2)
  693.                      IF (N2.EQ.NBASE1) GOTO 2121
  694.                      NBASE=NBASE+1
  695.                      JG=NBASE
  696.                      SEGADJ,LBASE
  697.                      LBASE.LECT(NBASE)=N2
  698.                      DO M=1,4
  699.                         IF (IPT2.NUM(M,N2).EQ.NREFA) NA=M
  700.                         IF (IPT2.NUM(M,N2).EQ.NREFB) NB=M
  701.                         IF (IPT2.NUM(M,N2).EQ.NREFD) NC=M
  702.                      ENDDO
  703.                      NREFD=IPT2.NUM(ND,N2)
  704.                      N1=N2
  705.                   ENDDO
  706. 2121              CONTINUE
  707. C              Cas ou NPI1 est confondu avec un des sommets
  708.                ELSEIF (NCP.EQ.3) THEN
  709.                   NBNN=IPT2.NUM(/1)
  710.                   NBELEM=NBELE2
  711.                   NBSOUS=0
  712.                   NBREF=0
  713.                   SEGADJ,IPT2
  714. C                 Erreur de triangulation : vérifiez qu'il n'y a pas de
  715. C                 points confondus
  716.                   CALL ERREUR(846)
  717.                   GOTO 2501
  718.                ENDIF
  719.             ENDIF
  720.          ENDIF
  721. C
  722. C------- 2.2 Construction de la cavite
  723. C
  724. C        2.2.1 Construction de la cavite (MELEME IPT3)
  725. C        par adjacence a partir de l'element NBASE1 de la base
  726. C
  727. C        LCAVI : liste des numeros des elements de IPT2 appartenant a
  728. C                la cavite, taille initialisee a NBELE2
  729.          JG=NBELE2
  730.          SEGINI,LCAVI
  731. C        Le premier element de la cavite est NBASE1
  732.          LCAVI.LECT(1)=NBASE1
  733. C        NCAVI : nombre d'elements dans la cavite (toujours a jour !)
  734.          NCAVI=1
  735. C        Initialisation de la cavite IPT3 a NBELE2
  736.          NBELEM=NBELE2
  737.          NBNN=IPT2.NUM(/1)
  738.          SEGINI,IPT3
  739.          IPT3.ITYPEL=IPT2.ITYPEL
  740.          DO K=1,IPT3.NUM(/1)
  741.             IPT3.NUM(K,NCAVI)=IPT2.NUM(K,NBASE1)
  742.          ENDDO
  743. C        Parcours des adjacents de maniere recursive jusqu'a stabilite
  744.          ITRA=1
  745.          DO J=1,NBELE2
  746. C           Element a traiter de la liste LCAVI
  747.             NEJ=LCAVI.LECT(ITRA)
  748.             IF ((NEJ.EQ.0).OR.(J.EQ.NBELE2)) GOTO 2213
  749. C           On prend les elements adjacents de NEJ
  750.             DO K=1,IDIMP1
  751.                NEAJ=ILEA1.LEAC(NEJ,K,1)
  752. C              On verifie qu'il y a bien un element adjacent
  753.                IF (NEAJ.EQ.0) GOTO 2212
  754. C              On verifie que NEAJ n'est pas deja dans LCAVI
  755.                DO L=1,LCAVI.LECT(/1)
  756.                   IF (LCAVI.LECT(L).EQ.0) GOTO 2211
  757.                   IF (LCAVI.LECT(L).EQ.NEAJ) GOTO 2212
  758.                ENDDO
  759. 2211           CONTINUE
  760. C              XRBJ    : rayon de la boule circonscrite a l'elem. NEAJ
  761. C              XCBJ(i) : ieme coordonnee de son centre
  762.                XRBJ=ITRC1.TRC(NEAJ,1)
  763.                DO L=1,IDIM
  764.                   XCBJ(L)=ITRC1.TRC(NEAJ,L+1)
  765.                ENDDO
  766. C              DIJ1 : distance de NPI1 au centre de la boule de
  767. C                     l'element considere
  768.                DIJ1=0.D0
  769.                DO L=1,IDIM
  770.                   DIJ1=DIJ1+((XPI1(L)-XCBJ(L))**2)
  771.                ENDDO
  772. C
  773.                IF (MPOVAL .NE. 0)
  774.      &         DIJ1 = DIJ1 - ((IRN1.RAYN(IPT1.NUM(1,I)))**2)
  775. C
  776.                DIJ1=(DIJ1/(XRBJ**2))-1.
  777. C              Si NEAJ est un nouvel element de la cavite,
  778.                IF (DIJ1.LT.ZERO) THEN
  779. C                 alors on l'ajoute a LCAVI et on incremente NCAVI
  780.                   NCAVI=NCAVI+1
  781.                   LCAVI.LECT(NCAVI)=NEAJ
  782. C                 et on remplis IPT3 avec ce nouvel element
  783.                   DO L=1,IPT3.NUM(/1)
  784.                      IPT3.NUM(L,NCAVI)=IPT2.NUM(L,NEAJ)
  785.                   ENDDO
  786.                ENDIF
  787. 2212           CONTINUE
  788.             ENDDO
  789.             ITRA=ITRA+1
  790.          ENDDO
  791. 2213     CONTINUE
  792. C        Ajustement de LCAVI et IPT3 a la bonne taille si besoin
  793.          IF (LCAVI.LECT(/1).GT.NCAVI) THEN
  794.             JG=NCAVI
  795.             SEGADJ,LCAVI
  796.             NBNN=IPT3.NUM(/1)
  797.             NBELEM=NCAVI
  798.             SEGADJ,IPT3
  799.          ENDIF
  800.          NBELE3=IPT3.NUM(/2)
  801. C
  802. C        2.2.2 Verification : la cavite doit etre connexe, chaque
  803. C        element doit avoir au moins un element adjacent, c'est a dire
  804. C        partageant un cote/face. De plus, un cote/face ne peut etre
  805. C        commun qu'a deux elements au plus.
  806. C        On determine egalement le contour/enveloppe, IPT4, de la cavite
  807.          NNBELE2=NBELE2
  808. 2221     CONTINUE
  809.          NBELE2=NNBELE2
  810.          NNCON=0
  811.          JG=IDIM
  812.          SEGINI,LREF1
  813.          NBNN=IPT3.NUM(/1)-1
  814.          NBELEM=NBELE3*3
  815.          SEGINI,IPT4
  816.          IF (IDIM.EQ.2) IPT4.ITYPEL=2
  817.          IF (IDIM.EQ.3) IPT4.ITYPEL=4
  818.          NBELE4=0
  819. C        On parcoure les elements de la cavite
  820.          DO I1=1,IPT3.NUM(/2)
  821. C           Nombre d'adjacents a l'element I1
  822. C           un element est adjacent s'il partage un cote/face
  823.             NVI1=0
  824. C           On parcoure les 3/4 cotes/faces de l'element I1
  825.             DO I2=1,IPT3.NUM(/1)
  826. C              Nombre d'adjacents a la face I2 (1 ou 0)
  827.                NVI2=0
  828.                IF (IDIM.EQ.2) THEN
  829.                   IF (I2.EQ.1) THEN
  830.                      LREF1.LECT(1)=IPT3.NUM(2,I1)
  831.                      LREF1.LECT(2)=IPT3.NUM(3,I1)
  832.                   ELSEIF (I2.EQ.2) THEN
  833.                      LREF1.LECT(1)=IPT3.NUM(1,I1)
  834.                      LREF1.LECT(2)=IPT3.NUM(3,I1)
  835.                   ELSEIF (I2.EQ.3) THEN
  836.                      LREF1.LECT(1)=IPT3.NUM(1,I1)
  837.                      LREF1.LECT(2)=IPT3.NUM(2,I1)
  838.                   ENDIF
  839.                ELSEIF (IDIM.EQ.3) THEN
  840.                   IF (I2.EQ.1) THEN
  841.                      LREF1.LECT(1)=IPT3.NUM(2,I1)
  842.                      LREF1.LECT(2)=IPT3.NUM(3,I1)
  843.                      LREF1.LECT(3)=IPT3.NUM(4,I1)
  844.                   ELSEIF (I2.EQ.2) THEN
  845.                      LREF1.LECT(1)=IPT3.NUM(1,I1)
  846.                      LREF1.LECT(2)=IPT3.NUM(3,I1)
  847.                      LREF1.LECT(3)=IPT3.NUM(4,I1)
  848.                   ELSEIF (I2.EQ.3) THEN
  849.                      LREF1.LECT(1)=IPT3.NUM(1,I1)
  850.                      LREF1.LECT(2)=IPT3.NUM(2,I1)
  851.                      LREF1.LECT(3)=IPT3.NUM(4,I1)
  852.                   ELSEIF (I2.EQ.4) THEN
  853.                      LREF1.LECT(1)=IPT3.NUM(1,I1)
  854.                      LREF1.LECT(2)=IPT3.NUM(2,I1)
  855.                      LREF1.LECT(3)=IPT3.NUM(3,I1)
  856.                   ENDIF
  857.                ENDIF
  858.                DO I11=1,IPT3.NUM(/2)
  859.                   NNCI11=0
  860.                   IF (I11.EQ.I1) GOTO 2222
  861.                   DO I22=1,IPT3.NUM(/1)
  862.                      NREFX=IPT3.NUM(I22,I11)
  863.                      DO K=1,IDIM
  864.                         IF (NREFX.EQ.LREF1.LECT(K)) THEN
  865.                            NNCI11=NNCI11+1
  866.                         ENDIF
  867.                      ENDDO
  868.                   ENDDO
  869.                   IF (NNCI11.EQ.IDIM) THEN
  870.                      NVI2=NVI2+1
  871.                      NVI1=NVI1+1
  872.                   ENDIF
  873. 2222              CONTINUE
  874.                ENDDO
  875. C              Cas ou I1 a plus de 1 voisin sur la face I2
  876.                IF (NVI2.GT.1) THEN
  877.                   NBNN=IPT2.NUM(/1)
  878.                   NBELEM=NBELE2
  879.                   NBSOUS=0
  880.                   NBREF=0
  881.                   SEGADJ,IPT2
  882. C                 Erreur de triangulation : vérifiez qu'il n'y a pas de
  883. C                 points confondus
  884.                   CALL ERREUR(846)
  885.                   GOTO 2501
  886.                ENDIF
  887. C              Si la face I2 n'est que sur l'element I1, alors on
  888. C              l'ajoute au maillage du contour/enceloppe de la cavite
  889.                IF (NVI2.EQ.0) THEN
  890.                   NBELE4=NBELE4+1
  891.                   IF (IPT4.NUM(/2).LT.NBELE4) THEN
  892.                      NBNN=IPT4.NUM(/1)
  893.                      NBELEM=IPT4.NUM(/2)+100
  894.                      SEGADJ,IPT4
  895.                   ENDIF
  896.                   DO K=1,IPT4.NUM(/1)
  897.                      IF (IDIM.EQ.2) THEN
  898.                         IF (I2.EQ.1) THEN
  899.                            IPT4.NUM(1,NBELE4)=IPT3.NUM(2,I1)
  900.                            IPT4.NUM(2,NBELE4)=IPT3.NUM(3,I1)
  901.                         ELSEIF (I2.EQ.2) THEN
  902.                            IPT4.NUM(1,NBELE4)=IPT3.NUM(1,I1)
  903.                            IPT4.NUM(2,NBELE4)=IPT3.NUM(3,I1)
  904.                         ELSEIF (I2.EQ.3) THEN
  905.                            IPT4.NUM(1,NBELE4)=IPT3.NUM(1,I1)
  906.                            IPT4.NUM(2,NBELE4)=IPT3.NUM(2,I1)
  907.                         ENDIF
  908.                      ELSEIF (IDIM.EQ.3) THEN
  909.                         IF (I2.EQ.1) THEN
  910.                            IPT4.NUM(1,NBELE4)=IPT3.NUM(2,I1)
  911.                            IPT4.NUM(2,NBELE4)=IPT3.NUM(3,I1)
  912.                            IPT4.NUM(3,NBELE4)=IPT3.NUM(4,I1)
  913.                         ELSEIF (I2.EQ.2) THEN
  914.                            IPT4.NUM(1,NBELE4)=IPT3.NUM(1,I1)
  915.                            IPT4.NUM(2,NBELE4)=IPT3.NUM(3,I1)
  916.                            IPT4.NUM(3,NBELE4)=IPT3.NUM(4,I1)
  917.                         ELSEIF (I2.EQ.3) THEN
  918.                            IPT4.NUM(1,NBELE4)=IPT3.NUM(1,I1)
  919.                            IPT4.NUM(2,NBELE4)=IPT3.NUM(2,I1)
  920.                            IPT4.NUM(3,NBELE4)=IPT3.NUM(4,I1)
  921.                         ELSEIF (I2.EQ.4) THEN
  922.                            IPT4.NUM(1,NBELE4)=IPT3.NUM(1,I1)
  923.                            IPT4.NUM(2,NBELE4)=IPT3.NUM(2,I1)
  924.                            IPT4.NUM(3,NBELE4)=IPT3.NUM(3,I1)
  925.                         ENDIF
  926.                      ENDIF
  927.                   ENDDO
  928.                ENDIF
  929.             ENDDO
  930. C           Cas ou I1 n'a aucun voisin => cavite non connexe, on retire
  931. C           l'element de la cavite
  932.             IF (NVI1.EQ.0) THEN
  933.                IF (NBELE3.NE.1) THEN
  934.                   NNCON=1
  935.                   NESUP=LCAVI.LECT(I1)
  936.                ENDIF
  937.             ENDIF
  938.          ENDDO
  939.          NBNN=IPT4.NUM(/1)
  940.          NBELEM=NBELE4
  941.          SEGADJ,IPT4
  942.          IF (NNCON.EQ.1) GOTO 2341
  943. C
  944. C        2.2.3 Verification : la cavite et son contour/enveloppe doivent
  945. C        avoir le meme nombre de noeuds
  946.          CALL ECROBJ('MAILLAGE',IPT3)
  947.          IF (IERR.NE.0) RETURN
  948.          CALL NBNO
  949.          SEGACT,IPT3
  950.          CALL LIRENT(NBNO3,1,IRETOU)
  951.          CALL ECROBJ('MAILLAGE',IPT4)
  952.          IF (IERR.NE.0) RETURN
  953.          CALL NBNO
  954.          SEGACT,IPT4
  955.          CALL LIRENT(NBNO4,1,IRETOU)
  956.          IF (NBNO3.NE.NBNO4) THEN
  957.             NBNN=IPT2.NUM(/1)
  958.             NBELEM=NBELE2
  959.             NBSOUS=0
  960.             NBREF=0
  961.             SEGADJ,IPT2
  962. C           Erreur de triangulation : vérifiez qu'il n'y a pas de
  963. C           points confondus
  964.             CALL ERREUR(846)
  965.             GOTO 2501
  966.          ENDIF
  967. C
  968. C------- 2.3 Construction de la boule
  969. C
  970. C        2.3.1 Construction de la liste des elements de la boule LBOUL
  971. C        La boule est formee par les elements construits en joignant
  972. C        NPI1 aux elements du contour de la cavite.
  973.          NBOUL=NBELE4
  974.          JG=NBOUL
  975.          SEGINI,LBOUL
  976. C        Nouveau nombre d'elements de la triangulation
  977.          NEW2=NBELE2-NCAVI+NBOUL
  978. C        On regarde si le nombre d'elements dans la cavite est inferieur
  979. C        au nombre d'elements de son contour/enveloppe car il faudra un
  980. C        traitement special lors de la creation de la boule
  981.          N43=NBELE4-NBELE3
  982. C        Cas general ou l'on va ajouter de nouveaux elements a IPT2
  983.          IF (N43.GE.0) THEN
  984.             DO J=1,NBOUL
  985. C              Les premiers elements de LBOUL sont ceux de LCAVI
  986.                IF (J.LE.NCAVI) THEN
  987.                   LBOUL.LECT(J)=LCAVI.LECT(J)
  988. C              et je complete LBOUL avec les nouveaux numeros d'elements
  989.                ELSE
  990.                   LBOUL.LECT(J)=NBELE2+1
  991.                   NBELE2=NBELE2+1
  992.                ENDIF
  993.             ENDDO
  994. C           J'ajuste IPT2, ITRC1 et ILEA1 si besoin
  995.             IF (NEW2.GT.IPT2.NUM(/2)) THEN
  996.                NBNN=IPT2.NUM(/1)
  997.                NBELEM=NBELE2+ICMIN
  998.                SEGADJ,IPT2
  999.                NBE1=NBELEM
  1000.                SEGADJ,ITRC1
  1001.                NBL1=NBELEM
  1002.                SEGADJ,ILEA1
  1003.             ENDIF
  1004. C        Cas particulier ou l'on va enlever des elements a IPT2
  1005.          ELSE
  1006.             NSUP=NCAVI-NBOUL
  1007.             JG=NSUP
  1008.             SEGINI,LSUP
  1009. C           LBOUL est fait des NBOUL premiers elements de LCAVI
  1010.             DO J=1,NCAVI
  1011.                NELJ=LCAVI.LECT(J)
  1012.                IF (J.LE.NBOUL) THEN
  1013.                   LBOUL.LECT(J)=NELJ
  1014.                ELSE
  1015.                   LSUP.LECT(J-NBOUL)=NELJ
  1016. C                 Mise a zero des numeros des noeuds de l'element NELJ
  1017.                   DO K=1,IPT2.NUM(/1)
  1018.                      IPT2.NUM(K,NELJ)=0
  1019.                   ENDDO
  1020.                ENDIF
  1021.             ENDDO
  1022.          ENDIF
  1023.          NBELE2=NEW2
  1024. C
  1025. C        2.3.2 Je recherche les elements exterieurs a la cavite
  1026. C        partageant une face de l'enveloppe IPT4, pour cela on regarde
  1027. C        les elements de IPT3 et leurs adjacents grace a la table de
  1028. C        connectivite
  1029.          JG=NBELE4
  1030.          SEGINI,LEADJ
  1031.          SEGINI,LNVADJ
  1032.          JG=IDIM
  1033.          SEGINI,LREF1
  1034. C        On parcoure les elements II de IPT4 du contour/enveloppe
  1035.          DO II=1,NBELE4
  1036. C           Numeros des noeuds JJ de l'element II
  1037.             DO JJ=1,IDIM
  1038.                LREF1.LECT(JJ)=IPT4.NUM(JJ,II)
  1039.             ENDDO
  1040. C           On parcoure les elements J de IPT3 de la cavite
  1041.             DO J=1,NBELE3
  1042.                NCOM=0
  1043.                NSOM=0
  1044.                DO K=1,IDIMP1
  1045.                   NREFX=IPT3.NUM(K,J)
  1046. C                 Test si le noeud K de l'element J est commun a II
  1047.                   DO L=1,IDIM
  1048.                      IF (NREFX.EQ.LREF1.LECT(L)) THEN
  1049.                         NCOM=NCOM+1
  1050.                         NSOM=NSOM+K
  1051.                      ENDIF
  1052.                   ENDDO
  1053. C                 Si l'on a trouve l'element appuye sur la face II
  1054.                   IF (NCOM.EQ.IDIM) THEN
  1055.                      NELK2=LCAVI.LECT(J)
  1056.                      GOTO 2321
  1057.                   ENDIF
  1058.                ENDDO
  1059.             ENDDO
  1060. 2321        CONTINUE
  1061. C           Determination du numero du noeud voyeur dans IPT2
  1062.             IF (IDIM.EQ.2) NNOK2=6-NSOM
  1063.             IF (IDIM.EQ.3) NNOK2=10-NSOM
  1064.             LEADJ.LECT(II)=ILEA1.LEAC(NELK2,NNOK2,1)
  1065.             LNVADJ.LECT(II)=ILEA1.LEAC(NELK2,NNOK2,2)
  1066.          ENDDO
  1067. C
  1068. C        2.3.3 Verification : le contour/enveloppe de la cavite doit
  1069. C        etre etoile par rapport au point NPI1
  1070.          NESUP=0
  1071.          IF (IDIM.EQ.3) THEN
  1072.             DO J=1,NBELE4
  1073.                NA=IPT4.NUM(1,J)
  1074.                NB=IPT4.NUM(2,J)
  1075.                NC=IPT4.NUM(3,J)
  1076.                XA=XCOOR((NA-1)*(IDIM+1)+1)
  1077.                YA=XCOOR((NA-1)*(IDIM+1)+2)
  1078.                ZA=XCOOR((NA-1)*(IDIM+1)+3)
  1079.                XB=XCOOR((NB-1)*(IDIM+1)+1)
  1080.                YB=XCOOR((NB-1)*(IDIM+1)+2)
  1081.                ZB=XCOOR((NB-1)*(IDIM+1)+3)
  1082.                XC=XCOOR((NC-1)*(IDIM+1)+1)
  1083.                YC=XCOOR((NC-1)*(IDIM+1)+2)
  1084.                ZC=XCOOR((NC-1)*(IDIM+1)+3)
  1085. C              Normale au plan ABC : n = AB^AC
  1086.                XN=((YB-YA)*(ZC-ZA))-((ZB-ZA)*(YC-YA))
  1087.                YN=((ZB-ZA)*(XC-XA))-((XB-XA)*(ZC-ZA))
  1088.                ZN=((XB-XA)*(YC-YA))-((YB-YA)*(XC-XA))
  1089.                XXN=SQRT((XN**2)+(YN**2)+(ZN**2))
  1090.                XN=XN/XXN
  1091.                YN=YN/XXN
  1092.                ZN=ZN/XXN
  1093. C              Test si NPI1 est dans le plan ABC
  1094.                PS1=(XN*(XPI1(1)-XA))+(YN*(XPI1(2)-YA))+(ZN*(XPI1(3)-ZA))
  1095.                XX0=SQRT(((XPI1(1)-XA)**2)+((XPI1(2)-YA)**2)+
  1096.      &              ((XPI1(3)-ZA)**2))
  1097.                PS11=PS1/XX0
  1098.                IF (ABS(PS11).LT.ZERO) THEN
  1099.                   NESUP=ILEA1.LEAC(LEADJ.LECT(J),LNVADJ.LECT(J),1)
  1100.                   GOTO 2341
  1101.                ENDIF
  1102. C              Point D voyant la face ABC mais situe hors de la cavite
  1103.                IF (LEADJ.LECT(J).EQ.0) GOTO 2331
  1104.                ND=IPT2.NUM(LNVADJ.LECT(J),LEADJ.LECT(J))
  1105.                XD=XCOOR((ND-1)*(IDIM+1)+1)
  1106.                YD=XCOOR((ND-1)*(IDIM+1)+2)
  1107.                ZD=XCOOR((ND-1)*(IDIM+1)+3)
  1108. C              Test si NPI1 est du meme cote que D par rapport au plan
  1109. C              ABC
  1110.                PS2=(XN*(XD-XA)) + (YN*(YD-YA)) + (ZN*(ZD-ZA))
  1111.                XX0=SQRT(((XD-XA)**2)+((YD-YA)**2)+((ZD-ZA)**2))
  1112.                PS22=PS2/XXN
  1113.                IF (ABS(PS22).LT.ZERO) PS2=0.D0
  1114.                IF ((SIGN(1D0,PS1)).EQ.(SIGN(1D0,PS2))) THEN
  1115.                   NESUP=ILEA1.LEAC(LEADJ.LECT(J),LNVADJ.LECT(J),1)
  1116.                   GOTO 2341
  1117.                ELSE
  1118.                   GOTO 2331
  1119.                ENDIF
  1120. 2331           CONTINUE
  1121.             ENDDO
  1122.          ENDIF
  1123. C
  1124. C        2.3.4 Suppression d'un element (NESUP) de la cavite
  1125. 2341     CONTINUE
  1126.          IF (NESUP.NE.0) THEN
  1127. C           Test si NESUP est dans la base, pas touche a la base !
  1128.             DO K=1,NBASE
  1129.                IF (LBASE.LECT(K).EQ.NESUP) THEN
  1130.                   NBNN=IPT2.NUM(/1)
  1131.                   NBELEM=NBELE2
  1132.                   NBSOUS=0
  1133.                   NBREF=0
  1134.                   SEGADJ,IPT2
  1135. C                 Erreur de triangulation : vérifiez qu'il n'y a pas de
  1136. C                 points confondus
  1137.                   CALL ERREUR(846)
  1138.                   GOTO 2501
  1139.                ENDIF
  1140.             ENDDO
  1141.             NBNN=IPT3.NUM(/1)
  1142.             NBELEM=NBELE3-1
  1143.             SEGINI,IPT5
  1144.             IPT5.ITYPEL=IPT3.ITYPEL
  1145. C           Suppression de l'element NESUP le la liste LCAVI
  1146.             NS=0
  1147.             IF (LCAVI.LECT(NCAVI).EQ.NESUP) THEN
  1148.                NCAVI=NCAVI-1
  1149.                NS=1
  1150.                DO JJ=1,NCAVI
  1151.                   DO KK=1,IPT3.NUM(/1)
  1152.                      IPT5.NUM(KK,JJ)=IPT3.NUM(KK,JJ)
  1153.                   ENDDO
  1154.                ENDDO
  1155.             ELSE
  1156.                DO K=1,NCAVI-1
  1157.                   IF (LCAVI.LECT(K).EQ.NESUP) THEN
  1158.                      NS=NS+1
  1159.                      NCAVI=NCAVI-1
  1160.                   ENDIF
  1161.                   LCAVI.LECT(K)=LCAVI.LECT(K+NS)
  1162.                   DO KK=1,IPT3.NUM(/1)
  1163.                      IPT5.NUM(KK,K)=IPT3.NUM(KK,K+NS)
  1164.                   ENDDO
  1165.                ENDDO
  1166.             ENDIF
  1167. C           Ajustement de la liste LCAVI
  1168.             JG=NCAVI
  1169.             SEGADJ,LCAVI
  1170.             SEGSUP,IPT3
  1171.             IPT3=IPT5
  1172.             NBELE3=NCAVI
  1173. C           On repart calculer le contour/enveloppe de la cavite
  1174.             GOTO 2221
  1175.          ENDIF
  1176. C
  1177. C        2.3.5 Je parcours les elements du contour/enveloppe pour
  1178. C        construire le maillage de la boule
  1179.          DO J=1,NBELE4
  1180.             JEL=LBOUL.LECT(J)
  1181. C           Attribution des premiers noeuds des nouveux elements
  1182. C           de IPT2 : on choisi ceux du contour IPT4
  1183.             DO K=1,IPT4.NUM(/1)
  1184.                IPT2.NUM(K,JEL)=IPT4.NUM(K,J)
  1185.             ENDDO
  1186. C           Attribution du dernier noeud des nouveaux elements de IPT2 :
  1187. C           il s'agit de NPI1 au entre de la cavite
  1188.             IPT2.NUM(IDIMP1,JEL)=NPI1
  1189. C           Mise a jour de la table de connectivite (partie 1 sur 2)
  1190. C           uniquement sur les derniers noeuds (IDIMP1) des elements de
  1191. C           la boule (adjacents extexrieurs a la boule).
  1192. C           Pour cela, on utilise LEADJ et LNVADJ calcules juste avant
  1193.             ILEA1.LEAC(JEL,IDIMP1,1)=LEADJ.LECT(J)
  1194.             ILEA1.LEAC(JEL,IDIMP1,2)=LNVADJ.LECT(J)
  1195.             IF (LEADJ.LECT(J).NE.0) THEN
  1196.                ILEA1.LEAC(LEADJ.LECT(J),LNVADJ.LECT(J),1)=JEL
  1197.                ILEA1.LEAC(LEADJ.LECT(J),LNVADJ.LECT(J),2)=IDIMP1
  1198.             ENDIF
  1199. C           et je rajoute au segment ITRC1 ce qu'il faut :
  1200. C
  1201.             CALL BCIRCO(IRN1,IPT2,JEL,XCBJ,XRBJ)
  1202. C
  1203.             ITRC1.TRC(JEL,1)=XRBJ
  1204.             DO K=1,IDIM
  1205.                ITRC1.TRC(JEL,K+1)=XCBJ(K)
  1206.             ENDDO
  1207.          ENDDO
  1208. C
  1209. C        2.3.6 Mise a jour de la table de connectivite (partie 2 sur 2)
  1210. C        On parcoure les elements de la boule et on determine toutes les
  1211. C        elements adjacents internes a la boule
  1212.          JG=IDIM-1
  1213.          SEGINI,LREF1
  1214.          DO K=1,NBOUL
  1215.             NELK=LBOUL.LECT(K)
  1216. C           On recherche l'adjacent de NELK vu par le noeud M
  1217.             DO M=1,IDIM
  1218. C              Numeros des noeuds de l'element NELK a observer
  1219.                IF (IDIM.EQ.2) THEN
  1220.                   IF (M.EQ.1) THEN
  1221.                      LREF1.LECT(1)=IPT2.NUM(2,NELK)
  1222.                   ELSEIF (M.EQ.2) THEN
  1223.                      LREF1.LECT(1)=IPT2.NUM(1,NELK)
  1224.                   ENDIF
  1225.                ENDIF
  1226.                IF (IDIM.EQ.3) THEN
  1227.                   IF (M.EQ.1) THEN
  1228.                      LREF1.LECT(1)=IPT2.NUM(2,NELK)
  1229.                      LREF1.LECT(2)=IPT2.NUM(3,NELK)
  1230.                   ELSEIF (M.EQ.2) THEN
  1231.                      LREF1.LECT(1)=IPT2.NUM(1,NELK)
  1232.                      LREF1.LECT(2)=IPT2.NUM(3,NELK)
  1233.                   ELSEIF (M.EQ.3) THEN
  1234.                      LREF1.LECT(1)=IPT2.NUM(1,NELK)
  1235.                      LREF1.LECT(2)=IPT2.NUM(2,NELK)
  1236.                   ENDIF
  1237.                ENDIF
  1238. C              On parcoure les autres elements de la boule
  1239.                NNOKK=0
  1240.                DO KK=1,NBOUL
  1241.                   NELKK=LBOUL.LECT(KK)
  1242.                   IF (NELKK.EQ.NELK) GOTO 2361
  1243.                   NCOM=1
  1244.                   NSOMKK=IDIMP1
  1245. C                 On parcoure les 2/3 premiers noeuds LL de NELKK
  1246.                   DO LL=1,IDIM
  1247.                      NREFX=IPT2.NUM(LL,NELKK)
  1248. C                    Test si le noeud LL de NELKK est commun a NELK
  1249.                      DO L=1,IDIM-1
  1250.                         IF (NREFX.EQ.LREF1.LECT(L)) THEN
  1251.                            NCOM=NCOM+1
  1252.                            NSOMKK=NSOMKK+LL
  1253.                         ENDIF
  1254.                      ENDDO
  1255. C                    Test si NELKK est bien un element adjacent a NELK
  1256.                      IF (NCOM.EQ.IDIM) GOTO 2362
  1257.                   ENDDO
  1258. 2361              CONTINUE
  1259.                ENDDO
  1260. 2362           CONTINUE
  1261. C              Determination du numero du noeud de voyeur de K dans IPT2
  1262.                IF (IDIM.EQ.2) NNOKK=6-NSOMKK
  1263.                IF (IDIM.EQ.3) NNOKK=10-NSOMKK
  1264.                ILEA1.LEAC(NELK,M,1)=NELKK
  1265.                ILEA1.LEAC(NELK,M,2)=NNOKK
  1266.             ENDDO
  1267.          ENDDO
  1268. C
  1269. C        2.3.7 Renumerotation et nettoyage des elements de IPT2 dans le
  1270. C        cas particulier ou le nombre d'elements a diminue
  1271.          IF (N43.LT.0) THEN
  1272.             DO J=1,IPT2.NUM(/2)
  1273. C              Determination du nombre d'elements a decaler
  1274.                NDEC=0
  1275.                DO K=1,NSUP
  1276. C                 Si l'element examine n'est plus dans la triangulation
  1277. C                 alors on ne fait rien
  1278.                   IF (J.EQ.LSUP.LECT(K)) GOTO 2371
  1279.                   IF (J.GE.LSUP.LECT(K)) NDEC=NDEC+1
  1280.                ENDDO
  1281. C              Nouveau numero que portera l'element J
  1282.                JJ=J-NDEC
  1283. C              Renumerotation dans IPT2
  1284.                DO K=1,IPT2.NUM(/1)
  1285.                   IPT2.NUM(K,JJ)=IPT2.NUM(K,J)
  1286.                ENDDO
  1287. C              Renumerotation dans ILEA1
  1288.                DO K=1,IDIMP1
  1289.                   IOLDVJK=ILEA1.LEAC(J,K,1)
  1290. C                 Determination du nombre d'elements a decaler
  1291.                   NDEC2=0
  1292.                   DO KK=1,NSUP
  1293.                      IF (IOLDVJK.EQ.LSUP.LECT(KK)) THEN
  1294.                         NBNN=IPT2.NUM(/1)
  1295.                         NBELEM=NBELE2
  1296.                         NBSOUS=0
  1297.                         NBREF=0
  1298.                         SEGADJ,IPT2
  1299. C                       Erreur de triangulation : vérifiez qu'il n'y a
  1300. C                       pas de points confondus
  1301.                         CALL ERREUR(846)
  1302.                         GOTO 2501
  1303.                      ENDIF
  1304.                      IF (IOLDVJK.GE.LSUP.LECT(KK)) NDEC2=NDEC2+1
  1305.                   ENDDO
  1306.                   NEWVJK=IOLDVJK-NDEC2
  1307.                   ILEA1.LEAC(JJ,K,1)=NEWVJK
  1308.                   ILEA1.LEAC(JJ,K,2)=ILEA1.LEAC(J,K,2)
  1309.                ENDDO
  1310. C              Renumerotation dans ITCR1
  1311.                ITRC1.TRC(JJ,1)=ITRC1.TRC(J,1)
  1312.                DO K=1,IDIM
  1313.                   ITRC1.TRC(JJ,K+1)=ITRC1.TRC(J,K+1)
  1314.                ENDDO
  1315. 2371           CONTINUE
  1316.             ENDDO
  1317.             NBNN=IPT2.NUM(/1)
  1318.             NBELEM=NBELE2
  1319.             SEGADJ,IPT2
  1320.             NBL1=NBELE2
  1321.             SEGADJ,ILEA1
  1322.             NBE1=NBELE2
  1323.             SEGADJ,ITRC1
  1324.          ENDIF
  1325.          SEGSUP,IPT3,IPT4
  1326. C
  1327. C        2.3.8 Cas particulier de la dimension 1, on ajoute le point
  1328. C        NPI1 au maillage
  1329. 2381     IF (IDIM.EQ.1) THEN
  1330. C           Nouveau nombre d'element dans IPT2
  1331.             NEW2=NBELE2+1
  1332. C           J'ajuste IPT2, ITRC1 et ILEA1 si besoin
  1333.             IF (NEW2.GT.IPT2.NUM(/2)) THEN
  1334.                NBNN=IPT2.NUM(/1)
  1335.                NBELEM=NBELE2+ICMIN
  1336.                SEGADJ,IPT2
  1337.                NBE1=NBELEM
  1338.                SEGADJ,ITRC1
  1339.                NBL1=NBELEM
  1340.                SEGADJ,ILEA1
  1341.             ENDIF
  1342.             NBELE2=NEW2
  1343. C           Modification de IPT2 pour ajouter un element
  1344.             IPT2.NUM(1,NBASE1)=NREFA
  1345.             IPT2.NUM(2,NBASE1)=NPI1
  1346.             IPT2.NUM(1,NBELE2)=NPI1
  1347.             IPT2.NUM(2,NBELE2)=NREFB
  1348. C           Mise a jour de la table de connectivite ILEA1
  1349.             NED=ILEA1.LEAC(NBASE1,1,1)
  1350.             ILEA1.LEAC(NBASE1,1,1)=NBELE2
  1351.             ILEA1.LEAC(NBASE1,1,2)=2
  1352.             ILEA1.LEAC(NBELE2,2,1)=NBASE1
  1353.             ILEA1.LEAC(NBELE2,2,2)=1
  1354.             IF (NED.EQ.0) THEN
  1355.                ILEA1.LEAC(NBELE2,1,1)=0
  1356.                ILEA1.LEAC(NBELE2,1,2)=0
  1357.             ELSE
  1358.                ILEA1.LEAC(NBELE2,1,1)=NED
  1359.                ILEA1.LEAC(NBELE2,1,2)=2
  1360.                ILEA1.LEAC(NED,2,1)=NBELE2
  1361.                ILEA1.LEAC(NED,2,2)=1
  1362.             ENDIF
  1363.          ENDIF
  1364.       ENDDO
  1365. C     Ajustement final de IPT2, ITRC1 et ILEA1
  1366.       IF (IDIM.NE.1) THEN
  1367.          NBNN=IPT2.NUM(/1)
  1368.          NBELEM=NBELE2
  1369.          SEGADJ,IPT2
  1370.          NBE1=NBELE2
  1371.          SEGADJ,ITRC1
  1372.          NBL1=NBELE2
  1373.          SEGADJ,ILEA1
  1374.       ENDIF
  1375. C
  1376. C
  1377. C------- 2.4 Retour du maillage dans le cas ou l'on garde la boite
  1378. C     Si IBOI egal a 1, on renvoie IPT2 (maillage avec boite)
  1379. 2501  CONTINUE
  1380.       IF (IERR.NE.0) GOTO 999
  1381.       IF (IBOI.EQ.1) THEN
  1382.          IF (IDIM.NE.1) THEN
  1383.             SEGSUP,LBASE,LCAVI,LBOUL,LREF1,LEADJ,LNVADJ
  1384.             SEGSUP,LSUP,LSCP
  1385.          ENDIF
  1386.          GOTO 999
  1387.       ENDIF
  1388. C
  1389. C
  1390. C
  1391. C----------------------------------------------------------------------C
  1392. C                               PARTIE 3                               C
  1393. C      Elimination des elements appuyes sur les coins de la boite      C
  1394. C----------------------------------------------------------------------C
  1395. C
  1396. C------- 3.1 Suppression des elements touchant les noeuds de la boite
  1397. C     Je parcours les noeuds de la boite en cherchant les elements du
  1398. C     maillage IPT2 qui en contiennent. Ceux qui n'en contiennent pas
  1399. C     sont recopies dans IPT3, qui sera le maillage de sortie :
  1400.       NBSOUS=0
  1401.       NBREF=0
  1402.       NBNN=IPT2.NUM(/1)
  1403.       NBELEM=NBELE2
  1404.       SEGINI,IPT3
  1405.       IPT3.ITYPEL=IPT2.ITYPEL
  1406.       NBELE3=0
  1407.       DO I=1,NBELE2
  1408.          DO J=1,IPT2.NUM(/1)
  1409.             NPI2=IPT2.NUM(J,I)
  1410.             DO K=1,NNBOIT
  1411.                IF (NPI2.EQ.LNBOIT(K)) GOTO 3101
  1412.             ENDDO
  1413.          ENDDO
  1414.          NBELE3=NBELE3+1
  1415.          DO J=1,IPT2.NUM(/1)
  1416.             IPT3.NUM(J,NBELE3)=IPT2.NUM(J,I)
  1417.          ENDDO
  1418. 3101     CONTINUE
  1419.       ENDDO
  1420.       NBNN=IPT3.NUM(/1)
  1421.       NBELEM=NBELE3
  1422.       SEGADJ,IPT3
  1423. C
  1424. C------- 3.2 Cas ou IPT3 ne contient pas d'elements.
  1425.       IF (IPT3.NUM(/2).EQ.0) THEN
  1426. C        En dimension 2, cela veut dire que les points de IPT1 sont
  1427. C        alignes. On va alors rendre le maillage des lignes ne touchant
  1428. C        pas les noeuds de la boite
  1429.          IF (IDIM.EQ.2) GOTO 3203
  1430. C        En dimension 3, cela veut dire que les points sont coplanaires.
  1431. C        On refait la meme operation pour rendre le maillage faces ne
  1432. C        touchant pas les noeuds de la boite
  1433.          JG=IDIM
  1434.          SEGINI,LREF1
  1435.          JG=IPT2.NUM(/2)
  1436.          SEGINI,LSUP
  1437.          NBSOUS=0
  1438.          NBREF=0
  1439.          NBNN=IDIM
  1440.          NBELEM=NBELE2
  1441.          SEGINI,IPT3
  1442. C        On va rendre un maillage de TRI3
  1443.          IPT3.ITYPEL=4
  1444.          NBELE3=0
  1445. C        Boucle sur les elements de la triangulation avec la boite
  1446.          DO I=1,NBELE2
  1447.             NNS=0
  1448.             NSOM=0
  1449.             NEL=0
  1450. C           On regarde si le noeud J de l'element I fait partie de ceux
  1451. C           de la boite
  1452.             DO J=1,IPT2.NUM(/1)
  1453.                NPI2=IPT2.NUM(J,I)
  1454.                DO K=1,NNBOIT
  1455.                   IF (NPI2.EQ.LNBOIT(K)) THEN
  1456.                      NNS=NNS+1
  1457.                      NSOM=J
  1458.                      NEL=ILEA1.LEAC(I,J,1)
  1459.                   ENDIF
  1460.                ENDDO
  1461.             ENDDO
  1462. C           Si l'element I n'a qu'un noeud faisant partie de ceux de la
  1463. C           boite, alors on cree l'element appuyes sur ses autres noeuds
  1464.             IF (NNS.EQ.1) THEN
  1465.                DO K=1,LSUP.LECT(/1)
  1466.                   IF (LSUP.LECT(K).EQ.0) GOTO 3201
  1467.                   IF (I.EQ.LSUP.LECT(K)) THEN
  1468.                      GOTO 3202
  1469.                   ENDIF
  1470.                ENDDO
  1471. 3201           CONTINUE
  1472.                NBELE3=NBELE3+1
  1473.                LSUP.LECT(NBELE3)=NEL
  1474.                IF (NSOM.EQ.1) THEN
  1475.                   LREF1.LECT(1)=2
  1476.                   LREF1.LECT(2)=3
  1477.                   LREF1.LECT(3)=4
  1478.                ELSEIF (NSOM.EQ.2) THEN
  1479.                   LREF1.LECT(1)=1
  1480.                   LREF1.LECT(2)=3
  1481.                   LREF1.LECT(3)=4
  1482.                ELSEIF (NSOM.EQ.3) THEN
  1483.                   LREF1.LECT(1)=1
  1484.                   LREF1.LECT(2)=2
  1485.                   LREF1.LECT(3)=4
  1486.                ELSEIF (NSOM.EQ.4) THEN
  1487.                   LREF1.LECT(1)=1
  1488.                   LREF1.LECT(2)=2
  1489.                   LREF1.LECT(3)=3
  1490.                ENDIF
  1491.                DO J=1,IPT3.NUM(/1)
  1492.                   IPT3.NUM(J,NBELE3)=IPT2.NUM(LREF1.LECT(J),I)
  1493.                ENDDO
  1494.             ENDIF
  1495. 3202        CONTINUE
  1496.          ENDDO
  1497.          NBNN=IPT3.NUM(/1)
  1498.          NBELEM=NBELE3
  1499.          SEGADJ,IPT3
  1500.       ENDIF
  1501. C
  1502. C     Cas ou les points de IPT1 sont alignes. On change le maillage en
  1503. C     lignes et ne garde que les lignes ne touchant pas les noeuds de la
  1504. C     boite
  1505. 3203  CONTINUE
  1506.       IF (IPT3.NUM(/2).EQ.0) THEN
  1507.          CALL ECROBJ('MAILLAGE',IPT2)
  1508.          CALL CHANLG
  1509.          CALL LIROBJ('MAILLAGE',IPT4,1,IRETOU)
  1510.          IF (IERR.NE.0) RETURN
  1511.          SEGACT,IPT4
  1512.          NBELE4=IPT4.NUM(/2)
  1513.          NBSOUS=0
  1514.          NBREF=0
  1515.          NBNN=IPT4.NUM(/1)
  1516.          NBELEM=NBELE4
  1517.          SEGINI,IPT3
  1518. C        On va rendre un maillage de SEG2
  1519.          IPT3.ITYPEL=IPT4.ITYPEL
  1520.          NBELE3=0
  1521. C        Boucle sur les lignes de la triangulation avec la boite
  1522.          DO I=1,NBELE4
  1523.             NNS=0
  1524. C           On regarde si les noeuds l'element I font partie de ceux de
  1525. C           la boite
  1526.             DO J=1,IPT4.NUM(/1)
  1527.                NPI4=IPT4.NUM(J,I)
  1528.                DO K=1,NNBOIT
  1529.                   IF (NPI4.EQ.LNBOIT(K)) THEN
  1530.                      NNS=NNS+1
  1531.                   ENDIF
  1532.                ENDDO
  1533.             ENDDO
  1534. C           Si l'element I n'a aucun noeud faisant partie de ceux de la
  1535. C           boite, alors on cree l'element appuyes sur ses noeuds
  1536.             IF (NNS.EQ.0) THEN
  1537.                NBELE3=NBELE3+1
  1538.                IPT3.NUM(1,NBELE3)=IPT4.NUM(1,I)
  1539.                IPT3.NUM(2,NBELE3)=IPT4.NUM(2,I)
  1540.             ENDIF
  1541.          ENDDO
  1542.          NBNN=IPT3.NUM(/1)
  1543.          NBELEM=NBELE3
  1544.          SEGADJ,IPT3
  1545.       ENDIF
  1546. C
  1547. C     Un peu de menage avant de quitter
  1548.       IF (IDIM.EQ.1) THEN
  1549.          SEGSUP,LSUP
  1550.       ELSE
  1551.          SEGSUP,IPT2,LBASE,LCAVI,LBOUL,LREF1,LEADJ,LNVADJ
  1552.          SEGSUP,LSUP,LSCP
  1553.       ENDIF
  1554.       IPT2=IPT3
  1555.       NBPTS=NBPTS0
  1556.       SEGADJ,MCOORD
  1557. C
  1558. 999   RETURN
  1559.       END
  1560.  
  1561.  
  1562.  

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