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

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