Sources Esope | Cast3M

tria1
C TRIA1     SOURCE    CB215821  20/11/25    13:41:21     10792          
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE TRIA1(MPOVAL,IPT1,IVC,XDEN,IPT2)
C
C Triangulation de Delaunay d un ensemble de points
C
C IPT1 : Maillage de POI1 a trianguler
C IVC  : Indicateur pour verifier (IVC=1) ou non (IVC=0) la convexite
C        de la triangulation
C XDEN : Critere sur la taille de maille maximale cible (le maillage est
C        raffine en ajoutant des points au milieu des aretes)
C IPT2 : Maillage de la triangulation
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C
-INC SMCHPOI
 
-INC PPARAM
-INC CCOPTIO
-INC SMELEME
-INC SMCOORD
-INC CCGEOME
      PARAMETER (ICMIN=1000)
      DIMENSION XA(3),XB(3),XPMIL(3)
      DIMENSION LNBOIT(8)
C
      SEGMENT,MCIRCONS
         REAL*8 TRC(NBE1,4)
      ENDSEGMENT
      POINTEUR ITRC1.MCIRCONS
C
      SEGMENT,MADJACEN
         INTEGER LEAC(NBL1,IDIM+1,2)
      ENDSEGMENT
      POINTEUR ILEA1.MADJACEN
C
C
      SEGINI,IPT4=IPT1
C     Taille de la boite de triangulation :
      IF (IDIM.EQ.1) XBOI=2.
      IF (IDIM.EQ.2) XBOI=200.
      IF (IDIM.EQ.3) XBOI=500.
C     IBOI=1  -->  on garde la boite de triangulation
C     IBOI=0  -->  on ne la garde pas (valeur par defaut)
      IBOI=0
C     Triangulation de DELAUNAY
      NBVC=0
100   CONTINUE
c
      CALL DELAUN(MPOVAL,IPT4,XBOI,IBOI,IPT2,LNBOIT,ITRC1,ILEA1)
c
      IF (IDIM.NE.1) SEGSUP,ITRC1,ILEA1
C     En cas d'erreur dans DELAUN
      IF (IERR.EQ.2) THEN
         SEGSUP,IPT4
         IERR=0
         GOTO 999
      ENDIF
      SEGACT,IPT2
C     Verification de la convexite de la premiere triangulation
      IF  (IDIM.EQ.1) GOTO 102
      IF ((IDIM.EQ.2).AND. (IPT2.ITYPEL.EQ.2)) GOTO 102
      IF ((IDIM.EQ.3).AND.((IPT2.ITYPEL.EQ.2).OR.(IPT2.ITYPEL.EQ.4)))
     &   GOTO 102
      IF ((IVC.EQ.1).AND.(NBVC.EQ.0)) THEN
         ICONV=0
         NBVC=NBVC+1
         CALL VERCON(IPT2,ICONV)
         IF (ICONV.EQ.0) THEN
            DO I=1,10
               XBOI=5.*XBOI
c
             CALL DELAUN(MPOVAL,IPT4,XBOI,IBOI,IPT2,LNBOIT,ITRC1,ILEA1)
c
               IF (IDIM.NE.1) SEGSUP,ITRC1,ILEA1
               NBVC=NBVC+1
               CALL VERCON(IPT2,ICONV)
               IF (ICONV.EQ.1) GOTO 101
            ENDDO
101         CONTINUE
         ENDIF
         NBVC=1
      ENDIF
102   CONTINUE
C     Re-triangulation en ajoutant des noeuds sur les aretes
      SEGACT,IPT2
      IF (XDEN.NE.0.) THEN
         segact mcoord*mod
         NBELE0=nbpts
         NBELE00=NBELE0
         NBELE4=IPT4.NUM(/2)
C        Maillage des lignes de la triangulation initiale
         IF (IPT2.ITYPEL.EQ.2) THEN
            SEGINI,IPT3=IPT2
         ELSE
            CALL ECROBJ('MAILLAGE',IPT2)
            CALL CHANLG
            CALL LIROBJ('MAILLAGE',IPT3,1,IRETOU)
            IF (IERR.NE.0) RETURN
            CALL ACTOBJ('MAILLAGE',IPT3,1)
         ENDIF
C        Boucle sur les lignes pour l'ajout de noeuds milieux
         DO I=1,IPT3.NUM(/2)
C           Calcul de la distance de la ligne AB
            NA=IPT3.NUM(1,I)
            NB=IPT3.NUM(2,I)
            DAB=0.
            DO J=1,IDIM
               XA(J)=XCOOR((NA-1)*(IDIM+1)+J)
               XB(J)=XCOOR((NB-1)*(IDIM+1)+J)
               DAB=DAB+((XA(J)-XB(J))**2)
            ENDDO
            DAB=DAB**0.5
            IF (DAB.GT.XDEN) THEN
C              Creation d'un nouveau noeud au milieu de AB
               NBELE4=NBELE4+1
               NBELE0=NBELE0+1
               NBPTS0=nbpts
C              Ajustement du segment MCOORD si besoin
               IF (NBPTS0.LT.NBELE0) THEN
                  NBPTS=NBPTS0+ICMIN
                  SEGADJ,MCOORD
               ENDIF
C              Ecriture des coordonnees du nouveau noeud dans XCOOR
               DO J=1,IDIM
                  XPMIL(J)=0.5D0*(XA(J)+XB(J))
                  XCOOR(((NBELE0-1)*(IDIM+1))+J)=XPMIL(J)
               ENDDO
C              et sa densite
               XCOOR(((NBELE0-1)*(IDIM+1))+(IDIM+1))=0.D0
C              Ajustement du segment IPT4 si besoin
               NBE4=IPT4.NUM(/2)
               IF (NBE4.LT.NBELE4) THEN
                  NBNN=1
                  NBELEM=NBE4+ICMIN
                  NBSOUS=0
                  NBREF=0
                  SEGADJ,IPT4
               ENDIF
C              Ajout de ce point dans le maillage de points IPT4
               IPT4.NUM(1,NBELE4)=NBELE0
            ENDIF
         ENDDO
C        Ajustement final de MCOORD et IPT4
         IF (NBELE0.NE.NBELE00) THEN
            NBPTS=NBELE0
            SEGADJ,MCOORD
            NBNN=1
            NBELEM=NBELE4
            NBSOUS=0
            NBREF=0
            SEGADJ,IPT4
C           On remonte faire la triangulation de IPT4
            GOTO 100
         ENDIF
         SEGSUP,IPT3
      ENDIF
      SEGSUP,IPT4
C     Sortie/fin
999   RETURN
      END