Sources Esope | Cast3M

voro
C VORO      SOURCE    CB215821  23/01/25    21:15:39     11573          
      SUBROUTINE VORO
C     ******************************************************************
C     INTERFACE CASTEM 2000
C
C
C     TAB1 = VORO MAIL1 MAIL2 ;
C
C
C     OBJET :
C     _______
C
C     L'OPERATEUR VORONOI CONSTRUIT LA PARTITION DE VORONOI D'UN
C     ENSEMBLE DE POINTS MAIL1 RESTREINTE PAR UN CONTOUR/ENVELOPPE MAIL2
C
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
-INC SMLENTI
-INC SMTABLE
C
      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        =======================
C     --- 1.LECTURE DES DONNEES ---
C        =======================
C
C     LECTURE DES OBJETS COURANTS (ENTREES)
C     =====================================
      CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU)
      IF (IERR.NE.0) THEN
C        ON A PAS TROUVE LE MAILLAGE
         CALL ERREUR(503)
         GOTO 999
      ENDIF
      SEGACT,IPT1
      NBSZ1=IPT1.LISOUS(/1)
      NTYP1=IPT1.ITYPEL
C
      CALL LIROBJ('MAILLAGE',IPT2,1,IRETOU)
      IF (IERR.NE.0) THEN
C        ON A PAS TROUVE LE MAILLAGE
         SEGDES,IPT1
         CALL ERREUR(503)
         GOTO 999
      ENDIF
      SEGACT,IPT2
      NBSZ2=IPT2.LISOUS(/1)
      NTYP2=IPT2.ITYPEL
C
C     Parametre optionnel : critere de distance pour considerer deux
C     centres de cercles circonscrits confondus
C     Par defaut, on le calcule selon la taille du nuage de points
      IDIMP1=IDIM+1
      DDMAX=0.
      IP1=IPT1.NUM(1,1)
 
      SEGACT,MCOORD
      XP1=XCOOR((IP1-1)*IDIMP1+1)
      YP1=XCOOR((IP1-1)*IDIMP1+2)
      ZP1=XCOOR((IP1-1)*IDIMP1+3)
      XP2=XP1
      YP2=YP1
      ZP2=ZP1
      DO I=2,IPT1.NUM(/2)
         IPI=IPT1.NUM(1,I)
         XPI=XCOOR((IPI-1)*IDIMP1+1)
         YPI=XCOOR((IPI-1)*IDIMP1+2)
         ZPI=XCOOR((IPI-1)*IDIMP1+3)
         IF(XPI.LT.XP1) XP1=XPI
         IF(YPI.LT.YP1) YP1=YPI
         IF(ZPI.LT.ZP1) ZP1=ZPI
         IF(XPI.GT.XP2) XP2=XPI
         IF(YPI.GT.YP2) YP2=YPI
         IF(ZPI.GT.ZP2) ZP2=ZPI
      ENDDO
      SEGDES,MCOORD
 
      DDMAX=(XP2-XP1)+(YP2-YP1)
      IF (IDIM.EQ.3) DDMAX=DDMAX+(ZP2-ZP1)
      DDMAX=DDMAX*0.333333
C     WRITE(6,*) 'DDMAX=', DDMAX
C
C---- SI MAILLAGE DE POI1 ET CONTOUR/ENVELOPPE : PARTITION DE VORONOI --
C
C     L'operateur n'est disponible qu'en dimension 2 ou 3
      IF ((IDIM.NE.2).AND.(IDIM.NE.3)) THEN
         INTERR(1)=IDIM
C        FONCTION INDISPONIBLE EN DIMENSION %I1
         SEGDES,IPT1,IPT2
         CALL ERREUR(709)
         GOTO 999
      ENDIF
C     Test sur les maillages d'entree
C     --> Le nombre de sous zones doit etre nul
      IF ((NBSZ1.NE.0).OR.(NBSZ2.NE.0)) THEN
