Esope source files

   1 : C TRIA1     SOURCE    SH236661  13/11/25    21:15:10     7869           
   2 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
   3 :       SUBROUTINE TRIA1(MPOVAL,IPT1,IVC,XDEN,IPT2)
   4 : C
   5 : C Triangulation de Delaunay d un ensemble de points
   6 : C
   7 : C IPT1 : Maillage de POI1 a trianguler
   8 : C IVC  : Indicateur pour verifier (IVC=1) ou non (IVC=0) la convexite
   9 : C        de la triangulation
  10 : C XDEN : Critere sur la taille de maille maximale cible (le maillage est
  11 : C        raffine en ajoutant des points au milieu des aretes)
  12 : C IPT2 : Maillage de la triangulation
  13 : C
  14 :       IMPLICIT INTEGER(I-N)
  15 :       IMPLICIT REAL*8 (A-H,O-Z)
  16 : C
  17 : -INC SMCHPOI
  18 : -INC CCOPTIO
  19 : -INC SMELEME
  20 : -INC SMCOORD
  21 : -INC CCGEOME
  22 :       PARAMETER (ICMIN=1000)
  23 :       DIMENSION XA(3),XB(3),XPMIL(3)
  24 :       DIMENSION LNBOIT(8)
  25 : C
  26 :       SEGMENT,MCIRCONS
  27 :          REAL*8 TRC(NBE1,4)
  28 :       ENDSEGMENT
  29 :       POINTEUR ITRC1.MCIRCONS
  30 : C
  31 :       SEGMENT,MADJACEN
  32 :          INTEGER LEAC(NBL1,IDIM+1,2)
  33 :       ENDSEGMENT
  34 :       POINTEUR ILEA1.MADJACEN
  35 : C
  36 : C
  37 :       SEGINI,IPT4=IPT1
  38 : C     Taille de la boite de triangulation :
  39 :       IF (IDIM.EQ.1) XBOI=2.
  40 :       IF (IDIM.EQ.2) XBOI=200.
  41 :       IF (IDIM.EQ.3) XBOI=500.
  42 : C     IBOI=1  -->  on garde la boite de triangulation
  43 : C     IBOI=0  -->  on ne la garde pas (valeur par defaut)
  44 :       IBOI=0
  45 : C     Triangulation de DELAUNAY
  46 :       NBVC=0
  47 : 100   CONTINUE
  48 : c
  48 :       CALL DELAUN(MPOVAL,IPT4,XBOI,IBOI,IPT2,LNBOIT,ITRC1,ILEA1)
  49 : c
  50 :       IF (IDIM.NE.1) SEGSUP,ITRC1,ILEA1
  51 : C     En cas d'erreur dans DELAUN
  52 :       IF (IERR.EQ.2) THEN
  53 :          SEGSUP,IPT4
  54 :          IERR=0
  55 :          GOTO 999
  56 :       ENDIF
  57 :       SEGACT,IPT2
  58 : C     Verification de la convexite de la premiere triangulation
  59 :       IF (IDIM.EQ.1) GOTO 102
  60 :       IF ((IDIM.EQ.2).AND.(IPT2.ITYPEL.EQ.2)) GOTO 102
  61 :       IF ((IDIM.EQ.3).AND.((IPT2.ITYPEL.EQ.2).OR.(IPT2.ITYPEL.EQ.4)))
  62 :      &   GOTO 102
  63 :       IF ((IVC.EQ.1).AND.(NBVC.EQ.0)) THEN
  64 :          ICONV=0
  65 :          NBVC=NBVC+1
  65 :          CALL VERCON(IPT2,ICONV)
  66 :          IF (ICONV.EQ.0) THEN
  67 :             DO I=1,10
  68 :                XBOI=5.*XBOI
  69 : c
  69 :              CALL DELAUN(MPOVAL,IPT4,XBOI,IBOI,IPT2,LNBOIT,ITRC1,ILEA1)
  70 : c
  71 :                IF (IDIM.NE.1) SEGSUP,ITRC1,ILEA1
  72 :                NBVC=NBVC+1
  72 :                CALL VERCON(IPT2,ICONV)
  73 :                IF (ICONV.EQ.1) GOTO 101
  74 :             ENDDO
  75 : 101         CONTINUE
  76 :          ENDIF
  77 :          NBVC=1
  78 :       ENDIF
  79 : 102   CONTINUE
  80 : C     Re-triangulation en ajoutant des noeuds sur les aretes
  81 :       SEGACT,IPT2
  82 :       IF (XDEN.NE.0.) THEN
  83 :          NBELE0=XCOOR(/1)/(IDIM+1)
  84 :          NBELE00=NBELE0
  85 :          NBELE4=IPT4.NUM(/2)
  86 : C        Maillage des lignes de la triangulation initiale
  87 :          IF (IPT2.ITYPEL.EQ.2) THEN
  88 :             SEGINI,IPT3=IPT2
  89 :          ELSE
  89 :             CALL ECROBJ('MAILLAGE',IPT2)
  89 :             CALL CHANLG
  89 :             CALL LIROBJ('MAILLAGE',IPT3,1,IRETOU)
  90 :             IF (IERR.NE.0) RETURN
  91 :             SEGACT,IPT3
  92 :          ENDIF
  93 : C        Boucle sur les lignes pour l'ajout de noeuds milieux
  94 :          DO I=1,IPT3.NUM(/2)
  95 : C           Calcul de la distance de la ligne AB
  96 :             NA=IPT3.NUM(1,I)
  97 :             NB=IPT3.NUM(2,I)
  98 :             DAB=0.
  99 :             DO J=1,IDIM
 100 :                XA(J)=XCOOR((NA-1)*(IDIM+1)+J)
 101 :                XB(J)=XCOOR((NB-1)*(IDIM+1)+J)
 102 :                DAB=DAB+((XA(J)-XB(J))**2)
 103 :             ENDDO
 104 :             DAB=DAB**0.5
 105 :             IF (DAB.GT.XDEN) THEN
 106 : C              Creation d'un nouveau noeud au milieu de AB
 107 :                NBELE4=NBELE4+1
 108 :                NBELE0=NBELE0+1
 109 :                NBPTS0=XCOOR(/1)/(IDIM+1)
 110 : C              Ajustement du segment MCOORD si besoin
 111 :                IF (NBPTS0.LT.NBELE0) THEN
 112 :                   NBPTS=NBPTS0+ICMIN
 113 :                   SEGADJ,MCOORD
 114 :                ENDIF
 115 : C              Ecriture des coordonnees du nouveau noeud dans XCOOR
 116 :                DO J=1,IDIM
 117 :                   XPMIL(J)=0.5*(XA(J)+XB(J))
 118 :                   XCOOR(((NBELE0-1)*(IDIM+1))+J)=XPMIL(J)
 119 :                ENDDO
 120 : C              et sa densite
 121 :                XCOOR(((NBELE0-1)*(IDIM+1))+(IDIM+1))=0.
 122 : C              Ajustement du segment IPT4 si besoin
 123 :                NBE4=IPT4.NUM(/2)
 124 :                IF (NBE4.LT.NBELE4) THEN
 125 :                   NBNN=1
 126 :                   NBELEM=NBE4+ICMIN
 127 :                   NBSOUS=0
 128 :                   NBREF=0
 129 :                   SEGADJ,IPT4
 130 :                ENDIF
 131 : C              Ajout de ce point dans le maillage de points IPT4
 132 :                IPT4.NUM(1,NBELE4)=NBELE0
 133 :             ENDIF
 134 :          ENDDO
 135 : C        Ajustement final de MCOORD et IPT4
 136 :          IF (NBELE0.NE.NBELE00) THEN
 137 :             NBPTS=NBELE0
 138 :             SEGADJ,MCOORD
 139 :             NBNN=1
 140 :             NBELEM=NBELE4
 141 :             NBSOUS=0
 142 :             NBREF=0
 143 :             SEGADJ,IPT4
 144 : C           On remonte faire la triangulation de IPT4
 145 :             GOTO 100
 146 :          ENDIF
 147 :          SEGSUP,IPT3
 148 :       ENDIF
 149 :       SEGSUP,IPT4
 150 : C     Sortie/fin
 151 : 999   RETURN
 152 :       END
 153 : 
 154 :  
 155 :  
 156 :