1 : C VORO SOURCE PASCAL 14/02/14 21:15:00 7963 2 : SUBROUTINE VORO 3 : C ****************************************************************** 4 : C INTERFACE CASTEM 2000 5 : C 6 : C 7 : C TAB1 = VORO MAIL1 MAIL2 ; 8 : C 9 : C 10 : C OBJET : 11 : C _______ 12 : C 13 : C L'OPERATEUR VORONOI CONSTRUIT LA PARTITION DE VORONOI D'UN 14 : C ENSEMBLE DE POINTS MAIL1 RESTREINTE PAR UN CONTOUR/ENVELOPPE MAIL2 15 : C 16 : C ****************************************************************** 17 : IMPLICIT INTEGER(I-N) 18 : IMPLICIT REAL*8 (A-H,O-Z) 19 : C 20 : -INC SMCHPOI 21 : -INC CCOPTIO 22 : -INC SMELEME 23 : -INC SMCOORD 24 : -INC CCGEOME 25 : -INC SMLENTI 26 : -INC SMTABLE 27 : C 28 : DIMENSION LNBOIT(8) 29 : C 30 : SEGMENT,MCIRCONS 31 : REAL*8 TRC(NBE1,4) 32 : ENDSEGMENT 33 : POINTEUR ITRC1.MCIRCONS 34 : C 35 : SEGMENT,MADJACEN 36 : INTEGER LEAC(NBL1,IDIM+1,2) 37 : ENDSEGMENT 38 : POINTEUR ILEA1.MADJACEN 39 : C 40 : C ======================= 41 : C --- 1.LECTURE DES DONNEES --- 42 : C ======================= 43 : C 44 : C LECTURE DES OBJETS COURANTS (ENTREES) 45 : C ===================================== 45 : CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU) 46 : IF (IERR.NE.0) THEN 47 : C ON A PAS TROUVE LE MAILLAGE 47 : CALL ERREUR(503) 48 : GOTO 999 49 : ENDIF 50 : SEGACT,IPT1 51 : NBSZ1=IPT1.LISOUS(/1) 52 : NTYP1=IPT1.ITYPEL 53 : C 53 : CALL LIROBJ('MAILLAGE',IPT2,1,IRETOU) 54 : IF (IERR.NE.0) THEN 55 : C ON A PAS TROUVE LE MAILLAGE 56 : SEGDES,IPT1 56 : CALL ERREUR(503) 57 : GOTO 999 58 : ENDIF 59 : SEGACT,IPT2 60 : NBSZ2=IPT2.LISOUS(/1) 61 : NTYP2=IPT2.ITYPEL 62 : C 63 : C Parametre optionnel : critere de distance pour considerer deux 64 : C centres de cercles circonscrits confondus 65 : C Par defaut, on le calcule selon la taille du nauge de points 66 : IDIMP1=IDIM+1 67 : DDMAX=0. 68 : IP1=IPT1.NUM(1,1) 69 : XP1=XCOOR((IP1-1)*IDIMP1+1) 70 : YP1=XCOOR((IP1-1)*IDIMP1+2) 71 : ZP1=XCOOR((IP1-1)*IDIMP1+3) 72 : XP2=XP1 73 : YP2=YP1 74 : ZP2=ZP1 75 : DO I=2,IPT1.NUM(/2) 76 : IPI=IPT1.NUM(1,I) 77 : XPI=XCOOR((IPI-1)*IDIMP1+1) 78 : YPI=XCOOR((IPI-1)*IDIMP1+2) 79 : ZPI=XCOOR((IPI-1)*IDIMP1+3) 80 : IF(XPI.LT.XP1) XP1=XPI 81 : IF(YPI.LT.YP1) YP1=YPI 82 : IF(ZPI.LT.ZP1) ZP1=ZPI 83 : IF(XPI.GT.XP2) XP2=XPI 84 : IF(YPI.GT.YP2) YP2=YPI 85 : IF(ZPI.GT.ZP2) ZP2=ZPI 86 : ENDDO 87 : DDMAX=(XP2-XP1)+(YP2-YP1) 88 : IF (IDIM.EQ.3) DDMAX=DDMAX+(ZP2-ZP1) 89 : DDMAX=DDMAX*0.333333 90 : C WRITE(6,*) 'DDMAX=', DDMAX 91 : C 92 : C---- SI MAILLAGE DE POI1 ET CONTOUR/ENVELOPPE : PARTITION DE VORONOI -- 93 : C 94 : C L'operateur n'est disponible qu'en dimension 2 ou 3 95 : IF ((IDIM.NE.2).AND.(IDIM.NE.3)) THEN 96 : INTERR(1)=IDIM 97 : C FONCTION INDISPONIBLE EN DIMENSION %I1 98 : SEGDES,IPT1,IPT2 98 : CALL ERREUR(709) 99 : GOTO 999 100 : ENDIF 101 : C Test sur les maillages d'entree 102 : C --> Le nombre de sous zones doit etre nul 103 : IF ((NBSZ1.NE.0).OR.(NBSZ2.NE.0)) THEN 104 : C MAILLAGE INCORRECT 105 : SEGDES,IPT1,IPT2 105 : CALL ERREUR(426) 106 : GOTO 999 107 : ENDIF 108 : C --> Le maillage MAIL1 doit etre constitue d'elements POI1 109 : IF (NTYP1.NE.1) THEN 110 : C TYPE D'ELEMENTS INCORRECT 111 : SEGDES,IPT1,IPT2 111 : CALL ERREUR(16) 112 : GOTO 999 113 : ENDIF 114 : C --> Le maillage MAIL2 doit etre constitue d'elements SEG2 ou TRI3 115 : IF (((IDIM.EQ.2).AND.(NTYP2.NE.2)).OR. 116 : & ((IDIM.EQ.3).AND.(NTYP2.NE.4))) THEN 117 : C TYPE D'ELEMENTS INCORRECT 118 : SEGDES,IPT1,IPT2 118 : CALL ERREUR(16) 119 : GOTO 999 120 : ENDIF 121 : C --> Les elements du maillage MAIL2 doivent etre orientes de la 122 : C meme maniere (appel a l'operateur VERSENS) 122 : CALL ECROBJ('MAILLAGE',IPT2) 122 : CALL VERSEN 122 : CALL LIROBJ('MAILLAGE',IPT2,1,IRETOU) 123 : IF (IERR.NE.0) THEN 124 : SEGDES,IPT1 125 : GOTO 999 126 : ENDIF 127 : C --> Le maillage contour/enveloppe MAIL2 doit etre connexe 127 : CALL ECROBJ('MAILLAGE',IPT2) 127 : CALL CCON 127 : CALL DIMENS 127 : CALL LIRENT(NCCON,1,IRETOU) 128 : IF (NCCON.NE.1) THEN 129 : INTERR(1)=NCCON 130 : C CET OBJET CONTIENT %1 COMPOSANTES CONNEXES 131 : SEGDES,IPT1,IPT2 131 : CALL ERREUR(754) 132 : GOTO 999 133 : ENDIF 134 : C 135 : MPOVAL = 0 136 : MCHPOI = 0 136 : CALL LIROBJ('CHPOINT ',MCHPOI,0,IRETOU) 137 : IF(IRETOU.EQ.1) THEN 138 : SEGACT MCHPOI 139 : MSOUPO=IPCHP(1) 140 : SEGACT MSOUPO 141 : MPOVAL = IPOVAL 142 : SEGACT MPOVAL 143 : ENDIF 144 : C 145 : C On recree le maillage de points avec une reference vers le 146 : C maillage du contour/enveloppe pour etre pris en compte dans DELAUN 147 : C lors du calcul de la taille de la boite de triangulation 148 : SEGINI,IPT3=IPT1 149 : NBNN=IPT3.NUM(/1) 150 : NBELEM=IPT3.NUM(/2) 151 : NBSOUS=0 152 : NBREF=1 153 : SEGADJ,IPT3 154 : IPT3.LISREF(1)=IPT2 155 : C Appel a DELAUN avec IBOI=1 pour le calcul de la triangulation 156 : C de Delaunay de IPT3 157 : XBOI=50. 158 : IBOI=1 158 : CALL DELAUN(MPOVAL,IPT3,XBOI,IBOI,IPT4,LNBOIT,ITRC1,ILEA1) 159 : IF (IERR.NE.0) GOTO 999 160 : SEGACT,IPT4 161 : C Appel a VORO1 pour le calcul de la partition de Voronoi non coupee 161 : CALL VORO1(IPT1,IPT4,LNBOIT,ILEA1,ITRC1,DDMAX,MTAB1) 162 : IF (IERR.NE.0) GOTO 999 163 : IF (MPOVAL.NE.0) SEGDES,MPOVAL,MSOUPO,MCHPOI 164 : C Decoupage de la partition de Voronoi selon le contour/enveloppe 165 : C avec VORO2 165 : CALL VORO2(IPT1,IPT2,DDMAX,MCHPOI,MTAB1) 166 : IF (IERR.NE.0) GOTO 999 167 : C Ecriture sortie/fin 168 : SEGSUP,IPT3,IPT4,ITRC1,ILEA1 169 : SEGDES,IPT1,IPT2 169 : IF (IERR.EQ.0) CALL ECROBJ('TABLE ',MTAB1) 170 : GOTO 999 171 : C 172 : C 173 : 999 RETURN 174 : END 175 :
© Cast3M 2003 - All rights reserved.
Disclaimer