C        MAILLAGE INCORRECT
         SEGDES,IPT1,IPT2
         CALL ERREUR(426)
         GOTO 999
      ENDIF
C     --> Le maillage MAIL1 doit etre constitue d'elements POI1
      IF (NTYP1.NE.1) THEN
C        TYPE D'ELEMENTS INCORRECT
         SEGDES,IPT1,IPT2
         CALL ERREUR(16)
         GOTO 999
      ENDIF
C     --> Le maillage MAIL2 doit etre constitue d'elements SEG2 ou TRI3
      IF (((IDIM.EQ.2).AND.(NTYP2.NE.2)).OR.
     &    ((IDIM.EQ.3).AND.(NTYP2.NE.4))) THEN
C        TYPE D'ELEMENTS INCORRECT
         SEGDES,IPT1,IPT2
         CALL ERREUR(16)
         GOTO 999
      ENDIF
C     --> Les elements du maillage MAIL2 doivent etre orientes de la
C     meme maniere (appel a l'operateur VERSENS)
      CALL ECROBJ('MAILLAGE',IPT2)
      CALL VERSEN
      CALL LIROBJ('MAILLAGE',IPT2,1,IRETOU)
      IF (IERR.NE.0) THEN
         SEGDES,IPT1
         GOTO 999
      ENDIF
C     --> Le maillage contour/enveloppe MAIL2 doit etre connexe
      CALL ECROBJ('MAILLAGE',IPT2)
      CALL CCON
      CALL DIMENS
      CALL LIRENT(NCCON,1,IRETOU)
      IF (NCCON.NE.1) THEN
         INTERR(1)=NCCON
C        CET OBJET CONTIENT %1 COMPOSANTES CONNEXES
         SEGDES,IPT1,IPT2
         CALL ERREUR(754)
         GOTO 999
      ENDIF
C
      MPOVAL = 0
      MCHPOI = 0
      CALL LIROBJ('CHPOINT ',MCHPOI,0,IRETOU)
      IF(IRETOU.EQ.1) THEN
         SEGACT MCHPOI
         MSOUPO=IPCHP(1)
         SEGACT MSOUPO
         MPOVAL = IPOVAL
         SEGACT MPOVAL
      ENDIF
C
C     On recree le maillage de points avec une reference vers le
C     maillage du contour/enveloppe pour etre pris en compte dans DELAUN
C     lors du calcul de la taille de la boite de triangulation
      SEGINI,IPT3=IPT1
      NBNN=IPT3.NUM(/1)
      NBELEM=IPT3.NUM(/2)
      NBSOUS=0
      NBREF=1
      SEGADJ,IPT3
      IPT3.LISREF(1)=IPT2
C     Appel a DELAUN avec IBOI=1 pour le calcul de la triangulation
C     de Delaunay de IPT3
      XBOI=50.
      IBOI=1
      CALL DELAUN(MPOVAL,IPT3,XBOI,IBOI,IPT4,LNBOIT,ITRC1,ILEA1)
      IF (IERR.NE.0) GOTO 999
      SEGACT,IPT4
C     Appel a VORO1 pour le calcul de la partition de Voronoi non coupee
      CALL VORO1(IPT1,IPT4,LNBOIT,ILEA1,ITRC1,DDMAX,MTAB1)
      IF (IERR.NE.0) GOTO 999
      IF (MPOVAL.NE.0) SEGDES,MPOVAL,MSOUPO,MCHPOI
C     Decoupage de la partition de Voronoi selon le contour/enveloppe
C     avec VORO2
      CALL VORO2(IPT1,IPT2,DDMAX,MCHPOI,MTAB1)
      IF (IERR.NE.0) GOTO 999
C     Ecriture sortie/fin
      SEGSUP,IPT3,IPT4,ITRC1,ILEA1
      SEGDES,IPT1,IPT2
      IF (IERR.EQ.0) CALL ECROBJ('TABLE   ',MTAB1)
      GOTO 999
C
C
  999 RETURN
      END