Télécharger tria2.eso

Retour à la liste

tria2
C TRIA2     SOURCE    PV        20/03/24    21:22:40     10554          
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE TRIA2(IPT1,XDEN,IPT5)
C
C Maillage de l'interieur d'un contour/enveloppe par triangulation de
C Delaunay puis ajout de points
C
C IPT1 : Maillage initial constitue :
C        - de SEG2 (contour ferme) ou de TRI3 (surface) en dimension 2
C        - de TRI3 (enveloppe fermee) ou de TET4 (volume) en dimension 3
C XDEN : Critere sur la taille de maille maximale cible (le maillage est
C        raffine en ajoutant des points au milieu des aretes)
C IPT5 : Maillage de sortie constitue :
C        - de TRI3 en dimension 2
C        - de TET4 en dimension 3
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C
 
-INC PPARAM
-INC CCOPTIO
-INC SMELEME
-INC SMCOORD
-INC CCGEOME
-INC SMLENTI
      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
      POINTEUR LNCONT.MLENTI
C
      SEGACT,IPT1
      NBSZ=IPT1.LISOUS(/1)
      NTYP=IPT1.ITYPEL
      I1=0
C     On change le maillage d'entree en POI1 --> IPT2
      SEGINI,IPT3=IPT1
      CALL CHANGE(IPT3,1)
      SEGINI,IPT2=IPT3
101   CONTINUE
C     Maillage englobant a utiliser dans INCL --> IPT3
      IF (I1.EQ.0) THEN
         IF ((IDIM.EQ.3).AND.(NTYP.EQ.4)) THEN
            CALL ECROBJ('MAILLAGE',IPT1)
            CALL VOLUME
            CALL LIROBJ('MAILLAGE',IPT3,1,IRETOU)
            IF (IERR.NE.0) RETURN
            SEGACT,IPT3
         ELSE
            SEGINI,IPT3=IPT1
         ENDIF
      ELSE
         SEGINI,IPT3=IPT5
      ENDIF
C     Calcul du contour/enveloppe de reference --> IPT4
      IF (IDIM.EQ.2) THEN
         IF (IPT3.ITYPEL.EQ.2) THEN
            SEGINI,IPT4=IPT1
         ELSE
            CALL ECROBJ('MAILLAGE',IPT3)
            CALL PRCONT
            CALL LIROBJ('MAILLAGE',IPT4,1,IRETOU)
            IF (IERR.NE.0) RETURN
            SEGACT,IPT4
         ENDIF
      ELSEIF (IDIM.EQ.3) THEN
         CALL ECROBJ('MAILLAGE',IPT3)
         CALL ENVVOL
         CALL LIROBJ('MAILLAGE',IPT4,1,IRETOU)
         IF (IERR.NE.0) RETURN
         SEGACT,IPT4
      ENDIF
C     Creation d'une liste des numeros de noeuds
      IF (I1.EQ.0) THEN
         segact mcoord *mod
         JG=nbpts
         SEGINI,LNCONT
C        et marquage des noeuds situes sur le contour/enveloppe IPT4
         DO I=1,IPT4.NUM(/2)
            DO J=1,IPT4.NUM(/1)
               LNCONT.LECT(IPT4.NUM(J,I))=1
            ENDDO
         ENDDO
      ENDIF
C     Triangulation de DELAUNAY de IPT2 --> IPT5
C     (sans ajout de points)
      IF (IDIM.EQ.2) XBOI=200.
      IF (IDIM.EQ.3) XBOI=500.
c
      MPOVAL=0
      CALL DELAUN(MPOVAL,IPT2,XBOI,0,IPT5,LNBOIT,ITRC1,ILEA1)
c
      IF (IDIM.NE.1) SEGSUP,ITRC1,ILEA1
      SEGACT,IPT2,IPT5
C     En cas d'erreur dans DELAUN
      IF (IERR.EQ.2) THEN
         IERR=0
C        On rend le maillage pour controle
         SEGSUP,IPT2,IPT3,LNCONT
         SEGDES,IPT1
         GOTO 999
      ENDIF
