1 : C DELAUN SOURCE PASCAL 14/02/17 21:15:00 7965 2 : SUBROUTINE DELAUN(MPOVAL,IPT1,XBOI,IBOI,IPT2,LNBOIT,ITRC1,ILEA1) 3 : C----------------------------------------------------------------------C 4 : C Triangulation d'un maillage de points (POI1), 1D, 2D ou 3D, C 5 : C suivant l'agorithme de DELAUNAY donne dans P.-L. George, C 6 : C Delaunay triangulation and meshing, Hermes, 1998, p. 41. C 7 : C C 8 : C IPT1 = pointeur sur le maillage de POI1 en entree ; C 9 : C suppose ACTIF en ENTREE et en SORTIE C 10 : C XBOI = rapport entre le "diametre" du nuage de points (IPT1) et C 11 : C la taille de la boite de triangulation ; C 12 : C IBOI = indicateur permettant d'indiquer si on veut "garder" la C 13 : C boite de triangulation : C 14 : C - IBOI = 0 : on ne la garde pas (defaut) ; C 15 : C - IBOI = 1 : on la garde ; C 16 : C C 17 : C IPT2 = pointeur sur le maillage triangule en sortie C 18 : C LNBOIT = liste des noeuds formant la boite de triangulation C 19 : C ITRC1 = table contenant les centres des cercles circonscrits C 20 : C ILEA1 = table des d'adjacence de la triangulation C 21 : C C 22 : C----------------------------------------------------------------------C 23 : C Version V0 : centre et rayon cercles circons. calcules chaque fois C 24 : C Version V1 : centre et rayon cercles circons. stockes dans SEGMENT C 25 : C MCIRCONS + elements adjacents stockes dans SEGMENT C 26 : C ILEA1 C 27 : C----------------------------------------------------------------------C 28 : IMPLICIT INTEGER(I-N) 29 : IMPLICIT REAL*8 (A-H,O-Z) 30 : c 31 : -INC SMCHPOI 32 : -INC CCOPTIO 33 : -INC SMELEME 34 : -INC SMCOORD 35 : -INC CCGEOME 36 : -INC SMLENTI 37 : -INC SMTABLE 38 : PARAMETER (ICMIN=1000) 39 : PARAMETER (ZERO=1D-14) 40 : DIMENSION LNBOIT(8) 41 : DIMENSION XCMIMA(2,3),X0(4,3),XCMIL(3),XPI1(3),XCBJ(3),XG0(3) 42 : LOGICAL BMOD,BHORS,BFF 43 : C 44 : POINTEUR LCAVI.MLENTI,LBOUL.MLENTI,LBASE.MLENTI,LREF1.MLENTI 45 : POINTEUR LEADJ.MLENTI,LNVADJ.MLENTI,LSUP.MLENTI,LSCP.MLENTI 46 : C 47 : SEGMENT,MCIRCONS 48 : REAL*8 TRC(NBE1,4) 49 : ENDSEGMENT 50 : POINTEUR ITRC1.MCIRCONS 51 : C 52 : SEGMENT,MADJACEN 53 : INTEGER LEAC(NBL1,IDIM+1,2) 54 : ENDSEGMENT 55 : POINTEUR ILEA1.MADJACEN 56 : C 57 : SEGMENT,MLIAIS 58 : INTEGER LIAS(NB1,NB1) 59 : ENDSEGMENT 60 : POINTEUR ITL1.MLIAIS 61 : C 62 : SEGMENT RAYONOE 63 : REAL*8 RAYN(NBNER) 64 : ENDSEGMENT 65 : POINTEUR IRN1.RAYONOE 66 : c 67 : C----------------------------------------------------------------------C 68 : C PARTIE 1 C 69 : C Construction de la boite englobant tous les points C 70 : C----------------------------------------------------------------------C 71 : C 72 : C---- 1.1 Je determine la taille du maillage 73 : SEGACT,MCOORD 74 : IDIMP1=IDIM+1 75 : IPT2=IPT1 76 : IF ((IBOI.EQ.1).AND.(IPT1.LISREF(/1).NE.0)) THEN 77 : IPT2=IPT1.LISREF(1) 78 : SEGACT,IPT2 79 : ENDIF 80 : IREF1=(IPT1.NUM(1,1)-1)*IDIMP1 81 : DDMAX=0. 82 : DO I=1,IDIM 83 : XCMIMA(1,I)=XCOOR(IREF1+I) 84 : XCMIMA(2,I)=XCOOR(IREF1+I) 85 : DO J=2,IPT1.NUM(/2) 86 : IREFJ=(IPT1.NUM(1,J)-1)*IDIMP1 87 : XCNJ1=XCOOR(IREFJ+I) 88 : IF (XCNJ1.LT.XCMIMA(1,I)) XCMIMA(1,I)=XCNJ1 89 : IF (XCNJ1.GT.XCMIMA(2,I)) XCMIMA(2,I)=XCNJ1 90 : ENDDO 91 : IF ((IBOI.EQ.1).AND.(IPT1.LISREF(/1).NE.0)) THEN 92 : DO J=1,IPT2.NUM(/2) 93 : IREFJ=(IPT2.NUM(1,J)-1)*IDIMP1 94 : XCNJ1=XCOOR(IREFJ+I) 95 : IF (XCNJ1.LT.XCMIMA(1,I)) XCMIMA(1,I)=XCNJ1 96 : IF (XCNJ1.GT.XCMIMA(2,I)) XCMIMA(2,I)=XCNJ1 97 : ENDDO 98 : ENDIF 99 : DDMAXI=0.5*(XCMIMA(2,I)-XCMIMA(1,I)) 100 : IF (DDMAXI.GT.DDMAX) DDMAX=DDMAXI 101 : ENDDO 102 : C 103 : C---- 1.2 On construit la boite 104 : C 1.2.1 Points coins de la boite 105 : DBOIT=XBOI*DDMAX 106 : NBPTS0=XCOOR(/1)/IDIMP1 107 : NBPTS=NBPTS0+(2**IDIM) 108 : SEGADJ,MCOORD 109 : IF (IDIM.EQ.1) THEN 110 : XCMIL(1)=0.5*(XCMIMA(1,1)+XCMIMA(2,1)) 111 : IREF1=(NBPTS-2)*IDIMP1 112 : IREF2=(NBPTS-1)*IDIMP1 113 : XCOOR(IREF1+1)=XCMIL(1)-DBOIT 114 : XCOOR(IREF2+1)=XCMIL(1)+DBOIT 115 : XCOOR(IREF1+2)=0. 116 : XCOOR(IREF2+2)=0. 117 : ELSEIF (IDIM.EQ.2) THEN 118 : DO I=1,IDIM 119 : XCMIL(I)=0.5*(XCMIMA(1,I)+XCMIMA(2,I)) 120 : ENDDO 121 : IREF1=(NBPTS-4)*IDIMP1 122 : IREF2=(NBPTS-3)*IDIMP1 123 : IREF3=(NBPTS-2)*IDIMP1 124 : IREF4=(NBPTS-1)*IDIMP1 125 : XCOOR(IREF1+1)=XCMIL(1)-DBOIT 126 : XCOOR(IREF2+1)=XCMIL(1)+DBOIT 127 : XCOOR(IREF3+1)=XCMIL(1)+DBOIT 128 : XCOOR(IREF4+1)=XCMIL(1)-DBOIT 129 : XCOOR(IREF1+2)=XCMIL(2)-DBOIT 130 : XCOOR(IREF2+2)=XCMIL(2)-DBOIT 131 : XCOOR(IREF3+2)=XCMIL(2)+DBOIT 132 : XCOOR(IREF4+2)=XCMIL(2)+DBOIT 133 : XCOOR(IREF1+3)=0. 134 : XCOOR(IREF2+3)=0. 135 : XCOOR(IREF3+3)=0. 136 : XCOOR(IREF4+3)=0. 137 : ELSEIF (IDIM.EQ.3) THEN 138 : DO I=1,IDIM 139 : XCMIL(I)=0.5*(XCMIMA(1,I)+XCMIMA(2,I)) 140 : ENDDO 141 : IREF1=(NBPTS-8)*IDIMP1 142 : IREF2=(NBPTS-7)*IDIMP1 143 : IREF3=(NBPTS-6)*IDIMP1 144 : IREF4=(NBPTS-5)*IDIMP1 145 : IREF5=(NBPTS-4)*IDIMP1 146 : IREF6=(NBPTS-3)*IDIMP1 147 : IREF7=(NBPTS-2)*IDIMP1 148 : IREF8=(NBPTS-1)*IDIMP1 149 : XCOOR(IREF1+1)=XCMIL(1)-DBOIT 150 : XCOOR(IREF2+1)=XCMIL(1)+DBOIT 151 : XCOOR(IREF3+1)=XCMIL(1)+DBOIT 152 : XCOOR(IREF4+1)=XCMIL(1)-DBOIT 153 : XCOOR(IREF5+1)=XCMIL(1)-DBOIT 154 : XCOOR(IREF6+1)=XCMIL(1)+DBOIT 155 : XCOOR(IREF7+1)=XCMIL(1)+DBOIT 156 : XCOOR(IREF8+1)=XCMIL(1)-DBOIT 157 : XCOOR(IREF1+2)=XCMIL(2)-DBOIT 158 : XCOOR(IREF2+2)=XCMIL(2)-DBOIT 159 : XCOOR(IREF3+2)=XCMIL(2)+DBOIT 160 : XCOOR(IREF4+2)=XCMIL(2)+DBOIT 161 : XCOOR(IREF5+2)=XCMIL(2)-DBOIT 162 : XCOOR(IREF6+2)=XCMIL(2)-DBOIT 163 : XCOOR(IREF7+2)=XCMIL(2)+DBOIT 164 : XCOOR(IREF8+2)=XCMIL(2)+DBOIT 165 : XCOOR(IREF1+3)=XCMIL(3)-DBOIT 166 : XCOOR(IREF2+3)=XCMIL(3)-DBOIT 167 : XCOOR(IREF3+3)=XCMIL(3)-DBOIT 168 : XCOOR(IREF4+3)=XCMIL(3)-DBOIT 169 : XCOOR(IREF5+3)=XCMIL(3)+DBOIT 170 : XCOOR(IREF6+3)=XCMIL(3)+DBOIT 171 : XCOOR(IREF7+3)=XCMIL(3)+DBOIT 172 : XCOOR(IREF8+3)=XCMIL(3)+DBOIT 173 : XCOOR(IREF1+4)=0. 174 : XCOOR(IREF2+4)=0. 175 : XCOOR(IREF3+4)=0. 176 : XCOOR(IREF4+4)=0. 177 : XCOOR(IREF5+4)=0. 178 : XCOOR(IREF6+4)=0. 179 : XCOOR(IREF7+4)=0. 180 : XCOOR(IREF8+4)=0. 181 : ENDIF 182 : C 183 : C 1.2.2 Maillage de la boite 184 : NBSOUS=0 185 : NBREF=0 186 : NNOE1=NBPTS-(2**IDIM)+1 187 : IF (IDIM.EQ.1) THEN 188 : C Creation des elements SEG2 en definissant les noeuds 189 : NBNN=2 190 : NBELEM=IPT1.NUM(/2)+1 191 : SEGINI,IPT2 192 : IPT2.ITYPEL=2 193 : NBELE2=1 194 : IPT2.NUM(1,1)=NNOE1 195 : IPT2.NUM(2,1)=NNOE1+1 196 : C Liste des numeros des noeuds de la boite 197 : NNBOIT=2 198 : LNBOIT(1)=NNOE1 199 : LNBOIT(2)=NNOE1+1 200 : C Initialisation du segment ILEA1 contenant la table de 201 : C connectivite des elements de IPT2 202 : NBL1=NBELEM 203 : SEGINI,ILEA1 204 : ILEA1.LEAC(1,1,1)= 0 205 : ILEA1.LEAC(1,1,2)= 0 206 : ILEA1.LEAC(1,2,1)= 0 207 : ILEA1.LEAC(1,2,2)= 0 208 : ELSEIF (IDIM.EQ.2) THEN 209 : C Creation des elements TRI3 en definissant les noeuds 210 : NBNN=3 211 : NBELEM=IPT1.NUM(/2)*3 212 : SEGINI,IPT2 213 : IPT2.ITYPEL=4 214 : NBELE2=2 215 : IPT2.NUM(1,1)=NNOE1 216 : IPT2.NUM(2,1)=NNOE1+1 217 : IPT2.NUM(3,1)=NNOE1+3 218 : IPT2.NUM(1,2)=NNOE1+1 219 : IPT2.NUM(2,2)=NNOE1+2 220 : IPT2.NUM(3,2)=NNOE1+3 221 : C Liste des numeros des noeuds de la boite 222 : NNBOIT=4 223 : LNBOIT(1)=NNOE1 224 : LNBOIT(2)=NNOE1+1 225 : LNBOIT(3)=NNOE1+2 226 : LNBOIT(4)=NNOE1+3 227 : C Initialisation du segment ILEA1 contenant la table de 228 : C connectivite des elements de IPT2 229 : NBL1=NBELEM 230 : SEGINI,ILEA1 231 : ILEA1.LEAC(1,1,1)= 2 232 : ILEA1.LEAC(1,1,2)= 2 233 : ILEA1.LEAC(1,2,1)= 0 234 : ILEA1.LEAC(1,2,2)= 0 235 : ILEA1.LEAC(1,3,1)= 0 236 : ILEA1.LEAC(1,3,2)= 0 237 : ILEA1.LEAC(2,1,1)= 0 238 : ILEA1.LEAC(2,1,2)= 0 239 : ILEA1.LEAC(2,2,1)= 1 240 : ILEA1.LEAC(2,2,2)= 1 241 : ILEA1.LEAC(2,3,1)= 0 242 : ILEA1.LEAC(2,3,2)= 0 243 : ELSEIF (IDIM.EQ.3) THEN 244 : C Creation des elements TET4 en definissant les noeuds 245 : NBNN=4 246 : NBELEM=IPT1.NUM(/2)*7 247 : SEGINI,IPT2 248 : IPT2.ITYPEL=23 249 : NBELE2=5 250 : IPT2.NUM(1,1)=NNOE1 251 : IPT2.NUM(2,1)=NNOE1+1 252 : IPT2.NUM(3,1)=NNOE1+3 253 : IPT2.NUM(4,1)=NNOE1+4 254 : IPT2.NUM(1,2)=NNOE1+1 255 : IPT2.NUM(2,2)=NNOE1+2 256 : IPT2.NUM(3,2)=NNOE1+3 257 : IPT2.NUM(4,2)=NNOE1+6 258 : IPT2.NUM(1,3)=NNOE1+1 259 : IPT2.NUM(2,3)=NNOE1+4 260 : IPT2.NUM(3,3)=NNOE1+5 261 : IPT2.NUM(4,3)=NNOE1+6 262 : IPT2.NUM(1,4)=NNOE1+3 263 : IPT2.NUM(2,4)=NNOE1+4 264 : IPT2.NUM(3,4)=NNOE1+6 265 : IPT2.NUM(4,4)=NNOE1+7 266 : IPT2.NUM(1,5)=NNOE1+1 267 : IPT2.NUM(2,5)=NNOE1+3 268 : IPT2.NUM(3,5)=NNOE1+4 269 : IPT2.NUM(4,5)=NNOE1+6 270 : C Liste des numeros des noeuds de la boite 271 : NNBOIT=8 272 : LNBOIT(1)=NNOE1 273 : LNBOIT(2)=NNOE1+1 274 : LNBOIT(3)=NNOE1+2 275 : LNBOIT(4)=NNOE1+3 276 : LNBOIT(5)=NNOE1+4 277 : LNBOIT(6)=NNOE1+5 278 : LNBOIT(7)=NNOE1+6 279 : LNBOIT(8)=NNOE1+7 280 : C Initialisation du segment ILEA1 contenant la table de 281 : C connectivite des elements de IPT2 282 : NBL1=NBELEM 283 : SEGINI,ILEA1 284 : ILEA1.LEAC(1,1,1)= 5 285 : ILEA1.LEAC(1,1,2)= 4 286 : ILEA1.LEAC(1,2,1)= 0 287 : ILEA1.LEAC(1,2,2)= 0 288 : ILEA1.LEAC(1,3,1)= 0 289 : ILEA1.LEAC(1,3,2)= 0 290 : ILEA1.LEAC(1,4,1)= 0 291 : ILEA1.LEAC(1,4,2)= 0 292 : ILEA1.LEAC(2,1,1)= 0 293 : ILEA1.LEAC(2,1,2)= 0 294 : ILEA1.LEAC(2,2,1)= 5 295 : ILEA1.LEAC(2,2,2)= 3 296 : ILEA1.LEAC(2,3,1)= 0 297 : ILEA1.LEAC(2,3,2)= 0 298 : ILEA1.LEAC(2,4,1)= 0 299 : ILEA1.LEAC(2,4,2)= 0 300 : ILEA1.LEAC(3,1,1)= 0 301 : ILEA1.LEAC(3,1,2)= 0 302 : ILEA1.LEAC(3,2,1)= 0 303 : ILEA1.LEAC(3,2,2)= 0 304 : ILEA1.LEAC(3,3,1)= 5 305 : ILEA1.LEAC(3,3,2)= 2 306 : ILEA1.LEAC(3,4,1)= 0 307 : ILEA1.LEAC(3,4,2)= 0 308 : ILEA1.LEAC(4,1,1)= 0 309 : ILEA1.LEAC(4,1,2)= 0 310 : ILEA1.LEAC(4,2,1)= 0 311 : ILEA1.LEAC(4,2,2)= 0 312 : ILEA1.LEAC(4,3,1)= 0 313 : ILEA1.LEAC(4,3,2)= 0 314 : ILEA1.LEAC(4,4,1)= 5 315 : ILEA1.LEAC(4,4,2)= 1 316 : ILEA1.LEAC(5,1,1)= 4 317 : ILEA1.LEAC(5,1,2)= 4 318 : ILEA1.LEAC(5,2,1)= 3 319 : ILEA1.LEAC(5,2,2)= 3 320 : ILEA1.LEAC(5,3,1)= 2 321 : ILEA1.LEAC(5,3,2)= 2 322 : ILEA1.LEAC(5,4,1)= 1 323 : ILEA1.LEAC(5,4,2)= 1 324 : ENDIF 325 : C 326 : C---- 1.3 Initialisation du segment contenant les rayons et centres 327 : C de la triangulation : 328 : IF (IDIM.NE.1) THEN 329 : NBE1=IPT2.NUM(/2) 330 : SEGINI,ITRC1 331 : NBNER = LNBOIT(NNBOIT) 332 : SEGINI,IRN1 333 : C 334 : if (MPOVAL .ne. 0) then 335 : DO I=1,IPT1.NUM(/2) 336 : IRN1.RAYN(IPT1.NUM(1,I))= MPOVAL.VPOCHA(I,1) 337 : END DO 338 : endif 339 : C 340 : DO I0=1,NBELE2 341 : C BCIRCO renvoie le centre et le rayon de l element considere 342 : C 342 : CALL BCIRCO(IRN1,IPT2,I0,XCBJ,XRBJ) 343 : C 344 : ITRC1.TRC(I0,1)=XRBJ 345 : DO J0=1,IDIM 346 : ITRC1.TRC(I0,J0+1)=XCBJ(J0) 347 : ENDDO 348 : ENDDO 349 : ENDIF 350 : C 351 : C 352 : C 353 : C 354 : C----------------------------------------------------------------------C 355 : C PARTIE 2 C 356 : C Triangulation incrementale C 357 : C----------------------------------------------------------------------C 358 : C 359 : JG=1 360 : SEGINI,LSUP 361 : C On boucle sur tous les points a trianguler 362 : NBSOUS=0 363 : NBREF=0 364 : NBNN=IPT2.NUM(/1) 365 : DO I=1,IPT1.NUM(/2) 366 : C WRITE(6,*) 367 : C WRITE(6,*) 'Iteration I=',I,' /',IPT1.NUM(/2) 368 : C 369 : C------- 2.1 Construction de la base 370 : C 371 : C Je tire un point a ajouter a la triangulation 372 : NPI1=IPT1.NUM(1,I) 373 : C XPI1 : vecteur contenant les coordonnees du point 374 : DO J=1,IDIM 375 : XPI1(J)=XCOOR((NPI1-1)*IDIMP1+J) 376 : ENDDO 377 : C 378 : C 2.1.1 Recherche d'un element de la base (NBASE1) 379 : C Je demarre avec le dernier elem. de IPT2 (NEL0) et cherche si 380 : C NPI1 est dedans. Si oui, c'est gagne, sinon : 381 : C je determine le voisin de NEL0 oriente vers NPI1 382 : C => nouveau NEL0, puis je recommence jusqu'a y arriver 383 : NEL0=NBELE2 384 : NBASE=1 385 : C Cas particulier de la dimension 1 386 : IF (IDIM.EQ.1) THEN 387 : DO J=1,NBELE2 388 : C Coordonnees des noeuds A et B de l'element NEL0 389 : C A est a gauche et B est a droite 390 : NREFA=IPT2.NUM(1,NEL0) 391 : NREFB=IPT2.NUM(2,NEL0) 392 : XA=XCOOR((NREFA-1)*IDIMP1+1) 393 : XB=XCOOR((NREFB-1)*IDIMP1+1) 394 : C Cas ou NPI1 est confondu avec A ou B 395 : X1=ABS(XPI1(1)) 396 : IF (X1.LT.ZERO) X1=1. 397 : DTEST1=ABS(XPI1(1)-XA)/X1 398 : DTEST2=ABS(XPI1(1)-XB)/X1 399 : IF ((DTEST1.LT.ZERO).OR.(DTEST2.LT.ZERO)) THEN 400 : NBNN=IPT2.NUM(/1) 401 : NBELEM=NBELE2 402 : NBSOUS=0 403 : NBREF=0 404 : SEGADJ,IPT2 405 : C Erreur de triangulation : vérifiez qu'il n'y a pas de 406 : C points confondus 406 : CALL ERREUR(846) 407 : GOTO 2501 408 : C Cas ou le NPT1 est entre A et B, on a trouve la base 409 : ELSEIF ((XPI1(1).GT.XA).AND.(XPI1(1).LT.XB)) THEN 410 : NBASE1=NEL0 411 : GOTO 2381 412 : C Si NPI1 est a gauche de A, on va chercher le voisin 413 : ELSEIF (XPI1(1).LT.XA) THEN 414 : NEL0=ILEA1.LEAC(NEL0,2,1) 415 : C Si NPI1 est a droite de B, on va chercher le voisin 416 : ELSEIF (XPI1(1).GT.XB) THEN 417 : NEL0=ILEA1.LEAC(NEL0,1,1) 418 : ENDIF 419 : ENDDO 420 : C Si l'on passe ici, c'est que l'on a pas trouve NBASE1 ! 421 : NBNN=IPT2.NUM(/1) 422 : NBELEM=NBELE2 423 : NBSOUS=0 424 : NBREF=0 425 : SEGADJ,IPT2 426 : C Erreur de triangulation : vérifiez qu'il n'y a pas de 427 : C points confondus 427 : CALL ERREUR(846) 428 : GOTO 2501 429 : ENDIF 430 : C Cas generaux en dimension 2 et 3 431 : JG=1 432 : SEGINI,LSCP 433 : NEL00=NEL0 434 : NEL000=NEL0 435 : DO J=1,NBELE2 436 : NCP=0 437 : C Coordonnees des noeuds de l'element NEL0 438 : DO K=1,IDIMP1 439 : NREF0=IPT2.NUM(K,NEL0) 440 : DO KK=1,IDIM 441 : X0(K,KK)=XCOOR((NREF0-1)*IDIMP1+KK) 442 : ENDDO 443 : ENDDO 444 : C Coordonnees du barycentre G de NEL0 445 : DO KK=1,IDIM 446 : XSOM=0. 447 : DO K=1,IDIMP1 448 : XSOM=XSOM+X0(K,KK) 449 : ENDDO 450 : XG0(KK)=XSOM/(IDIMP1) 451 : ENDDO 452 : IK1=0 453 : IK2=0 454 : IK3=0 455 : C En dimension 2 456 : IF (IDIM.EQ.2) THEN 457 : C On effectue le meme test sur chaque cote 458 : DO K=1,3 459 : IF (K.EQ.1) THEN 460 : NA=2 461 : NB=3 462 : ELSEIF (K.EQ.2) THEN 463 : NA=1 464 : NB=3 465 : ELSEIF (K.EQ.3) THEN 466 : NA=1 467 : NB=2 468 : ENDIF 469 : C Normale au cote AB du triangle 470 : XN=X0(NA,2)-X0(NB,2) 471 : YN=X0(NB,1)-X0(NA,1) 472 : XXN=SQRT((XN**2)+(YN**2)) 473 : XN=XN/XXN 474 : YN=YN/XXN 475 : C On teste si NPI1 est du meme cote que G par rapport a 476 : C la droite (AB) 477 : PS1=(XN*(XPI1(1)-X0(NA,1)))+(YN*(XPI1(2)-X0(NA,2))) 478 : XX0=SQRT(((XPI1(1)-X0(NA,1))**2)+ 479 : & ((XPI1(2)-X0(NA,2))**2)) 480 : PS11=PS1/XX0 481 : IF (ABS(PS11).LT.ZERO) THEN 482 : NSCP=K 483 : NCP=1 484 : NBASE=2 485 : GOTO 2111 486 : ENDIF 487 : PS2=(XN*(XG0(1)-X0(NA,1)))+(YN*(XG0(2)-X0(NA,2))) 488 : IF ((SIGN(1D0,PS1)).EQ.(SIGN(1D0,PS2))) THEN 489 : GOTO 2111 490 : ELSE 491 : NADJ0=K 492 : GOTO 2112 493 : ENDIF 494 : 2111 CONTINUE 495 : IF (K.EQ.1) IK1=1 496 : IF (K.EQ.2) IK2=1 497 : IF (K.EQ.3) IK3=1 498 : ENDDO 499 : IF ((IK1.EQ.1).AND.(IK2.EQ.1).AND.(IK3.EQ.1)) THEN 500 : C NPI1 est dans NEL0, on a trouve un element de la base 501 : NBASE1=NEL0 502 : GOTO 2115 503 : ENDIF 504 : 2112 CONTINUE 505 : C NPI1 et G ne sont pas du meme cote par rapport au cote 506 : C vu par le noeud NADJ0 de NEL0 => on passe a l'element 507 : C voisin NEL1 508 : IF (J.GE.3) NEL000=NEL00 509 : IF (J.GE.2) NEL00=NEL0 510 : NEL0=ILEA1.LEAC(NEL0,NADJ0,1) 511 : IF (NEL0.EQ.NEL00) THEN 512 : C PRINT*,'Erreur dans la table des elements adjacents' 513 : C PRINT*,'L''element',NEL0,'est voisin de lui meme !' 513 : CALL ERREUR(846) 514 : GOTO 2501 515 : ENDIF 516 : IF (NEL0.EQ.NEL000) THEN 517 : C PRINT*,'Erreur dans la construction du maillage' 518 : C PRINT*,'Les elements',NEL0,NEL00,'se chevauchent !' 518 : CALL ERREUR(846) 519 : GOTO 2501 520 : ENDIF 521 : ENDIF 522 : C En dimension 3 523 : IF (IDIM.EQ.3) THEN 524 : IK4=0 525 : C On effectue le test sur chaque face 526 : DO K=1,4 527 : IF (K.EQ.1) THEN 528 : NA=2 529 : NB=3 530 : NC=4 531 : ELSEIF (K.EQ.2) THEN 532 : NA=1 533 : NB=3 534 : NC=4 535 : ELSEIF (K.EQ.3) THEN 536 : NA=1 537 : NB=2 538 : NC=4 539 : ELSEIF (K.EQ.4) THEN 540 : NA=1 541 : NB=2 542 : NC=3 543 : ENDIF 544 : C Normale au plan ABC du triangle (n = AB^AC) 545 : XN=((X0(NB,2)-X0(NA,2))*(X0(NC,3)-X0(NA,3))) - 546 : & ((X0(NB,3)-X0(NA,3))*(X0(NC,2)-X0(NA,2))) 547 : YN=((X0(NB,3)-X0(NA,3))*(X0(NC,1)-X0(NA,1))) - 548 : & ((X0(NB,1)-X0(NA,1))*(X0(NC,3)-X0(NA,3))) 549 : ZN=((X0(NB,1)-X0(NA,1))*(X0(NC,2)-X0(NA,2))) - 550 : & ((X0(NB,2)-X0(NA,2))*(X0(NC,1)-X0(NA,1))) 551 : XXN=SQRT((XN**2)+(YN**2)+(ZN**2)) 552 : XN=XN/XXN 553 : YN=YN/XXN 554 : ZN=ZN/XXN 555 : C On teste si NPI1 est du meme cote que G par rapport au 556 : C plan ABC 557 : PS1=(XN*(XPI1(1)-X0(NA,1))) + 558 : & (YN*(XPI1(2)-X0(NA,2))) + 559 : & (ZN*(XPI1(3)-X0(NA,3))) 560 : XX0=SQRT(((XPI1(1)-X0(NA,1))**2) + 561 : & ((XPI1(2)-X0(NA,2))**2)+ 562 : & ((XPI1(3)-X0(NA,3))**2)) 563 : PS11=PS1/XX0 564 : IF (ABS(PS11).LT.ZERO) THEN 565 : NCP=NCP+1 566 : JG=NCP 567 : SEGADJ,LSCP 568 : LSCP.LECT(NCP)=K 569 : GOTO 2113 570 : ENDIF 571 : PS2=(XN*(XG0(1)-X0(NA,1))) + 572 : & (YN*(XG0(2)-X0(NA,2))) + 573 : & (ZN*(XG0(3)-X0(NA,3))) 574 : IF ((SIGN(1D0,PS1)).EQ.(SIGN(1D0,PS2))) THEN 575 : GOTO 2113 576 : ELSE 577 : GOTO 2114 578 : ENDIF 579 : 2113 CONTINUE 580 : IF (K.EQ.1) IK1=1 581 : IF (K.EQ.2) IK2=1 582 : IF (K.EQ.3) IK3=1 583 : IF (K.EQ.4) IK4=1 584 : ENDDO 585 : IF ((IK1.EQ.1).AND.(IK2.EQ.1).AND. 586 : & (IK3.EQ.1).AND.(IK4.EQ.1)) THEN 587 : C NPI1 est bien dans NEL0, on a donc trouve la base 588 : NBASE1=NEL0 589 : GOTO 2115 590 : ENDIF 591 : 2114 CONTINUE 592 : C NPI1 et G ne sont pas du meme cote par rapport au plan 593 : C vu par le noeud NADJ0 de NEL0 => on passe a l'element 594 : C adjacent a NEL0 vu par le noeud K 595 : IF (J.GE.3) NEL000=NEL00 596 : IF (J.GE.2) NEL00=NEL0 597 : NEL0=ILEA1.LEAC(NEL0,K,1) 598 : IF (NEL0.EQ.NEL00) THEN 599 : C PRINT*,'Erreur dans la table des elements adjacents' 600 : C PRINT*,'L''element',NEL0,'est voisin de lui meme !' 600 : CALL ERREUR(846) 601 : GOTO 2501 602 : ENDIF 603 : IF (NEL0.EQ.NEL000) THEN 604 : C PRINT*,'Erreur dans la construction du maillage' 605 : C PRINT*,'Les elements',NEL0,NEL00,'se chevauchent !' 605 : CALL ERREUR(846) 606 : GOTO 2501 607 : ENDIF 608 : ENDIF 609 : ENDDO 610 : C Si l'on passe ici, c'est que l'on a pas trouve NBASE1 ! 611 : NBNN=IPT2.NUM(/1) 612 : NBELEM=NBELE2 613 : NBSOUS=0 614 : NBREF=0 615 : SEGADJ,IPT2 616 : C Erreur de triangulation : vérifiez qu'il n'y a pas de points 617 : C confondus 617 : CALL ERREUR(846) 618 : GOTO 2501 619 : 2115 CONTINUE 620 : C Test si NPI1 n'est pas confondu avec un des noeuds de NBASE1 621 : DO J=1,IDIMP1 622 : NREFJ=IPT2.NUM(J,NBASE1) 623 : NC=0 624 : DO JJ=1,IDIM 625 : X1=ABS(XPI1(JJ)) 626 : IF (X1.LT.ZERO) X1=1. 627 : DTEST=ABS(XPI1(JJ)-XCOOR((NREFJ-1)*IDIMP1+JJ))/X1 628 : IF (DTEST.LT.ZERO) NC=NC+1 629 : ENDDO 630 : IF (NC.EQ.IDIM) THEN 631 : NBNN=IPT2.NUM(/1) 632 : NBELEM=NBELE2 633 : NBSOUS=0 634 : NBREF=0 635 : SEGADJ,IPT2 636 : C Erreur de triangulation : vérifiez qu'il n'y a pas de 637 : C points confondus 637 : CALL ERREUR(846) 638 : GOTO 2501 639 : ENDIF 640 : ENDDO 641 : C 642 : C 2.1.2 Ajout eventuel d'autres elements a la base 643 : C LBASE : liste des numeros des elements de IPT2 appartenant a 644 : C la base, taille initialisee a NBASE 645 : JG=NBASE 646 : SEGINI,LBASE 647 : LBASE.LECT(1)=NBASE1 648 : C Ajout des elments supplementaires a LBASE si necessaire 649 : IF (NCP.NE.0) THEN 650 : C En dimension 2, un seul adjacent est ajoute a LBASE, celui 651 : C partageant le cote 652 : IF (IDIM.EQ.2) THEN 653 : LBASE.LECT(2)=ILEA1.LEAC(NBASE1,NSCP,1) 654 : C En dimension 3 plusieurs cas sont possibles 655 : ELSEIF (IDIM.EQ.3) THEN 656 : C Cas ou NPI1 est sur un des 4 plans 657 : IF (NCP.EQ.1) THEN 658 : C On recupere l'adjacent partageant le plan 659 : NBASE=2 660 : JG=NBASE 661 : SEGADJ,LBASE 662 : LBASE.LECT(2)=ILEA1.LEAC(NBASE1,LSCP.LECT(1),1) 663 : C Cas ou NPI1 est sur une des aretes 664 : ELSEIF (NCP.EQ.2) THEN 665 : C On recupere tous les elements partageant l'arete 666 : NA=1 667 : NB=2 668 : NC=LSCP.LECT(1) 669 : ND=LSCP.LECT(2) 670 : NREFC=IPT2.NUM(NC,NBASE1) 671 : NREFD=IPT2.NUM(ND,NBASE1) 672 : DO L=1,4 673 : IF ((NA.EQ.NC).OR.(NA.EQ.ND).OR.(NA.EQ.NB)) NA=NA+1 674 : IF ((NB.EQ.NC).OR.(NB.EQ.ND).OR.(NB.EQ.NA)) NB=NB+1 675 : ENDDO 676 : NREFA=IPT2.NUM(NA,NBASE1) 677 : NREFB=IPT2.NUM(NB,NBASE1) 678 : N1=NBASE1 679 : DO L=1,100 680 : N2=ILEA1.LEAC(N1,NC,1) 681 : ND=ILEA1.LEAC(N1,NC,2) 682 : IF (N2.EQ.NBASE1) GOTO 2121 683 : NBASE=NBASE+1 684 : JG=NBASE 685 : SEGADJ,LBASE 686 : LBASE.LECT(NBASE)=N2 687 : DO M=1,4 688 : IF (IPT2.NUM(M,N2).EQ.NREFA) NA=M 689 : IF (IPT2.NUM(M,N2).EQ.NREFB) NB=M 690 : IF (IPT2.NUM(M,N2).EQ.NREFD) NC=M 691 : ENDDO 692 : NREFD=IPT2.NUM(ND,N2) 693 : N1=N2 694 : ENDDO 695 : 2121 CONTINUE 696 : C Cas ou NPI1 est confondu avec un des sommets 697 : ELSEIF (NCP.EQ.3) THEN 698 : NBNN=IPT2.NUM(/1) 699 : NBELEM=NBELE2 700 : NBSOUS=0 701 : NBREF=0 702 : SEGADJ,IPT2 703 : C Erreur de triangulation : vérifiez qu'il n'y a pas de 704 : C points confondus 704 : CALL ERREUR(846) 705 : GOTO 2501 706 : ENDIF 707 : ENDIF 708 : ENDIF 709 : C 710 : C------- 2.2 Construction de la cavite 711 : C 712 : C 2.2.1 Construction de la cavite (MELEME IPT3) 713 : C par adjacence a partir de l'element NBASE1 de la base 714 : C 715 : C LCAVI : liste des numeros des elements de IPT2 appartenant a 716 : C la cavite, taille initialisee a NBELE2 717 : JG=NBELE2 718 : SEGINI,LCAVI 719 : C Le premier element de la cavite est NBASE1 720 : LCAVI.LECT(1)=NBASE1 721 : C NCAVI : nombre d'elements dans la cavite (toujours a jour !) 722 : NCAVI=1 723 : C Initialisation de la cavite IPT3 a NBELE2 724 : NBELEM=NBELE2 725 : NBNN=IPT2.NUM(/1) 726 : SEGINI,IPT3 727 : IPT3.ITYPEL=IPT2.ITYPEL 728 : DO K=1,IPT3.NUM(/1) 729 : IPT3.NUM(K,NCAVI)=IPT2.NUM(K,NBASE1) 730 : ENDDO 731 : C Parcours des adjacents de maniere recursive jusqu'a stabilite 732 : ITRA=1 733 : DO J=1,NBELE2 734 : C Element a traiter de la liste LCAVI 735 : NEJ=LCAVI.LECT(ITRA) 736 : IF ((NEJ.EQ.0).OR.(J.EQ.NBELE2)) GOTO 2213 737 : C On prend les elements adjacents de NEJ 738 : DO K=1,IDIMP1 739 : NEAJ=ILEA1.LEAC(NEJ,K,1) 740 : C On verifie qu'il y a bien un element adjacent 741 : IF (NEAJ.EQ.0) GOTO 2212 742 : C On verifie que NEAJ n'est pas deja dans LCAVI 743 : DO L=1,LCAVI.LECT(/1) 744 : IF (LCAVI.LECT(L).EQ.0) GOTO 2211 745 : IF (LCAVI.LECT(L).EQ.NEAJ) GOTO 2212 746 : ENDDO 747 : 2211 CONTINUE 748 : C XRBJ : rayon de la boule circonscrite a l'elem. NEAJ 749 : C XCBJ(i) : ieme coordonnee de son centre 750 : XRBJ=ITRC1.TRC(NEAJ,1) 751 : DO L=1,IDIM 752 : XCBJ(L)=ITRC1.TRC(NEAJ,L+1) 753 : ENDDO 754 : C DIJ1 : distance de NPI1 au centre de la boule de 755 : C l'element considere 756 : DIJ1=0. 757 : DO L=1,IDIM 758 : DIJ1=DIJ1+((XPI1(L)-XCBJ(L))**2) 759 : ENDDO 760 : C 761 : IF (MPOVAL .NE. 0) 762 : & DIJ1 = DIJ1 - ((IRN1.RAYN(IPT1.NUM(1,I)))**2) 763 : C 764 : DIJ1=(DIJ1/(XRBJ**2))-1. 765 : C Si NEAJ est un nouvel element de la cavite, 766 : IF (DIJ1.LT.ZERO) THEN 767 : C alors on l'ajoute a LCAVI et on incremente NCAVI 768 : NCAVI=NCAVI+1 769 : LCAVI.LECT(NCAVI)=NEAJ 770 : C et on remplis IPT3 avec ce nouvel element 771 : DO L=1,IPT3.NUM(/1) 772 : IPT3.NUM(L,NCAVI)=IPT2.NUM(L,NEAJ) 773 : ENDDO 774 : ENDIF 775 : 2212 CONTINUE 776 : ENDDO 777 : ITRA=ITRA+1 778 : ENDDO 779 : 2213 CONTINUE 780 : C Ajustement de LCAVI et IPT3 a la bonne taille si besoin 781 : IF (LCAVI.LECT(/1).GT.NCAVI) THEN 782 : JG=NCAVI 783 : SEGADJ,LCAVI 784 : NBNN=IPT3.NUM(/1) 785 : NBELEM=NCAVI 786 : SEGADJ,IPT3 787 : ENDIF 788 : NBELE3=IPT3.NUM(/2) 789 : C 790 : C 2.2.2 Verification : la cavite doit etre connexe, chaque 791 : C element doit avoir au moins un element adjacent, c'est a dire 792 : C partageant un cote/face. De plus, un cote/face ne peut etre 793 : C commun qu'a deux elements au plus. 794 : C On determine egalement le contour/enveloppe, IPT4, de la cavite 795 : NNBELE2=NBELE2 796 : 2221 CONTINUE 797 : NBELE2=NNBELE2 798 : NNCON=0 799 : JG=IDIM 800 : SEGINI,LREF1 801 : NBNN=IPT3.NUM(/1)-1 802 : NBELEM=NBELE3*3 803 : SEGINI,IPT4 804 : IF (IDIM.EQ.2) IPT4.ITYPEL=2 805 : IF (IDIM.EQ.3) IPT4.ITYPEL=4 806 : NBELE4=0 807 : C On parcoure les elements de la cavite 808 : DO I1=1,IPT3.NUM(/2) 809 : C Nombre d'adjacents a l'element I1 810 : C un element est adjacent s'il partage un cote/face 811 : NVI1=0 812 : C On parcoure les 3/4 cotes/faces de l'element I1 813 : DO I2=1,IPT3.NUM(/1) 814 : C Nombre d'adjacents a la face I2 (1 ou 0) 815 : NVI2=0 816 : IF (IDIM.EQ.2) THEN 817 : IF (I2.EQ.1) THEN 818 : LREF1.LECT(1)=IPT3.NUM(2,I1) 819 : LREF1.LECT(2)=IPT3.NUM(3,I1) 820 : ELSEIF (I2.EQ.2) THEN 821 : LREF1.LECT(1)=IPT3.NUM(1,I1) 822 : LREF1.LECT(2)=IPT3.NUM(3,I1) 823 : ELSEIF (I2.EQ.3) THEN 824 : LREF1.LECT(1)=IPT3.NUM(1,I1) 825 : LREF1.LECT(2)=IPT3.NUM(2,I1) 826 : ENDIF 827 : ELSEIF (IDIM.EQ.3) THEN 828 : IF (I2.EQ.1) THEN 829 : LREF1.LECT(1)=IPT3.NUM(2,I1) 830 : LREF1.LECT(2)=IPT3.NUM(3,I1) 831 : LREF1.LECT(3)=IPT3.NUM(4,I1) 832 : ELSEIF (I2.EQ.2) THEN 833 : LREF1.LECT(1)=IPT3.NUM(1,I1) 834 : LREF1.LECT(2)=IPT3.NUM(3,I1) 835 : LREF1.LECT(3)=IPT3.NUM(4,I1) 836 : ELSEIF (I2.EQ.3) THEN 837 : LREF1.LECT(1)=IPT3.NUM(1,I1) 838 : LREF1.LECT(2)=IPT3.NUM(2,I1) 839 : LREF1.LECT(3)=IPT3.NUM(4,I1) 840 : ELSEIF (I2.EQ.4) THEN 841 : LREF1.LECT(1)=IPT3.NUM(1,I1) 842 : LREF1.LECT(2)=IPT3.NUM(2,I1) 843 : LREF1.LECT(3)=IPT3.NUM(3,I1) 844 : ENDIF 845 : ENDIF 846 : DO I11=1,IPT3.NUM(/2) 847 : NNCI11=0 848 : IF (I11.EQ.I1) GOTO 2222 849 : DO I22=1,IPT3.NUM(/1) 850 : NREFX=IPT3.NUM(I22,I11) 851 : DO K=1,IDIM 852 : IF (NREFX.EQ.LREF1.LECT(K)) THEN 853 : NNCI11=NNCI11+1 854 : ENDIF 855 : ENDDO 856 : ENDDO 857 : IF (NNCI11.EQ.IDIM) THEN 858 : NVI2=NVI2+1 859 : NVI1=NVI1+1 860 : ENDIF 861 : 2222 CONTINUE 862 : ENDDO 863 : C Cas ou I1 a plus de 1 voisin sur la face I2 864 : IF (NVI2.GT.1) THEN 865 : NBNN=IPT2.NUM(/1) 866 : NBELEM=NBELE2 867 : NBSOUS=0 868 : NBREF=0 869 : SEGADJ,IPT2 870 : C Erreur de triangulation : vérifiez qu'il n'y a pas de 871 : C points confondus 871 : CALL ERREUR(846) 872 : GOTO 2501 873 : ENDIF 874 : C Si la face I2 n'est que sur l'element I1, alors on 875 : C l'ajoute au maillage du contour/enceloppe de la cavite 876 : IF (NVI2.EQ.0) THEN 877 : NBELE4=NBELE4+1 878 : IF (IPT4.NUM(/2).LT.NBELE4) THEN 879 : NBNN=IPT4.NUM(/1) 880 : NBELEM=IPT4.NUM(/2)+100 881 : SEGADJ,IPT4 882 : ENDIF 883 : DO K=1,IPT4.NUM(/1) 884 : IF (IDIM.EQ.2) THEN 885 : IF (I2.EQ.1) THEN 886 : IPT4.NUM(1,NBELE4)=IPT3.NUM(2,I1) 887 : IPT4.NUM(2,NBELE4)=IPT3.NUM(3,I1) 888 : ELSEIF (I2.EQ.2) THEN 889 : IPT4.NUM(1,NBELE4)=IPT3.NUM(1,I1) 890 : IPT4.NUM(2,NBELE4)=IPT3.NUM(3,I1) 891 : ELSEIF (I2.EQ.3) THEN 892 : IPT4.NUM(1,NBELE4)=IPT3.NUM(1,I1) 893 : IPT4.NUM(2,NBELE4)=IPT3.NUM(2,I1) 894 : ENDIF 895 : ELSEIF (IDIM.EQ.3) THEN 896 : IF (I2.EQ.1) THEN 897 : IPT4.NUM(1,NBELE4)=IPT3.NUM(2,I1) 898 : IPT4.NUM(2,NBELE4)=IPT3.NUM(3,I1) 899 : IPT4.NUM(3,NBELE4)=IPT3.NUM(4,I1) 900 : ELSEIF (I2.EQ.2) THEN 901 : IPT4.NUM(1,NBELE4)=IPT3.NUM(1,I1) 902 : IPT4.NUM(2,NBELE4)=IPT3.NUM(3,I1) 903 : IPT4.NUM(3,NBELE4)=IPT3.NUM(4,I1) 904 : ELSEIF (I2.EQ.3) THEN 905 : IPT4.NUM(1,NBELE4)=IPT3.NUM(1,I1) 906 : IPT4.NUM(2,NBELE4)=IPT3.NUM(2,I1) 907 : IPT4.NUM(3,NBELE4)=IPT3.NUM(4,I1) 908 : ELSEIF (I2.EQ.4) THEN 909 : IPT4.NUM(1,NBELE4)=IPT3.NUM(1,I1) 910 : IPT4.NUM(2,NBELE4)=IPT3.NUM(2,I1) 911 : IPT4.NUM(3,NBELE4)=IPT3.NUM(3,I1) 912 : ENDIF 913 : ENDIF 914 : ENDDO 915 : ENDIF 916 : ENDDO 917 : C Cas ou I1 n'a aucun voisin => cavite non connexe, on retire 918 : C l'element de la cavite 919 : IF (NVI1.EQ.0) THEN 920 : IF (NBELE3.NE.1) THEN 921 : NNCON=1 922 : NESUP=LCAVI.LECT(I1) 923 : ENDIF 924 : ENDIF 925 : ENDDO 926 : NBNN=IPT4.NUM(/1) 927 : NBELEM=NBELE4 928 : SEGADJ,IPT4 929 : IF (NNCON.EQ.1) GOTO 2341 930 : C 931 : C 2.2.3 Verification : la cavite et son contour/enveloppe doivent 932 : C avoir le meme nombre de noeuds 932 : CALL ECROBJ('MAILLAGE',IPT3) 933 : IF (IERR.NE.0) RETURN 933 : CALL NBNO 934 : SEGACT,IPT3 934 : CALL LIRENT(NBNO3,1,IRETOU) 934 : CALL ECROBJ('MAILLAGE',IPT4) 935 : IF (IERR.NE.0) RETURN 935 : CALL NBNO 936 : SEGACT,IPT4 936 : CALL LIRENT(NBNO4,1,IRETOU) 937 : IF (NBNO3.NE.NBNO4) THEN 938 : NBNN=IPT2.NUM(/1) 939 : NBELEM=NBELE2 940 : NBSOUS=0 941 : NBREF=0 942 : SEGADJ,IPT2 943 : C Erreur de triangulation : vérifiez qu'il n'y a pas de 944 : C points confondus 944 : CALL ERREUR(846) 945 : GOTO 2501 946 : ENDIF 947 : C 948 : C------- 2.3 Construction de la boule 949 : C 950 : C 2.3.1 Construction de la liste des elements de la boule LBOUL 951 : C La boule est formee par les elements construits en joignant 952 : C NPI1 aux elements du contour de la cavite. 953 : NBOUL=NBELE4 954 : JG=NBOUL 955 : SEGINI,LBOUL 956 : C Nouveau nombre d'elements de la triangulation 957 : NEW2=NBELE2-NCAVI+NBOUL 958 : C On regarde si le nombre d'elements dans la cavite est inferieur 959 : C au nombre d'elements de son contour/enveloppe car il faudra un 960 : C traitement special lors de la creation de la boule 961 : N43=NBELE4-NBELE3 962 : C Cas general ou l'on va ajouter de nouveaux elements a IPT2 963 : IF (N43.GE.0) THEN 964 : DO J=1,NBOUL 965 : C Les premiers elements de LBOUL sont ceux de LCAVI 966 : IF (J.LE.NCAVI) THEN 967 : LBOUL.LECT(J)=LCAVI.LECT(J) 968 : C et je complete LBOUL avec les nouveaux numeros d'elements 969 : ELSE 970 : LBOUL.LECT(J)=NBELE2+1 971 : NBELE2=NBELE2+1 972 : ENDIF 973 : ENDDO 974 : C J'ajuste IPT2, ITRC1 et ILEA1 si besoin 975 : IF (NEW2.GT.IPT2.NUM(/2)) THEN 976 : NBNN=IPT2.NUM(/1) 977 : NBELEM=NBELE2+ICMIN 978 : SEGADJ,IPT2 979 : NBE1=NBELEM 980 : SEGADJ,ITRC1 981 : NBL1=NBELEM 982 : SEGADJ,ILEA1 983 : ENDIF 984 : C Cas particulier ou l'on va enlever des elements a IPT2 985 : ELSE 986 : NSUP=NCAVI-NBOUL 987 : JG=NSUP 988 : SEGINI,LSUP 989 : C LBOUL est fait des NBOUL premiers elements de LCAVI 990 : DO J=1,NCAVI 991 : NELJ=LCAVI.LECT(J) 992 : IF (J.LE.NBOUL) THEN 993 : LBOUL.LECT(J)=NELJ 994 : ELSE 995 : LSUP.LECT(J-NBOUL)=NELJ 996 : C Mise a zero des numeros des noeuds de l'element NELJ 997 : DO K=1,IPT2.NUM(/1) 998 : IPT2.NUM(K,NELJ)=0 999 : ENDDO 1000 : ENDIF 1001 : ENDDO 1002 : ENDIF 1003 : NBELE2=NEW2 1004 : C 1005 : C 2.3.2 Je recherche les elements exterieurs a la cavite 1006 : C partageant une face de l'enveloppe IPT4, pour cela on regarde 1007 : C les elements de IPT3 et leurs adjacents grace a la table de 1008 : C connectivite 1009 : JG=NBELE4 1010 : SEGINI,LEADJ 1011 : SEGINI,LNVADJ 1012 : JG=IDIM 1013 : SEGINI,LREF1 1014 : C On parcoure les elements II de IPT4 du contour/enveloppe 1015 : DO II=1,NBELE4 1016 : C Numeros des noeuds JJ de l'element II 1017 : DO JJ=1,IDIM 1018 : LREF1.LECT(JJ)=IPT4.NUM(JJ,II) 1019 : ENDDO 1020 : C On parcoure les elements J de IPT3 de la cavite 1021 : DO J=1,NBELE3 1022 : NCOM=0 1023 : NSOM=0 1024 : DO K=1,IDIMP1 1025 : NREFX=IPT3.NUM(K,J) 1026 : C Test si le noeud K de l'element J est commun a II 1027 : DO L=1,IDIM 1028 : IF (NREFX.EQ.LREF1.LECT(L)) THEN 1029 : NCOM=NCOM+1 1030 : NSOM=NSOM+K 1031 : ENDIF 1032 : ENDDO 1033 : C Si l'on a trouve l'element appuye sur la face II 1034 : IF (NCOM.EQ.IDIM) THEN 1035 : NELK2=LCAVI.LECT(J) 1036 : GOTO 2321 1037 : ENDIF 1038 : ENDDO 1039 : ENDDO 1040 : 2321 CONTINUE 1041 : C Determination du numero du noeud voyeur dans IPT2 1042 : IF (IDIM.EQ.2) NNOK2=6-NSOM 1043 : IF (IDIM.EQ.3) NNOK2=10-NSOM 1044 : LEADJ.LECT(II)=ILEA1.LEAC(NELK2,NNOK2,1) 1045 : LNVADJ.LECT(II)=ILEA1.LEAC(NELK2,NNOK2,2) 1046 : ENDDO 1047 : C 1048 : C 2.3.3 Verification : le contour/enveloppe de la cavite doit 1049 : C etre etoile par rapport au point NPI1 1050 : NESUP=0 1051 : IF (IDIM.EQ.3) THEN 1052 : DO J=1,NBELE4 1053 : NA=IPT4.NUM(1,J) 1054 : NB=IPT4.NUM(2,J) 1055 : NC=IPT4.NUM(3,J) 1056 : XA=XCOOR((NA-1)*(IDIM+1)+1) 1057 : YA=XCOOR((NA-1)*(IDIM+1)+2) 1058 : ZA=XCOOR((NA-1)*(IDIM+1)+3) 1059 : XB=XCOOR((NB-1)*(IDIM+1)+1) 1060 : YB=XCOOR((NB-1)*(IDIM+1)+2) 1061 : ZB=XCOOR((NB-1)*(IDIM+1)+3) 1062 : XC=XCOOR((NC-1)*(IDIM+1)+1) 1063 : YC=XCOOR((NC-1)*(IDIM+1)+2) 1064 : ZC=XCOOR((NC-1)*(IDIM+1)+3) 1065 : C Normale au plan ABC : n = AB^AC 1066 : XN=((YB-YA)*(ZC-ZA))-((ZB-ZA)*(YC-YA)) 1067 : YN=((ZB-ZA)*(XC-XA))-((XB-XA)*(ZC-ZA)) 1068 : ZN=((XB-XA)*(YC-YA))-((YB-YA)*(XC-XA)) 1069 : XXN=SQRT((XN**2)+(YN**2)+(ZN**2)) 1070 : XN=XN/XXN 1071 : YN=YN/XXN 1072 : ZN=ZN/XXN 1073 : C Test si NPI1 est dans le plan ABC 1074 : PS1=(XN*(XPI1(1)-XA))+(YN*(XPI1(2)-YA))+(ZN*(XPI1(3)-ZA)) 1075 : XX0=SQRT(((XPI1(1)-XA)**2)+((XPI1(2)-YA)**2)+ 1076 : & ((XPI1(3)-ZA)**2)) 1077 : PS11=PS1/XX0 1078 : IF (ABS(PS11).LT.ZERO) THEN 1079 : NESUP=ILEA1.LEAC(LEADJ.LECT(J),LNVADJ.LECT(J),1) 1080 : GOTO 2341 1081 : ENDIF 1082 : C Point D voyant la face ABC mais situe hors de la cavite 1083 : IF (LEADJ.LECT(J).EQ.0) GOTO 2331 1084 : ND=IPT2.NUM(LNVADJ.LECT(J),LEADJ.LECT(J)) 1085 : XD=XCOOR((ND-1)*(IDIM+1)+1) 1086 : YD=XCOOR((ND-1)*(IDIM+1)+2) 1087 : ZD=XCOOR((ND-1)*(IDIM+1)+3) 1088 : C Test si NPI1 est du meme cote que D par rapport au plan 1089 : C ABC 1090 : PS2=(XN*(XD-XA)) + (YN*(YD-YA)) + (ZN*(ZD-ZA)) 1091 : XX0=SQRT(((XD-XA)**2)+((YD-YA)**2)+((ZD-ZA)**2)) 1092 : PS22=PS2/XXN 1093 : IF (ABS(PS22).LT.ZERO) PS2=0. 1094 : IF ((SIGN(1D0,PS1)).EQ.(SIGN(1D0,PS2))) THEN 1095 : NESUP=ILEA1.LEAC(LEADJ.LECT(J),LNVADJ.LECT(J),1) 1096 : GOTO 2341 1097 : ELSE 1098 : GOTO 2331 1099 : ENDIF 1100 : 2331 CONTINUE 1101 : ENDDO 1102 : ENDIF 1103 : C 1104 : C 2.3.4 Suppression d'un element (NESUP) de la cavite 1105 : 2341 CONTINUE 1106 : IF (NESUP.NE.0) THEN 1107 : C Test si NESUP est dans la base, pas touche a la base ! 1108 : DO K=1,NBASE 1109 : IF (LBASE.LECT(K).EQ.NESUP) THEN 1110 : NBNN=IPT2.NUM(/1) 1111 : NBELEM=NBELE2 1112 : NBSOUS=0 1113 : NBREF=0 1114 : SEGADJ,IPT2 1115 : C Erreur de triangulation : vérifiez qu'il n'y a pas de 1116 : C points confondus 1116 : CALL ERREUR(846) 1117 : GOTO 2501 1118 : ENDIF 1119 : ENDDO 1120 : NBNN=IPT3.NUM(/1) 1121 : NBELEM=NBELE3-1 1122 : SEGINI,IPT5 1123 : IPT5.ITYPEL=IPT3.ITYPEL 1124 : C Suppression de l'element NESUP le la liste LCAVI 1125 : NS=0 1126 : IF (LCAVI.LECT(NCAVI).EQ.NESUP) THEN 1127 : NCAVI=NCAVI-1 1128 : NS=1 1129 : DO JJ=1,NCAVI 1130 : DO KK=1,IPT3.NUM(/1) 1131 : IPT5.NUM(KK,JJ)=IPT3.NUM(KK,JJ) 1132 : ENDDO 1133 : ENDDO 1134 : ELSE 1135 : DO K=1,NCAVI-1 1136 : IF (LCAVI.LECT(K).EQ.NESUP) THEN 1137 : NS=NS+1 1138 : NCAVI=NCAVI-1 1139 : ENDIF 1140 : LCAVI.LECT(K)=LCAVI.LECT(K+NS) 1141 : DO KK=1,IPT3.NUM(/1) 1142 : IPT5.NUM(KK,K)=IPT3.NUM(KK,K+NS) 1143 : ENDDO 1144 : ENDDO 1145 : ENDIF 1146 : C Ajustement de la liste LCAVI 1147 : JG=NCAVI 1148 : SEGADJ,LCAVI 1149 : SEGSUP,IPT3 1150 : IPT3=IPT5 1151 : NBELE3=NCAVI 1152 : C On repart calculer le contour/enveloppe de la cavite 1153 : GOTO 2221 1154 : ENDIF 1155 : C 1156 : C 2.3.5 Je parcours les elements du contour/enveloppe pour 1157 : C construire le maillage de la boule 1158 : DO J=1,NBELE4 1159 : JEL=LBOUL.LECT(J) 1160 : C Attribution des premiers noeuds des nouveux elements 1161 : C de IPT2 : on choisi ceux du contour IPT4 1162 : DO K=1,IPT4.NUM(/1) 1163 : IPT2.NUM(K,JEL)=IPT4.NUM(K,J) 1164 : ENDDO 1165 : C Attribution du dernier noeud des nouveaux elements de IPT2 : 1166 : C il s'agit de NPI1 au entre de la cavite 1167 : IPT2.NUM(IDIMP1,JEL)=NPI1 1168 : C Mise a jour de la table de connectivite (partie 1 sur 2) 1169 : C uniquement sur les derniers noeuds (IDIMP1) des elements de 1170 : C la boule (adjacents extexrieurs a la boule). 1171 : C Pour cela, on utilise LEADJ et LNVADJ calcules juste avant 1172 : ILEA1.LEAC(JEL,IDIMP1,1)=LEADJ.LECT(J) 1173 : ILEA1.LEAC(JEL,IDIMP1,2)=LNVADJ.LECT(J) 1174 : IF (LEADJ.LECT(J).NE.0) THEN 1175 : ILEA1.LEAC(LEADJ.LECT(J),LNVADJ.LECT(J),1)=JEL 1176 : ILEA1.LEAC(LEADJ.LECT(J),LNVADJ.LECT(J),2)=IDIMP1 1177 : ENDIF 1178 : C et je rajoute au segment ITRC1 ce qu'il faut : 1179 : C 1179 : CALL BCIRCO(IRN1,IPT2,JEL,XCBJ,XRBJ) 1180 : C 1181 : ITRC1.TRC(JEL,1)=XRBJ 1182 : DO K=1,IDIM 1183 : ITRC1.TRC(JEL,K+1)=XCBJ(K) 1184 : ENDDO 1185 : ENDDO 1186 : C 1187 : C 2.3.6 Mise a jour de la table de connectivite (partie 2 sur 2) 1188 : C On parcoure les elements de la boule et on determine toutes les 1189 : C elements adjacents internes a la boule 1190 : JG=IDIM-1 1191 : SEGINI,LREF1 1192 : DO K=1,NBOUL 1193 : NELK=LBOUL.LECT(K) 1194 : C On recherche l'adjacent de NELK vu par le noeud M 1195 : DO M=1,IDIM 1196 : C Numeros des noeuds de l'element NELK a observer 1197 : IF (IDIM.EQ.2) THEN 1198 : IF (M.EQ.1) THEN 1199 : LREF1.LECT(1)=IPT2.NUM(2,NELK) 1200 : ELSEIF (M.EQ.2) THEN 1201 : LREF1.LECT(1)=IPT2.NUM(1,NELK) 1202 : ENDIF 1203 : ENDIF 1204 : IF (IDIM.EQ.3) THEN 1205 : IF (M.EQ.1) THEN 1206 : LREF1.LECT(1)=IPT2.NUM(2,NELK) 1207 : LREF1.LECT(2)=IPT2.NUM(3,NELK) 1208 : ELSEIF (M.EQ.2) THEN 1209 : LREF1.LECT(1)=IPT2.NUM(1,NELK) 1210 : LREF1.LECT(2)=IPT2.NUM(3,NELK) 1211 : ELSEIF (M.EQ.3) THEN 1212 : LREF1.LECT(1)=IPT2.NUM(1,NELK) 1213 : LREF1.LECT(2)=IPT2.NUM(2,NELK) 1214 : ENDIF 1215 : ENDIF 1216 : C On parcoure les autres elements de la boule 1217 : NNOKK=0 1218 : DO KK=1,NBOUL 1219 : NELKK=LBOUL.LECT(KK) 1220 : IF (NELKK.EQ.NELK) GOTO 2361 1221 : NCOM=1 1222 : NSOMKK=IDIMP1 1223 : C On parcoure les 2/3 premiers noeuds LL de NELKK 1224 : DO LL=1,IDIM 1225 : NREFX=IPT2.NUM(LL,NELKK) 1226 : C Test si le noeud LL de NELKK est commun a NELK 1227 : DO L=1,IDIM-1 1228 : IF (NREFX.EQ.LREF1.LECT(L)) THEN 1229 : NCOM=NCOM+1 1230 : NSOMKK=NSOMKK+LL 1231 : ENDIF 1232 : ENDDO 1233 : C Test si NELKK est bien un element adjacent a NELK 1234 : IF (NCOM.EQ.IDIM) GOTO 2362 1235 : ENDDO 1236 : 2361 CONTINUE 1237 : ENDDO 1238 : 2362 CONTINUE 1239 : C Determination du numero du noeud de voyeur de K dans IPT2 1240 : IF (IDIM.EQ.2) NNOKK=6-NSOMKK 1241 : IF (IDIM.EQ.3) NNOKK=10-NSOMKK 1242 : ILEA1.LEAC(NELK,M,1)=NELKK 1243 : ILEA1.LEAC(NELK,M,2)=NNOKK 1244 : ENDDO 1245 : ENDDO 1246 : C 1247 : C 2.3.7 Renumerotation et nettoyage des elements de IPT2 dans le 1248 : C cas particulier ou le nombre d'elements a diminue 1249 : IF (N43.LT.0) THEN 1250 : DO J=1,IPT2.NUM(/2) 1251 : C Determination du nombre d'elements a decaler 1252 : NDEC=0 1253 : DO K=1,NSUP 1254 : C Si l'element examine n'est plus dans la triangulation 1255 : C alors on ne fait rien 1256 : IF (J.EQ.LSUP.LECT(K)) GOTO 2371 1257 : IF (J.GE.LSUP.LECT(K)) NDEC=NDEC+1 1258 : ENDDO 1259 : C Nouveau numero que portera l'element J 1260 : JJ=J-NDEC 1261 : C Renumerotation dans IPT2 1262 : DO K=1,IPT2.NUM(/1) 1263 : IPT2.NUM(K,JJ)=IPT2.NUM(K,J) 1264 : ENDDO 1265 : C Renumerotation dans ILEA1 1266 : DO K=1,IDIMP1 1267 : IOLDVJK=ILEA1.LEAC(J,K,1) 1268 : C Determination du nombre d'elements a decaler 1269 : NDEC2=0 1270 : DO KK=1,NSUP 1271 : IF (IOLDVJK.EQ.LSUP.LECT(KK)) THEN 1272 : NBNN=IPT2.NUM(/1) 1273 : NBELEM=NBELE2 1274 : NBSOUS=0 1275 : NBREF=0 1276 : SEGADJ,IPT2 1277 : C Erreur de triangulation : vérifiez qu'il n'y a 1278 : C pas de points confondus 1278 : CALL ERREUR(846) 1279 : GOTO 2501 1280 : ENDIF 1281 : IF (IOLDVJK.GE.LSUP.LECT(KK)) NDEC2=NDEC2+1 1282 : ENDDO 1283 : NEWVJK=IOLDVJK-NDEC2 1284 : ILEA1.LEAC(JJ,K,1)=NEWVJK 1285 : ILEA1.LEAC(JJ,K,2)=ILEA1.LEAC(J,K,2) 1286 : ENDDO 1287 : C Renumerotation dans ITCR1 1288 : ITRC1.TRC(JJ,1)=ITRC1.TRC(J,1) 1289 : DO K=1,IDIM 1290 : ITRC1.TRC(JJ,K+1)=ITRC1.TRC(J,K+1) 1291 : ENDDO 1292 : 2371 CONTINUE 1293 : ENDDO 1294 : NBNN=IPT2.NUM(/1) 1295 : NBELEM=NBELE2 1296 : SEGADJ,IPT2 1297 : NBL1=NBELE2 1298 : SEGADJ,ILEA1 1299 : NBE1=NBELE2 1300 : SEGADJ,ITRC1 1301 : ENDIF 1302 : SEGSUP,IPT3,IPT4 1303 : C 1304 : C 2.3.8 Cas particulier de la dimension 1, on ajoute le point 1305 : C NPI1 au maillage 1306 : 2381 IF (IDIM.EQ.1) THEN 1307 : C Nouveau nombre d'element dans IPT2 1308 : NEW2=NBELE2+1 1309 : C J'ajuste IPT2, ITRC1 et ILEA1 si besoin 1310 : IF (NEW2.GT.IPT2.NUM(/2)) THEN 1311 : NBNN=IPT2.NUM(/1) 1312 : NBELEM=NBELE2+ICMIN 1313 : SEGADJ,IPT2 1314 : NBE1=NBELEM 1315 : SEGADJ,ITRC1 1316 : NBL1=NBELEM 1317 : SEGADJ,ILEA1 1318 : ENDIF 1319 : NBELE2=NEW2 1320 : C Modification de IPT2 pour ajouter un element 1321 : IPT2.NUM(1,NBASE1)=NREFA 1322 : IPT2.NUM(2,NBASE1)=NPI1 1323 : IPT2.NUM(1,NBELE2)=NPI1 1324 : IPT2.NUM(2,NBELE2)=NREFB 1325 : C Mise a jour de la table de connectivite ILEA1 1326 : NED=ILEA1.LEAC(NBASE1,1,1) 1327 : ILEA1.LEAC(NBASE1,1,1)=NBELE2 1328 : ILEA1.LEAC(NBASE1,1,2)=2 1329 : ILEA1.LEAC(NBELE2,2,1)=NBASE1 1330 : ILEA1.LEAC(NBELE2,2,2)=1 1331 : IF (NED.EQ.0) THEN 1332 : ILEA1.LEAC(NBELE2,1,1)=0 1333 : ILEA1.LEAC(NBELE2,1,2)=0 1334 : ELSE 1335 : ILEA1.LEAC(NBELE2,1,1)=NED 1336 : ILEA1.LEAC(NBELE2,1,2)=2 1337 : ILEA1.LEAC(NED,2,1)=NBELE2 1338 : ILEA1.LEAC(NED,2,2)=1 1339 : ENDIF 1340 : ENDIF 1341 : ENDDO 1342 : C Ajustement final de IPT2, ITRC1 et ILEA1 1343 : IF (IDIM.NE.1) THEN 1344 : NBNN=IPT2.NUM(/1) 1345 : NBELEM=NBELE2 1346 : SEGADJ,IPT2 1347 : NBE1=NBELE2 1348 : SEGADJ,ITRC1 1349 : NBL1=NBELE2 1350 : SEGADJ,ILEA1 1351 : ENDIF 1352 : C 1353 : C 1354 : C------- 2.4 Retour du maillage dans le cas ou l'on garde la boite 1355 : C Si IBOI egal a 1, on renvoie IPT2 (maillage avec boite) 1356 : 2501 CONTINUE 1357 : IF (IERR.NE.0) GOTO 999 1358 : IF (IBOI.EQ.1) THEN 1359 : IF (IDIM.NE.1) THEN 1360 : SEGSUP,LBASE,LCAVI,LBOUL,LREF1,LEADJ,LNVADJ 1361 : SEGSUP,LSUP,LSCP 1362 : SEGDES,IPT2 1363 : ENDIF 1364 : GOTO 999 1365 : ENDIF 1366 : C 1367 : C 1368 : C 1369 : C----------------------------------------------------------------------C 1370 : C PARTIE 3 C 1371 : C Elimination des elements appuyes sur les coins de la boite C 1372 : C----------------------------------------------------------------------C 1373 : C 1374 : C------- 3.1 Suppression des elements touchant les noeuds de la boite 1375 : C Je parcours les noeuds de la boite en cherchant les elements du 1376 : C maillage IPT2 qui en contiennent. Ceux qui n'en contiennent pas 1377 : C sont recopies dans IPT3, qui sera le maillage de sortie : 1378 : NBSOUS=0 1379 : NBREF=0 1380 : NBNN=IPT2.NUM(/1) 1381 : NBELEM=NBELE2 1382 : SEGINI,IPT3 1383 : IPT3.ITYPEL=IPT2.ITYPEL 1384 : NBELE3=0 1385 : DO I=1,NBELE2 1386 : DO J=1,IPT2.NUM(/1) 1387 : NPI2=IPT2.NUM(J,I) 1388 : DO K=1,NNBOIT 1389 : IF (NPI2.EQ.LNBOIT(K)) GOTO 3101 1390 : ENDDO 1391 : ENDDO 1392 : NBELE3=NBELE3+1 1393 : DO J=1,IPT2.NUM(/1) 1394 : IPT3.NUM(J,NBELE3)=IPT2.NUM(J,I) 1395 : ENDDO 1396 : 3101 CONTINUE 1397 : ENDDO 1398 : NBNN=IPT3.NUM(/1) 1399 : NBELEM=NBELE3 1400 : SEGADJ,IPT3 1401 : C 1402 : C------- 3.2 Cas ou IPT3 ne contient pas d'elements. 1403 : IF (IPT3.NUM(/2).EQ.0) THEN 1404 : C En dimension 2, cela veut dire que les points de IPT1 sont 1405 : C alignes. On va alors rendre le maillage des lignes ne touchant 1406 : C pas les noeuds de la boite 1407 : IF (IDIM.EQ.2) GOTO 3203 1408 : C En dimension 3, cela veut dire que les points sont coplanaires. 1409 : C On refait la meme operation pour rendre le maillage faces ne 1410 : C touchant pas les noeuds de la boite 1411 : JG=IDIM 1412 : SEGINI,LREF1 1413 : JG=IPT2.NUM(/2) 1414 : SEGINI,LSUP 1415 : NBSOUS=0 1416 : NBREF=0 1417 : NBNN=IDIM 1418 : NBELEM=NBELE2 1419 : SEGINI,IPT3 1420 : C On va rendre un maillage de TRI3 1421 : IPT3.ITYPEL=4 1422 : NBELE3=0 1423 : C Boucle sur les elements de la triangulation avec la boite 1424 : DO I=1,NBELE2 1425 : NNS=0 1426 : NSOM=0 1427 : NEL=0 1428 : C On regarde si le noeud J de l'element I fait partie de ceux 1429 : C de la boite 1430 : DO J=1,IPT2.NUM(/1) 1431 : NPI2=IPT2.NUM(J,I) 1432 : DO K=1,NNBOIT 1433 : IF (NPI2.EQ.LNBOIT(K)) THEN 1434 : NNS=NNS+1 1435 : NSOM=J 1436 : NEL=ILEA1.LEAC(I,J,1) 1437 : ENDIF 1438 : ENDDO 1439 : ENDDO 1440 : C Si l'element I n'a qu'un noeud faisant partie de ceux de la 1441 : C boite, alors on cree l'element appuyes sur ses autres noeuds 1442 : IF (NNS.EQ.1) THEN 1443 : DO K=1,LSUP.LECT(/1) 1444 : IF (LSUP.LECT(K).EQ.0) GOTO 3201 1445 : IF (I.EQ.LSUP.LECT(K)) THEN 1446 : GOTO 3202 1447 : ENDIF 1448 : ENDDO 1449 : 3201 CONTINUE 1450 : NBELE3=NBELE3+1 1451 : LSUP.LECT(NBELE3)=NEL 1452 : IF (NSOM.EQ.1) THEN 1453 : LREF1.LECT(1)=2 1454 : LREF1.LECT(2)=3 1455 : LREF1.LECT(3)=4 1456 : ELSEIF (NSOM.EQ.2) THEN 1457 : LREF1.LECT(1)=1 1458 : LREF1.LECT(2)=3 1459 : LREF1.LECT(3)=4 1460 : ELSEIF (NSOM.EQ.3) THEN 1461 : LREF1.LECT(1)=1 1462 : LREF1.LECT(2)=2 1463 : LREF1.LECT(3)=4 1464 : ELSEIF (NSOM.EQ.4) THEN 1465 : LREF1.LECT(1)=1 1466 : LREF1.LECT(2)=2 1467 : LREF1.LECT(3)=3 1468 : ENDIF 1469 : DO J=1,IPT3.NUM(/1) 1470 : IPT3.NUM(J,NBELE3)=IPT2.NUM(LREF1.LECT(J),I) 1471 : ENDDO 1472 : ENDIF 1473 : 3202 CONTINUE 1474 : ENDDO 1475 : NBNN=IPT3.NUM(/1) 1476 : NBELEM=NBELE3 1477 : SEGADJ,IPT3 1478 : ENDIF 1479 : C 1480 : C Cas ou les points de IPT1 sont alignes. On change le maillage en 1481 : C lignes et ne garde que les lignes ne touchant pas les noeuds de la 1482 : C boite 1483 : 3203 CONTINUE 1484 : IF (IPT3.NUM(/2).EQ.0) THEN 1484 : CALL ECROBJ('MAILLAGE',IPT2) 1484 : CALL CHANLG 1484 : CALL LIROBJ('MAILLAGE',IPT4,1,IRETOU) 1485 : IF (IERR.NE.0) RETURN 1486 : SEGACT,IPT4 1487 : NBELE4=IPT4.NUM(/2) 1488 : NBSOUS=0 1489 : NBREF=0 1490 : NBNN=IPT4.NUM(/1) 1491 : NBELEM=NBELE4 1492 : SEGINI,IPT3 1493 : C On va rendre un maillage de SEG2 1494 : IPT3.ITYPEL=IPT4.ITYPEL 1495 : NBELE3=0 1496 : C Boucle sur les lignes de la triangulation avec la boite 1497 : DO I=1,NBELE4 1498 : NNS=0 1499 : C On regarde si les noeuds l'element I font partie de ceux de 1500 : C la boite 1501 : DO J=1,IPT4.NUM(/1) 1502 : NPI4=IPT4.NUM(J,I) 1503 : DO K=1,NNBOIT 1504 : IF (NPI4.EQ.LNBOIT(K)) THEN 1505 : NNS=NNS+1 1506 : ENDIF 1507 : ENDDO 1508 : ENDDO 1509 : C Si l'element I n'a aucun noeud faisant partie de ceux de la 1510 : C boite, alors on cree l'element appuyes sur ses noeuds 1511 : IF (NNS.EQ.0) THEN 1512 : NBELE3=NBELE3+1 1513 : IPT3.NUM(1,NBELE3)=IPT4.NUM(1,I) 1514 : IPT3.NUM(2,NBELE3)=IPT4.NUM(2,I) 1515 : ENDIF 1516 : ENDDO 1517 : NBNN=IPT3.NUM(/1) 1518 : NBELEM=NBELE3 1519 : SEGADJ,IPT3 1520 : SEGDES,IPT4 1521 : ENDIF 1522 : C 1523 : C Un peu de menage avant de quitter 1524 : IF (IDIM.EQ.1) THEN 1525 : SEGSUP,LSUP 1526 : SEGDES,IPT2 1527 : ELSE 1528 : SEGDES,IPT3 1529 : SEGSUP,IPT2,LBASE,LCAVI,LBOUL,LREF1,LEADJ,LNVADJ 1530 : SEGSUP,LSUP,LSCP 1531 : ENDIF 1532 : IPT2=IPT3 1533 : NBPTS=NBPTS0 1534 : SEGADJ,MCOORD 1535 : C 1536 : 999 RETURN 1537 : END 1538 : 1539 : 1540 : 1541 : 1542 :
© Cast3M 2003 - All rights reserved.
Disclaimer