1 : C TETRA SOURCE PV 07/10/11 21:15:13 5907 2 : C---------------------------------------------------------------------| 3 : C | 4 : SUBROUTINE TETRA(II,JJ,IF1,IF2,IGAGNE,IBOUT) 5 : C | 6 : C CETTE SUBROUTINE TENTE DE CREER UN TETRAEDRE A PARTIR | 7 : C DES 2 TRIANGLES IF1 ET IF2. | 8 : C | 9 : C -1- RECHERCHE DES FACETTES IF3 ET IF4 | 10 : C -2- CREATION DES FACETTES INEXISTENTES | 11 : C -3- MISE A JOUR DU MAILLAGE DE SURFACE | 12 : C -4- TEST DU VOLUME ELEMENTAIRE CREE | 13 : C -5- MEMORISATION DU VOLUME CREE | 14 : C | 15 : C - IGAGNE=1 EN CAS DE SUCCES | 16 : C - IGAGNE=0 EN CAS D'ECHEC | 17 : C | 18 : C---------------------------------------------------------------------| 19 : C 20 : IMPLICIT INTEGER(I-N) 21 : IMPLICIT REAL*8(A-H,O-Z) 22 : -INC TDEMAIT 23 : -INC CCOPTIO 24 : LOGICAL REPONS,FACET,SOLTET,DIAGO,IN2,VERDIV,IN 25 : dimension ipts(10),jpts(10),qualti(10),qualtj(10) 26 : C 27 : inouv=0 28 : idepl=0 29 : if (nfc(4,if1).ne.0) return 30 : if (nfc(4,if2).ne.0) return 31 : C 32 : * on se refuse a mailler un tetraedre trop ouvert 33 : * if (TETA(IF1,IF2,II,JJ).lt.-100) then 34 : * write (6,*) 'if1 if2 ii jj TETA ',IF1,IF2,TETA(IF1,IF2,II,JJ) 35 : * goto 155 36 : * endif 37 : nfcini=nfcmax 38 : nptini=nptmax 39 : ip1=0 40 : ip2=0 41 : N3=0 42 : N4=0 43 : ICTF=0 44 : ICTV=0 45 : ipas=0 46 : IP=IPRED(IF1,II) 47 : JP=ISUCC(IF2,II) 48 : * 49 : if (TETA(IF1,IF2,II,JJ).lt.-100D0) goto 250 50 : C 51 : C -1- RECHERCHE DES FACETTES IF3 ET IF4 52 : C ------------------------------------- 53 : C RECHERCHE D'UNE 3EME FACETTE IF3 54 : C -------------------------------- 55 : C 56 : IF3=IFACE3(IP,II,JP) 57 : IF (IF3.NE.0) THEN 58 : * si if3 dans le mauvais sens rien a faire 59 : if (isucc(if3,ii).ne.ip) then 60 : write (6,*) ' tetra face a l''envers if3 ' 61 : return 62 : endif 63 : * SI IF3 QUADRANGULAIRE RENVOYER SUR PYRA1 64 : IF (NFC(4,IF3).NE.0) THEN 64 : CALL PYRA1(IP,II,IF1,IF3,IGAGNE,IBOUT) 65 : RETURN 66 : ENDIF 67 : N3=IF3 68 : * WRITE(6,1010)IF3 69 : 1010 FORMAT(' ** IF3=',I3) 70 : ENDIF 71 : C 72 : C 73 : C 74 : C RECHERCHE D'UNE 4EME FACETTE IF4 75 : C -------------------------------- 76 : C 77 : IF4=IFACE3(IP,JJ,JP) 78 : IF (IF4.NE.0) THEN 79 : * si if4 dans le mauvais sens rien a faire 80 : if (isucc(if4,jj).ne.jp) then 81 : write (6,*) ' tetra face a l''envers if4 ' 82 : return 83 : endif 84 : IF (NFC(4,IF4).NE.0) THEN 85 : * SI IF4 QUADRANGULAIRE RENVOYER SUR PYRA1 85 : CALL PYRA1(JJ,IP,IF1,IF4,IGAGNE,IBOUT) 86 : RETURN 87 : ENDIF 88 : N4=IF4 89 : * WRITE(6,1020)IF4 90 : 1020 FORMAT(' ** IF4=',I3) 91 : ENDIF 92 : C 93 : C -2- CREATION DES FACETTES INEXISTANTES 94 : C -------------------------------------- 95 : IF (N3.EQ.0) THEN 96 : C 97 : C CREATION DE LA FACETTE IF3 98 : C -------------------------- 99 : NFCMAX=NFCMAX+1 100 : IF3=NFCMAX 101 : ICTF=ICTF+1 102 : C 103 : NFC(1,IF3)=II 104 : NFC(2,IF3)=JP 105 : NFC(3,IF3)=IP 106 : NFC(4,IF3)=0 107 : C 108 : ENDIF 109 : C 110 : IF (N4.EQ.0) THEN 111 : C 112 : C CREATION DE LA FACETTE IF4 113 : C -------------------------- 114 : NFCMAX=NFCMAX+1 115 : IF4=NFCMAX 116 : ICTF=ICTF+1 117 : C 118 : NFC(1,IF4)=JJ 119 : NFC(2,IF4)=IP 120 : NFC(3,IF4)=JP 121 : NFC(4,IF4)=0 122 : C 123 : ENDIF 124 : C 125 : C VERIFICATION VALIDITE ARETE CREE 126 : IF (n3.EQ.0.AND.n4.EQ.0) THEN 127 : IF (DIAGO(IP,JP,diacri)) then 128 : iff=ifdiag 129 : lpts=0 130 : nfcmax=nfcini 131 : if (iff.ne.0) then 132 : if (nfc(1,iff).eq.ip.or.nfc(2,iff).eq.ip.or.nfc(3,iff).eq.ip) 133 : > then 134 : ipt1=isucc(iff,ip) 135 : if (ipt1.ne.ii.and.ipt1.ne.jj) then 136 : lpts=lpts+1 137 : ipts(lpts)=ipt1 138 : endif 139 : ipt1=ipred(iff,ip) 140 : if (ipt1.ne.ii.and.ipt1.ne.jj) then 141 : lpts=lpts+1 142 : ipts(lpts)=ipt1 143 : endif 144 : endif 145 : if (nfc(1,iff).eq.jp.or.nfc(2,iff).eq.jp.or.nfc(3,iff).eq.jp) 146 : > then 147 : ipt1=isucc(iff,jp) 148 : if (ipt1.ne.ii.and.ipt1.ne.jj) then 149 : lpts=lpts+1 150 : ipts(lpts)=ipt1 151 : endif 152 : ipt1=ipred(iff,jp) 153 : if (ipt1.ne.ii.and.ipt1.ne.jj) then 154 : lpts=lpts+1 155 : ipts(lpts)=ipt1 156 : endif 157 : endif 158 : lpts1=0 159 : do 700 lpt=1,lpts 160 : do 701 lpt1=1,lpts1 161 : if (ipts(lpt).eq.jpts(lpt1)) goto 700 162 : 701 continue 163 : qualti(lpt)=min(qualt(ii,jj,ip,ipts(lpt)), 164 : > qualt(jj,ii,jp,ipts(lpt))) 165 : if (qualti(lpt).le.0.001d0) goto 700 166 : lpts1=lpts1+1 167 : jpts(lpts1)=ipts(lpt) 168 : qualtj(lpts1)=qualti(lpt) 169 : 700 continue 170 : lpts=lpts1 171 : lpts1=0 172 : do 702 lpt=1,lpts 173 : do 703 l=1,lpts1 174 : if (qualtj(lpt).gt.qualti(l)) goto 704 175 : 703 continue 176 : l=lpts1+1 177 : 704 continue 178 : lpts1=lpts1+1 179 : do 705 l1=lpts1,l+1,-1 180 : qualti(l1)=qualti(l1-1) 181 : ipts(l1)=ipts(l1-1) 182 : 705 continue 183 : qualti(l)=qualtj(lpt) 184 : ipts(l)=jpts(lpt) 185 : 702 continue 186 : do 706 lpt=1,lpts 187 : * write (6,*) 'tetra appel tetra2 type 1 avec ', 188 : * > ipts(lpt),iff,qualti(lpt) 188 : call tetra2(ii,jj,if1,if2,igagne,ipts(lpt)) 189 : if (igagne.eq.1) return 190 : 706 continue 191 : if (npf(5,ii).ne.0.or.npf(5,jj).ne.0) goto 250 192 : * write (6,*) ' tetra diago 4-4 => point milieu',ip,jp 193 : goto 350 194 : endif 195 : * write (6,*) ' tetra echec 1 lpts',lpts 196 : return 197 : ENDIF 198 : * test in2 199 : * IF (.NOT.IN2(ip,jp,nfcini)) THEN 200 : * write (6,*) 'tetra test in2 echoue avec ',ip,jp 201 : * nfcmax=nfcini 202 : * nptmax=nptini 203 : * return 204 : * ENDIF 205 : ENDIF 206 : C 207 : C -3- MISE A JOUR DU MAILLAGE DE SURFACE 208 : C -------------------------------------- 208 : CALL REPSUB(IF1) 208 : CALL REPSUB(IF2) 208 : CALL REPSUB(IF3) 208 : CALL REPSUB(IF4) 209 : C 210 : C -4- TEST DU VOLUME ELEMENTAIRE CREE 211 : C ----------------------------------- 212 : IF (N3.EQ.0) THEN 213 : IF (.NOT.FACET(IF3)) then 214 : * write(6,*) ' tetra facet if3 invalide' 215 : GOTO 150 216 : endif 217 : ENDIF 218 : IF (N4.EQ.0) THEN 219 : IF (.NOT.FACET(IF4)) then 220 : * write(6,*) ' tetra facet if4 invalide' 221 : GOTO 150 222 : endif 223 : ENDIF 224 : * 225 : * verification longueur de l'arete a creer 226 : * 227 : DNORM=(XYZ(1,IP)-XYZ(1,JP))**2 228 : # +(XYZ(2,IP)-XYZ(2,JP))**2 229 : # +(XYZ(3,IP)-XYZ(3,JP))**2 230 : DTEST=tetrl*XYZ(4,IP)*XYZ(4,JP) 231 : * write (6,*) ' tetra dnorm dtest ',dnorm,dtest 232 : IF (DNORM.GT.DTEST) THEN 232 : CALL REPSUB(IF4) 232 : CALL REPSUB(IF3) 232 : CALL REPSUB(IF2) 232 : CALL REPSUB(IF1) 233 : NFCMAX=NFCini 234 : GOTO 250 235 : ENDIF 236 : * verification validite du tetraedre 237 : IF (.NOT.SOLTET(IF1,IF2,IF3,IF4,IPTT)) GOTO 151 238 : C 239 : C LE VOLUME CREE EST VALIDE. 240 : C -------------------------- 241 : C -5- MEMORISATION DU VOLUME CREE 242 : C ------------------------------- 243 : * write (6,*) ' tetra creation tetraedre ' 244 : NVOL=NVOL+1 245 : IF (NFV(1,IF1).EQ.0) NFV(1,IF1)=NVOL 246 : IF (NFV(1,IF1).NE.NVOL) NFV(2,IF1)=NVOL 247 : IF (NFV(1,IF2).EQ.0) NFV(1,IF2)=NVOL 248 : IF (NFV(1,IF2).NE.NVOL) NFV(2,IF2)=NVOL 249 : IF (NFV(1,IF3).EQ.0) NFV(1,IF3)=NVOL 250 : IF (NFV(1,IF3).NE.NVOL) NFV(2,IF3)=NVOL 251 : IF (NFV(1,IF4).EQ.0) NFV(1,IF4)=NVOL 252 : IF (NFV(1,IF4).NE.NVOL) NFV(2,IF4)=NVOL 253 : IVOL(9,NVOL)=25 254 : DO 130 I=1,3 255 : IVOL(I,NVOL)=NFC(I,IF1) 256 : 130 CONTINUE 257 : IVOL(4,NVOL)=JP 258 : C 259 : * WRITE(6,1070)NVOL,(IVOL(I,NVOL),I=1,9) 260 : *1070 FORMAT(I3,4X,14I4) 261 : do 961 npfi=40,1,-1 262 : if (npf(npfi,ii).ne.0) goto 962 263 : 961 continue 264 : 962 continue 265 : do 963 npfj=40,1,-1 266 : if (npf(npfj,jj).ne.0) goto 964 267 : 963 continue 268 : 964 continue 269 : qual=qualt(ii,jj,ip,jp) 270 : if (iimpi.eq.1) write (6,1100) nvol,nfcmax,nfacet,qual, 271 : > (ivol(i,nvol),i=1,4) 272 : 1100 format (' TETRA ',i5,2i6,f8.4,4i6) 273 : C 274 : * DO 140 J=1,NPTMAX 275 : * WRITE(6,1080)J,(NPF(I,J),I=1,40) 276 : *1080 FORMAT(I4,4X,40I3) 277 : * 140 CONTINUE 278 : C 279 : IGAGNE=1 280 : C 281 : RETURN 282 : C 283 : 350 CONTINUE 284 : * IP JP est invalide (diago) 285 : * IP JP II JJ ont 4 facettes 286 : * on va essayer de mettre un point au centre de 287 : * gravite de ii jj ip jp et des 2 noeuds supplementaires 288 : * des facettes touchant IP JP II JJ 289 : IF3=NOISIN(IP,II,IF1) 290 : if (nfc(4,if3).ne.0) return 291 : IP1=ISUCC(IF3,IP) 292 : IF5=NOISIN(II,JP,IF2) 293 : if (nfc(4,if5).ne.0) return 294 : IP3=ISUCC(IF5,II) 295 : IF4=NOISIN(JP,JJ,IF2) 296 : if (nfc(4,if4).ne.0) return 297 : IP2=ISUCC(IF4,JP) 298 : IF6=NOISIN(JJ,IP,IF1) 299 : if (nfc(4,if6).ne.0) return 300 : IP4=ISUCC(IF6,JJ) 301 : * write (6,*) ' II JJ IP JP IP1 IP2 IP3 IP4', 302 : * > II,JJ,IP,JP,IP1,IP2,IP3,IP4 303 : NFCMAX=NFCini 304 : idepl=3 305 : inouv=1 306 : DO 351 I=1,4 307 : XYZ(I,NPTMAX+1)=(XYZ(I,II)+XYZ(I,JJ)+XYZ(I,IP)+ 308 : > XYZ(I,JP)+(XYZ(I,IP1)+XYZ(I,IP3))/2+ 309 : > (XYZ(I,IP2)+XYZ(I,IP4))/2)/6 310 : 351 CONTINUE 310 : * call vcrit(nptmax+1) 311 : 356 CONTINUE 312 : * verif de validite de la position du point 313 : iechec=0 314 : v1=0. 315 : v2=0. 316 : v3=0. 317 : v4=0. 318 : v5=0. 319 : v6=0. 320 : v7=0. 321 : v8=0. 322 : v9=0. 323 : v10=0. 324 : v11=0. 325 : v12=0. 326 : v1=vol(nptmax+1,ii,ip,ip1) 327 : if (v1.gt.0.) then 328 : iechec=1 329 : * write (6,*) ' tetra vol 1 positif ',v1 330 : endif 331 : v2=vol(nptmax+1,ii,ip1,jp) 332 : if (v2.gt.0.) then 333 : iechec=1 334 : * write (6,*) ' tetra vol 2 positif ',v2 335 : endif 336 : if (ip1.ne.ip3) then 337 : v3=vol(nptmax+1,ii,ip,ip3) 338 : if (v3.gt.0.) then 339 : iechec=1 340 : * write (6,*) ' tetra vol 3 positif ',v3 341 : endif 342 : v4=vol(nptmax+1,ii,ip3,jp) 343 : if (v4.gt.0.) then 344 : iechec=1 345 : * write (6,*) ' tetra vol 4 positif ',v4 346 : endif 347 : endif 348 : v5=vol(nptmax+1,ip2,jj,jp) 349 : if (v5.gt.0.) then 350 : iechec=1 351 : * write (6,*) ' tetra vol 5 positif ',v5 352 : endif 353 : v6=vol(nptmax+1,ip2,ip,jj) 354 : if (v6.gt.0.) then 355 : iechec=1 356 : * write (6,*) ' tetra vol 6 positif ',v6 357 : endif 358 : if (ip2.ne.ip4) then 359 : v7=vol(nptmax+1,ip4,jj,jp) 360 : if (v7.gt.0.) then 361 : iechec=1 362 : * write (6,*) ' tetra vol 7 positif ',v7 363 : endif 364 : v8=vol(nptmax+1,ip4,ip,jj) 365 : if (v8.gt.0.) then 366 : iechec=1 367 : * write (6,*) ' tetra vol 8 positif ',v8 368 : endif 369 : endif 370 : if (ip2.eq.ip4.and.ip1.eq.ip3.and.ip1.ne.ip2) then 371 : v9=vol(nptmax+1,ip2,ip1,ip) 372 : if (v9.gt.0.) then 373 : iechec=1 374 : * write (6,*) ' tetra vol 9 positif ',v9 375 : endif 376 : v10=vol(nptmax+1,ip1,ip2,jp) 377 : if (v10.gt.0.) then 378 : iechec=1 379 : * write (6,*) ' tetra vol 10 positif ',v10 380 : endif 381 : endif 382 : v11=vol(nptmax+1,ii,jj,ip) 383 : if (v11.gt.0.) then 384 : iechec=1 385 : * write (6,*) ' tetra vol 11 positif ',v11 386 : endif 387 : v12=vol(nptmax+1,jj,ii,jp) 388 : if (v12.gt.0.) then 389 : iechec=1 390 : * write (6,*) ' tetra vol 12 positif ',v12 391 : endif 392 : vtot=v1+v2+v3+v4+v5+v6+v7+v8+v9+v10+v11+v12 393 : if (iechec.eq.1) then 394 : * deplacement du point 395 : if (idepl.ne.0) then 396 : xyz(1,nptmax+1)=1.10*xyz(1,nptmax+1)-0.10*(xyz(1,ii)+xyz(1,jj))/2 397 : xyz(2,nptmax+1)=1.10*xyz(2,nptmax+1)-0.10*(xyz(2,ii)+xyz(2,jj))/2 398 : xyz(3,nptmax+1)=1.10*xyz(3,nptmax+1)-0.10*(xyz(3,ii)+xyz(3,jj))/2 399 : xyz(4,nptmax+1)=1.10*xyz(4,nptmax+1)-0.10*(xyz(4,ii)+xyz(4,jj))/2 399 : * call vcrit(nptmax+1) 400 : * write (6,*) ' tetra 1 deplacement du point ',nptmax+1 401 : idepl=idepl-1 402 : goto 356 403 : endif 404 : * goto 250 405 : * write (6,*) ' tetra echec 3 pt plus deplacable' 406 : return 407 : endif 408 : * write (6,*) ' tetra v* ',vtot,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10 409 : if (v1.gt.vtot/1000) return 410 : if (v2.gt.vtot/1000) return 411 : if (ip1.ne.ip3.and.v3.gt.vtot/1000) return 412 : if (ip1.ne.ip3.and.v4.gt.vtot/1000) return 413 : if (v5.gt.vtot/1000) return 414 : if (v6.gt.vtot/1000) return 415 : if (ip2.ne.ip4.and.v7.gt.vtot/1000) return 416 : if (ip2.ne.ip4.and.v8.gt.vtot/1000) return 417 : if (ip2.eq.ip4.and.ip1.eq.ip3.and.ip1.ne.ip2.and. 418 : > v9.gt.vtot/1000) return 419 : if (ip2.eq.ip4.and.ip1.eq.ip3.and.ip1.ne.ip2.and. 420 : > v10.gt.vtot/1000) return 421 : if (v11.gt.vtot/1000) return 422 : if (v12.gt.vtot/1000) return 423 : NPTMAX=NPTMAX+1 424 : GOTO 352 425 : 250 CONTINUE 426 : * IL FAUT FAIRE DEUX ELEMENTS 427 : idepl=3 428 : inouv=1 429 : NFCMAX=NFCini 430 : NPTMAX=NPTMAX+1 431 : * CREATION POINT 432 : xyz(4,nptmax)=(xyz(4,ii)+xyz(4,jj))/2 433 : * deplacement du point pour l'eloigner de ii jj 434 : xn1=(xyz(2,jj)-xyz(2,ii))*(xyz(3,ip)-xyz(3,ii))- 435 : > (xyz(3,jj)-xyz(3,ii))*(xyz(2,ip)-xyz(2,ii)) 436 : yn1=(xyz(3,jj)-xyz(3,ii))*(xyz(1,ip)-xyz(1,ii))- 437 : > (xyz(1,jj)-xyz(1,ii))*(xyz(3,ip)-xyz(3,ii)) 438 : zn1=(xyz(1,jj)-xyz(1,ii))*(xyz(2,ip)-xyz(2,ii))- 439 : > (xyz(2,jj)-xyz(2,ii))*(xyz(1,ip)-xyz(1,ii)) 440 : sn1=sqrt(xn1**2+yn1**2+zn1**2) 441 : xn1=xn1/sn1 442 : yn1=yn1/sn1 443 : zn1=zn1/sn1 444 : xn2=(xyz(2,jp)-xyz(2,ii))*(xyz(3,jj)-xyz(3,ii))- 445 : > (xyz(3,jp)-xyz(3,ii))*(xyz(2,jj)-xyz(2,ii)) 446 : yn2=(xyz(3,jp)-xyz(3,ii))*(xyz(1,jj)-xyz(1,ii))- 447 : > (xyz(1,jp)-xyz(1,ii))*(xyz(3,jj)-xyz(3,ii)) 448 : zn2=(xyz(1,jp)-xyz(1,ii))*(xyz(2,jj)-xyz(2,ii))- 449 : > (xyz(2,jp)-xyz(2,ii))*(xyz(1,jj)-xyz(1,ii)) 450 : sn2=sqrt(xn2**2+yn2**2+zn2**2) 451 : xn2=xn2/sn2 452 : yn2=yn2/sn2 453 : zn2=zn2/sn2 454 : * xn=(xn1+xn2)/2 455 : * yn=(yn1+yn2)/2 456 : * zn=(zn1+zn2)/2 457 : xn=(xn1+xn2)/2 458 : yn=(yn1+yn2)/2 459 : zn=(zn1+zn2)/2 460 : sn=sqrt(xn**2+yn**2+zn**2) 461 : xn=xn/sn 462 : yn=yn/sn 463 : zn=zn/sn 464 : * xmil=(xyz(1,ii)+xyz(1,jj))/2 465 : * ymil=(xyz(2,ii)+xyz(2,jj))/2 466 : * zmil=(xyz(3,ii)+xyz(3,jj))/2 467 : xv=xyz(1,jj)-xyz(1,ii) 468 : yv=xyz(2,jj)-xyz(2,ii) 469 : zv=xyz(3,jj)-xyz(3,ii) 470 : sv=sqrt(xv**2+yv**2+zv**2) 471 : xv=xv/sv 472 : yv=yv/sv 473 : zv=zv/sv 474 : xli=xv*(xyz(1,ip)-xyz(1,ii))+yv*(xyz(2,ip)-xyz(2,ii))+ 475 : > zv*(xyz(3,ip)-xyz(3,ii)) 476 : xlj=xv*(xyz(1,jp)-xyz(1,ii))+yv*(xyz(2,jp)-xyz(2,ii))+ 477 : > zv*(xyz(3,jp)-xyz(3,ii)) 478 : * xl=(xli+xlj+2*sv+2*0)/6 479 : xl=0.5*sv 480 : * xl=(xli+xlj)/2 481 : xmil=xyz(1,ii)+xl*xv 482 : ymil=xyz(2,ii)+xl*yv 483 : zmil=xyz(3,ii)+xl*zv 484 : expf = xyz(4,nptmax) 485 : xyz(1,nptmax)=xmil+xn*expf*expfac 486 : xyz(2,nptmax)=ymil+yn*expf*expfac 487 : xyz(3,nptmax)=zmil+zn*expf*expfac 487 : * call vcrit(nptmax) 488 : IP1=0 489 : IP2=0 490 : IP3=0 491 : IP4=0 492 : goto 352 493 : 355 CONTINUE 494 : * deplacement du point 495 : xyz(1,nptmax)=0.70*xyz(1,nptmax)+0.30*(xyz(1,ii)+xyz(1,jj))/2 496 : xyz(2,nptmax)=0.70*xyz(2,nptmax)+0.30*(xyz(2,ii)+xyz(2,jj))/2 497 : xyz(3,nptmax)=0.70*xyz(3,nptmax)+0.30*(xyz(3,ii)+xyz(3,jj))/2 498 : xyz(4,nptmax)=0.70*xyz(4,nptmax)+0.30*(xyz(4,ii)+xyz(4,jj))/2 498 : * call vcrit(nptmax) 499 : * write (6,*) ' tetra 2 deplacement du point ',nptmax 500 : idepl=idepl-1 501 : 352 CONTINUE 502 : IPTT=NPTMAX 502 : CALL DIST(nptmax,nptaux,GL,IOK,ii,jj,ip,jp,ip1,ip2,ip3,ip4,0,0) 503 : IF (IOK.EQ.0) THEN 504 : * WRITE (6,*) ' tetra DIST ',nptaux 505 : if (nptaux.eq.0) then 506 : if (idepl.ne.0) goto 355 507 : NPTMAX=nptini 508 : NFCMAX=NFCini 509 : * write (6,*) ' tetra echec 4 dist et pt non deplacable' 510 : return 511 : endif 512 : if (idepl.ne.0) goto 355 513 : NPTMAX=nptini 514 : NFCMAX=NFCini 515 : iptt=nptaux 516 : idepl=0 517 : inouv=0 518 : * return 519 : ELSEIF (gl.lt.xyz(4,nptmax)/4) then 520 : * WRITE (6,*) ' tetra GL ' 521 : if (idepl.ne.0) goto 355 522 : NPTMAX=nptini 523 : NFCMAX=NFCini 524 : * write (6,*) ' tetra echec 5 gl et pt non deplacable' 525 : return 526 : ENDIF 527 : 354 continue 528 : * verification diago avant l'appel a tetra2 529 : iff=0 530 : lpts=0 531 : if (diago(iptt,ii,diacrd)) iff=ifdiag 532 : if (diago(iptt,jj,diacrd)) iff=ifdiag 533 : if (diago(iptt,ip,diacrd)) iff=ifdiag 534 : if (diago(iptt,jp,diacrd)) iff=ifdiag 535 : * write (6,*) ' tetra iff ',iff 536 : if (iff.ne.0) then 537 : nptmaa=nptmax 538 : nptmax=nptini 539 : if (nfc(4,iff).ne.0) return 540 : if (nfc(1,iff).eq.ip.or.nfc(2,iff).eq.ip.or.nfc(3,iff).eq.ip) 541 : > then 542 : ipt1=isucc(iff,ip) 543 : if (ipt1.ne.ii.and.ipt1.ne.jj.and. 544 : > ipt1.ne.jp.and.ipt1.ne.iptt) then 545 : lpts=lpts+1 546 : ipts(lpts)=ipt1 547 : endif 548 : ipt1=ipred(iff,ip) 549 : if (ipt1.ne.ii.and.ipt1.ne.jj.and. 550 : > ipt1.ne.jp.and.ipt1.ne.iptt) then 551 : lpts=lpts+1 552 : ipts(lpts)=ipt1 553 : endif 554 : endif 555 : if (nfc(1,iff).eq.jp.or.nfc(2,iff).eq.jp.or.nfc(3,iff).eq.jp) 556 : > then 557 : ipt1=isucc(iff,jp) 558 : if (ipt1.ne.ii.and.ipt1.ne.jj.and. 559 : > ipt1.ne.ip.and.ipt1.ne.iptt) then 560 : lpts=lpts+1 561 : ipts(lpts)=ipt1 562 : endif 563 : ipt1=ipred(iff,jp) 564 : if (ipt1.ne.ii.and.ipt1.ne.jj.and. 565 : > ipt1.ne.ip.and.ipt1.ne.iptt) then 566 : lpts=lpts+1 567 : ipts(lpts)=ipt1 568 : endif 569 : endif 570 : if (nfc(1,iff).eq.ii.or.nfc(2,iff).eq.ii.or.nfc(3,iff).eq.ii) 571 : > then 572 : ipt1=isucc(iff,ii) 573 : if (ipt1.ne.ip.and.ipt1.ne.jp.and. 574 : > ipt1.ne.jj.and.ipt1.ne.iptt) then 575 : lpts=lpts+1 576 : ipts(lpts)=ipt1 577 : endif 578 : ipt1=ipred(iff,ii) 579 : if (ipt1.ne.ip.and.ipt1.ne.jp.and. 580 : > ipt1.ne.jj.and.ipt1.ne.iptt) then 581 : lpts=lpts+1 582 : ipts(lpts)=ipt1 583 : endif 584 : endif 585 : if (nfc(1,iff).eq.jj.or.nfc(2,iff).eq.jj.or.nfc(3,iff).eq.jj) 586 : > then 587 : ipt1=isucc(iff,jj) 588 : if (ipt1.ne.ip.and.ipt1.ne.jp.and. 589 : > ipt1.ne.ii.and.ipt1.ne.iptt) then 590 : lpts=lpts+1 591 : ipts(lpts)=ipt1 592 : endif 593 : ipt1=ipred(iff,jj) 594 : if (ipt1.ne.ip.and.ipt1.ne.jp.and. 595 : > ipt1.ne.ii.and.ipt1.ne.iptt) then 596 : lpts=lpts+1 597 : ipts(lpts)=ipt1 598 : endif 599 : endif 600 : if (nfc(1,iff).eq.iptt.or.nfc(2,iff).eq.iptt.or. 601 : > nfc(3,iff).eq.iptt) then 602 : ipt1=isucc(iff,iptt) 603 : if (ipt1.ne.ip.and.ipt1.ne.jp.and. 604 : > ipt1.ne.ii.and.ipt1.ne.jj) then 605 : lpts=lpts+1 606 : ipts(lpts)=ipt1 607 : endif 608 : ipt1=ipred(iff,iptt) 609 : if (ipt1.ne.ip.and.ipt1.ne.jp.and. 610 : > ipt1.ne.ii.and.ipt1.ne.jj) then 611 : lpts=lpts+1 612 : ipts(lpts)=ipt1 613 : endif 614 : endif 615 : lpts1=0 616 : lpts2=lpts 617 : do 710 lpt=1,lpts 618 : do 711 lpt1=1,lpts1 619 : if (ipts(lpt).eq.jpts(lpt1)) goto 710 620 : 711 continue 621 : qualti(lpt)=min(qualt(ii,jj,ip,ipts(lpt)), 622 : > qualt(jj,ii,jp,ipts(lpt))) 623 : if (qualti(lpt).le.0.001d0) goto 710 624 : lpts1=lpts1+1 625 : jpts(lpts1)=ipts(lpt) 626 : qualtj(lpts1)=qualti(lpt) 627 : 710 continue 628 : lpts=lpts1 629 : lpts1=0 630 : do 712 lpt=1,lpts 631 : do 713 l=1,lpts1 632 : if (qualtj(lpt).gt.qualti(l)) goto 714 633 : 713 continue 634 : l=lpts1+1 635 : 714 continue 636 : lpts1=lpts1+1 637 : do 715 l1=lpts1,l+1,-1 638 : qualti(l1)=qualti(l1-1) 639 : ipts(l1)=ipts(l1-1) 640 : 715 continue 641 : qualti(l)=qualtj(lpt) 642 : ipts(l)=jpts(lpt) 643 : 712 continue 644 : do 716 lpt=1,lpts 645 : * write (6,*) 'tetra appel tetra2 type 2 avec ', 646 : * > ipts(lpt),iff,qualti(lpt) 646 : call tetra2(ii,jj,if1,if2,igagne,ipts(lpt)) 647 : if (igagne.eq.1) return 648 : 716 continue 649 : * write (6,*) ' tetra echec 6 lpts2 lpts',lpts2,lpts,ip,jp 650 : * return 651 : nptmax=nptmaa 652 : endif 653 : * test supplementaire 654 : if (inouv.eq.1) then 655 : xyz(1,nptmax+1)=1.90*xyz(1,nptmax)-0.90*(xyz(1,ii)+xyz(1,jj))/2 656 : xyz(2,nptmax+1)=1.90*xyz(2,nptmax)-0.90*(xyz(2,ii)+xyz(2,jj))/2 657 : xyz(3,nptmax+1)=1.90*xyz(3,nptmax)-0.90*(xyz(3,ii)+xyz(3,jj))/2 658 : xyz(4,nptmax+1)=1.90*xyz(4,nptmax)-0.90*(xyz(4,ii)+xyz(4,jj))/2 658 : * call vcrit(nptmax+1) 659 : IF (.NOT.IN2(ii,nptmax+1,nfcini)) THEN 660 : * deplacement du point 661 : if (idepl.ne.0) goto 355 662 : NPTMAX=nptini 663 : NFCMAX=NFCini 664 : * write (6,*) ' tetra 2eme test in2 echoue ',ii,jj,iptt 665 : RETURN 666 : ENDIF 667 : endif 668 : quala=qualt(ii,jj,ip,iptt) 669 : qualb=qualt(jj,ii,jp,iptt) 670 : * write (6,*) ' tetra quala qualb ',quala,qualb 671 : if (quala.gt.0.001.and.qualb.gt.0.001) 671 : > call tetra2(ii,jj,if1,if2,igagne,iptt) 672 : if (igagne.eq.1) return 673 : * deplacement du point 674 : if (idepl.ne.0) goto 355 675 : nptmax=nptini 676 : idjf=1 677 : goto 152 678 : 150 CONTINUE 679 : idjf=0 680 : C 681 : C LE VOLUME CREE EST INVALIDE: IL FAUT DONC DETRUIRE LES FACETTES 682 : C CREEES. ------------------------------------------------------- 682 : CALL REPSUB(IF4) 682 : CALL REPSUB(IF3) 682 : CALL REPSUB(IF2) 682 : CALL REPSUB(IF1) 683 : C 684 : NFCMAX=NFCini 685 : 152 CONTINUE 686 : C 687 : * on essaye en cas d'intersection de facette d'utiliser cette information 687 : * call facetx(k1,k2,k3,ifr) 688 : k1=iirfac 689 : k2=j1rfac 690 : k3=j2rfac 691 : if (k1.gt.nptmax) k1=0 692 : if (k2.gt.nptmax) k2=0 693 : if (k3.gt.nptmax) k3=0 694 : ifr=ifrfac 695 : * write (6,*) ' tetra retour de facetx ',k1,k2,k3,ifr 696 : if (k1.eq.0) then 697 : * write (6,*) ' tetra echec 7 facetx => 0 ' 698 : if (ifdiag.eq.0) return 699 : k1=nfc(1,ifdiag) 700 : k2=nfc(2,ifdiag) 701 : k3=nfc(3,ifdiag) 702 : endif 703 : if (ip.eq.k1.or.jp.eq.k1) then 704 : if (ifr.eq.if3.or.ifr.eq.0) then 705 : iptt=0 706 : if (jj.eq.k2) iptt=k3 707 : if (jj.eq.k3) iptt=k2 707 : if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt) 708 : if (igagne.eq.1) return 709 : endif 710 : if (ifr.eq.if4.or.ifr.eq.0) then 711 : iptt=0 712 : if (ii.eq.k2) iptt=k3 713 : if (ii.eq.k3) iptt=k2 713 : if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt) 714 : if (igagne.eq.1) return 715 : endif 716 : endif 717 : if (ip.eq.k2.or.jp.eq.k2) then 718 : if (ifr.eq.if3.or.ifr.eq.0) then 719 : iptt=0 720 : if (jj.eq.k1) iptt=k3 721 : if (jj.eq.k3) iptt=k1 721 : if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt) 722 : if (igagne.eq.1) return 723 : endif 724 : if (ifr.eq.if4.or.ifr.eq.0) then 725 : iptt=0 726 : if (ii.eq.k1) iptt=k3 727 : if (ii.eq.k3) iptt=k1 727 : if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt) 728 : if (igagne.eq.1) return 729 : endif 730 : endif 731 : if (ip.eq.k3.or.jp.eq.k3) then 732 : if (ifr.eq.if3.or.ifr.eq.0) then 733 : iptt=0 734 : if (jj.eq.k1) iptt=k2 735 : if (jj.eq.k2) iptt=k1 735 : if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt) 736 : if (igagne.eq.1) return 737 : endif 738 : if (ifr.eq.if4.or.ifr.eq.0) then 739 : iptt=0 740 : if (ii.eq.k1) iptt=k2 741 : if (ii.eq.k2) iptt=k1 741 : if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt) 742 : if (igagne.eq.1) return 743 : endif 744 : endif 745 : if (idjf.eq.1) then 746 : * write (6,*) ' tetra echec 8 apres facetx',k1,k2,k3,ifr 747 : return 748 : endif 749 : if (npf(5,ii).ne.0.or.npf(5,jj).ne.0) goto 250 750 : * write (6,*) ' tetra diago 4-4 => point milieu',ip,jp 751 : goto 350 752 : RETURN 753 : 151 CONTINUE 753 : CALL REPSUB(IF4) 753 : CALL REPSUB(IF3) 753 : CALL REPSUB(IF2) 753 : CALL REPSUB(IF1) 754 : NFCMAX=NFCini 754 : if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt) 755 : if (igagne.eq.1) return 756 : if (npf(5,ii).ne.0.or.npf(5,jj).ne.0) goto 250 757 : * write (6,*) ' tetra diago 4-4 => point milieu',ip,jp 758 : goto 350 759 : return 760 : 155 CONTINUE 761 : * write (6,*) ' tetra probleme avec nouvelle arrete ' 761 : * CALL COM33(II,JJ,IF1,IF2,IGAGNE) 762 : C 763 : * write (6,*) ' tetra echec 9' 764 : RETURN 765 : END 766 : 767 : 768 :
© Cast3M 2003 - All rights reserved.
Disclaimer