1 : C VERCON SOURCE PASCAL 12/01/24 21:15:18 7263 2 : SUBROUTINE VERCON(IPT1,ICONV) 3 : C----------------------------------------------------------------------- 4 : C 5 : C VERIFICATION DE LA CONVEXITE D'UN MAILLAGE (2D, 3D) 6 : C 7 : C----------------------------------------------------------------------- 8 : C 9 : C Appelé par VERMAI 10 : C 11 : C Entrée : 12 : C IPT1 : pointeur sur le maillage (MELEME) a tester. 13 : C Ce maillage doit respecter certaines conditions: 14 : C - en dimension 2, le contour doit etre fait de SEG2 15 : C - en dimension 3, l'enveloppe doit etre faite de TRI3 16 : C 17 : C Sortie : 18 : C ICONV : entier indiquant si le maillage est convexe 19 : C = 1 si le maillage est convexe 20 : C = 0 sinon 21 : C 22 : C Un message est emis si le contour, ou de l'enveloppe, ne 23 : C respecte pas les conditions requises 24 : C 25 : C 26 : C Note : en dimension 2, on ne supprime pas le MELEME IPT2 (contour) 27 : C car lors de l'appel a PRCONT, le contour IPT2 est reference 28 : C dans le maillage IPT1 29 : 30 : C----------------------------------------------------------------------- 31 : IMPLICIT INTEGER(I-N) 32 : IMPLICIT REAL*8(A-H,O-Z) 33 : -INC SMCOORD 34 : -INC CCOPTIO 35 : -INC SMELEME 36 : -INC SMLENTI 37 : C 38 : SEGMENT,MADJACEN 39 : INTEGER LEAC(NBL1,IDIM,2) 40 : ENDSEGMENT 41 : POINTEUR ILEA1.MADJACEN 42 : C 43 : POINTEUR LNPT1.MLENTI 44 : C 45 : C Tolerance pour la planeite des element adjacents 46 : ZERO=1.E-14 47 : C 48 : C 1. CONSTRUCTION DU CONTOUR OU DE L'ENVELOPPE DU MAILLAGE 49 : C 50 : C J'ecris le maillage dans la pile 50 : CALL ECROBJ('MAILLAGE',IPT1) 51 : IF (IERR.NE.0) RETURN 52 : C puis je fais appel a l'operateur ad hoc pour calculer : 53 : C - soit le contour, en dimension 2 54 : C - soit l'enveloppe, en dimension 3 55 : IF (IDIM.EQ.2) THEN 55 : CALL PRCONT 56 : ELSEIF (IDIM.EQ.3) THEN 56 : CALL ENVVOL 57 : ENDIF 58 : C Je lis le contour/enveloppe dans la pile, je le designe par IPT2 58 : CALL LIROBJ('MAILLAGE',IPT2,1,IRETOU) 59 : IF (IERR.NE.0) RETURN 60 : SEGACT,IPT2 61 : C 62 : C 2. TESTS SUR LE MAILLAGE DU CONTOUR/ENVELOPPE 63 : C 64 : C Test si le type d'elements est bien celui attendu, a savoir: 65 : C des TRI3 en dimension 3 ou bien des SEG2 en dimension 2 66 : IF (IDIM.EQ.2) THEN 67 : IF ((IPT2.ITYPEL).NE.2) THEN 68 : ICONV=0 69 : MOTERR(1:8) = 'utilise ' 69 : CALL ERREUR(926) 70 : RETURN 71 : ENDIF 72 : ELSEIF (IDIM.EQ.3) THEN 73 : IF ((IPT2.ITYPEL).NE.4) THEN 74 : ICONV=0 75 : MOTERR(1:8) = 'utilise ' 75 : CALL ERREUR(926) 76 : RETURN 77 : ENDIF 78 : ENDIF 79 : C 80 : C 3. VERIFICATION DE LA CONVEXITE 81 : C 82 : C 3.1 EN DIMENSION 2 83 : C On parcoure le contour de proche en proche et on verifie la 84 : C convexite: le contour doit tourner toujours dans le meme sens 85 : C 86 : IF (IDIM.EQ.2) THEN 87 : I1=0 88 : JG=IPT2.NUM(/2) 89 : SEGINI,LNPT1 90 : C On part du premier element du contour 91 : NEL1=1 92 : NNOA1=1 93 : NNOB1=2 94 : C Boucle sur les elements (indice I) 95 : DO I=1,JG 96 : LNPT1.LECT(I)=1 97 : C Coordonnees des points A et B, extremites de l'element I 98 : NREFA=IPT2.NUM(NNOA1,NEL1) 99 : XA=XCOOR((NREFA-1)*(IDIM+1)+1) 100 : YA=XCOOR((NREFA-1)*(IDIM+1)+2) 101 : NREFB=IPT2.NUM(NNOB1,NEL1) 102 : XB=XCOOR((NREFB-1)*(IDIM+1)+1) 103 : YB=XCOOR((NREFB-1)*(IDIM+1)+2) 104 : C Le voisin du dernier element est le premier 105 : IF (I.EQ.JG) THEN 106 : NEL2=1 107 : NNOB2=1 108 : NNOC2=2 109 : GOTO 10 110 : ENDIF 111 : C Boucle sur les autres elements du contour 112 : C on recherche le voisin de I partageant le point B 113 : DO J=1,JG 114 : C Si on a deja traite l'element J, alors on itere 115 : IF (LNPT1.LECT(J).EQ.1) GOTO 20 116 : C Boucle sur les noeuds de l'element J 117 : DO NNOB2=1,IDIM 118 : C On teste si le noeud NNOB2 correspond au noeud B 119 : IF ((IPT2.NUM(NNOB2,J)).EQ.NREFB) THEN 120 : NEL2=J 121 : IF (NNOB2.EQ.1) NNOC2=2 122 : IF (NNOB2.EQ.2) NNOC2=1 123 : GOTO 10 124 : ENDIF 125 : ENDDO 126 : 20 CONTINUE 127 : ENDDO 128 : 10 CONTINUE 129 : C Coordonnees du point C 130 : NREFC=IPT2.NUM(NNOC2,NEL2) 131 : XC=XCOOR((NREFC-1)*(IDIM+1)+1) 132 : YC=XCOOR((NREFC-1)*(IDIM+1)+2) 133 : C Calcul de la composante Z du produit vectoriel AB^BC 134 : C => on obtient ainsi le sens de rotation 135 : ZPV=((XB-XA)*(YC-YB))-((YB-YA)*(XC-XB)) 136 : X0A=((XA**2)+(YA**2))**0.5 137 : X0B=((XB**2)+(YC**2))**0.5 138 : X0C=((XC**2)+(YC**2))**0.5 139 : ZPV1=ZPV/(X0A*X0B*X0C) 140 : C Determination du sens de rotation initial 141 : IF (I1.EQ.0) THEN 142 : C Si le produit vectoriel est nul, on ne peut pas connaitre le 143 : C sens initial, alors on passe a l'element suivant 144 : IF (ABS(ZPV1).LT.ZERO) THEN 145 : GOTO 30 146 : ELSE 147 : SIG_0=SIGN(1D0,ZPV) 148 : I1=1 149 : ENDIF 150 : ELSE 151 : C Si le sens de rotation change, par rapport a l'initial, 152 : C alors le maillage n'est pas convexe => on sort 153 : IF (ABS(ZPV1).LT.ZERO) THEN 154 : GOTO 30 155 : ELSE 156 : SIG_I=SIGN(1D0,ZPV) 157 : IF (SIG_I.NE.SIG_0) THEN 158 : ICONV=0 159 : SEGSUP,LNPT1 160 : RETURN 161 : ENDIF 162 : ENDIF 163 : ENDIF 164 : 30 CONTINUE 165 : C On passe a l'element suivant, B devient A 166 : NEL1=NEL2 167 : NNOA1=NNOB2 168 : IF (NNOA1.EQ.1) NNOB1=2 169 : IF (NNOA1.EQ.2) NNOB1=1 170 : ENDDO 171 : ICONV=1 172 : SEGSUP,LNPT1 173 : RETURN 174 : C 175 : C 3.2 EN DIMENSION 3 176 : C On parcoure toutes les facettes de l'enveloppe et on recherche 177 : C les 3 elements adjacents. Une fois trouves, on teste la 178 : C convexite: les noeuds des elements adjacents doivent tous etre 179 : C du meme cote du plan definit par la facette 180 : C 181 : ELSEIF (IDIM.EQ.3) THEN 182 : I1=0 183 : NBL1=IPT2.NUM(/2) 184 : SEGINI,ILEA1 185 : C Boucle sur les elements (indice I) 186 : DO I=1,NBL1 187 : C Numeros des noeuds A, B et C de l'element I 188 : NREFA=IPT2.NUM(1,I) 189 : NREFB=IPT2.NUM(2,I) 190 : NREFC=IPT2.NUM(3,I) 191 : C Numeros des adjacents et des noeuds D, E et F 192 : NELD=0 193 : NELE=0 194 : NELF=0 195 : NREFD=0 196 : NREFE=0 197 : NREFF=0 198 : C On part a la recherche des adjacents de l'element I 199 : C Si des indices de la table de connectivite existent deja, 200 : C alors on va chercher les numeros des elements adjacents et 201 : C ceux des noeuds D, E ou F 202 : IF ((ILEA1.LEAC(I,1,1)).NE.0) THEN 203 : NELD=ILEA1.LEAC(I,1,1) 204 : NREFD=IPT2.NUM((ILEA1.LEAC(I,1,2)),NELD) 205 : ENDIF 206 : IF ((ILEA1.LEAC(I,2,1)).NE.0) THEN 207 : NELE=ILEA1.LEAC(I,2,1) 208 : NREFE=IPT2.NUM((ILEA1.LEAC(I,2,2)),NELE) 209 : ENDIF 210 : IF ((ILEA1.LEAC(I,3,1)).NE.0) THEN 211 : NELF=ILEA1.LEAC(I,3,1) 212 : NREFF=IPT2.NUM((ILEA1.LEAC(I,3,2)),NELF) 213 : ENDIF 214 : C Boucle sur les elements d'indice superieurs a I (indice II) 215 : DO II=I+1,NBL1 216 : C Si l'element II est un des adjacents connus, on itere 217 : IF ((II.EQ.NELD).OR.(II.EQ.NELE).OR.(II.EQ.NELF)) GOTO 40 218 : C Indicateurs si l'un des noeuds de II correspond a A, B ou C 219 : IA=0 220 : IB=0 221 : IC=0 222 : C Boucle sur les noeuds (indice JJ) de l'element II 223 : DO JJ=1,IDIM 224 : NREFX=IPT2.NUM(JJ,II) 225 : C On teste si le noeud JJ correspond aux noeuds A, B ou C 226 : IF (NREFX.EQ.NREFA) IA=JJ 227 : IF (NREFX.EQ.NREFB) IB=JJ 228 : IF (NREFX.EQ.NREFC) IC=JJ 229 : C A t on trouve l'adjacent vu par A ? 230 : IF ((IB.NE.0).AND.(IC.NE.0)) THEN 231 : IF ((IB+IC).EQ.3) NNOJJ=3 232 : IF ((IB+IC).EQ.4) NNOJJ=2 233 : IF ((IB+IC).EQ.5) NNOJJ=1 234 : NREFD=IPT2.NUM(NNOJJ,II) 235 : ILEA1.LEAC(I,1,1)=II 236 : ILEA1.LEAC(I,1,2)=NNOJJ 237 : ILEA1.LEAC(II,NNOJJ,1)=I 238 : ILEA1.LEAC(II,NNOJJ,2)=1 239 : ENDIF 240 : C A t on trouve l'adjacent vu par B ? 241 : IF ((IA.NE.0).AND.(IC.NE.0)) THEN 242 : IF ((IA+IC).EQ.3) NNOJJ=3 243 : IF ((IA+IC).EQ.4) NNOJJ=2 244 : IF ((IA+IC).EQ.5) NNOJJ=1 245 : NREFE=IPT2.NUM(NNOJJ,II) 246 : ILEA1.LEAC(I,2,1)=II 247 : ILEA1.LEAC(I,2,2)=NNOJJ 248 : ILEA1.LEAC(II,NNOJJ,1)=I 249 : ILEA1.LEAC(II,NNOJJ,2)=2 250 : ENDIF 251 : C A t on trouve l'adjacent vu par C ? 252 : IF ((IA.NE.0).AND.(IB.NE.0)) THEN 253 : IF ((IA+IB).EQ.3) NNOJJ=3 254 : IF ((IA+IB).EQ.4) NNOJJ=2 255 : IF ((IA+IB).EQ.5) NNOJJ=1 256 : NREFF=IPT2.NUM(NNOJJ,II) 257 : ILEA1.LEAC(I,3,1)=II 258 : ILEA1.LEAC(I,3,2)=NNOJJ 259 : ILEA1.LEAC(II,NNOJJ,1)=I 260 : ILEA1.LEAC(II,NNOJJ,2)=3 261 : ENDIF 262 : ENDDO 263 : 40 CONTINUE 264 : C A t on trouve tous les adjacents ? 265 : IF ((NREFD*NREFE*NREFF).NE.0) GOTO 50 266 : ENDDO 267 : C On a trouve les 3 voisins de l'element I 268 : 50 CONTINUE 269 : C Coordonnees des noeuds A, B, C, D, E et F 270 : XA=XCOOR((NREFA-1)*(IDIM+1)+1) 271 : YA=XCOOR((NREFA-1)*(IDIM+1)+2) 272 : ZA=XCOOR((NREFA-1)*(IDIM+1)+3) 273 : XB=XCOOR((NREFB-1)*(IDIM+1)+1) 274 : YB=XCOOR((NREFB-1)*(IDIM+1)+2) 275 : ZB=XCOOR((NREFB-1)*(IDIM+1)+3) 276 : XC=XCOOR((NREFC-1)*(IDIM+1)+1) 277 : YC=XCOOR((NREFC-1)*(IDIM+1)+2) 278 : ZC=XCOOR((NREFC-1)*(IDIM+1)+3) 279 : XD=XCOOR((NREFD-1)*(IDIM+1)+1) 280 : YD=XCOOR((NREFD-1)*(IDIM+1)+2) 281 : ZD=XCOOR((NREFD-1)*(IDIM+1)+3) 282 : XE=XCOOR((NREFE-1)*(IDIM+1)+1) 283 : YE=XCOOR((NREFE-1)*(IDIM+1)+2) 284 : ZE=XCOOR((NREFE-1)*(IDIM+1)+3) 285 : XF=XCOOR((NREFF-1)*(IDIM+1)+1) 286 : YF=XCOOR((NREFF-1)*(IDIM+1)+2) 287 : ZF=XCOOR((NREFF-1)*(IDIM+1)+3) 288 : C Calcul d'un vecteur normal au plan ABC, on prend le produit 289 : C vectoriel n = AB^AC 290 : XN=((YB-YA)*(ZC-ZA))-((ZB-ZA)*(YC-YA)) 291 : YN=((ZB-ZA)*(XC-XA))-((XB-XA)*(ZC-ZA)) 292 : ZN=((XB-XA)*(YC-YA))-((YB-YA)*(XC-XA)) 293 : C Calcul des produits scalaires n.AD n.AE et n.AF pour connaitre 294 : C la position de D, E et F par rapport au plan ABC 295 : PSD=(XN*(XD-XA))+(YN*(YD-YA))+(ZN*(ZD-ZA)) 296 : PSE=(XN*(XE-XA))+(YN*(YE-YA))+(ZN*(ZE-ZA)) 297 : PSF=(XN*(XF-XA))+(YN*(YF-YA))+(ZN*(ZF-ZA)) 298 : X0N=((XN**2)+(YN**2)+(ZN**2))**0.5 299 : X0AD=(((XD-XA)**2)+((YD-YA)**2)+((ZD-ZA)**2))**0.5 300 : X0AE=(((XE-XA)**2)+((YE-YA)**2)+((ZE-ZA)**2))**0.5 301 : X0AF=(((XF-XA)**2)+((YF-YA)**2)+((ZF-ZA)**2))**0.5 302 : PSD1=PSD/(X0N*X0AD) 303 : PSE1=PSE/(X0N*X0AE) 304 : PSF1=PSF/(X0N*X0AF) 305 : C Cas particulier ou le point D est dans le plan ABC 306 : IF (ABS(PSD1).LT.ZERO) THEN 307 : SIGD=0. 308 : ELSE 309 : SIGD=SIGN(1D0,PSD) 310 : ENDIF 311 : C Cas particulier ou le point E est dans le plan ABC 312 : IF (ABS(PSE1).LT.ZERO) THEN 313 : SIGE=0. 314 : ELSE 315 : SIGE=SIGN(1D0,PSE) 316 : ENDIF 317 : C Cas particulier ou le point F est dans le plan ABC 318 : IF (ABS(PSF1).LT.ZERO) THEN 319 : SIGF=0. 320 : ELSE 321 : SIGF=SIGN(1D0,PSF) 322 : ENDIF 323 : C Si les 3 points ne sont pas tous situe du meme cote du plan 324 : C ABC, alors le maillage n'est pas convexe => on sort 325 : C - si D seulement est dans le plan ABC, on teste juste E et F 326 : IF ((SIGD.EQ.0.).AND.(SIGE.NE.0.).AND.(SIGF.NE.0.)) THEN 327 : IF (SIGE.NE.SIGF) THEN 328 : ICONV=0 329 : SEGSUP,ILEA1,IPT2 330 : RETURN 331 : ENDIF 332 : ENDIF 333 : C - si E seulement est dans le plan ABC, on teste juste D et F 334 : IF ((SIGE.EQ.0.).AND.(SIGD.NE.0.).AND.(SIGF.NE.0.)) THEN 335 : IF (SIGD.NE.SIGF) THEN 336 : ICONV=0 337 : SEGSUP,ILEA1,IPT2 338 : RETURN 339 : ENDIF 340 : ENDIF 341 : C - si F seulement est dans le plan ABC, on teste juste D et E 342 : IF ((SIGF.EQ.0.).AND.(SIGD.NE.0.).AND.(SIGE.NE.0.)) THEN 343 : IF (SIGD.NE.SIGE) THEN 344 : ICONV=0 345 : SEGSUP,ILEA1,IPT2 346 : RETURN 347 : ENDIF 348 : ENDIF 349 : C - cas general (D, E et F hors du plan ABC) 350 : IF ((SIGD.NE.0.).AND.(SIGE.NE.0.).AND.(SIGF.NE.0.)) THEN 351 : IF ((SIGD.NE.SIGE).OR.(SIGE.NE.SIGF)) THEN 352 : ICONV=0 353 : SEGSUP,ILEA1,IPT2 354 : RETURN 355 : ENDIF 356 : ENDIF 357 : C - dans tous les autres cas (2 ou 3 des points D, E, F sont 358 : C dans le plan ABC), on itere 359 : ENDDO 360 : ICONV=1 361 : SEGSUP,ILEA1,IPT2 362 : ENDIF 363 : RETURN 364 : END 365 : 366 :
© Cast3M 2003 - All rights reserved.
Disclaimer