C     On isole les elements de IPT5 inclus dans IPT3 par appel a
C     INCLUS option 'BARY'
      CALL ECRCHA('BARY')
      IF (IDIM.EQ.3) CALL ECRCHA('VOLU')
      CALL ECROBJ('MAILLAGE',IPT3)
      CALL ECROBJ('MAILLAGE',IPT5)
      CALL INCLUS
C     Les elements inclus sont ranges dans IPT5 (on ecrase le precedent)
      SEGSUP,IPT5
      CALL LIROBJ('MAILLAGE',IPT5,1,IRETOU)
      IF (IERR.NE.0) RETURN
      SEGACT,IPT5
C     Creation du contour/enveloppe de IPT5 --> IPT6
      CALL ECROBJ('MAILLAGE',IPT5)
      IF (IDIM.EQ.2) THEN
         CALL PRCONT
      ELSEIF (IDIM.EQ.3) THEN
         CALL ENVVOL
      ENDIF
      CALL LIROBJ('MAILLAGE',IPT6,1,IRETOU)
      IF (IERR.NE.0) RETURN
      SEGACT,IPT6
C     Test si les noeuds du contour IPT6 sont marques
      DO I=1,IPT6.NUM(/2)
         DO J=1,IPT6.NUM(/1)
            IF (LNCONT.LECT(IPT6.NUM(J,I)).NE.1) THEN
C              Maillage incorrect
               CALL ERREUR(426)
               SEGSUP,IPT2,IPT3,LNCONT
               SEGDES,IPT1
               GOTO 999
            ENDIF
         ENDDO
      ENDDO
      IF (XDEN.NE.0.) THEN
C        Maillage des lignes de la triangulation --> IPT7
         NBELE0=nbpts
         NBELE00=NBELE0
         NBELE2=IPT2.NUM(/2)
         CALL ECROBJ('MAILLAGE',IPT5)
         CALL CHANLG
         CALL LIROBJ('MAILLAGE',IPT7,1,IRETOU)
         IF (IERR.NE.0) RETURN
         SEGACT,IPT7
C        On anticipe l'ajout de noeuds sur le contour/enveloppe de
C        reference, il faudra les marquer dans LNCONT
         NLI=IPT7.NUM(/2)
         LNC0=LNCONT.LECT(/1)
         JG=LNC0+NLI
         SEGADJ,LNCONT
C        Boucle sur les lignes pour l'ajout de noeuds milieux
         DO I=1,IPT7.NUM(/2)
C           Calcul de la distance de la ligne AB
            NA=IPT7.NUM(1,I)
            NB=IPT7.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
               NBELE2=NBELE2+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.5*(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.
C              Ajustement du segment IPT2 si besoin
               NBE2=IPT2.NUM(/2)
               IF (NBE2.LT.NBELE2) THEN
                  NBNN=1
                  NBELEM=NBE2+ICMIN
                  NBSOUS=0
                  NBREF=0
                  SEGADJ,IPT2
               ENDIF
C              Ajout de ce point dans le maillage de points IPT2
               IPT2.NUM(1,NBELE2)=NBELE0
C              Si le point est sur le bord, on le marque dans LNCONT
               IF ((LNCONT.LECT(NA).EQ.1).AND.(LNCONT.LECT(NB).EQ.1))
     &         THEN
                  LNCONT.LECT(NBELE0)=1
               ENDIF
            ENDIF
         ENDDO
C        Ajustement de LNCONT
         JG=LNC0+(NBELE0-NBELE00)
         SEGADJ,LNCONT
C        Ajustement final de MCOORD et IPT2
         IF (NBELE0.NE.NBELE00) THEN
            NBPTS=NBELE0
            SEGADJ,MCOORD
            NBNN=1
            NBELEM=NBELE2
            NBSOUS=0
            NBREF=0
            SEGADJ,IPT2
C           Suppression des maillages temporaires utilises
            SEGSUP,IPT3,IPT7
C           On remonte pour refaire la triangulation
            I1=1
            GOTO 101
         ENDIF
         SEGSUP,IPT7
      ENDIF
C     Ecriture sortie/fin
      SEGSUP,IPT2,IPT3,LNCONT
      SEGDES,IPT1
C     Sortie/fin
999   RETURN
      END
 
 
 
 
 
 
 

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