OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fvmesh.F File Reference
#include "implicit_f.inc"
#include "com01_c.inc"
#include "units_c.inc"
#include "sysunit.inc"
#include "task_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine fvmesh1 (ibuf, elem, x, ivolu, bric, xb, rvolu, nel, nbric, tbric, sfac, dxm, nba, nela, nna, tba, tfaca, tagels, ibufa, elema, tagela, ixs, nnf)
subroutine itribox (tri, box, norm, nverts, poly, nvmax)
subroutine polclip (polyin, npin, a, n, polyout, npout)
subroutine facepoly (quad, tri, ntri, ipoly, rpoly, n, normf, normt, nsmax, nnp, nrpmax, nb, nv, lmin, npolmax, nppmax, info, ielnod, xxx, elema, ibuf, nela, ibric, ifac, mvid, ilvout, ibufa, tagela, xxxa)
subroutine polyhedr (ipoly, rpoly, polb, npolb, polh, npolh, nrpmax, nphmax, ibric, lmin, info, npolhmax, nppmax, tole2)
subroutine coorloc (vpx, vpy, vpz, v1x, v1y, v1z, v2x, v2y, v2z, ksi, eta)
subroutine gpolcut (ipoly, rpoly, ipoly_old, rpoly_old, inedge, rnedge, nbnedge, nx, ny, nz, x0, y0, z0, ins, rns, nn, nhol, inz, iz, nns3, npoly, ns, ielnod, ielnod_old)
subroutine horiedge (ipoly, rpoly, nx, ny, nz, nbnedge, inedge, rnedge, x0, y0, z0, inz, nns3, nref, aref, nnsp)
subroutine horipoly (inedge, rnedge, ledge, nedge, ipoly, rpoly, iz, ielnod, npoly, nx, ny, nz, inz, ibric)

Function/Subroutine Documentation

◆ coorloc()

subroutine coorloc ( vpx,
vpy,
vpz,
v1x,
v1y,
v1z,
v2x,
v2y,
v2z,
ksi,
eta )

Definition at line 4275 of file fvmesh.F.

4278C-----------------------------------------------
4279C I m p l i c i t T y p e s
4280C-----------------------------------------------
4281#include "implicit_f.inc"
4282C-----------------------------------------------
4283C D u m m y A r g u m e n t s
4284C-----------------------------------------------
4285 my_real
4286 . vpx, vpy, vpz, v1x, v1y,
4287 . v1z, v2x, v2y, v2z, ksi,
4288 . eta
4289C-----------------------------------------------
4290C L o c a l V a r i a b l e s
4291C-----------------------------------------------
4292 my_real
4293 . sp1, sp2, s11, s22, s12, d
4294C
4295 sp1=vpx*v1x+vpy*v1y+vpz*v1z
4296 sp2=vpx*v2x+vpy*v2y+vpz*v2z
4297 s11=v1x*v1x+v1y*v1y+v1z*v1z
4298 s22=v2x*v2x+v2y*v2y+v2z*v2z
4299 s12=v1x*v2x+v1y*v2y+v1z*v2z
4300C
4301 d=s11*s22-s12*s12
4302C
4303 ksi=(sp1*s22-sp2*s12)/d
4304 eta=(sp2*s11-sp1*s12)/d
4305C
4306 RETURN
my_real
#define my_real
Definition cppsort.cpp:32

◆ facepoly()

subroutine facepoly ( quad,
integer, dimension(nsmax,*) tri,
integer ntri,
integer, dimension(6+nppmax+1+npolmax,*) ipoly,
rpoly,
n,
normf,
normt,
integer nsmax,
integer nnp,
integer nrpmax,
integer nb,
integer nv,
lmin,
integer npolmax,
integer nppmax,
integer info,
integer, dimension(nppmax,*) ielnod,
xxx,
integer, dimension(3,*) elema,
integer, dimension(*) ibuf,
integer nela,
integer ibric,
integer ifac,
integer mvid,
integer ilvout,
integer, dimension(*) ibufa,
integer, dimension(*) tagela,
xxxa )

Definition at line 3251 of file fvmesh.F.

3257 USE message_mod
3258C-----------------------------------------------
3259C I m p l i c i t T y p e s
3260C-----------------------------------------------
3261#include "implicit_f.inc"
3262C-----------------------------------------------
3263C D u m m y A r g u m e n t s
3264C-----------------------------------------------
3265 INTEGER TRI(NSMAX,*), NTRI, NPPMAX, IPOLY(6+NPPMAX+1+NPOLMAX,*),
3266 . NSMAX, NRPMAX, NNP, NB, NV, NPOLMAX, INFO,
3267 . IELNOD(NPPMAX,*), ELEMA(3,*), IBUF(*), NELA, IBRIC, IFAC,
3268 . MVID, ILVOUT, IBUFA(*), TAGELA(*)
3269 my_real
3270 . quad(3,*), rpoly(nrpmax+3*npolmax,*), n(*), normf(3,*),
3271 . normt(3,*), lmin, xxx(3,*), xxxa(3,*)
3272C-----------------------------------------------
3273C L o c a l V a r i a b l e s
3274C-----------------------------------------------
3275 INTEGER I, NINTER(4), NP, J, JJ, IDBL, NPI, NSEG, REDIR(NTRI), K,
3276 . KK, JJ1, NNO, ID1, ID2, ISEG(3,5*NTRI), NSEGF, ISTOP,
3277 . ISSEG, IPSEG, ITAGSEG(5*NTRI), ITAG(5*NTRI*2), ICLOSE,
3278 . INEXT, NN, NPOLY, IADPOLY(NTRI), LENPOLY(NTRI), IAD, II,
3279 . POLY(2,5*NTRI+1), I0, JJJ, REDIR_TMP(NTRI), NPI_TMP,
3280 . ELSEG(5*NTRI,2), ELINTER(4,NTRI), ELNODE(5*NTRI*2),
3281 . NNS, LENMAX, NEDGE, J1, J2, IEDGE, K1, K2,
3282 . TEDGE(3*NTRI,3), EDGE(3*NTRI,3), N1, N2, IPOUT, M, MM,
3283 . NHOL, IADHOL(NPOLMAX), ISEG3_OLD(5*NTRI), IDIFF, NSEC,
3284 . KSMIN
3285 my_real
3286 . tole, tole2, a(3), x1, y1, z1, x2, y2, z2, ss1, ss2,
3287 . x12, y12, z12, xa1, ya1, za1, alpha, xm, ym, zm,
3288 . stmp(4,3), seg(5*ntri,3,2), nx, ny, nz, xa2, ya2, za2,
3289 . xa12, ya12, za12, xa1p1, ya1p1, za1p1, beta,
3290 . pinter(4,ntri,4), xl(ntri), xlref, xp1, yp1, zp1, xp2,
3291 . yp2, zp2, node(3,5*ntri*2), xx1, yy1, zz1, dd1, dd2,
3292 . xlc, ll, tolei, t1, t2, dd, xa1p2, ya1p2, za1p2,
3293 . xedge(3*ntri,4), xx2, yy2, zz2, nr1, nr2, vx1, vy1, vz1,
3294 . vx2, vy2, vz2, ss, vvx, vvy, vvz, xx, yy, zz,
3295 . phol(3,npolmax), xloc(2,nppmax), xsec(nppmax), ylmin,
3296 . ylmax, ylsec, xsmin1, xsmin2, xs, ys
3297C
3298C
3299* TOLE2=EM5
3300 tole2 = zero
3301 tole2 = epsilon(tole2)
3302 tolei=tole2*ten
3303 tole=tole2*lmin
3304C
3305 a(1)=quad(1,1)
3306 a(2)=quad(2,1)
3307 a(3)=quad(3,1)
3308C Recensement des aretes de triangles
3309 nedge=0
3310 DO i=1,ntri
3311 DO j=1,3
3312 jj=j+1
3313 IF (j==3) jj=1
3314 j1=tri(i,j)
3315 j2=tri(i,jj)
3316 iedge=0
3317 DO k=1,nedge
3318 k1=edge(k,1)
3319 k2=edge(k,2)
3320 IF ((k1==j1.AND.k2==j2).OR.
3321 . (k2==j1.AND.k1==j2)) iedge=k
3322 ENDDO
3323 IF (iedge==0) THEN
3324 nedge=nedge+1
3325 tedge(i,j)=nedge
3326 edge(nedge,1)=j1
3327 edge(nedge,2)=j2
3328 ELSE
3329 tedge(i,j)=iedge
3330 ENDIF
3331 ENDDO
3332 ENDDO
3333C Calcul des intersections des aretes avec le plan de la facette
3334 DO i=1,nedge
3335 n1=edge(i,1)
3336 n2=edge(i,2)
3337 IF (n1>0) THEN
3338 x1=xxx(1,n1)
3339 y1=xxx(2,n1)
3340 z1=xxx(3,n1)
3341 ELSEIF (n1<0) THEN
3342 x1=xxxa(1,-n1)
3343 y1=xxxa(2,-n1)
3344 z1=xxxa(3,-n1)
3345 ENDIF
3346 IF (n2>0) THEN
3347 x2=xxx(1,n2)
3348 y2=xxx(2,n2)
3349 z2=xxx(3,n2)
3350 ELSEIF (n2<0) THEN
3351 x2=xxxa(1,-n2)
3352 y2=xxxa(2,-n2)
3353 z2=xxxa(3,-n2)
3354 ENDIF
3355 ss1=(x1-a(1))*n(1)+(y1-a(2))*n(2)+(z1-a(3))*n(3)
3356 ss2=(x2-a(1))*n(1)+(y2-a(2))*n(2)+(z2-a(3))*n(3)
3357 edge(i,3)=0
3358 IF (abs(ss1)<=tole.OR.abs(ss2)<=tole) THEN
3359 IF (abs(ss1)<=tole.AND.abs(ss2)<=tole) cycle
3360 ELSEIF (ss1<zero.AND.ss2<zero) THEN
3361 cycle
3362 ELSEIF (ss1>=zero.AND.ss2>=zero) THEN
3363 cycle
3364 ENDIF
3365 x12=x2-x1
3366 y12=y2-y1
3367 z12=z2-z1
3368 xa1=x1-a(1)
3369 ya1=y1-a(2)
3370 za1=z1-a(3)
3371 alpha=-(xa1*n(1)+ya1*n(2)+za1*n(3))
3372 . /(x12*n(1)+y12*n(2)+z12*n(3))
3373 xm=x1+alpha*x12
3374 ym=y1+alpha*y12
3375 zm=z1+alpha*z12
3376C
3377 edge(i,3)=1
3378 xedge(i,1)=xm
3379 xedge(i,2)=ym
3380 xedge(i,3)=zm
3381 xedge(i,4)=ss1*ss2
3382 ENDDO
3383C Calcul des intersections des triangles avec le plan de la facette (segments)
3384 DO i=1,4
3385 ninter(i)=0
3386 ENDDO
3387 DO i=1,ntri
3388 np=0
3389 elseg(i,1)=i
3390 elseg(i,2)=i
3391 DO j=1,3
3392 iedge=tedge(i,j)
3393 IF (edge(iedge,3)==0) cycle
3394 np=np+1
3395 stmp(1,np)=xedge(iedge,1)
3396 stmp(2,np)=xedge(iedge,2)
3397 stmp(3,np)=xedge(iedge,3)
3398 stmp(4,np)=xedge(iedge,4)
3399 ENDDO
3400 IF (np==0) THEN
3401 seg(i,1,1)=zero
3402 seg(i,2,1)=zero
3403 seg(i,3,1)=zero
3404 seg(i,1,2)=zero
3405 seg(i,2,2)=zero
3406 seg(i,3,2)=zero
3407 cycle
3408 ELSEIF (np==3) THEN
3409 np=0
3410 idbl=0
3411 DO j=1,3
3412 IF (stmp(4,j)<zero) THEN
3413 np=np+1
3414 seg(i,1,np)=stmp(1,j)
3415 seg(i,2,np)=stmp(2,j)
3416 seg(i,3,np)=stmp(3,j)
3417 ELSEIF (stmp(4,j)==zero) THEN
3418 IF (idbl==0) THEN
3419 np=np+1
3420 seg(i,1,np)=stmp(1,j)
3421 seg(i,2,np)=stmp(2,j)
3422 seg(i,3,np)=stmp(3,j)
3423 idbl=1
3424 ENDIF
3425 ENDIF
3426 ENDDO
3427 ELSEIF (np==2) THEN
3428 seg(i,1,1)=stmp(1,1)
3429 seg(i,2,1)=stmp(2,1)
3430 seg(i,3,1)=stmp(3,1)
3431 seg(i,1,2)=stmp(1,2)
3432 seg(i,2,2)=stmp(2,2)
3433 seg(i,3,2)=stmp(3,2)
3434 ENDIF
3435C Intersection (eventuelle) de chaque segment avec les aretes de la facette
3436 DO j=1,4
3437 npi=ninter(j)
3438 IF (j/=4) THEN
3439 jj=j+1
3440 ELSE
3441 jj=1
3442 ENDIF
3443 nx=normf(1,j)
3444 ny=normf(2,j)
3445 nz=normf(3,j)
3446 x1=seg(i,1,1)
3447 y1=seg(i,2,1)
3448 z1=seg(i,3,1)
3449 x2=seg(i,1,2)
3450 y2=seg(i,2,2)
3451 z2=seg(i,3,2)
3452 x12=x2-x1
3453 y12=y2-y1
3454 z12=z2-z1
3455 xa1=quad(1,j)
3456 ya1=quad(2,j)
3457 za1=quad(3,j)
3458 xa2=quad(1,jj)
3459 ya2=quad(2,jj)
3460 za2=quad(3,jj)
3461 xa12=xa2-xa1
3462 ya12=ya2-ya1
3463 za12=za2-za1
3464 xa1p1=x1-xa1
3465 ya1p1=y1-ya1
3466 za1p1=z1-za1
3467 xa1p2=x2-xa1
3468 ya1p2=y2-ya1
3469 za1p2=z2-za1
3470C
3471 ss1=xa1p1*nx+ya1p1*ny+za1p1*nz
3472 ss2=xa1p2*nx+ya1p2*ny+za1p2*nz
3473 IF (ss1>zero.AND.ss2>zero) THEN
3474 seg(i,1,1)=zero
3475 seg(i,2,1)=zero
3476 seg(i,3,1)=zero
3477 seg(i,1,2)=zero
3478 seg(i,2,2)=zero
3479 seg(i,3,2)=zero
3480 cycle
3481 ENDIF
3482C
3483 ss2=x12*nx+y12*ny+z12*nz
3484 IF (ss2==zero) cycle
3485 alpha=-ss1/ss2
3486 IF (alpha<-tolei.OR.alpha>one+tolei) cycle
3487 xm=x1+alpha*x12
3488 ym=y1+alpha*y12
3489 zm=z1+alpha*z12
3490 beta=((xm-xa1)*xa12+(ym-ya1)*ya12+(zm-za1)*za12)
3491 . /(xa12**2+ya12**2+za12**2)
3492 IF (beta<-tolei.OR.beta>one+tolei) cycle
3493 ss1=(x1-xa1)*nx+(y1-ya1)*ny+(z1-za1)*nz
3494 ss2=(x2-xa1)*nx+(y2-ya1)*ny+(z2-za1)*nz
3495C
3496 IF (ss1*ss2>zero) THEN
3497 IF (abs(ss1)<=abs(ss2)) THEN
3498 seg(i,1,1)=xm
3499 seg(i,2,1)=ym
3500 seg(i,3,1)=zm
3501 ELSE
3502 seg(i,1,2)=xm
3503 seg(i,2,2)=ym
3504 seg(i,3,2)=zm
3505 ENDIF
3506 ELSEIF (ss1*ss2<zero) THEN
3507 IF (ss1<zero) THEN
3508 seg(i,1,2)=xm
3509 seg(i,2,2)=ym
3510 seg(i,3,2)=zm
3511 ELSEIF (ss2<zero) THEN
3512 seg(i,1,1)=xm
3513 seg(i,2,1)=ym
3514 seg(i,3,1)=zm
3515 ENDIF
3516 ELSE
3517 IF (ss1==zero) THEN
3518 IF (ss2<zero) THEN
3519 seg(i,1,1)=xm
3520 seg(i,2,1)=ym
3521 seg(i,3,1)=zm
3522 ELSEIF (ss2>=zero) THEN
3523 seg(i,1,2)=xm
3524 seg(i,2,2)=ym
3525 seg(i,3,2)=zm
3526 ENDIF
3527 ELSEIF (ss2==zero) THEN
3528 IF (ss1<zero) THEN
3529 seg(i,1,2)=xm
3530 seg(i,2,2)=ym
3531 seg(i,3,2)=zm
3532 ELSEIF (ss1>=zero) THEN
3533 seg(i,1,1)=xm
3534 seg(i,2,1)=ym
3535 seg(i,3,1)=zm
3536 ENDIF
3537 ENDIF
3538 ENDIF
3539C
3540 npi=npi+1
3541 ninter(j)=npi
3542 pinter(j,npi,1)=xm
3543 pinter(j,npi,2)=ym
3544 pinter(j,npi,3)=zm
3545 elinter(j,npi)=i
3546 nx=normt(1,i)
3547 ny=normt(2,i)
3548 nz=normt(3,i)
3549 ss1=xa12*nx+ya12*ny+za12*nz
3550 IF (ss1<=zero) THEN
3551 pinter(j,npi,4)=one
3552 ELSEIF (ss1>=zero) THEN
3553 pinter(j,npi,4)=-one
3554 ENDIF
3555 ENDDO
3556 ENDDO
3557C
3558 nseg=ntri
3559 DO i=1,4
3560 IF (i/=4) THEN
3561 ii=i+1
3562 ELSE
3563 ii=1
3564 ENDIF
3565 xa1=quad(1,i)
3566 ya1=quad(2,i)
3567 za1=quad(3,i)
3568 xa2=quad(1,ii)
3569 ya2=quad(2,ii)
3570 za2=quad(3,ii)
3571 xa12=xa2-xa1
3572 ya12=ya2-ya1
3573 za12=za2-za1
3574 IF (ninter(i)==0) cycle
3575C Tri des points sur les bords de l'arete
3576 npi=ninter(i)
3577 DO j=1,npi
3578 xm=pinter(i,j,1)
3579 ym=pinter(i,j,2)
3580 zm=pinter(i,j,3)
3581 xl(j)=((xm-xa1)*xa12+(ym-ya1)*ya12+(zm-za1)*za12)
3582 . /(xa12**2+ya12**2+za12**2)
3583 ENDDO
3584 DO j=1,npi
3585 redir(j)=j
3586 ENDDO
3587 DO j=1,npi
3588 xlref=xl(redir(j))
3589 DO k=j+1,npi
3590 xlc=xl(redir(k))
3591 IF (xlc<=xlref) THEN
3592 kk=redir(k)
3593 redir(k)=redir(j)
3594 redir(j)=kk
3595 xlref=xlc
3596 ENDIF
3597 ENDDO
3598 ENDDO
3599C Identification des doubles
3600 DO j=1,npi
3601 jj=redir(j)
3602 IF (jj>0) THEN
3603 x1=pinter(i,jj,1)
3604 y1=pinter(i,jj,2)
3605 z1=pinter(i,jj,3)
3606 t1=pinter(i,jj,4)
3607 DO k=j+1,npi
3608 kk=redir(k)
3609 IF (kk>0) THEN
3610 x2=pinter(i,kk,1)
3611 y2=pinter(i,kk,2)
3612 z2=pinter(i,kk,3)
3613 t2=pinter(i,kk,4)
3614 dd=(x2-x1)**2+(y2-y1)**2+(z2-z1)**2+(t2-t1)**2
3615 IF (dd<=tole) redir(k)=0
3616 ENDIF
3617 ENDDO
3618 ENDIF
3619 ENDDO
3620C Elimination des segments de longueur nulle
3621 DO j=1,npi
3622 jj=redir(j)
3623 IF (jj>0) THEN
3624 x1=pinter(i,jj,1)
3625 y1=pinter(i,jj,2)
3626 z1=pinter(i,jj,3)
3627 DO k=j+1,npi
3628 kk=redir(k)
3629 IF (kk>0) THEN
3630 x2=pinter(i,kk,1)
3631 y2=pinter(i,kk,2)
3632 z2=pinter(i,kk,3)
3633 dd=(x2-x1)**2+(y2-y1)**2+(z2-z1)**2
3634 IF (dd<=tole) THEN
3635 redir(j)=0
3636 redir(k)=0
3637 ENDIF
3638 ENDIF
3639 ENDDO
3640 ENDIF
3641 ENDDO
3642C Condensation du tableau REDIR
3643 DO j=1,npi
3644 redir_tmp(j)=redir(j)
3645 ENDDO
3646 npi_tmp=npi
3647 npi=0
3648 DO j=1,npi_tmp
3649 IF (redir_tmp(j)>0) THEN
3650 npi=npi+1
3651 redir(npi)=redir_tmp(j)
3652 ENDIF
3653 ENDDO
3654 IF (npi==0) cycle
3655C Creation des segments sur les aretes
3656 IF (pinter(i,redir(1),4)==-one) THEN
3657 nseg=nseg+1
3658 seg(nseg,1,1)=xa1
3659 seg(nseg,2,1)=ya1
3660 seg(nseg,3,1)=za1
3661 elseg(nseg,1)=0
3662 seg(nseg,1,2)=pinter(i,redir(1),1)
3663 seg(nseg,2,2)=pinter(i,redir(1),2)
3664 seg(nseg,3,2)=pinter(i,redir(1),3)
3665 elseg(nseg,2)=elinter(i,redir(1))
3666 ENDIF
3667 DO j=1,npi-1
3668 jj=redir(j)
3669 jj1=redir(j+1)
3670 IF (pinter(i,jj,4)==-one) cycle
3671 xp1=pinter(i,jj,1)
3672 yp1=pinter(i,jj,2)
3673 zp1=pinter(i,jj,3)
3674 xp2=pinter(i,jj1,1)
3675 yp2=pinter(i,jj1,2)
3676 zp2=pinter(i,jj1,3)
3677 nseg=nseg+1
3678 seg(nseg,1,1)=xp1
3679 seg(nseg,2,1)=yp1
3680 seg(nseg,3,1)=zp1
3681 elseg(nseg,1)=elinter(i,jj)
3682 seg(nseg,1,2)=xp2
3683 seg(nseg,2,2)=yp2
3684 seg(nseg,3,2)=zp2
3685 elseg(nseg,2)=elinter(i,jj1)
3686 ENDDO
3687 IF (pinter(i,redir(npi),4)==one) THEN
3688 nseg=nseg+1
3689 seg(nseg,1,1)=pinter(i,redir(npi),1)
3690 seg(nseg,2,1)=pinter(i,redir(npi),2)
3691 seg(nseg,3,1)=pinter(i,redir(npi),3)
3692 elseg(nseg,1)=elinter(i,redir(npi))
3693 seg(nseg,1,2)=xa2
3694 seg(nseg,2,2)=ya2
3695 seg(nseg,3,2)=za2
3696 elseg(nseg,2)=0
3697 ENDIF
3698 ENDDO
3699C Recensement des noeuds
3700 nno=0
3701 DO i=1,nseg
3702 x1=seg(i,1,1)
3703 y1=seg(i,2,1)
3704 z1=seg(i,3,1)
3705 x2=seg(i,1,2)
3706 y2=seg(i,2,2)
3707 z2=seg(i,3,2)
3708 IF (x1==zero.AND.y1==zero.AND.z1==zero.AND.
3709 . x2==zero.AND.y2==zero.AND.z2==zero) THEN
3710 iseg(1,i)=0
3711 iseg(2,i)=0
3712 iseg(3,i)=1
3713 cycle
3714 ENDIF
3715 id1=0
3716 id2=0
3717 DO j=1,nno
3718 xx1=node(1,j)
3719 yy1=node(2,j)
3720 zz1=node(3,j)
3721 dd1=(xx1-x1)**2+(yy1-y1)**2+(zz1-z1)**2
3722 dd2=(xx1-x2)**2+(yy1-y2)**2+(zz1-z2)**2
3723 IF (dd1<=tole) id1=j
3724 IF (dd2<=tole) id2=j
3725 ENDDO
3726 ll=(x2-x1)**2+(y2-y1)**2+(z2-z1)**2
3727 IF (id1==0) THEN
3728 nno=nno+1
3729 node(1,nno)=x1
3730 node(2,nno)=y1
3731 node(3,nno)=z1
3732 elnode(nno)=elseg(i,1)
3733 id1=nno
3734 ENDIF
3735 IF (ll<=tole) id2=id1
3736 IF (id2==0) THEN
3737 nno=nno+1
3738 node(1,nno)=x2
3739 node(2,nno)=y2
3740 node(3,nno)=z2
3741 elnode(nno)=elseg(i,2)
3742 id2=nno
3743 ENDIF
3744 iseg(1,i)=id1
3745 iseg(2,i)=id2
3746 iseg(3,i)=0
3747 IF (id1==id2) iseg(3,i)=1
3748 ENDDO
3749C Elimination des segments doubles
3750 DO i=1,nseg
3751 IF (iseg(3,i)==1) cycle
3752 id1=iseg(1,i)
3753 id2=iseg(2,i)
3754 DO j=i+1,nseg
3755 IF ((iseg(1,j)==id1.AND.iseg(2,j)==id2).OR.
3756 . (iseg(1,j)==id2.AND.iseg(2,j)==id1)) iseg(3,j)=1
3757 ENDDO
3758 ENDDO
3759C
3760 nsegf=0
3761 DO i=1,nseg
3762 IF (iseg(3,i)==0) nsegf=nsegf+1
3763 ENDDO
3764 IF (nsegf<=1) RETURN
3765C
3766C Reconstruction des polygones
3767 istop=0
3768 npoly=0
3769 np=0
3770C
3771 DO i=1,nseg
3772 iseg3_old(i)=iseg(3,i)
3773 ENDDO
3774C
3775 DO WHILE (istop==0)
3776 i=1
3777 DO WHILE (iseg(3,i)==1.AND.i<=nseg)
3778 i=i+1
3779 ENDDO
3780 IF (i==nseg+1) THEN
3781 istop=1
3782 cycle
3783 ENDIF
3784 isseg=i
3785 ipseg=2
3786 DO j=1,nno
3787 itag(j)=0
3788 ENDDO
3789 DO j=1,nseg
3790 itagseg(j)=0
3791 ENDDO
3792 itag(iseg(1,isseg))=i
3793 itagseg(isseg)=i
3794 iclose=0
3795 i0=i-1
3796 DO WHILE (iclose==0.AND.i<nseg)
3797 i=i+1
3798 inext=0
3799 DO j=1,nseg
3800 IF (j==isseg.OR.iseg(3,j)==1.OR.inext/=0) cycle
3801 id1=iseg(1,j)
3802 id2=iseg(2,j)
3803 IF (id1==iseg(ipseg,isseg)) THEN
3804 inext=1
3805 isseg=j
3806 ipseg=2
3807 itag(iseg(1,isseg))=i
3808 itagseg(j)=i
3809 ELSEIF (id2==iseg(ipseg,isseg)) THEN
3810 inext=1
3811 isseg=j
3812 ipseg=1
3813 itag(iseg(2,isseg))=i
3814 itagseg(j)=-i
3815 ENDIF
3816 ENDDO
3817 IF (inext==0) THEN
3818 iseg(3,isseg)=1
3819 i=nseg
3820 ELSE
3821 nn=itag(iseg(ipseg,isseg))
3822 IF (nn/=0) iclose=nn
3823 ENDIF
3824 ENDDO
3825 IF (iclose/=0) THEN
3826 npoly=npoly+1
3827 iadpoly(npoly)=np
3828 iad=np
3829 DO j=1,nseg
3830 jj=itagseg(j)
3831 IF (abs(jj)>=iclose) THEN
3832 np=np+1
3833 poly(1,iad+abs(jj)-i0)=j
3834 IF (jj>0) THEN
3835 poly(2,iad+abs(jj)-i0)=2
3836 ELSEIF (jj<0) THEN
3837 poly(2,iad+abs(jj)-i0)=1
3838 ENDIF
3839 iseg(3,j)=1
3840 ENDIF
3841 ENDDO
3842 lenpoly(npoly)=np-iadpoly(npoly)
3843 ENDIF
3844C
3845 idiff=0
3846 DO j=1,nseg
3847 idiff=max(idiff,abs(iseg3_old(j)-iseg(3,j)))
3848 ENDDO
3849 IF (idiff==0) THEN
3850 CALL ancmsg(msgid=9,anmode=aninfo,
3851 . i1=mvid)
3852 IF (ilvout/=0) THEN
3853 CALL ancmsg(msgid=10,anmode=aninfo_blind,
3854 . i1=ibric,i2=ifac)
3855 END IF
3856 CALL arret(2)
3857 ENDIF
3858 ENDDO
3859C
3860 info=0
3861 IF (npoly>npolmax) THEN
3862 npolmax=npoly
3863 info=1
3864 ENDIF
3865 lenmax=0
3866 DO i=1,npoly
3867 lenmax=max(lenmax,lenpoly(i))
3868 ENDDO
3869 IF (lenmax>nppmax) THEN
3870 nppmax=lenmax
3871 info=1
3872 ENDIF
3873 IF (info==1) RETURN
3874 nnp=npoly
3875C
3876 nns=0
3877 DO i=1,npoly
3878 iad=iadpoly(i)
3879 ipoly(1,i)=2
3880 ipoly(2,i)=lenpoly(i)
3881 ipoly(3,i)=nb
3882 ipoly(4,i)=nv
3883 ipoly(5,i)=0
3884 ipoly(6,i)=0
3885 rpoly(1,i)=zero
3886 rpoly(2,i)=n(1)
3887 rpoly(3,i)=n(2)
3888 rpoly(4,i)=n(3)
3889 DO j=1,lenpoly(i)
3890 jj=poly(1,iad+j)
3891 jjj=poly(2,iad+j)
3892 id1=iseg(jjj,jj)
3893 nns=nns+1
3894 ipoly(6+j,i)=nns
3895 ielnod(j,i)=elnode(id1)
3896 rpoly(4+3*(j-1)+1,i)=node(1,id1)
3897 rpoly(4+3*(j-1)+2,i)=node(2,id1)
3898 rpoly(4+3*(j-1)+3,i)=node(3,id1)
3899 ENDDO
3900 ENDDO
3901C
3902C Recherche des trous
3903 info=0
3904 DO i=1,npoly
3905C
3906 nhol=0
3907 DO j=1,npoly
3908 iadhol(j)=0
3909 ENDDO
3910C
3911 DO j=1,npoly
3912 IF (i==j) cycle
3913 alpha=zero
3914 x1=rpoly(5,j)
3915 y1=rpoly(6,j)
3916 z1=rpoly(7,j)
3917 DO k=1,lenpoly(i)
3918 kk=k+1
3919 IF (k==lenpoly(i)) kk=1
3920 xx1=rpoly(4+3*(k-1)+1,i)
3921 yy1=rpoly(4+3*(k-1)+2,i)
3922 zz1=rpoly(4+3*(k-1)+3,i)
3923 xx2=rpoly(4+3*(kk-1)+1,i)
3924 yy2=rpoly(4+3*(kk-1)+2,i)
3925 zz2=rpoly(4+3*(kk-1)+3,i)
3926 vx1=xx1-x1
3927 vy1=yy1-y1
3928 vz1=zz1-z1
3929 vx2=xx2-x1
3930 vy2=yy2-y1
3931 vz2=zz2-z1
3932 nr1=sqrt(vx1**2+vy1**2+vz1**2)
3933 nr2=sqrt(vx2**2+vy2**2+vz2**2)
3934 vx1=vx1/nr1
3935 vy1=vy1/nr1
3936 vz1=vz1/nr1
3937 vx2=vx2/nr2
3938 vy2=vy2/nr2
3939 vz2=vz2/nr2
3940 ss=vx1*vx2+vy1*vy2+vz1*vz2
3941 vvx=vy1*vz2-vz1*vy2
3942 vvy=vz1*vx2-vx1*vz2
3943 vvz=vx1*vy2-vy1*vx2
3944 ss1=n(1)*vvx+n(2)*vvy+n(3)*vvz
3945 IF (ss1>=zero) THEN
3946 alpha=alpha+acos(ss)
3947 ELSE
3948 alpha=alpha-acos(ss)
3949 ENDIF
3950 ENDDO
3951C
3952 IF (abs(alpha)>=one) THEN
3953C Le premier point du polygone i est a l'interieur du polygone j
3954C On teste tous les autres
3955 ipout=0
3956 DO k=2,lenpoly(j)
3957 x2=rpoly(4+3*(k-1)+1,j)
3958 y2=rpoly(4+3*(k-1)+2,j)
3959 z2=rpoly(4+3*(k-1)+3,j)
3960 alpha=zero
3961 DO m=1,lenpoly(i)
3962 mm=m+1
3963 IF (m==lenpoly(i)) mm=1
3964 xx1=rpoly(4+3*(m-1)+1,i)
3965 yy1=rpoly(4+3*(m-1)+2,i)
3966 zz1=rpoly(4+3*(m-1)+3,i)
3967 xx2=rpoly(4+3*(mm-1)+1,i)
3968 yy2=rpoly(4+3*(mm-1)+2,i)
3969 zz2=rpoly(4+3*(mm-1)+3,i)
3970 vx1=xx1-x1
3971 vy1=yy1-y1
3972 vz1=zz1-z1
3973 vx2=xx2-x1
3974 vy2=yy2-y1
3975 vz2=zz2-z1
3976 nr1=sqrt(vx1**2+vy1**2+vz1**2)
3977 nr2=sqrt(vx2**2+vy2**2+vz2**2)
3978 vx1=vx1/nr1
3979 vy1=vy1/nr1
3980 vz1=vz1/nr1
3981 vx2=vx2/nr2
3982 vy2=vy2/nr2
3983 vz2=vz2/nr2
3984 ss=vx1*vx2+vy1*vy2+vz1*vz2
3985 vvx=vy1*vz2-vz1*vy2
3986 vvy=vz1*vx2-vx1*vz2
3987 vvz=vx1*vy2-vy1*vx2
3988 ss1=n(1)*vvx+n(2)*vvy+n(3)*vvz
3989 IF (ss1>=zero) THEN
3990 alpha=alpha+acos(ss)
3991 ELSE
3992 alpha=alpha-acos(ss)
3993 ENDIF
3994 ENDDO
3995 IF (abs(alpha)<one) ipout=1
3996 ENDDO
3997C
3998 IF (ipout==1) cycle
3999C Le polygone j est un trou dans le polygone i
4000 ipoly(1,j)=-1
4001 nhol=nhol+1
4002 iadhol(nhol)=lenpoly(i)
4003 IF (lenpoly(i)+lenpoly(j)>nppmax) THEN
4004 nppmax=lenpoly(i)+lenpoly(j)
4005 info=1
4006 ELSE
4007 DO k=1,lenpoly(j)
4008 ipoly(6+iadhol(nhol)+k,i)=ipoly(6+k,j)
4009 ielnod(iadhol(nhol)+k,i)=ielnod(k,j)
4010 rpoly(4+3*iadhol(nhol)+3*(k-1)+1,i)=
4011 . rpoly(4+3*(k-1)+1,j)
4012 rpoly(4+3*iadhol(nhol)+3*(k-1)+2,i)=
4013 . rpoly(4+3*(k-1)+2,j)
4014 rpoly(4+3*iadhol(nhol)+3*(k-1)+3,i)=
4015 . rpoly(4+3*(k-1)+3,j)
4016 ENDDO
4017 ENDIF
4018 lenpoly(i)=lenpoly(i)+lenpoly(j)
4019C Point interieur polygone j
4020 vx1=quad(1,2)-quad(1,1)
4021 vy1=quad(2,2)-quad(2,1)
4022 vz1=quad(3,2)-quad(3,1)
4023 ss=sqrt(vx1**2+vy1**2+vz1**2)
4024 vx1=vx1/ss
4025 vy1=vy1/ss
4026 vz1=vz1/ss
4027 vx2=n(2)*vz1-n(3)*vy1
4028 vy2=n(3)*vx1-n(1)*vz1
4029 vz2=n(1)*vy1-n(2)*vx1
4030 x1=rpoly(5,j)
4031 y1=rpoly(6,j)
4032 z1=rpoly(7,j)
4033 xloc(1,1)=zero
4034 xloc(2,1)=zero
4035 ylmin=ep30
4036 ylmax=-ep30
4037 DO k=2,lenpoly(j)
4038 xx=rpoly(4+3*(k-1)+1,j)
4039 yy=rpoly(4+3*(k-1)+2,j)
4040 zz=rpoly(4+3*(k-1)+3,j)
4041 vvx=xx-x1
4042 vvy=yy-y1
4043 vvz=zz-z1
4044 xloc(1,k)=vvx*vx1+vvy*vy1+vvz*vz1
4045 xloc(2,k)=vvx*vx2+vvy*vy2+vvz*vz2
4046 IF (xloc(2,k)<ylmin) ylmin=xloc(2,k)
4047 IF (xloc(2,k)>ylmax) ylmax=xloc(2,k)
4048 ENDDO
4049 ylsec=half*(ylmin+ylmax)
4050C
4051 nsec=0
4052 DO k=1,lenpoly(j)
4053 kk=k+1
4054 IF (k==lenpoly(j)) kk=1
4055 x1=xloc(1,k)
4056 y1=xloc(2,k)
4057 x2=xloc(1,kk)
4058 y2=xloc(2,kk)
4059 IF (y1-y2/=zero) THEN
4060 alpha=(ylsec-y2)/(y1-y2)
4061 IF (alpha>=zero.AND.alpha<=one) THEN
4062 nsec=nsec+1
4063 xsec(nsec)=alpha*x1+(one-alpha)*x2
4064 ENDIF
4065 ELSE
4066 IF (y1==ylsec) THEN
4067 nsec=nsec+1
4068 xsec(nsec)=x1
4069 nsec=nsec+1
4070 xsec(nsec)=x2
4071 ENDIF
4072 ENDIF
4073 ENDDO
4074C
4075 xsmin1=ep30
4076 DO k=1,nsec
4077 IF (xsec(k)<xsmin1) THEN
4078 xsmin1=xsec(k)
4079 ksmin=k
4080 ENDIF
4081 ENDDO
4082 xsmin2=ep30
4083 DO k=1,nsec
4084 IF (k==ksmin) cycle
4085 IF (xsec(k)<xsmin2) xsmin2=xsec(k)
4086 ENDDO
4087C
4088 xs=half*(xsmin1+xsmin2)
4089 ys=ylsec
4090 phol(1,nhol)=rpoly(5,j)+xs*vx1+ys*vx2
4091 phol(2,nhol)=rpoly(6,j)+xs*vy1+ys*vy2
4092 phol(3,nhol)=rpoly(7,j)+xs*vz1+ys*vz2
4093 ENDIF
4094 ENDDO
4095C
4096 IF (info==0) THEN
4097 ipoly(2,i)=lenpoly(i)
4098 ipoly(6+lenpoly(i)+1,i)=nhol
4099 DO j=1,nhol
4100 ipoly(6+lenpoly(i)+1+j,i)=iadhol(j)
4101 rpoly(4+3*lenpoly(i)+3*(j-1)+1,i)=phol(1,j)
4102 rpoly(4+3*lenpoly(i)+3*(j-1)+2,i)=phol(2,j)
4103 rpoly(4+3*lenpoly(i)+3*(j-1)+3,i)=phol(3,j)
4104 ENDDO
4105 ENDIF
4106 ENDDO
4107C
4108 RETURN
alpha
#define alpha
Definition eval.h:35
max
#define max(a, b)
Definition macros.h:21
message_mod
Definition message_mod.F:1249
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:889
arret
subroutine arret(nn)
Definition arret.F:87

◆ fvmesh1()

subroutine fvmesh1 ( integer, dimension(*) ibuf,
integer, dimension(3,*) elem,
x,
integer, dimension(*) ivolu,
integer, dimension(8,*) bric,
xb,
rvolu,
integer nel,
integer nbric,
integer, dimension(13,*) tbric,
sfac,
dxm,
integer nba,
integer nela,
integer nna,
integer, dimension(2,*) tba,
integer, dimension(12,*) tfaca,
integer, dimension(*) tagels,
integer, dimension(*) ibufa,
integer, dimension(3,*) elema,
integer, dimension(*) tagela,
integer, dimension(nixs,*) ixs,
integer nnf )

Definition at line 44 of file fvmesh.F.

49C-----------------------------------------------
50C M o d u l e s
51C-----------------------------------------------
52 USE fvbag_mod
53C-----------------------------------------------
54C I m p l i c i t T y p e s
55C-----------------------------------------------
56#include "implicit_f.inc"
57C-----------------------------------------------
58C C o m m o n B l o c k s
59C-----------------------------------------------
60#include "com01_c.inc"
61#include "units_c.inc"
62#include "sysunit.inc"
63#include "task_c.inc"
64C-----------------------------------------------
65C D u m m y A r g u m e n t s
66C-----------------------------------------------
67 INTEGER IBUF(*), ELEM(3,*), IVOLU(*), BRIC(8,*),
68 . NEL, NBRIC, TBRIC(13,*), NBA, NELA, NNA, TBA(2,*),
69 . TFACA(12,*), TAGELS(*), IBUFA(*), ELEMA(3,*), TAGELA(*),
70 . IXS(NIXS,*), NNF
71 my_real
72 . x(3,*), xb(3,*), rvolu(*), sfac(6,4,*), dxm
73C-----------------------------------------------
74C L o c a l V a r i a b l e s
75C-----------------------------------------------
76 INTEGER ILVOUT, NLAYER, NFACMAX, NPPMAX, I, J, IEL, N1, N2, N3,
77 . NN1, NN2, NN3, GRBRIC, IAD, NPOLY, NNS, ITAGB(NBRIC),
78 . NVMAX, NG, INFO, ICUT, NVERTS, NF1, NF2, NF3, NF4, NS,
79 . NSMAX, NV, K, KK, NNP, NPOLMAX, NRPMAX, N, M, MM, NPOLH,
80 . NPOLB, ITY, NN, NPHMAX, NPOLHMAX, JJ, JJJ, II,
81 . NNS_OLD, NNP_OLD, NPA, JMAX, IMAX, JJMAX, NP, NNTR,
82 . NPOLY_OLD, IPSURF, IC1, IC2, NHOL, NELP, NPOLH_OLD,
83 . L, LL, IFV, NNS2, NPOLY2, NPL, POLC, NSPOLY, NCPOLY,
84 . NPOLHF, NNB, FILEN, LENP, LENH, IP, JMIN, LENIMAX,
85 . LENRMAX, KKK, NSEG, IMIN, NFAC, N4, NN4, IB, IFOUND,
86 . ITAGBA(NBA), IBSA(NBA), NALL, III,
87 . NNS_ANIM, NBZ, NBNEDGE, NNS3, NPOLY_NEW, NNSP, NEDGE,
88 . ITYZ, INZ, J0, NNS1, K0, I1, I2, IDEP, NSTR, NCTR, IIZ,
89 . NNSA
90 parameter(nvmax=12)
91 my_real
92 . volg, x1, y1, z1, x2, y2, z2, x3, y3, z3, x12, y12, z12,
93 . x13, y13, z13, nrx, nry, nrz, area2, elarea(nel),
94 . norm(3,nel), xbmax, ybmax, zbmax, xbmin, ybmin, zbmin,
95 . xc, yc, zc, xx, yy, zz, pp(3,3), xl7(3), bcenter(3),
96 . bhalfsize(3), xtmax, ytmax, ztmax, xtmin, ytmin, ztmin,
97 . tverts(9), tmptri(3,3), tmpbox(3,8), tmpnorm(3,6), tole,
98 . xg, yg, zg, fv0(3), fv1(3), fv2(3), fu0(3), fu1(3),
99 . fu2(3), quad(3,4), nr(3), area, nx, ny, nz, nnx,
100 . nny, nnz, x0, y0, z0, lmax2, tole2, dd, vm, volumin,
101 . areamax, volph, areap, areael, vpx, vpy, vpz,
102 . v1x, v1y, v1z, v2x, v2y, v2z, nrm1, vx, vy, vz, ss,
103 . normf(3,4), ksi, eta, vx1, vy1, vz1, vx2, vy2, vz2,
104 . vmin, vtmp, x4, y4, z4, x14, y14, z14, norma(3,nela),
105 . dd2, xn(3), zlmin, zlmax, zl, vnorm, vx3, vy3, vz3, lz,
106 . dz, alpha, gamai, cpai, cpbi, cpci, rmwi, pini, ti, cpi,
107 . cvi, rhoi, zl1, zl2, zl3, zlc, val, xxx(3,nnf),
108 . xxxa(3,nna)
109 CHARACTER CHFVB*7, CHMESH*5, FILNAM*116
110C
111 INTEGER, ALLOCATABLE :: FACET(:,:), IPOLY(:,:),
112 . IELNOD(:,:), POLH(:,:), IPOLY_F(:,:),
113 . POLH_F(:,:), IFVNOD(:), IFVNOD_OLD(:),
114 . REDIR_POLY(:), PSEG(:,:), REDIR(:),
115 . PTRI(:,:), REDIR_POLH(:), POLB(:),
116 . IPOLY_OLD(:), POLH_NEW(:,:), ITAGT(:),
117 . TRI(:,:), ADR(:), ADR_OLD(:), IMERGED(:),
118 . IPOLY_F_OLD(:,:), INEDGE(:,:), NREF(:,:),
119 . IZ(:,:), LEDGE(:), IFVNOD2(:,:), ITAGN(:)
120 my_real
121 . , ALLOCATABLE :: normt(:,:), poly(:,:), rpoly(:,:),
122 . rpoly_f(:,:), volu(:), pnodes(:,:),
123 . pholes(:,:), rpoly_old(:), volusort(:),
124 . volu_old(:), rpoly_f_old(:,:),
125 . rfvnod_old(:,:), xnew(:,:), rnedge(:,:),
126 . aref(:,:), rfvnod2(:,:), xns(:,:),
127 . xns2(:,:), xxxsa(:,:)
128C
129 INTEGER FAC(4,6), FACNOR(4,6), FAC4(3,4)
130 DATA fac /1,4,3,2,
131 . 5,6,7,8,
132 . 1,2,6,5,
133 . 2,3,7,6,
134 . 3,4,8,7,
135 . 4,1,5,8/
136 DATA facnor /3,4,5,6,
137 . 3,4,5,6,
138 . 1,4,2,6,
139 . 1,5,2,3,
140 . 1,6,2,4,
141 . 1,3,2,5/
142 DATA fac4 /1,5,3,
143 . 3,5,6,
144 . 6,5,1,
145 . 1,3,6/
146C liste chainee polygones
147 TYPE polygone
148 INTEGER IZ(3), NNSP
149 INTEGER, DIMENSION(:), POINTER :: IPOLY, IELNOD
150 INTEGER, DIMENSION(:,:), POINTER :: NREF
151 my_real
152 . , DIMENSION(:), POINTER :: rpoly
153 my_real
154 . , DIMENSION(:,:), POINTER :: aref
155 TYPE(POLYGONE), POINTER :: PTR
156 END TYPE polygone
157 TYPE(POLYGONE), POINTER :: FIRST, PTR_PREC, PTR_CUR, PTR_TMP,
158 . PTR_OLD, FIRST2, PTR_CUR2, PTR_PREC2,
159 . PTR_TMP2
160C liste chainee polyhedres
161 TYPE polyhedre
162 INTEGER RANK
163 INTEGER, DIMENSION(:), POINTER :: POLH
164 TYPE(POLYHEDRE), POINTER :: PTR
165 END TYPE polyhedre
166 TYPE(POLYHEDRE), POINTER :: PHFIRST, PH_PREC, PH_CUR, PH_TMP
167C
168C Structure de donnees maillage volumes finis
169C - NNS : Nombre de nouveaux noeuds ( = IVOLU(46) )
170C - NNTR : Nombre de nouveaux triangles ( = IVOLU(47) )
171C - NPOLY : Nombre de polygones ( = IVOLU(48) )
172C - NPOLH : Nombre de polyhedres ( volumes finis) ( = IVOLU(49) )
173C
174C - IFVNOD(NNS) : Triangle airbag de reference pour le nouveau noeud
175C - RFVNOD(2,NNS) : Coordonnees locales du noeud dans le triangle
176C
177C - IFVTRI(6,NNTR) : Connectivite nouveaux triangles (1->3)
178C Element de reference si appuye sur surface (4)
179C Polyhedres connectes si interne (5->6)
180C
181C - IFVPOLY(NNTR) : Attribution triangles-polygones
182C - IFVTADR(NPOLY+1) : Adresses dans le tableau precedent
183C
184C - IFVPOLH(NPOLY) : Attribution polygones-polyhedres
185C - IFVPADR(NPOLH+1) : Adresses dans le tableau precedent
186C
187 ifv=ivolu(45)
188 nnsa=fvspmd(ifv)%NNSA
189 ALLOCATE(xxxsa(3,nnsa))
190 IF (nspmd == 1) THEN
191C traitement necessaire pour rester p/on
192 DO i=1,fvspmd(ifv)%NN_L
193 i1=fvspmd(ifv)%IBUF_L(1,i)
194 i2=fvspmd(ifv)%IBUF_L(2,i)
195 xxx(1,i1)=x(1,i2)
196 xxx(2,i1)=x(2,i2)
197 xxx(3,i1)=x(3,i2)
198 ENDDO
199 DO i=1,fvspmd(ifv)%NNA_L
200 i1=fvspmd(ifv)%IBUFA_L(1,i)
201 i2=fvspmd(ifv)%IBUFA_L(2,i)
202 xxxa(1,i1)=x(1,i2)
203 xxxa(2,i1)=x(2,i2)
204 xxxa(3,i1)=x(3,i2)
205 ENDDO
206 DO i=1,fvspmd(ifv)%NNSA_L
207 i1=fvspmd(ifv)%IBUFSA_L(1,i)
208 i2=fvspmd(ifv)%IBUFSA_L(2,i)
209 xxxsa(1,i1)=x(1,i2)
210 xxxsa(2,i1)=x(2,i2)
211 xxxsa(3,i1)=x(3,i2)
212 ENDDO
213 ELSE
214 CALL spmd_fvb_gath(ifv, x, xxx, xxxa, xxxsa,
215 . 2 )
216 IF (ispmd/=fvspmd(ifv)%PMAIN-1) RETURN
217 ENDIF
218 NULLIFY(first)
219 NULLIFY(phfirst)
220C
221 ilvout=ivolu(44)
222C
223 nlayer=ivolu(40)
224 nfacmax=ivolu(41)
225 nppmax=ivolu(42)
226C
227 DO i=1,nbric
228 tbric(7,i)=0
229 DO j=1,6
230 tbric(7+j,i)=0
231 ENDDO
232 ENDDO
233C---------------------------------------------------
234C Normales elementaires et volume de l'airbag
235C---------------------------------------------------
236 volg=zero
237 DO iel=1,nel
238 n1=elem(1,iel)
239 n2=elem(2,iel)
240 n3=elem(3,iel)
241 x1=xxx(1,n1)
242 x2=xxx(1,n2)
243 x3=xxx(1,n3)
244 y1=xxx(2,n1)
245 y2=xxx(2,n2)
246 y3=xxx(2,n3)
247 z1=xxx(3,n1)
248 z2=xxx(3,n2)
249 z3=xxx(3,n3)
250 x12=x2-x1
251 y12=y2-y1
252 z12=z2-z1
253 x13=x3-x1
254 y13=y3-y1
255 z13=z3-z1
256 nrx=y12*z13-z12*y13
257 nry=z12*x13-x12*z13
258 nrz=x12*y13-y12*x13
259 area2=sqrt(nrx**2+nry**2+nrz**2)
260 elarea(iel)=half*area2
261 norm(1,iel)=nrx/area2
262 norm(2,iel)=nry/area2
263 norm(3,iel)=nrz/area2
264C
265 volg=volg+one_over_6*(x1*nrx+y1*nry+z1*nrz)
266 ENDDO
267C
268 DO iel=1,nela
269 n1=elema(1,iel)
270 n2=elema(2,iel)
271 n3=elema(3,iel)
272 IF (tagela(iel)>0) THEN
273 x1=xxx(1,n1)
274 y1=xxx(2,n1)
275 z1=xxx(3,n1)
276 x2=xxx(1,n2)
277 y2=xxx(2,n2)
278 z2=xxx(3,n2)
279 x3=xxx(1,n3)
280 y3=xxx(2,n3)
281 z3=xxx(3,n3)
282 ELSEIF (tagela(iel)<0) THEN
283 x1=xxxa(1,n1)
284 y1=xxxa(2,n1)
285 z1=xxxa(3,n1)
286 x2=xxxa(1,n2)
287 y2=xxxa(2,n2)
288 z2=xxxa(3,n2)
289 x3=xxxa(1,n3)
290 y3=xxxa(2,n3)
291 z3=xxxa(3,n3)
292 ENDIF
293 x12=x2-x1
294 y12=y2-y1
295 z12=z2-z1
296 x13=x3-x1
297 y13=y3-y1
298 z13=z3-z1
299 nrx=y12*z13-z12*y13
300 nry=z12*x13-x12*z13
301 nrz=x12*y13-y12*x13
302 area2=sqrt(nrx**2+nry**2+nrz**2)
303 norma(1,iel)=nrx/area2
304 norma(2,iel)=nry/area2
305 norma(3,iel)=nrz/area2
306 ENDDO
307C---------------------------------------------------
308C Creation de tous les polygones
309C---------------------------------------------------
310 npoly=0
311 nns=0
312 npoly2=0
313 nns2=0
314 IF (nela==0) THEN
315 npoly=0
316 npolh=0
317 GOTO 370
318 ENDIF
319C
320 DO i=1,nbric
321 itagb(i)=0
322 ENDDO
323C 1. Facettes airbag coupees
324 IF (ilvout/=0) THEN
325 WRITE(istdo,*)
326 WRITE(istdo,'(A25,I10,A23)')
327 . ' ** MONITORED VOLUME ID: ',ivolu(1),' - BUILDING POLYGONS **'
328 ENDIF
329C
330C nullify sinon test sur associated vrai second appel fvmesh1
331 NULLIFY(first2)
332 DO i=1,nbric
333C IF (ILVOUT<0) CALL PROGBAR_C(I,NBRIC)
334C
335 ptr_old => ptr_cur
336 npoly_old=npoly
337 nns_old=nns
338 100 CONTINUE
339 IF (ALLOCATED(facet)) DEALLOCATE(facet)
340 IF (ALLOCATED(tri)) DEALLOCATE(tri)
341 IF (ALLOCATED(normt)) DEALLOCATE(normt)
342 IF (ALLOCATED(poly)) DEALLOCATE(poly)
343 ALLOCATE(facet(6,1+nfacmax), tri(nfacmax,4),
344 . normt(3,nfacmax), poly(3,nvmax))
345 DO j=1,6
346 facet(j,1)=0
347 ENDDO
348C
349 xbmax=-ep30
350 ybmax=-ep30
351 zbmax=-ep30
352 xbmin=ep30
353 ybmin=ep30
354 zbmin=ep30
355 xc=zero
356 yc=zero
357 zc=zero
358 DO j=1,8
359 xx=xb(1,bric(j,i))
360 yy=xb(2,bric(j,i))
361 zz=xb(3,bric(j,i))
362 xbmax=max(xbmax,xx)
363 ybmax=max(ybmax,yy)
364 zbmax=max(zbmax,zz)
365 xbmin=min(xbmin,xx)
366 ybmin=min(ybmin,yy)
367 zbmin=min(zbmin,zz)
368 xc=xc+one_over_8*xx
369 yc=yc+one_over_8*yy
370 zc=zc+one_over_8*zz
371 ENDDO
372C Passage repere local brique
373 pp(1,1)=sfac(4,2,i)
374 pp(2,1)=sfac(4,3,i)
375 pp(3,1)=sfac(4,4,i)
376 pp(1,2)=sfac(5,2,i)
377 pp(2,2)=sfac(5,3,i)
378 pp(3,2)=sfac(5,4,i)
379 pp(1,3)=sfac(2,2,i)
380 pp(2,3)=sfac(2,3,i)
381 pp(3,3)=sfac(2,4,i)
382C
383 xx=xb(1,bric(7,i))
384 yy=xb(2,bric(7,i))
385 zz=xb(3,bric(7,i))
386 xl7(1)=pp(1,1)*(xx-xc)+pp(2,1)*(yy-yc)+pp(3,1)*(zz-zc)
387 xl7(2)=pp(1,2)*(xx-xc)+pp(2,2)*(yy-yc)+pp(3,2)*(zz-zc)
388 xl7(3)=pp(1,3)*(xx-xc)+pp(2,3)*(yy-yc)+pp(3,3)*(zz-zc)
389 bcenter(1)=zero
390 bcenter(2)=zero
391 bcenter(3)=zero
392 bhalfsize(1)=xl7(1)
393 bhalfsize(2)=xl7(2)
394 bhalfsize(3)=xl7(3)
395 info=0
396 DO iel=1,nela
397 n1=elema(1,iel)
398 n2=elema(2,iel)
399 n3=elema(3,iel)
400 IF (tagela(iel)>0) THEN
401 x1=xxx(1,n1)
402 y1=xxx(2,n1)
403 z1=xxx(3,n1)
404 x2=xxx(1,n2)
405 y2=xxx(2,n2)
406 z2=xxx(3,n2)
407 x3=xxx(1,n3)
408 y3=xxx(2,n3)
409 z3=xxx(3,n3)
410 ELSEIF (tagela(iel)<0) THEN
411 x1=xxxa(1,n1)
412 y1=xxxa(2,n1)
413 z1=xxxa(3,n1)
414 x2=xxxa(1,n2)
415 y2=xxxa(2,n2)
416 z2=xxxa(3,n2)
417 x3=xxxa(1,n3)
418 y3=xxxa(2,n3)
419 z3=xxxa(3,n3)
420 ENDIF
421 xtmax=max(x1,x2,x3)
422 ytmax=max(y1,y2,y3)
423 ztmax=max(z1,z2,z3)
424 xtmin=min(x1,x2,x3)
425 ytmin=min(y1,y2,y3)
426 ztmin=min(z1,z2,z3)
427 IF (xbmax<xtmin.OR.ybmax<ytmin.OR.zbmax<ztmin.OR.
428 . xbmin>xtmax.OR.ybmin>ytmax.OR.zbmin>ztmax)
429 . cycle
430C
431 tverts(1)=pp(1,1)*(x1-xc)+pp(2,1)*(y1-yc)+pp(3,1)*(z1-zc)
432 tverts(2)=pp(1,2)*(x1-xc)+pp(2,2)*(y1-yc)+pp(3,2)*(z1-zc)
433 tverts(3)=pp(1,3)*(x1-xc)+pp(2,3)*(y1-yc)+pp(3,3)*(z1-zc)
434 tverts(4)=pp(1,1)*(x2-xc)+pp(2,1)*(y2-yc)+pp(3,1)*(z2-zc)
435 tverts(5)=pp(1,2)*(x2-xc)+pp(2,2)*(y2-yc)+pp(3,2)*(z2-zc)
436 tverts(6)=pp(1,3)*(x2-xc)+pp(2,3)*(y2-yc)+pp(3,3)*(z2-zc)
437 tverts(7)=pp(1,1)*(x3-xc)+pp(2,1)*(y3-yc)+pp(3,1)*(z3-zc)
438 tverts(8)=pp(1,2)*(x3-xc)+pp(2,2)*(y3-yc)+pp(3,2)*(z3-zc)
439 tverts(9)=pp(1,3)*(x3-xc)+pp(2,3)*(y3-yc)+pp(3,3)*(z3-zc)
440C
441 CALL tribox3(icut, bcenter, bhalfsize, tverts)
442C Liste elargie pour faces laterales
443 IF (icut==0) THEN
444 tole=em5
445* TOLE=ZERO
446 xg=third*(x1+x2+x3)
447 yg=third*(y1+y2+y3)
448 zg=third*(z1+z2+z3)
449 x1=xg+(one+tole)*(x1-xg)
450 y1=yg+(one+tole)*(y1-yg)
451 z1=zg+(one+tole)*(z1-zg)
452 x2=xg+(one+tole)*(x2-xg)
453 y2=yg+(one+tole)*(y2-yg)
454 z2=zg+(one+tole)*(z2-zg)
455 x3=xg+(one+tole)*(x3-xg)
456 y3=yg+(one+tole)*(y3-yg)
457 z3=zg+(one+tole)*(z3-zg)
458C
459 tverts(1)=pp(1,1)*(x1-xc)+pp(2,1)*(y1-yc)+pp(3,1)*(z1-zc)
460 tverts(2)=pp(1,2)*(x1-xc)+pp(2,2)*(y1-yc)+pp(3,2)*(z1-zc)
461 tverts(3)=pp(1,3)*(x1-xc)+pp(2,3)*(y1-yc)+pp(3,3)*(z1-zc)
462 tverts(4)=pp(1,1)*(x2-xc)+pp(2,1)*(y2-yc)+pp(3,1)*(z2-zc)
463 tverts(5)=pp(1,2)*(x2-xc)+pp(2,2)*(y2-yc)+pp(3,2)*(z2-zc)
464 tverts(6)=pp(1,3)*(x2-xc)+pp(2,3)*(y2-yc)+pp(3,3)*(z2-zc)
465 tverts(7)=pp(1,1)*(x3-xc)+pp(2,1)*(y3-yc)+pp(3,1)*(z3-zc)
466 tverts(8)=pp(1,2)*(x3-xc)+pp(2,2)*(y3-yc)+pp(3,2)*(z3-zc)
467 tverts(9)=pp(1,3)*(x3-xc)+pp(2,3)*(y3-yc)+pp(3,3)*(z3-zc)
468C
469 CALL tribox3(icut, bcenter, bhalfsize, tverts)
470 IF (icut==1) icut=2
471 ENDIF
472C
473 IF (icut>=1) THEN
474 IF (icut==1) THEN
475 tbric(7,i)=2
476 tmptri(1,1)=x1
477 tmptri(2,1)=y1
478 tmptri(3,1)=z1
479 tmptri(1,2)=x2
480 tmptri(2,2)=y2
481 tmptri(3,2)=z2
482 tmptri(1,3)=x3
483 tmptri(2,3)=y3
484 tmptri(3,3)=z3
485 DO j=1,8
486 tmpbox(1,j)=xb(1,bric(j,i))
487 tmpbox(2,j)=xb(2,bric(j,i))
488 tmpbox(3,j)=xb(3,bric(j,i))
489 ENDDO
490 DO j=1,6
491 tmpnorm(1,j)=-sfac(j,2,i)
492 tmpnorm(2,j)=-sfac(j,3,i)
493 tmpnorm(3,j)=-sfac(j,4,i)
494 ENDDO
495 CALL itribox(tmptri, tmpbox, tmpnorm, nverts, poly,
496 . nvmax )
497 IF (nverts>0) THEN
498 npoly=npoly+1
499 IF (.NOT.ASSOCIATED(first)) THEN
500 ALLOCATE(first)
501 ptr_cur => first
502 ELSE
503 ALLOCATE(ptr_cur)
504 ptr_prec%PTR => ptr_cur
505 ENDIF
506 ALLOCATE(ptr_cur%IPOLY(6+nverts),
507 . ptr_cur%IELNOD(nverts),
508 . ptr_cur%RPOLY(4+3*nverts))
509 ptr_cur%IPOLY(1)=1
510 ptr_cur%IPOLY(2)=nverts
511 ptr_cur%IPOLY(3)=iel
512 ptr_cur%IPOLY(4)=i
513 ptr_cur%IPOLY(5)=0
514 ptr_cur%IPOLY(6)=0
515 ptr_cur%RPOLY(1)=zero
516 ptr_cur%RPOLY(2)=norma(1,iel)
517 ptr_cur%RPOLY(3)=norma(2,iel)
518 ptr_cur%RPOLY(4)=norma(3,iel)
519 DO j=1,nverts
520 nns=nns+1
521 ptr_cur%IPOLY(6+j)=nns
522 ptr_cur%IELNOD(j)=iel
523 ptr_cur%RPOLY(4+3*(j-1)+1)=poly(1,j)
524 ptr_cur%RPOLY(4+3*(j-1)+2)=poly(2,j)
525 ptr_cur%RPOLY(4+3*(j-1)+3)=poly(3,j)
526 ENDDO
527 ptr_prec => ptr_cur
528 ENDIF
529 ENDIF
530C Facettes coupees (1 quadrangle=2 triangles)
531 tole=em5
532* TOLE=ZERO
533 xg=third*(x1+x2+x3)
534 yg=third*(y1+y2+y3)
535 zg=third*(z1+z2+z3)
536 x1=xg+(one+tole)*(x1-xg)
537 y1=yg+(one+tole)*(y1-yg)
538 z1=zg+(one+tole)*(z1-zg)
539 x2=xg+(one+tole)*(x2-xg)
540 y2=yg+(one+tole)*(y2-yg)
541 z2=zg+(one+tole)*(z2-zg)
542 x3=xg+(one+tole)*(x3-xg)
543 y3=yg+(one+tole)*(y3-yg)
544 z3=zg+(one+tole)*(z3-zg)
545C
546 fv0(1)=x1
547 fv0(2)=y1
548 fv0(3)=z1
549 fv1(1)=x2
550 fv1(2)=y2
551 fv1(3)=z2
552 fv2(1)=x3
553 fv2(2)=y3
554 fv2(3)=z3
555 DO j=1,6
556 nf1=bric(fac(1,j),i)
557 nf2=bric(fac(2,j),i)
558 nf3=bric(fac(3,j),i)
559 nf4=bric(fac(4,j),i)
560 fu0(1)=xb(1,nf1)
561 fu0(2)=xb(2,nf1)
562 fu0(3)=xb(3,nf1)
563 fu1(1)=xb(1,nf2)
564 fu1(2)=xb(2,nf2)
565 fu1(3)=xb(3,nf2)
566 fu2(1)=xb(1,nf3)
567 fu2(2)=xb(2,nf3)
568 fu2(3)=xb(3,nf3)
569C
570 CALL tritri3(icut, fv0, fv1, fv2, fu0, fu1, fu2)
571 IF (icut==1) THEN
572 tbric(7+j,i)=2
573 ns=facet(j,1)
574 ns=ns+1
575 facet(j,1)=ns
576 IF (ns>nfacmax) THEN
577 info=1
578 ELSE
579 facet(j,1+ns)=iel
580 ENDIF
581 cycle
582 ENDIF
583C
584 fu1(1)=xb(1,nf3)
585 fu1(2)=xb(2,nf3)
586 fu1(3)=xb(3,nf3)
587 fu2(1)=xb(1,nf4)
588 fu2(2)=xb(2,nf4)
589 fu2(3)=xb(3,nf4)
590C
591 CALL tritri3(icut, fv0, fv1, fv2, fu0, fu1, fu2)
592 IF (icut==1) THEN
593 tbric(7+j,i)=2
594 ns=facet(j,1)
595 ns=ns+1
596 facet(j,1)=ns
597 IF (ns>nfacmax) THEN
598 info=1
599 ELSE
600 facet(j,1+ns)=iel
601 ENDIF
602 ENDIF
603 ENDDO
604 ENDIF
605 ENDDO
606C
607 IF (info==1) THEN
608 nsmax=0
609 DO j=1,6
610 nsmax=max(nsmax,facet(j,1))
611 ENDDO
612 nfacmax=nsmax
613C Destruction des polygones crees pour rien
614 IF (npoly_old>0) THEN
615 ptr_cur => ptr_old%PTR
616 ELSE
617 ptr_cur => first
618 ENDIF
619 DO j=1,npoly-npoly_old
620 ptr_tmp => ptr_cur%PTR
621 DEALLOCATE(ptr_cur%IPOLY, ptr_cur%RPOLY,
622 . ptr_cur%IELNOD)
623 DEALLOCATE(ptr_cur)
624 ptr_cur => ptr_tmp
625 ENDDO
626 IF (npoly_old>0) THEN
627 ptr_prec => ptr_old
628 ELSE
629 NULLIFY(first)
630 ENDIF
631 npoly=npoly_old
632 nns=nns_old
633C
634 GOTO 100
635 ENDIF
636C 2. Faces laterales
637 IF (tbric(7,i)<=1) THEN
638 DEALLOCATE(facet, tri, normt, poly)
639 cycle
640 END IF
641 DO j=1,6
642 nv=tbric(j,i)
643 IF (nv==0) cycle
644 IF (itagb(nv)==1) cycle
645 IF (tbric(7+j,i)==2) THEN
646 DO k=1,4
647 kk=fac(k,j)
648 quad(1,k)=xb(1,bric(kk,i))
649 quad(2,k)=xb(2,bric(kk,i))
650 quad(3,k)=xb(3,bric(kk,i))
651 ENDDO
652 ns=facet(j,1)
653 DO k=1,ns
654 iel=facet(j,1+k)
655 IF (tagela(iel)>0) THEN
656 DO l=1,3
657 ll=elema(l,iel)
658 tri(k,l)=ll
659 ENDDO
660 ELSEIF (tagela(iel)<0) THEN
661 DO l=1,3
662 ll=-elema(l,iel)
663 tri(k,l)=ll
664 ENDDO
665 ENDIF
666 tri(k,4)=iel
667 normt(1,k)=norma(1,iel)
668 normt(2,k)=norma(2,iel)
669 normt(3,k)=norma(3,iel)
670 ENDDO
671 DO k=1,4
672 normf(1,k)=sfac(facnor(k,j),2,i)
673 normf(2,k)=sfac(facnor(k,j),3,i)
674 normf(3,k)=sfac(facnor(k,j),4,i)
675 ENDDO
676 nr(1)=sfac(j,2,i)
677 nr(2)=sfac(j,3,i)
678 nr(3)=sfac(j,4,i)
679C
680 nnp=0
681 npolmax=nlayer
682C
683 200 CONTINUE
684 nrpmax=nppmax*3+4
685 IF (ALLOCATED(ipoly)) DEALLOCATE(ipoly)
686 IF (ALLOCATED(rpoly)) DEALLOCATE(rpoly)
687 IF (ALLOCATED(ielnod)) DEALLOCATE(ielnod)
688 ALLOCATE(ipoly(6+nppmax+1+npolmax,npolmax),
689 . rpoly(nrpmax+3*npolmax,npolmax),
690 . ielnod(nppmax,npolmax))
691C
692 CALL facepoly(quad, tri, ns, ipoly, rpoly,
693 . nr, normf, normt, nfacmax, nnp,
694 . nrpmax, i, nv, dxm , npolmax,
695 . nppmax, info, ielnod, xxx, elema,
696 . ibuf, nela, i, j, ivolu(1),
697 . ilvout, ibufa, tagela, xxxa )
698 IF (info==1) GOTO 200
699C
700 DO n=1,nnp
701 IF (ipoly(1,n)==-1) cycle
702 npoly2=npoly2+1
703 IF (.NOT.ASSOCIATED(first2)) THEN
704 ALLOCATE(first2)
705 ptr_cur2 => first2
706 ELSE
707 ALLOCATE(ptr_cur2)
708 ptr_prec2%PTR => ptr_cur2
709 ENDIF
710 nn=ipoly(2,n)
711 nhol=ipoly(6+nn+1,n)
712 ALLOCATE(ptr_cur2%IPOLY(6+nn+1+nhol),
713 . ptr_cur2%IELNOD(nn),
714 . ptr_cur2%RPOLY(4+3*nn+3*nhol))
715 DO m=1,6
716 ptr_cur2%IPOLY(m)=ipoly(m,n)
717 ENDDO
718 DO m=1,4
719 ptr_cur2%RPOLY(m)=rpoly(m,n)
720 ENDDO
721 DO m=1,ipoly(2,n)
722 nns2=nns2+1
723 ptr_cur2%IPOLY(6+m)=nns2
724 mm=ielnod(m,n)
725 ptr_cur2%IELNOD(m)=facet(j,1+mm)
726 ptr_cur2%RPOLY(4+3*(m-1)+1)=rpoly(4+3*(m-1)+1,n)
727 ptr_cur2%RPOLY(4+3*(m-1)+2)=rpoly(4+3*(m-1)+2,n)
728 ptr_cur2%RPOLY(4+3*(m-1)+3)=rpoly(4+3*(m-1)+3,n)
729 ENDDO
730 ptr_cur2%IPOLY(6+nn+1)=nhol
731 DO m=1,nhol
732 ptr_cur2%IPOLY(6+nn+1+m)=ipoly(6+nn+1+m,n)
733 ptr_cur2%RPOLY(4+3*nn+3*(m-1)+1)=
734 . rpoly(4+3*nn+3*(m-1)+1,n)
735 ptr_cur2%RPOLY(4+3*nn+3*(m-1)+2)=
736 . rpoly(4+3*nn+3*(m-1)+2,n)
737 ptr_cur2%RPOLY(4+3*nn+3*(m-1)+3)=
738 . rpoly(4+3*nn+3*(m-1)+3,n)
739 ENDDO
740 ptr_prec2 => ptr_cur2
741 ENDDO
742C
743 DEALLOCATE(ipoly, rpoly, ielnod)
744 ENDIF
745 ENDDO
746 itagb(i)=1
747C
748 DEALLOCATE(facet, tri, normt, poly)
749 ENDDO
750C Concatenation des deux listes
751 ptr_cur2 => first2
752 ptr_prec => ptr_cur
753 DO i=1,npoly2
754 npoly=npoly+1
755 ALLOCATE(ptr_cur)
756 ptr_prec%PTR => ptr_cur
757 nn=ptr_cur2%IPOLY(2)
758 nhol=ptr_cur2%IPOLY(6+nn+1)
759 ALLOCATE(ptr_cur%IPOLY(6+nn+1+nhol),
760 . ptr_cur%IELNOD(nn),
761 . ptr_cur%RPOLY(4+3*nn+3*nhol))
762 DO j=1,6
763 ptr_cur%IPOLY(j)=ptr_cur2%IPOLY(j)
764 ENDDO
765 DO j=1,4
766 ptr_cur%RPOLY(j)=ptr_cur2%RPOLY(j)
767 ENDDO
768 DO j=1,nn
769 ptr_cur%IPOLY(6+j)=nns+ptr_cur2%IPOLY(6+j)
770 ptr_cur%IELNOD(j)=ptr_cur2%IELNOD(j)
771 ptr_cur%RPOLY(4+3*(j-1)+1)=ptr_cur2%RPOLY(4+3*(j-1)+1)
772 ptr_cur%RPOLY(4+3*(j-1)+2)=ptr_cur2%RPOLY(4+3*(j-1)+2)
773 ptr_cur%RPOLY(4+3*(j-1)+3)=ptr_cur2%RPOLY(4+3*(j-1)+3)
774 ENDDO
775 ptr_cur%IPOLY(6+nn+1)=nhol
776 DO j=1,nhol
777 ptr_cur%IPOLY(6+nn+1+j)=ptr_cur2%IPOLY(6+nn+1+j)
778 ptr_cur%RPOLY(4+3*nn+3*(j-1)+1)=
779 . ptr_cur2%RPOLY(4+3*nn+3*(j-1)+1)
780 ptr_cur%RPOLY(4+3*nn+3*(j-1)+2)=
781 . ptr_cur2%RPOLY(4+3*nn+3*(j-1)+2)
782 ptr_cur%RPOLY(4+3*nn+3*(j-1)+3)=
783 . ptr_cur2%RPOLY(4+3*nn+3*(j-1)+3)
784 ENDDO
785 ptr_tmp2 => ptr_cur2%PTR
786 DEALLOCATE(ptr_cur2%IPOLY, ptr_cur2%RPOLY, ptr_cur2%IELNOD)
787 DEALLOCATE(ptr_cur2)
788 ptr_cur2 => ptr_tmp2
789 ptr_prec => ptr_cur
790 ENDDO
791C
792 nns=nns+nns2
793C---------------------------------------------------
794C Decoupage des polygones selon la verticale locale
795C---------------------------------------------------
796 nbz=ivolu(65)
797 IF (nbz>1) THEN
798 ptr_cur => first
799 DO i=1,npoly
800 ptr_cur%NNSP=0
801 ptr_cur => ptr_cur%PTR
802 ENDDO
803C
804 vx3=rvolu(35)
805 vy3=rvolu(36)
806 vz3=rvolu(37)
807 vx1=rvolu(38)
808 vy1=rvolu(39)
809 vz1=rvolu(40)
810C Repere local
811 vnorm=sqrt(vx3**2+vy3**2+vz3**2)
812 vx3=vx3/vnorm
813 vy3=vy3/vnorm
814 vz3=vz3/vnorm
815 ss=vx3*vx1+vy3*vy1+vz3*vz1
816 vx1=vx1-ss*vx3
817 vy1=vy1-ss*vy3
818 vz1=vz1-ss*vz3
819 vnorm=sqrt(vx1**2+vy1**2+vz1**2)
820 vx1=vx1/vnorm
821 vy1=vy1/vnorm
822 vz1=vz1/vnorm
823 vx2=vy3*vz1-vz3*vy1
824 vy2=vz3*vx1-vx3*vz1
825 vz2=vx3*vy1-vy3*vx1
826C
827 x0=rvolu(41)
828 y0=rvolu(42)
829 z0=rvolu(43)
830 lz=rvolu(53)
831 dz=two*lz/nbz
832 dxm=min(dxm,dz)
833C
834 zbmin=ep30
835 zbmax=-ep30
836 DO iel=1,nel
837 n1=elem(1,iel)
838 n2=elem(2,iel)
839 n3=elem(3,iel)
840 x1=xxx(1,n1)
841 x2=xxx(1,n2)
842 x3=xxx(1,n3)
843 y1=xxx(2,n1)
844 y2=xxx(2,n2)
845 y3=xxx(2,n3)
846 z1=xxx(3,n1)
847 z2=xxx(3,n2)
848 z3=xxx(3,n3)
849 x12=x2-x1
850 y12=y2-y1
851 z12=z2-z1
852 x13=x3-x1
853 y13=y3-y1
854 z13=z3-z1
855 zl1=(x1-x0)*vx3+(y1-y0)*vy3+(z1-z0)*vz3
856 zl2=(x2-x0)*vx3+(y2-y0)*vy3+(z2-z0)*vz3
857 zl3=(x3-x0)*vx3+(y3-y0)*vy3+(z3-z0)*vz3
858 zbmin=min(zbmin,zl1)
859 zbmin=min(zbmin,zl2)
860 zbmin=min(zbmin,zl3)
861 zbmax=max(zbmax,zl1)
862 zbmax=max(zbmax,zl2)
863 zbmax=max(zbmax,zl3)
864 ENDDO
865C Clipping des polygones
866 x0=x0-lz*vx3
867 y0=y0-lz*vy3
868 z0=z0-lz*vz3
869 zlc=-lz
870 ALLOCATE(inedge(6,npoly*nppmax), rnedge(6,npoly*nppmax))
871 nbnedge=0
872 nns3=0
873 npoly_new=npoly
874 iiz=0
875 DO i=1,nbz+1
876 IF (zlc<zbmin.OR.zlc>zbmax) THEN
877 x0=x0+dz*vx3
878 y0=y0+dz*vy3
879 z0=z0+dz*vz3
880 zlc=zlc+dz
881 cycle
882 ENDIF
883 iiz=iiz+1
884 ptr_cur => first
885 DO j=1,npoly
886 zlmin=ep30
887 zlmax=-ep30
888 DO k=1,ptr_cur%IPOLY(2)
889 xx=ptr_cur%RPOLY(4+3*(k-1)+1)
890 yy=ptr_cur%RPOLY(4+3*(k-1)+2)
891 zz=ptr_cur%RPOLY(4+3*(k-1)+3)
892 zl=(xx-x0)*vx3+(yy-y0)*vy3+(zz-z0)*vz3
893 zlmin=min(zlmin,zl)
894 zlmax=max(zlmax,zl)
895 ENDDO
896 IF (zlmin*zlmax<zero) THEN
897 ity=ptr_cur%IPOLY(1)
898 nn=ptr_cur%IPOLY(2)
899 nhol=0
900 IF (ity==2) nhol=ptr_cur%IPOLY(6+nn+1)
901 ALLOCATE(ipoly(6+2*nn+1+nhol,nn),
902 . rpoly(4+3*2*nn+3*nhol,nn),
903 . ielnod(2*nn,nn),
904 . nref(2,nn), aref(4,nn),
905 . iz(3,nn))
906 CALL gpolcut(
907 . ipoly, rpoly, ptr_cur%IPOLY, ptr_cur%RPOLY, inedge,
908 . rnedge, nbnedge, vx3, vy3, vz3,
909 . x0, y0, z0, nref, aref,
910 . nn, nhol, iiz, iz, nns3,
911 . nnp , nnsp, ielnod, ptr_cur%IELNOD)
912C
913 ptr_tmp => ptr_cur%PTR
914 npoly_new=npoly_new-1
915 DO n=1,nnp
916 npoly_new=npoly_new+1
917 IF (n==1) THEN
918 DEALLOCATE(ptr_cur%IPOLY, ptr_cur%RPOLY,
919 . ptr_cur%IELNOD)
920 ELSE
921 ALLOCATE(ptr_cur)
922 ptr_prec%PTR => ptr_cur
923 ptr_cur%NNSP=0
924 ENDIF
925 nn=ipoly(2,n)
926 nhol=ipoly(6+nn+1,n)
927 ALLOCATE(ptr_cur%IPOLY(6+nn+1+nhol),
928 . ptr_cur%IELNOD(nn),
929 . ptr_cur%RPOLY(4+3*nn+3*nhol))
930 DO m=1,6
931 ptr_cur%IPOLY(m)=ipoly(m,n)
932 ENDDO
933 DO m=1,4
934 ptr_cur%RPOLY(m)=rpoly(m,n)
935 ENDDO
936 DO m=1,nn
937 mm=ipoly(6+m,n)
938 ptr_cur%IPOLY(6+m)=mm
939 ptr_cur%IELNOD(m)=ielnod(m,n)
940 ptr_cur%RPOLY(4+3*(m-1)+1)=rpoly(4+3*(m-1)+1,n)
941 ptr_cur%RPOLY(4+3*(m-1)+2)=rpoly(4+3*(m-1)+2,n)
942 ptr_cur%RPOLY(4+3*(m-1)+3)=rpoly(4+3*(m-1)+3,n)
943 ENDDO
944 nhol=ipoly(6+nn+1,n)
945 ptr_cur%IPOLY(6+nn+1)=nhol
946 DO m=1,nhol
947 ptr_cur%IPOLY(6+nn+1+m)=ipoly(6+nn+1+m,n)
948 ptr_cur%RPOLY(4+3*nn+3*(m-1)+1)=
949 . rpoly(4+3*nn+3*(m-1)+1,n)
950 ptr_cur%RPOLY(4+3*nn+3*(m-1)+2)=
951 . rpoly(4+3*nn+3*(m-1)+2,n)
952 ptr_cur%RPOLY(4+3*nn+3*(m-1)+3)=
953 . rpoly(4+3*nn+3*(m-1)+3,n)
954 ENDDO
955 ptr_cur%IZ(1)=1
956 ptr_cur%IZ(2)=iz(2,n)
957 IF (n==nnp) THEN
958 ptr_cur%NNSP=nnsp
959 ALLOCATE(ptr_cur%NREF(3,nnsp),
960 . ptr_cur%AREF(4,nnsp))
961 DO m=1,nnsp
962 nns3=nns3+1
963 ptr_cur%NREF(1,m)=nns3
964 ptr_cur%NREF(2,m)=nref(1,m)
965 ptr_cur%NREF(3,m)=nref(2,m)
966 ptr_cur%AREF(1,m)=aref(1,m)
967 ptr_cur%AREF(2,m)=aref(2,m)
968 ptr_cur%AREF(3,m)=aref(3,m)
969 ptr_cur%AREF(4,m)=aref(4,m)
970 ENDDO
971 ENDIF
972 ptr_prec => ptr_cur
973 IF (n==nnp) ptr_cur%PTR => ptr_tmp
974 ENDDO
975C
976 DEALLOCATE(ipoly, rpoly, ielnod, nref, aref, iz)
977 ELSE
978 IF (zlmin==zero) THEN
979 nn=ptr_cur%IPOLY(2)
980 ALLOCATE(nref(2,2*nn), aref(4,2*nn))
981 CALL horiedge(
982 . ptr_cur%IPOLY, ptr_cur%RPOLY, vx3, vy3, vz3,
983 . nbnedge, inedge, rnedge, x0, y0,
984 . z0, iiz , nns3, nref, aref,
985 . nnsp )
986C
987 IF (nnsp>0) THEN
988 ALLOCATE(ptr_cur%NREF(3,nnsp),
989 . ptr_cur%AREF(4,nnsp))
990 ptr_cur%NNSP=nnsp
991 DO n=1,nnsp
992 nns3=nns3+1
993 ptr_cur%NREF(1,n)=nns3
994 ptr_cur%NREF(2,n)=nref(1,n)
995 ptr_cur%NREF(3,n)=nref(2,n)
996 ptr_cur%AREF(1,n)=aref(1,n)
997 ptr_cur%AREF(2,n)=aref(2,n)
998 ptr_cur%AREF(3,n)=aref(3,n)
999 ptr_cur%AREF(4,n)=aref(4,n)
1000 ENDDO
1001 ENDIF
1002 DEALLOCATE(nref, aref)
1003 ENDIF
1004 IF (zlmin>=zero) THEN
1005 ptr_cur%IZ(1)=1
1006 ptr_cur%IZ(2)=iiz+1
1007 ELSEIF (iiz==1.AND.zlmax<=zero) THEN
1008 ptr_cur%IZ(1)=1
1009 ptr_cur%IZ(2)=1
1010 ENDIF
1011 ptr_prec => ptr_cur
1012 ENDIF
1013 ptr_cur => ptr_cur%PTR
1014 ENDDO
1015 npoly=npoly_new
1016 x0=x0+dz*vx3
1017 y0=y0+dz*vy3
1018 z0=z0+dz*vz3
1019 zlc=zlc+dz
1020 ENDDO
1021C
1022 ALLOCATE(redir(nns))
1023 ptr_cur => first
1024 nns=0
1025 DO i=1,npoly
1026 nn=ptr_cur%IPOLY(2)
1027 DO j=1,nn
1028 jj=ptr_cur%IPOLY(6+j)
1029 IF (jj>0) THEN
1030 nns=nns+1
1031 redir(jj)=nns
1032 ENDIF
1033 ENDDO
1034 ptr_cur => ptr_cur%PTR
1035 ENDDO
1036C
1037 ptr_cur => first
1038 DO i=1,npoly
1039 nnsp=ptr_cur%NNSP
1040 IF (nnsp>0) THEN
1041 DO j=1,nnsp
1042 n1=ptr_cur%NREF(2,j)
1043 n2=ptr_cur%NREF(3,j)
1044 IF (n1>0) ptr_cur%NREF(2,j)=redir(n1)
1045 IF (n2>0) ptr_cur%NREF(3,j)=redir(n2)
1046 ENDDO
1047 ENDIF
1048 ptr_cur => ptr_cur%PTR
1049 ENDDO
1050 DEALLOCATE(redir)
1051C Creation des nouveaux polygones horizontaux
1052 ptr_cur => first
1053 DO i=1,npoly
1054 ptr_prec => ptr_cur
1055 ptr_cur => ptr_cur%PTR
1056 ENDDO
1057 ALLOCATE(ledge(nbnedge))
1058 DO i=1,iiz
1059 DO j=1,nbric
1060 nedge=0
1061 DO k=1,nbnedge
1062 IF (inedge(6,k)/=i) cycle
1063 IF (inedge(1,k)==1.AND.inedge(5,k)/=j) cycle
1064 IF (inedge(1,k)==2.AND.inedge(4,k)/=j.AND.
1065 . inedge(5,k)/=j) cycle
1066 nedge=nedge+1
1067 ledge(nedge)=k
1068 ENDDO
1069 IF (nedge==0) cycle
1070 ALLOCATE(ipoly(6+2*nedge+1+nedge,nedge),
1071 . rpoly(4+6*nedge+3*nedge,nedge),
1072 . iz(3,nedge), ielnod(nedge,nedge))
1073C
1074 CALL horipoly(inedge, rnedge, ledge, nedge, ipoly,
1075 . rpoly, iz, ielnod, nnp, vx3,
1076 . vy3, vz3, i , j )
1077C
1078 DO n=1,nnp
1079 IF (ipoly(1,n)==-1) cycle
1080 npoly=npoly+1
1081 ALLOCATE(ptr_cur)
1082 ptr_prec%PTR => ptr_cur
1083 nn=ipoly(2,n)
1084 nhol=ipoly(6+nn+1,n)
1085 ALLOCATE(ptr_cur%IPOLY(6+nn+1+nhol),
1086 . ptr_cur%RPOLY(4+3*nn+3*nhol),
1087 . ptr_cur%IELNOD(nn))
1088 ptr_cur%NNSP=0
1089 DO m=1,6
1090 ptr_cur%IPOLY(m)=ipoly(m,n)
1091 ENDDO
1092 DO m=1,4
1093 ptr_cur%RPOLY(m)=rpoly(m,n)
1094 ENDDO
1095 DO m=1,nn
1096 ptr_cur%IPOLY(6+m)=ipoly(6+m,n)
1097 ptr_cur%IELNOD(m)=ielnod(m,n)
1098 ptr_cur%RPOLY(4+3*(m-1)+1)=rpoly(4+3*(m-1)+1,n)
1099 ptr_cur%RPOLY(4+3*(m-1)+2)=rpoly(4+3*(m-1)+2,n)
1100 ptr_cur%RPOLY(4+3*(m-1)+3)=rpoly(4+3*(m-1)+3,n)
1101 ENDDO
1102 nhol=ipoly(6+nn+1,n)
1103 ptr_cur%IPOLY(6+nn+1)=nhol
1104 DO m=1,nhol
1105 ptr_cur%IPOLY(6+nn+1+m)=ipoly(6+nn+1+m,n)
1106 ptr_cur%RPOLY(4+3*nn+3*(m-1)+1)=
1107 . rpoly(4+3*nn+3*(m-1)+1,n)
1108 ptr_cur%RPOLY(4+3*nn+3*(m-1)+2)=
1109 . rpoly(4+3*nn+3*(m-1)+2,n)
1110 ptr_cur%RPOLY(4+3*nn+3*(m-1)+3)=
1111 . rpoly(4+3*nn+3*(m-1)+3,n)
1112 ENDDO
1113 DO m=1,3
1114 ptr_cur%IZ(m)=iz(m,n)
1115 ENDDO
1116 ptr_prec => ptr_cur
1117 ENDDO
1118C
1119 DEALLOCATE(ipoly, rpoly, iz, ielnod)
1120 ENDDO
1121 ENDDO
1122 DEALLOCATE(ledge)
1123 ELSE
1124 ptr_cur => first
1125 DO i=1,npoly
1126 ptr_cur%NNSP=0
1127 ptr_cur%IZ(1)=1
1128 ptr_cur%IZ(2)=1
1129 ptr_cur => ptr_cur%PTR
1130 iiz=0
1131 ENDDO
1132 ENDIF
1133C---------------------------------------------------
1134C Construction des polyhedres
1135C---------------------------------------------------
1136 npolh=0
1137 IF (ilvout/=0) WRITE(istdo,'(A25,I10,A24)')
1138 . ' ** MONITORED VOLUME ID: ',ivolu(1),' - BUILDING POLYHEDRA **'
1139C
1140 DO i=1,nbric
1141 IF (ilvout<0) CALL progbar_c(i,nbric)
1142C
1143 IF (tbric(7,i)/=2) cycle
1144C Polygones de la brique
1145 npolb=0
1146 nppmax=0
1147 ptr_cur => first
1148 DO j=1,npoly
1149 ity=ptr_cur%IPOLY(1)
1150 IF ((ity==1.AND.ptr_cur%IPOLY(4)==i).OR.
1151 . (ity==2.AND.
1152 . (ptr_cur%IPOLY(3)==i.OR.ptr_cur%IPOLY(4)==i))) THEN
1153 npolb=npolb+1
1154 nppmax=max(nppmax,ptr_cur%IPOLY(2))
1155 ENDIF
1156 ptr_cur => ptr_cur%PTR
1157 ENDDO
1158 IF (npolb==0) cycle
1159C
1160 nrpmax=4+3*nppmax
1161 ALLOCATE(polb(npolb), ipoly(6+nppmax,npolb),
1162 . rpoly(nrpmax,npolb), redir(npolb))
1163 DO inz=1,iiz+1
1164 npolb=0
1165 ptr_cur => first
1166 DO j=1,npoly
1167 ity=ptr_cur%IPOLY(1)
1168 ityz=ptr_cur%IZ(1)
1169 IF (((ity==1.AND.ptr_cur%IPOLY(4)==i).OR.
1170 . (ity==2.AND.(ptr_cur%IPOLY(3)==i.OR.
1171 . ptr_cur%IPOLY(4)==i)))
1172 . .AND.
1173 . ((ityz==1.AND.ptr_cur%IZ(2)==inz).OR.
1174 . (ityz==2.AND.(ptr_cur%IZ(2)==inz.OR.
1175 . ptr_cur%IZ(3)==inz)))) THEN
1176 npolb=npolb+1
1177 redir(npolb)=j
1178 nn=ptr_cur%IPOLY(2)
1179 DO k=1,6+nn
1180 ipoly(k,npolb)=ptr_cur%IPOLY(k)
1181 ENDDO
1182 DO k=1,4+3*nn
1183 rpoly(k,npolb)=ptr_cur%RPOLY(k)
1184 ENDDO
1185 polb(npolb)=npolb
1186 ENDIF
1187 ptr_cur => ptr_cur%PTR
1188 ENDDO
1189 IF (npolb==0) cycle
1190C
1191 nphmax=npolb
1192 npolhmax=nlayer
1193 300 CONTINUE
1194 IF (ALLOCATED(polh)) DEALLOCATE(polh)
1195 ALLOCATE(polh(nphmax+2,npolhmax))
1196C
1197 CALL polyhedr(ipoly, rpoly , polb, npolb, polh,
1198 . nnp, nrpmax , nphmax, i, dxm ,
1199 . info , npolhmax, nppmax, rvolu(50))
1200 IF (info==1) GOTO 300
1201C Modifications des polygones
1202 ptr_cur => first
1203 npl=1
1204 polc=redir(polb(npl))
1205 DO j=1,npoly
1206 IF (j==polc) THEN
1207 ity=ptr_cur%IPOLY(1)
1208 IF (ity==1) THEN
1209 ptr_cur%IPOLY(5)=npolh+ipoly(5,npl)
1210 ELSEIF (ity==2) THEN
1211 IF (ptr_cur%IPOLY(5)==0) THEN
1212 ptr_cur%IPOLY(5)=npolh+ipoly(5,npl)
1213 ELSE
1214 ptr_cur%IPOLY(6)=npolh+ipoly(6,npl)
1215 ENDIF
1216 ENDIF
1217 IF (npl==npolb) GOTO 350
1218 npl=npl+1
1219 polc=redir(polb(npl))
1220 ENDIF
1221 ptr_cur => ptr_cur%PTR
1222 ENDDO
1223 350 CONTINUE
1224C
1225 DO n=1,nnp
1226 npolh=npolh+1
1227 IF (.NOT.ASSOCIATED(phfirst)) THEN
1228 ALLOCATE(phfirst)
1229 ph_cur => phfirst
1230 ELSE
1231 ALLOCATE(ph_cur)
1232 ph_prec%PTR => ph_cur
1233 ENDIF
1234 nn=polh(1,n)
1235 ALLOCATE(ph_cur%POLH(2+nn))
1236 ph_cur%RANK=npolh
1237 ph_cur%POLH(1)=polh(1,n)
1238 ph_cur%POLH(2)=polh(2,n)
1239 DO m=1,nn
1240 ph_cur%POLH(2+m)=redir(polh(2+m,n))
1241 ENDDO
1242 ph_prec => ph_cur
1243 ENDDO
1244C
1245 ENDDO
1246 DEALLOCATE(ipoly, rpoly, polb, polh, redir)
1247 ENDDO
1248C---------------------------------------------------
1249C Passage listes chainees -> tableaux classiques pour la suite des traitements
1250C---------------------------------------------------
1251 370 CONTINUE
1252C
1253 ptr_cur => first
1254 lenimax=0
1255 lenrmax=0
1256 nns=0
1257 nns2=0
1258 DO i=1,npoly
1259 nn=ptr_cur%IPOLY(2)
1260 nns=nns+nn
1261 nns2=nns2+ptr_cur%NNSP
1262 IF (ptr_cur%IPOLY(1)==1) THEN
1263 lenimax=max(lenimax,6+nn+1)
1264 lenrmax=max(lenrmax,4+3*nn)
1265 ELSEIF (ptr_cur%IPOLY(1)==2) THEN
1266 nhol=ptr_cur%IPOLY(6+nn+1)
1267 lenimax=max(lenimax,6+nn+1+nhol)
1268 lenrmax=max(lenrmax,4+3*nn+3*nhol)
1269 ENDIF
1270 ptr_cur => ptr_cur%PTR
1271 ENDDO
1272 ph_cur => phfirst
1273 nphmax=0
1274 DO i=1,npolh
1275 nn=ph_cur%POLH(1)
1276 nphmax=max(nphmax,nn)
1277 ph_cur => ph_cur%PTR
1278 ENDDO
1279C
1280 npoly_old=npoly
1281 npolh_old=npolh
1282 IF (nba>0) THEN
1283c IF (NSPMD == 1) THEN
1284c ALLOCATE(ITAGN(NUMNOD))
1285c NPOLH=NPOLH+NBA
1286c DO I=1,NBA
1287c ITAGBA(I)=0
1288c ENDDO
1289cC
1290c DO I=1,NUMNOD
1291c ITAGN(I)=0
1292c ENDDO
1293c DO IEL=1,NEL
1294c N1=IBUF(ELEM(1,IEL))
1295c N2=IBUF(ELEM(2,IEL))
1296c N3=IBUF(ELEM(3,IEL))
1297c ITAGN(N1)=1
1298c ITAGN(N2)=1
1299c ITAGN(N3)=1
1300c ENDDO
1301cC
1302c DO I=1,NBA
1303c NFAC=TBA(2,I)
1304c NNP=0
1305c DO J=1,NFAC
1306c ITY=TFACA(2*(J-1)+1,I)
1307c IF (ITY==1) THEN
1308cC Face interne brique/brique
1309c NV=TFACA(2*(J-1)+2,I)
1310c IF (ITAGBA(NV)==0) THEN
1311c LENIMAX=MAX(LENIMAX,6+3+1)
1312c LENRMAX=MAX(LENRMAX,6+3*3)
1313c IF (NFAC==4) THEN
1314c NPOLY=NPOLY+1
1315c NNS=NNS+3
1316c ELSEIF (NFAC==6) THEN
1317c NPOLY=NPOLY+2
1318c NNS=NNS+6
1319c ENDIF
1320c ENDIF
1321c IF (NFAC==4) THEN
1322c NNP=NNP+1
1323c ELSEIF (NFAC==6) THEN
1324c NNP=NNP+2
1325c ENDIF
1326c ELSEIF (ITY==2) THEN
1327cC Face brique appuyee sur surface
1328c IF (NFAC==4) THEN
1329c NPOLY=NPOLY+1
1330c NNS=NNS+3
1331c NNP=NNP+1
1332c ELSEIF (NFAC==6) THEN
1333c NPOLY=NPOLY+2
1334c NNS=NNS+6
1335c NNP=NNP+2
1336c ENDIF
1337c ELSEIF (ITY==3) THEN
1338cC Face interne brique/polyhedres
1339c PTR_CUR => FIRST
1340c DO K=1,NPOLY_OLD
1341c ITY=PTR_CUR%IPOLY(1)
1342c IF (ITY==1) THEN
1343c IEL=PTR_CUR%IPOLY(3)
1344c IF (-TAGELA(IEL)==I) NNP=NNP+1
1345c ENDIF
1346c PTR_CUR => PTR_CUR%PTR
1347c ENDDO
1348c ENDIF
1349c ENDDO
1350c ITAGBA(I)=1
1351c NPHMAX=MAX(NPHMAX,NNP)
1352cC
1353c II=TBA(1,I)
1354c NALL=1
1355c IF (NFAC==4) THEN
1356c DO J=1,4
1357c JJ=IXS(1+2*(J-1)+1,II)
1358c NALL=NALL*ITAGN(JJ)
1359c ENDDO
1360c ELSEIF (NFAC==6) THEN
1361c DO J=1,8
1362c JJ=IXS(1+J,II)
1363c NALL=NALL*ITAGN(JJ)
1364c ENDDO
1365c ENDIF
1366c IBSA(I)=NALL
1367c ENDDO
1368c DEALLOCATE(ITAGN)
1369c ELSE
1370 ALLOCATE(itagn(nnsa))
1371 npolh=npolh+nba
1372 DO i=1,nba
1373 itagba(i)=0
1374 ENDDO
1375C
1376 DO i=1,nnsa
1377 itagn(i)=0
1378 ENDDO
1379 DO iel=1,nel
1380 n1=fvspmd(ifv)%ELEMSA(1,iel)
1381 n2=fvspmd(ifv)%ELEMSA(2,iel)
1382 n3=fvspmd(ifv)%ELEMSA(3,iel)
1383 IF (n1>0) itagn(n1)=1
1384 IF (n1>0) itagn(n2)=1
1385 IF (n1>0) itagn(n3)=1
1386 ENDDO
1387C
1388 DO i=1,nba
1389 nfac=tba(2,i)
1390 nnp=0
1391 DO j=1,nfac
1392 ity=tfaca(2*(j-1)+1,i)
1393 IF (ity==1) THEN
1394C Face interne brique/brique
1395 nv=tfaca(2*(j-1)+2,i)
1396 IF (itagba(nv)==0) THEN
1397 lenimax=max(lenimax,6+3+1)
1398 lenrmax=max(lenrmax,6+3*3)
1399 IF (nfac==4) THEN
1400 npoly=npoly+1
1401 nns=nns+3
1402 ELSEIF (nfac==6) THEN
1403 npoly=npoly+2
1404 nns=nns+6
1405 ENDIF
1406 ENDIF
1407 IF (nfac==4) THEN
1408 nnp=nnp+1
1409 ELSEIF (nfac==6) THEN
1410 nnp=nnp+2
1411 ENDIF
1412 ELSEIF (ity==2) THEN
1413C Face brique appuyee sur surface
1414 IF (nfac==4) THEN
1415 npoly=npoly+1
1416 nns=nns+3
1417 nnp=nnp+1
1418 ELSEIF (nfac==6) THEN
1419 npoly=npoly+2
1420 nns=nns+6
1421 nnp=nnp+2
1422 ENDIF
1423 ELSEIF (ity==3) THEN
1424C Face interne brique/polyhedres
1425 ptr_cur => first
1426 DO k=1,npoly_old
1427 ity=ptr_cur%IPOLY(1)
1428 IF (ity==1) THEN
1429 iel=ptr_cur%IPOLY(3)
1430 IF (-tagela(iel)==i) nnp=nnp+1
1431 ENDIF
1432 ptr_cur => ptr_cur%PTR
1433 ENDDO
1434 ENDIF
1435 ENDDO
1436 itagba(i)=1
1437 nphmax=max(nphmax,nnp)
1438C
1439 nall=1
1440 IF (nfac==4) THEN
1441 DO j=1,4
1442 jj=fvspmd(ifv)%IXSA(2*(j-1)+1,i)
1443 nall=nall*itagn(jj)
1444 ENDDO
1445 ELSEIF (nfac==6) THEN
1446 DO j=1,8
1447 jj=fvspmd(ifv)%IXSA(j,i)
1448 nall=nall*itagn(jj)
1449 ENDDO
1450 ENDIF
1451 ibsa(i)=nall
1452 ENDDO
1453c ENDIF
1454 ENDIF
1455C
1456 ALLOCATE(ipoly_f(lenimax,npoly), rpoly_f(lenrmax,npoly),
1457 . polh_f(2+nphmax,npolh), ifvnod(nns), volu(npolh),
1458 . ifvnod2(2,nns2), rfvnod2(4,nns2), xns(3,nns),
1459 . xns2(3,nns2))
1460C
1461 npoly=npoly_old
1462 npolh=npolh_old
1463 ptr_cur => first
1464 nns=0
1465 DO i=1,npoly
1466 nn=ptr_cur%IPOLY(2)
1467 DO j=1,nn
1468 jj=ptr_cur%IPOLY(6+j)
1469 IF (jj>0) THEN
1470 nns=nns+1
1471 xns(1,nns)=ptr_cur%RPOLY(4+3*(j-1)+1)
1472 xns(2,nns)=ptr_cur%RPOLY(4+3*(j-1)+2)
1473 xns(3,nns)=ptr_cur%RPOLY(4+3*(j-1)+3)
1474 ENDIF
1475 ENDDO
1476 ptr_cur => ptr_cur%PTR
1477 ENDDO
1478C
1479 nns=0
1480 ptr_cur => first
1481 DO i=1,npoly
1482 nn=ptr_cur%IPOLY(2)
1483 DO j=1,6
1484 ipoly_f(j,i)=ptr_cur%IPOLY(j)
1485 ENDDO
1486 DO j=1,4+3*nn
1487 rpoly_f(j,i)=ptr_cur%RPOLY(j)
1488 ENDDO
1489 DO j=1,nn
1490 jj=ptr_cur%IPOLY(6+j)
1491 IF (jj>0) THEN
1492 nns=nns+1
1493 ipoly_f(6+j,i)=nns
1494 ifvnod(nns)=ptr_cur%IELNOD(j)
1495 ELSE
1496 ipoly_f(6+j,i)=jj
1497 ENDIF
1498 ENDDO
1499 IF (ptr_cur%IPOLY(1)==1) THEN
1500 ipoly_f(4,i)=ipoly_f(5,i)
1501 ipoly_f(5,i)=0
1502 ELSEIF (ptr_cur%IPOLY(1)==2) THEN
1503 nhol=ptr_cur%IPOLY(6+nn+1)
1504 ipoly_f(6+nn+1,i)=nhol
1505 DO j=1,nhol
1506 ipoly_f(6+nn+1+j,i)=ptr_cur%IPOLY(6+nn+1+j)
1507 rpoly_f(4+3*nn+3*(j-1)+1,i)=
1508 . ptr_cur%RPOLY(4+3*nn+3*(j-1)+1)
1509 rpoly_f(4+3*nn+3*(j-1)+2,i)=
1510 . ptr_cur%RPOLY(4+3*nn+3*(j-1)+2)
1511 rpoly_f(4+3*nn+3*(j-1)+3,i)=
1512 . ptr_cur%RPOLY(4+3*nn+3*(j-1)+3)
1513 ENDDO
1514 ENDIF
1515 nnsp=ptr_cur%NNSP
1516 IF (nnsp>0) THEN
1517 DO j=1,nnsp
1518 jj=ptr_cur%NREF(1,j)
1519 ifvnod2(1,jj)=ptr_cur%NREF(2,j)
1520 ifvnod2(2,jj)=ptr_cur%NREF(3,j)
1521 rfvnod2(1,jj)=ptr_cur%AREF(1,j)
1522 rfvnod2(2,jj)=ptr_cur%AREF(2,j)
1523 rfvnod2(3,jj)=ptr_cur%AREF(3,j)
1524 rfvnod2(4,jj)=ptr_cur%AREF(4,j)
1525 ENDDO
1526 ENDIF
1527 ptr_tmp => ptr_cur%PTR
1528 DEALLOCATE(ptr_cur%IPOLY, ptr_cur%RPOLY, ptr_cur%IELNOD)
1529 IF (nnsp>0) DEALLOCATE(ptr_cur%NREF, ptr_cur%AREF)
1530 DEALLOCATE(ptr_cur)
1531 ptr_cur => ptr_tmp
1532 ENDDO
1533C
1534 DO i=1,nns2
1535 n1=ifvnod2(1,i)
1536 n2=ifvnod2(2,i)
1537 idep=0
1538 IF (n1<0) THEN
1539 ii=-n1
1540 IF (ifvnod2(1,ii)/=n2.AND.ifvnod2(2,ii)/=n2) THEN
1541 WRITE(istdo,*) 'PROBLEM DEPENDANT NODE ',i
1542 CALL arret(2)
1543 ENDIF
1544 n1=ifvnod2(1,ii)
1545 n2=ifvnod2(2,ii)
1546 idep=1
1547 ELSEIF (n2<0) THEN
1548 ii=-n2
1549 IF (ifvnod2(1,ii)/=n1.AND.ifvnod2(2,ii)/=n1) THEN
1550 WRITE(istdo,*) 'PROBLEM DEPENDANT NODE ',i
1551 CALL arret(2)
1552 ENDIF
1553 n1=ifvnod2(1,ii)
1554 n2=ifvnod2(2,ii)
1555 idep=1
1556 ENDIF
1557 ifvnod2(1,i)=n1
1558 ifvnod2(2,i)=n2
1559 IF (idep==1) THEN
1560C Recalcul des coordonnees barycentriques du point dependant
1561 ii=1
1562 val=abs(xns(1,n1)-xns(1,n2))
1563 DO j=2,3
1564 IF (abs(xns(j,n1)-xns(j,n2))>val) THEN
1565 ii=j
1566 val=abs(xns(j,n1)-xns(j,n2))
1567 ENDIF
1568 ENDDO
1569 rfvnod2(1,i)=(rfvnod2(ii+1,i)-xns(ii,n2))/
1570 . (xns(ii,n1)-xns(ii,n2))
1571 ENDIF
1572 ENDDO
1573C
1574 ph_cur => phfirst
1575 DO i=1,npolh
1576 nn=ph_cur%POLH(1)
1577 DO j=1,2+nn
1578 polh_f(j,i)=ph_cur%POLH(j)
1579 ENDDO
1580 ph_tmp => ph_cur%PTR
1581 DEALLOCATE(ph_cur%POLH)
1582 DEALLOCATE(ph_cur)
1583 ph_cur => ph_tmp
1584 ENDDO
1585C
1586 IF (nba>0) THEN
1587 DO i=1,nba
1588 itagba(i)=0
1589 ENDDO
1590C
1591 npoly_old=npoly
1592 DO i=1,nba
1593 ii=tba(1,i)
1594 nfac=tba(2,i)
1595 nnp=0
1596 DO j=1,nfac
1597 ity=tfaca(2*(j-1)+1,i)
1598 IF (ity==1) THEN
1599C Face interne brique/brique
1600 nv=tfaca(2*(j-1)+2,i)
1601 IF (itagba(nv)==0) THEN
1602 IF (nfac==4) THEN
1603 npoly=npoly+1
1604 n1=fac4(1,j)
1605 n2=fac4(2,j)
1606 n3=fac4(3,j)
1607 ipoly_f(1,npoly)=2
1608 ipoly_f(2,npoly)=3
1609 ipoly_f(3,npoly)=0
1610 ipoly_f(4,npoly)=0
1611 ipoly_f(5,npoly)=npolh+i
1612 ipoly_f(6,npoly)=npolh+nv
1613 ipoly_f(6+1,npoly)=nns+1
1614 ipoly_f(6+2,npoly)=nns+2
1615 ipoly_f(6+3,npoly)=nns+3
1616 ipoly_f(6+3+1,npoly)=0
1617C
1618c IF (NSPMD == 1) THEN
1619c NN1=IXS(1+N1,II)
1620c NN2=IXS(1+N2,II)
1621c NN3=IXS(1+N3,II)
1622c IFVNOD(NNS+1)=-NN1
1623c IFVNOD(NNS+2)=-NN2
1624c IFVNOD(NNS+3)=-NN3
1625c NNS=NNS+3
1626c X1=X(1,NN1)
1627c Y1=X(2,NN1)
1628c Z1=X(3,NN1)
1629c X2=X(1,NN2)
1630c Y2=X(2,NN2)
1631c Z2=X(3,NN2)
1632c X3=X(1,NN3)
1633c Y3=X(2,NN3)
1634c Z3=X(3,NN3)
1635c ELSE
1636 nn1=fvspmd(ifv)%IXSA(n1,i)
1637 nn2=fvspmd(ifv)%IXSA(n2,i)
1638 nn3=fvspmd(ifv)%IXSA(n3,i)
1639 ifvnod(nns+1)=-nn1
1640 ifvnod(nns+2)=-nn2
1641 ifvnod(nns+3)=-nn3
1642 nns=nns+3
1643 x1=xxxsa(1,nn1)
1644 y1=xxxsa(2,nn1)
1645 z1=xxxsa(3,nn1)
1646 x2=xxxsa(1,nn2)
1647 y2=xxxsa(2,nn2)
1648 z2=xxxsa(3,nn2)
1649 x3=xxxsa(1,nn3)
1650 y3=xxxsa(2,nn3)
1651 z3=xxxsa(3,nn3)
1652c ENDIF
1653 x12=x2-x1
1654 y12=y2-y1
1655 z12=z2-z1
1656 x13=x3-x1
1657 y13=y3-y1
1658 z13=z3-z1
1659 nrx=y12*z13-z12*y13
1660 nry=z12*x13-x12*z13
1661 nrz=x12*y13-y12*x13
1662 area2=sqrt(nrx**2+nry**2+nrz**2)
1663 rpoly_f(2,npoly)=nrx/area2
1664 rpoly_f(3,npoly)=nry/area2
1665 rpoly_f(4,npoly)=nrz/area2
1666 rpoly_f(4+1,npoly)=x1
1667 rpoly_f(4+2,npoly)=y1
1668 rpoly_f(4+3,npoly)=z1
1669 rpoly_f(4+4,npoly)=x2
1670 rpoly_f(4+5,npoly)=y2
1671 rpoly_f(4+6,npoly)=z2
1672 rpoly_f(4+7,npoly)=x3
1673 rpoly_f(4+8,npoly)=y3
1674 rpoly_f(4+9,npoly)=z3
1675C
1676 nnp=nnp+1
1677 polh_f(2+nnp,npolh+i)=npoly
1678 ELSEIF (nfac==6) THEN
1679 n1=fac(1,j)
1680 n2=fac(2,j)
1681 n3=fac(3,j)
1682 n4=fac(4,j)
1683c IF (NSPMD == 1) THEN
1684c NN1=IXS(1+N1,II)
1685c NN2=IXS(1+N2,II)
1686c NN3=IXS(1+N3,II)
1687c NN4=IXS(1+N4,II)
1688c X1=X(1,NN1)
1689c Y1=X(2,NN1)
1690c Z1=X(3,NN1)
1691c X2=X(1,NN2)
1692c Y2=X(2,NN2)
1693c Z2=X(3,NN2)
1694c X3=X(1,NN3)
1695c Y3=X(2,NN3)
1696c Z3=X(3,NN3)
1697c X4=X(1,NN4)
1698c Y4=X(2,NN4)
1699c Z4=X(3,NN4)
1700c ELSE
1701 nn1=fvspmd(ifv)%IXSA(n1,i)
1702 nn2=fvspmd(ifv)%IXSA(n2,i)
1703 nn3=fvspmd(ifv)%IXSA(n3,i)
1704 nn4=fvspmd(ifv)%IXSA(n4,i)
1705 x1=xxxsa(1,nn1)
1706 y1=xxxsa(2,nn1)
1707 z1=xxxsa(3,nn1)
1708 x2=xxxsa(1,nn2)
1709 y2=xxxsa(2,nn2)
1710 z2=xxxsa(3,nn2)
1711 x3=xxxsa(1,nn3)
1712 y3=xxxsa(2,nn3)
1713 z3=xxxsa(3,nn3)
1714 x4=xxxsa(1,nn4)
1715 y4=xxxsa(2,nn4)
1716 z4=xxxsa(3,nn4)
1717c ENDIF
1718C
1719 npoly=npoly+1
1720 ipoly_f(1,npoly)=2
1721 ipoly_f(2,npoly)=3
1722 ipoly_f(3,npoly)=0
1723 ipoly_f(4,npoly)=0
1724 ipoly_f(5,npoly)=npolh+i
1725 ipoly_f(6,npoly)=npolh+nv
1726 ipoly_f(6+1,npoly)=nns+1
1727 ipoly_f(6+2,npoly)=nns+2
1728 ipoly_f(6+3,npoly)=nns+3
1729 ipoly_f(6+3+1,npoly)=0
1730 ifvnod(nns+1)=-nn1
1731 ifvnod(nns+2)=-nn2
1732 ifvnod(nns+3)=-nn3
1733 nns=nns+3
1734 x12=x2-x1
1735 y12=y2-y1
1736 z12=z2-z1
1737 x13=x3-x1
1738 y13=y3-y1
1739 z13=z3-z1
1740 nrx=y12*z13-z12*y13
1741 nry=z12*x13-x12*z13
1742 nrz=x12*y13-y12*x13
1743 area2=sqrt(nrx**2+nry**2+nrz**2)
1744 rpoly_f(2,npoly)=nrx/area2
1745 rpoly_f(3,npoly)=nry/area2
1746 rpoly_f(4,npoly)=nrz/area2
1747 rpoly_f(4+1,npoly)=x1
1748 rpoly_f(4+2,npoly)=y1
1749 rpoly_f(4+3,npoly)=z1
1750 rpoly_f(4+4,npoly)=x2
1751 rpoly_f(4+5,npoly)=y2
1752 rpoly_f(4+6,npoly)=z2
1753 rpoly_f(4+7,npoly)=x3
1754 rpoly_f(4+8,npoly)=y3
1755 rpoly_f(4+9,npoly)=z3
1756C
1757 nnp=nnp+1
1758 polh_f(2+nnp,npolh+i)=npoly
1759C
1760 npoly=npoly+1
1761 ipoly_f(1,npoly)=2
1762 ipoly_f(2,npoly)=3
1763 ipoly_f(3,npoly)=0
1764 ipoly_f(4,npoly)=0
1765 ipoly_f(5,npoly)=npolh+i
1766 ipoly_f(6,npoly)=npolh+nv
1767 ipoly_f(6+1,npoly)=nns+1
1768 ipoly_f(6+2,npoly)=nns+2
1769 ipoly_f(6+3,npoly)=nns+3
1770 ipoly_f(6+3+1,npoly)=0
1771 ifvnod(nns+1)=-nn1
1772 ifvnod(nns+2)=-nn3
1773 ifvnod(nns+3)=-nn4
1774 nns=nns+3
1775 x13=x3-x1
1776 y13=y3-y1
1777 z13=z3-z1
1778 x14=x4-x1
1779 y14=y4-y1
1780 z14=z4-z1
1781 nrx=y13*z14-z13*y14
1782 nry=z13*x14-x13*z14
1783 nrz=x13*y14-y13*x14
1784 area2=sqrt(nrx**2+nry**2+nrz**2)
1785 rpoly_f(2,npoly)=nrx/area2
1786 rpoly_f(3,npoly)=nry/area2
1787 rpoly_f(4,npoly)=nrz/area2
1788 rpoly_f(4+1,npoly)=x1
1789 rpoly_f(4+2,npoly)=y1
1790 rpoly_f(4+3,npoly)=z1
1791 rpoly_f(4+4,npoly)=x3
1792 rpoly_f(4+5,npoly)=y3
1793 rpoly_f(4+6,npoly)=z3
1794 rpoly_f(4+7,npoly)=x4
1795 rpoly_f(4+8,npoly)=y4
1796 rpoly_f(4+9,npoly)=z4
1797C
1798 nnp=nnp+1
1799 polh_f(2+nnp,npolh+i)=npoly
1800 ENDIF
1801 ELSE
1802 DO k=1,polh_f(1,npolh+nv)
1803 kk=polh_f(2+k,npolh+nv)
1804 IF (ipoly_f(1,kk)==2.AND.
1805 . ipoly_f(6,kk)==npolh+i) THEN
1806 nnp=nnp+1
1807 polh_f(2+nnp,npolh+i)=kk
1808 ENDIF
1809 ENDDO
1810 ENDIF
1811 ELSEIF (ity==2) THEN
1812C Face brique appuyee sur surface (traitee plus loin)
1813 ELSEIF (ity==3) THEN
1814C Face interne brique/polyhedres
1815 DO k=1,npoly_old
1816 ity=ipoly_f(1,k)
1817 IF (ity==1) THEN
1818 iel=ipoly_f(3,k)
1819 IF (-tagela(iel)==i) THEN
1820 nnp=nnp+1
1821 polh_f(2+nnp,npolh+i)=k
1822 ipoly_f(1,k)=2
1823 ipoly_f(3,k)=0
1824 ip=ipoly_f(4,k)
1825 ipoly_f(4,k)=0
1826 ipoly_f(5,k)=ip
1827 ipoly_f(6,k)=npolh+i
1828 ENDIF
1829 ENDIF
1830 ENDDO
1831 ENDIF
1832 ENDDO
1833 polh_f(1,npolh+i)=nnp
1834 polh_f(2,npolh+i)=-i
1835 itagba(i)=1
1836 ENDDO
1837C
1838 npolh=npolh+nba
1839C
1840 DO i=1,npoly_old
1841 IF (ipoly_f(1,i)==1) THEN
1842 iel=ipoly_f(3,i)
1843 IF (tagela(iel)>0) ipoly_f(3,i)=tagela(iel)
1844 ENDIF
1845 ENDDO
1846C
1847 DO iel=1,nel
1848 IF (tagels(iel)>0) THEN
1849 npoly=npoly+1
1850 ipoly_f(1,npoly)=1
1851 ipoly_f(2,npoly)=3
1852 ipoly_f(3,npoly)=iel
1853 ipoly_f(4,npoly)=npolh_old+tagels(iel)
1854 ipoly_f(5,npoly)=0
1855 ipoly_f(6,npoly)=0
1856 ipoly_f(7,npoly)=nns+1
1857 ipoly_f(8,npoly)=nns+2
1858 ipoly_f(9,npoly)=nns+3
1859c IF (NSPMD == 1) THEN
1860c N1=ELEM(1,IEL)
1861c N2=ELEM(2,IEL)
1862c N3=ELEM(3,IEL)
1863c NN1=IBUF(N1)
1864c NN2=IBUF(N2)
1865c NN3=IBUF(N3)
1866c IFVNOD(NNS+1)=-NN1
1867c IFVNOD(NNS+2)=-NN2
1868c IFVNOD(NNS+3)=-NN3
1869c NNS=NNS+3
1870c RPOLY_F(1,NPOLY)=ELAREA(IEL)
1871c RPOLY_F(2,NPOLY)=NORM(1,IEL)
1872c RPOLY_F(3,NPOLY)=NORM(2,IEL)
1873c RPOLY_F(4,NPOLY)=NORM(3,IEL)
1874c X1=X(1,NN1)
1875c Y1=X(2,NN1)
1876c Z1=X(3,NN1)
1877c X2=X(1,NN2)
1878c Y2=X(2,NN2)
1879c Z2=X(3,NN2)
1880c X3=X(1,NN3)
1881c Y3=X(2,NN3)
1882c Z3=X(3,NN3)
1883c ELSE
1884 nn1=fvspmd(ifv)%ELEMSA(1,iel)
1885 nn2=fvspmd(ifv)%ELEMSA(2,iel)
1886 nn3=fvspmd(ifv)%ELEMSA(3,iel)
1887 ifvnod(nns+1)=-nn1
1888 ifvnod(nns+2)=-nn2
1889 ifvnod(nns+3)=-nn3
1890 nns=nns+3
1891 rpoly_f(1,npoly)=elarea(iel)
1892 rpoly_f(2,npoly)=norm(1,iel)
1893 rpoly_f(3,npoly)=norm(2,iel)
1894 rpoly_f(4,npoly)=norm(3,iel)
1895 x1=xxxsa(1,nn1)
1896 y1=xxxsa(2,nn1)
1897 z1=xxxsa(3,nn1)
1898 x2=xxxsa(1,nn2)
1899 y2=xxxsa(2,nn2)
1900 z2=xxxsa(3,nn2)
1901 x3=xxxsa(1,nn3)
1902 y3=xxxsa(2,nn3)
1903 z3=xxxsa(3,nn3)
1904c ENDIF
1905 rpoly_f(5,npoly)=x1
1906 rpoly_f(6,npoly)=y1
1907 rpoly_f(7,npoly)=z1
1908 rpoly_f(8,npoly)=x2
1909 rpoly_f(9,npoly)=y2
1910 rpoly_f(10,npoly)=z2
1911 rpoly_f(11,npoly)=x3
1912 rpoly_f(12,npoly)=y3
1913 rpoly_f(13,npoly)=z3
1914C
1915 nnp=polh_f(1,npolh_old+tagels(iel))
1916 nnp=nnp+1
1917 polh_f(1,npolh_old+tagels(iel))=nnp
1918 polh_f(2+nnp,npolh_old+tagels(iel))=npoly
1919 ENDIF
1920 ENDDO
1921 ENDIF
1922C---------------------------------------------------
1923C Surface des polygones
1924C---------------------------------------------------
1925 DO i=1,npoly
1926 ity=ipoly_f(1,i)
1927 nn=ipoly_f(2,i)
1928 nhol=0
1929 IF (ity==2) nhol=ipoly_f(6+nn+1,i)
1930 area=zero
1931 nx=rpoly_f(2,i)
1932 ny=rpoly_f(3,i)
1933 nz=rpoly_f(4,i)
1934 IF (nhol==0) THEN
1935 x1=rpoly_f(5,i)
1936 y1=rpoly_f(6,i)
1937 z1=rpoly_f(7,i)
1938 DO j=1,nn-2
1939 jj=j+1
1940 jjj=jj+1
1941 x2=rpoly_f(4+3*(jj-1)+1,i)
1942 y2=rpoly_f(4+3*(jj-1)+2,i)
1943 z2=rpoly_f(4+3*(jj-1)+3,i)
1944 x3=rpoly_f(4+3*(jjj-1)+1,i)
1945 y3=rpoly_f(4+3*(jjj-1)+2,i)
1946 z3=rpoly_f(4+3*(jjj-1)+3,i)
1947 x12=x2-x1
1948 y12=y2-y1
1949 z12=z2-z1
1950 x13=x3-x1
1951 y13=y3-y1
1952 z13=z3-z1
1953 nnx=y12*z13-z12*y13
1954 nny=z12*x13-x12*z13
1955 nnz=x12*y13-y12*x13
1956 area=area+(nnx*nx+nny*ny+nnz*nz)
1957 ENDDO
1958 ELSE
1959 ALLOCATE(adr(nhol))
1960 DO j=1,nhol
1961 adr(j)=ipoly_f(6+nn+1+j,i)
1962 ENDDO
1963 adr(nhol+1)=nn
1964 x1=rpoly_f(5,i)
1965 y1=rpoly_f(6,i)
1966 z1=rpoly_f(7,i)
1967 DO j=1,adr(1)-2
1968 jj=j+1
1969 jjj=jj+1
1970 x2=rpoly_f(4+3*(jj-1)+1,i)
1971 y2=rpoly_f(4+3*(jj-1)+2,i)
1972 z2=rpoly_f(4+3*(jj-1)+3,i)
1973 x3=rpoly_f(4+3*(jjj-1)+1,i)
1974 y3=rpoly_f(4+3*(jjj-1)+2,i)
1975 z3=rpoly_f(4+3*(jjj-1)+3,i)
1976 x12=x2-x1
1977 y12=y2-y1
1978 z12=z2-z1
1979 x13=x3-x1
1980 y13=y3-y1
1981 z13=z3-z1
1982 nnx=y12*z13-z12*y13
1983 nny=z12*x13-x12*z13
1984 nnz=x12*y13-y12*x13
1985 area=area+(nnx*nx+nny*ny+nnz*nz)
1986 ENDDO
1987 area=abs(area)
1988C On retranche la surface des trous
1989 DO j=1,nhol
1990 x1=rpoly_f(4+3*adr(j)+1,i)
1991 y1=rpoly_f(4+3*adr(j)+2,i)
1992 z1=rpoly_f(4+3*adr(j)+3,i)
1993 area2=zero
1994 DO k=adr(j)+1,adr(j+1)-2
1995 kk=k+1
1996 kkk=kk+1
1997 x2=rpoly_f(4+3*(kk-1)+1,i)
1998 y2=rpoly_f(4+3*(kk-1)+2,i)
1999 z2=rpoly_f(4+3*(kk-1)+3,i)
2000 x3=rpoly_f(4+3*(kkk-1)+1,i)
2001 y3=rpoly_f(4+3*(kkk-1)+2,i)
2002 z3=rpoly_f(4+3*(kkk-1)+3,i)
2003 x12=x2-x1
2004 y12=y2-y1
2005 z12=z2-z1
2006 x13=x3-x1
2007 y13=y3-y1
2008 z13=z3-z1
2009 nnx=y12*z13-z12*y13
2010 nny=z12*x13-x12*z13
2011 nnz=x12*y13-y12*x13
2012 area2=area2+(nnx*nx+nny*ny+nnz*nz)
2013 ENDDO
2014 area=area-abs(area2)
2015 ENDDO
2016 DEALLOCATE(adr)
2017 ENDIF
2018 rpoly_f(1,i)=half*abs(area)
2019 ENDDO
2020C---------------------------------------------------
2021C Volume des polyhedres
2022C---------------------------------------------------
2023 DO i=1,npolh
2024 volu(i)=zero
2025 DO j=1,polh_f(1,i)
2026 jj=polh_f(2+j,i)
2027 area=rpoly_f(1,jj)
2028 ity=ipoly_f(1,jj)
2029 IF (ity==1) THEN
2030 nx=rpoly_f(2,jj)
2031 ny=rpoly_f(3,jj)
2032 nz=rpoly_f(4,jj)
2033 ELSEIF (ity==2) THEN
2034 IF (ipoly_f(5,jj)==i) THEN
2035 nx=rpoly_f(2,jj)
2036 ny=rpoly_f(3,jj)
2037 nz=rpoly_f(4,jj)
2038 ELSEIF (ipoly_f(6,jj)==i) THEN
2039 nx=-rpoly_f(2,jj)
2040 ny=-rpoly_f(3,jj)
2041 nz=-rpoly_f(4,jj)
2042 ENDIF
2043 ENDIF
2044 x1=rpoly_f(5,jj)
2045 y1=rpoly_f(6,jj)
2046 z1=rpoly_f(7,jj)
2047 volu(i)=volu(i)+third*area*(x1*nx+y1*ny+z1*nz)
2048 ENDDO
2049 ENDDO
2050C---------------------------------------------------
2051C Elimination des segments de longueur nulle dans les polygones
2052C---------------------------------------------------
2053 nns_old=nns
2054 ALLOCATE(ifvnod_old(nns_old), redir(nns_old))
2055 DO i=1,nns_old
2056 ifvnod_old(i)=ifvnod(i)
2057 ENDDO
2058 nns_old=0
2059 nns=0
2060C
2061 tole=rvolu(50)
2062 DO i=1,npoly
2063 nnp_old=ipoly_f(2,i)
2064C
2065 lmax2=zero
2066 x0=rpoly_f(5,i)
2067 y0=rpoly_f(6,i)
2068 z0=rpoly_f(7,i)
2069 DO j=2,nnp_old
2070 x1=rpoly_f(4+3*(j-1)+1,i)
2071 y1=rpoly_f(4+3*(j-1)+2,i)
2072 z1=rpoly_f(4+3*(j-1)+3,i)
2073 lmax2=max(lmax2,(x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
2074 x0=x1
2075 y0=y1
2076 z0=z1
2077 ENDDO
2078 x1=rpoly_f(5,i)
2079 y1=rpoly_f(6,i)
2080 z1=rpoly_f(7,i)
2081 lmax2=max(lmax2,(x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
2082 tole2=tole**2*lmax2
2083C
2084 nhol=0
2085 IF (ipoly_f(1,i)==2) nhol=ipoly_f(6+nnp_old+1,i)
2086 ALLOCATE(ipoly_old(nnp_old), rpoly_old(3*nnp_old),
2087 . adr_old(nhol+2), adr(nhol+2))
2088 DO j=1,nnp_old
2089 ipoly_old(j)=ipoly_f(6+j,i)
2090 ENDDO
2091 DO j=1,3*nnp_old
2092 rpoly_old(j)=rpoly_f(4+j,i)
2093 ENDDO
2094 nnp=0
2095 IF (nhol==0) THEN
2096 x0=rpoly_f(5,i)
2097 y0=rpoly_f(6,i)
2098 z0=rpoly_f(7,i)
2099 IF (ipoly_old(1)>0) THEN
2100 nns_old=nns_old+1
2101 nns=nns+1
2102 redir(nns_old)=nns
2103 ipoly_f(7,i)=nns
2104 ifvnod(nns)=ifvnod_old(nns_old)
2105 ELSE
2106 ipoly_f(7,i)=ipoly_old(1)
2107 ENDIF
2108 nnp=nnp+1
2109 j0=1
2110 nns1=nns
2111 DO j=2,nnp_old
2112 IF (ipoly_old(j)>0) nns_old=nns_old+1
2113 x1=rpoly_old(3*(j-1)+1)
2114 y1=rpoly_old(3*(j-1)+2)
2115 z1=rpoly_old(3*(j-1)+3)
2116 dd=(x1-x0)**2+(y1-y0)**2+(z1-z0)**2
2117 IF (dd>tole2) THEN
2118 nnp=nnp+1
2119 IF (ipoly_old(j)>0) THEN
2120 nns=nns+1
2121 ifvnod(nns)=ifvnod_old(nns_old)
2122 rpoly_f(4+3*(nnp-1)+1,i)=x1
2123 rpoly_f(4+3*(nnp-1)+2,i)=y1
2124 rpoly_f(4+3*(nnp-1)+3,i)=z1
2125 ipoly_f(6+nnp,i)=nns
2126 ELSE
2127 rpoly_f(4+3*(nnp-1)+1,i)=x1
2128 rpoly_f(4+3*(nnp-1)+2,i)=y1
2129 rpoly_f(4+3*(nnp-1)+3,i)=z1
2130 ipoly_f(6+nnp,i)=ipoly_old(j)
2131 ENDIF
2132 x0=x1
2133 y0=y1
2134 z0=z1
2135 ELSEIF (ipoly_old(j)>0.AND.
2136 . ipoly_old(j0)<0) THEN
2137 nns=nns+1
2138 ifvnod(nns)=ifvnod_old(nns_old)
2139 rpoly_f(4+3*(nnp-1)+1,i)=x1
2140 rpoly_f(4+3*(nnp-1)+2,i)=y1
2141 rpoly_f(4+3*(nnp-1)+3,i)=z1
2142 ipoly_f(6+nnp,i)=nns
2143 ENDIF
2144 IF (ipoly_old(j)>0) redir(nns_old)=nns
2145 j0=j
2146 ENDDO
2147 x1=rpoly_f(5,i)
2148 y1=rpoly_f(6,i)
2149 z1=rpoly_f(7,i)
2150 dd=(x1-x0)**2+(y1-y0)**2+(z1-z0)**2
2151 IF (dd<=tole2.AND.ipoly_old(1)>0) THEN
2152 redir(nns_old)=nns1
2153 nns=nns-1
2154 nnp=nnp-1
2155 ELSEIF (dd<=tole2.AND.ipoly_old(1)<0) THEN
2156 nnp=nnp-1
2157 rpoly_f(5,i)=x0
2158 rpoly_f(6,i)=y0
2159 rpoly_f(7,i)=z0
2160 ipoly_f(7,i)=nns
2161 ENDIF
2162 ipoly_f(6+nnp+1,i)=0
2163 ELSE
2164 adr_old(1)=0
2165 adr(1)=0
2166 DO j=1,nhol
2167 adr_old(j+1)=ipoly_f(6+nnp_old+1+j,i)
2168 ENDDO
2169 adr_old(nhol+2)=nnp_old
2170 DO j=1,nhol+1
2171 x0=rpoly_f(4+3*adr(j)+1,i)
2172 y0=rpoly_f(4+3*adr(j)+2,i)
2173 z0=rpoly_f(4+3*adr(j)+3,i)
2174 IF (ipoly_old(adr_old(j)+1)>0) THEN
2175 nns_old=nns_old+1
2176 nns=nns+1
2177 redir(nns_old)=nns
2178 ipoly_f(6+adr(j)+1,i)=nns
2179 ifvnod(nns)=ifvnod_old(nns_old)
2180 ELSE
2181 ipoly_f(6+adr(j)+1,i)=ipoly_old(adr_old(j)+1)
2182 ENDIF
2183 nnp=nnp+1
2184 k0=1
2185 nns1=nns
2186 DO k=adr_old(j)+2,adr_old(j+1)
2187 nns_old=nns_old+1
2188 x1=rpoly_old(3*(k-1)+1)
2189 y1=rpoly_old(3*(k-1)+2)
2190 z1=rpoly_old(3*(k-1)+3)
2191 dd=(x1-x0)**2+(y1-y0)**2+(z1-z0)**2
2192 IF (dd>tole2) THEN
2193 nnp=nnp+1
2194 IF (ipoly_old(k)>0) THEN
2195 nns=nns+1
2196 ifvnod(nns)=ifvnod_old(nns_old)
2197 rpoly_f(4+3*(nnp-1)+1,i)=x1
2198 rpoly_f(4+3*(nnp-1)+2,i)=y1
2199 rpoly_f(4+3*(nnp-1)+3,i)=z1
2200 ipoly_f(6+nnp,i)=nns
2201 ELSE
2202 rpoly_f(4+3*(nnp-1)+1,i)=x1
2203 rpoly_f(4+3*(nnp-1)+2,i)=y1
2204 rpoly_f(4+3*(nnp-1)+3,i)=z1
2205 ipoly_f(6+nnp,i)=ipoly_old(k)
2206 ENDIF
2207 x0=x1
2208 y0=y1
2209 z0=z1
2210 ELSEIF (ipoly_old(k)>0.AND.
2211 . ipoly_old(k0)<0) THEN
2212 nns=nns+1
2213 ifvnod(nns)=ifvnod_old(nns_old)
2214 rpoly_f(4+3*(nnp-1)+1,i)=x1
2215 rpoly_f(4+3*(nnp-1)+2,i)=y1
2216 rpoly_f(4+3*(nnp-1)+3,i)=z1
2217 ipoly_f(6+nnp,i)=nns
2218 ENDIF
2219 IF (ipoly_old(k)>0) redir(nns_old)=nns
2220 k0=k
2221 ENDDO
2222 x1=rpoly_f(4+3*(adr(j)-1)+1,i)
2223 y1=rpoly_f(4+3*(adr(j)-1)+2,i)
2224 z1=rpoly_f(4+3*(adr(j)-1)+3,i)
2225 dd=(x1-x0)**2+(y1-y0)**2+(z1-z0)**2
2226 IF (dd<=tole2.AND.ipoly_old(adr_old(j)+1)>0) THEN
2227 redir(nns_old)=nns1
2228 nns=nns-1
2229 nnp=nnp-1
2230 ELSEIF (dd<=tole2.AND.
2231 . ipoly_old(adr_old(j)+1)<0) THEN
2232 nnp=nnp-1
2233 rpoly_f(4+3*(adr(j)-1)+1,i)=x0
2234 rpoly_f(4+3*(adr(j)-1)+2,i)=y0
2235 rpoly_f(4+3*(adr(j)-1)+3,i)=z0
2236 ipoly_f(6+adr(j)+1,i)=nns
2237 ENDIF
2238 adr(j+1)=nnp
2239 ENDDO
2240 ipoly_f(6+nnp+1,i)=nhol
2241 DO j=1,nhol
2242 ipoly_f(6+nnp+1+j,i)=adr(j+1)
2243 rpoly_f(4+3*nnp+3*(j-1)+1,i)=
2244 . rpoly_f(4+3*nnp_old+3*(j-1)+1,i)
2245 rpoly_f(4+3*nnp+3*(j-1)+2,i)=
2246 . rpoly_f(4+3*nnp_old+3*(j-1)+2,i)
2247 rpoly_f(4+3*nnp+3*(j-1)+3,i)=
2248 . rpoly_f(4+3*nnp_old+3*(j-1)+3,i)
2249 ENDDO
2250 ENDIF
2251 ipoly_f(2,i)=nnp
2252C
2253 IF (nnp<nnp_old) THEN
2254C Nouvelle aire du polygone
2255 area=zero
2256 nx=rpoly_f(2,i)
2257 ny=rpoly_f(3,i)
2258 nz=rpoly_f(4,i)
2259 IF (nhol==0) THEN
2260 x1=rpoly_f(5,i)
2261 y1=rpoly_f(6,i)
2262 z1=rpoly_f(7,i)
2263 DO j=1,ipoly_f(2,i)-2
2264 jj=j+1
2265 jjj=jj+1
2266 x2=rpoly_f(4+3*(jj-1)+1,i)
2267 y2=rpoly_f(4+3*(jj-1)+2,i)
2268 z2=rpoly_f(4+3*(jj-1)+3,i)
2269 x3=rpoly_f(4+3*(jjj-1)+1,i)
2270 y3=rpoly_f(4+3*(jjj-1)+2,i)
2271 z3=rpoly_f(4+3*(jjj-1)+3,i)
2272 x12=x2-x1
2273 y12=y2-y1
2274 z12=z2-z1
2275 x13=x3-x1
2276 y13=y3-y1
2277 z13=z3-z1
2278 nnx=y12*z13-z12*y13
2279 nny=z12*x13-x12*z13
2280 nnz=x12*y13-y12*x13
2281 area=area+(nnx*nx+nny*ny+nnz*nz)
2282 ENDDO
2283 ELSE
2284 x1=rpoly_f(5,i)
2285 y1=rpoly_f(6,i)
2286 z1=rpoly_f(7,i)
2287 DO j=1,adr(1)-2
2288 jj=j+1
2289 jjj=jj+1
2290 x2=rpoly_f(4+3*(jj-1)+1,i)
2291 y2=rpoly_f(4+3*(jj-1)+2,i)
2292 z2=rpoly_f(4+3*(jj-1)+3,i)
2293 x3=rpoly_f(4+3*(jjj-1)+1,i)
2294 y3=rpoly_f(4+3*(jjj-1)+2,i)
2295 z3=rpoly_f(4+3*(jjj-1)+3,i)
2296 x12=x2-x1
2297 y12=y2-y1
2298 z12=z2-z1
2299 x13=x3-x1
2300 y13=y3-y1
2301 z13=z3-z1
2302 nnx=y12*z13-z12*y13
2303 nny=z12*x13-x12*z13
2304 nnz=x12*y13-y12*x13
2305 area=area+(nnx*nx+nny*ny+nnz*nz)
2306 ENDDO
2307 ENDIF
2308C
2309 DO j=1,nhol
2310 x1=rpoly_f(4+3*adr(j+1)+1,i)
2311 y1=rpoly_f(4+3*adr(j+1)+2,i)
2312 z1=rpoly_f(4+3*adr(j+1)+3,i)
2313 DO k=adr(j+1)+1,adr(j+2)
2314 kk=k+1
2315 kkk=kk+1
2316 x2=rpoly_f(4+3*(kk-1)+1,i)
2317 y2=rpoly_f(4+3*(kk-1)+2,i)
2318 z2=rpoly_f(4+3*(kk-1)+3,i)
2319 x3=rpoly_f(4+3*(kkk-1)+1,i)
2320 y3=rpoly_f(4+3*(kkk-1)+2,i)
2321 z3=rpoly_f(4+3*(kkk-1)+3,i)
2322 x12=x2-x1
2323 y12=y2-y1
2324 z12=z2-z1
2325 x13=x3-x1
2326 y13=y3-y1
2327 z13=z3-z1
2328 nnx=y12*z13-z12*y13
2329 nny=z12*x13-x12*z13
2330 nnz=x12*y13-y12*x13
2331 area=area-(nnx*nx+nny*ny+nnz*nz)
2332 ENDDO
2333 ENDDO
2334 rpoly_f(1,i)=half*abs(area)
2335 ENDIF
2336C
2337 DEALLOCATE(ipoly_old, rpoly_old, adr_old, adr)
2338 ENDDO
2339C
2340C
2341 DO i=1,nns2
2342 i1=ifvnod2(1,i)
2343 i2=ifvnod2(2,i)
2344 ifvnod2(1,i)=redir(i1)
2345 ifvnod2(2,i)=redir(i2)
2346 ENDDO
2347C
2348 DEALLOCATE(ifvnod_old, redir)
2349C---------------------------------------------------
2350C Fusion des polyhedres a volume infinitesimal a leur voisin
2351C---------------------------------------------------
2352 nphmax=2*nphmax
2353 info=0
2354 ALLOCATE(volu_old(npolh), rpoly_f_old(lenrmax,npoly),
2355 . ipoly_f_old(lenimax,npoly))
2356 DO i=1,npolh
2357 volu_old(i)=volu(i)
2358 ENDDO
2359 DO i=1,npoly
2360 DO j=1,lenimax
2361 ipoly_f_old(j,i)=ipoly_f(j,i)
2362 ENDDO
2363 DO j=1,lenrmax
2364 rpoly_f_old(j,i)=rpoly_f(j,i)
2365 ENDDO
2366 ENDDO
2367 400 IF (ALLOCATED(polh_new)) DEALLOCATE(polh_new)
2368 IF (ALLOCATED(imerged)) DEALLOCATE(imerged)
2369 ALLOCATE(polh_new(2+nphmax,npolh), imerged(npolh))
2370C
2371 IF (info==1) THEN
2372 DO i=1,npolh
2373 volu(i)=volu_old(i)
2374 ENDDO
2375 DO i=1,npoly
2376 DO j=1,lenimax
2377 ipoly_f(j,i)=ipoly_f_old(j,i)
2378 ENDDO
2379 DO j=1,lenrmax
2380 rpoly_f(j,i)=rpoly_f_old(j,i)
2381 ENDDO
2382 ENDDO
2383 ENDIF
2384C
2385 DO i=1,npolh
2386 imerged(i)=0
2387 DO j=1,2+polh_f(1,i)
2388 polh_new(j,i)=polh_f(j,i)
2389 ENDDO
2390 ENDDO
2391C
2392 DO i=1,npolh
2393 IF (polh_f(2,i)<0) cycle
2394 IF (volu(i)<zero) THEN
2395 DO j=1,polh_f(1,i)
2396 jj=polh_f(2+j,i)
2397 rpoly_f(1,jj)=-one
2398 ENDDO
2399 volu(i)=zero
2400 ENDIF
2401 ENDDO
2402C
2403C Volume moyen
2404 vm=zero
2405 npa=0
2406 DO i=1,npolh
2407 IF (volu(i)>zero) THEN
2408 vm=vm+volu(i)
2409 npa=npa+1
2410 ENDIF
2411 ENDDO
2412 vm=vm/npa
2413C
2414 rvolu(33)=vm
2415 volumin=vm*rvolu(31)
2416 info=0
2417 DO i=1,npolh
2418 IF (polh_new(2,i)<0) THEN
2419 ii=-polh_new(2,i)
2420C On ne merge pas les briques additionnelles qui ont des noeuds a l'interieur de l'airbag
2421 IF(ibsa(ii)==0) cycle
2422 ENDIF
2423C
2424 IF (volu(i)<=volumin) THEN
2425 areamax=zero
2426 jmax=0
2427 imax=0
2428 DO j=1,polh_new(1,i)
2429 jj=polh_new(2+j,i)
2430 ity=ipoly_f(1,jj)
2431 IF (ity==1) cycle
2432 area=rpoly_f(1,jj)
2433 IF (area>areamax) THEN
2434 IF (ipoly_f(5,jj)==i) THEN
2435 ii=ipoly_f(6,jj)
2436 IF (polh_new(2,ii)<0) THEN
2437 iii=-polh_new(2,ii)
2438 IF (ibsa(iii)==0) cycle
2439 ENDIF
2440C
2441 jmax=j
2442 imax=ii
2443 ELSEIF (ipoly_f(6,jj)==i) THEN
2444 ii=ipoly_f(5,jj)
2445 IF (polh_new(2,ii)<0) THEN
2446 iii=-polh_new(2,ii)
2447 IF (ibsa(iii)==0) cycle
2448 ENDIF
2449C
2450 jmax=j
2451 imax=ii
2452 ENDIF
2453 areamax=area
2454 ENDIF
2455 ENDDO
2456 IF (jmax/=0) THEN
2457 jjmax=polh_new(2+jmax,i)
2458 rpoly_f(1,jjmax)=-one
2459 ELSE
2460 jmax=1
2461 jjmax=polh_new(2+jmax,i)
2462 ity=ipoly_f(1,jjmax)
2463 IF (ity==2) rpoly_f(1,jjmax)=-one
2464C On merge tous les polyhedres de volume nul avec le polyhedre 1
2465 imax=1
2466 ENDIF
2467 IF (imerged(imax)==1) imax=polh_new(2,imax)
2468 np=polh_new(1,imax)
2469 volu(imax)=volu(imax)+volu(i)
2470 volu(i)=-one
2471 imerged(i)=1
2472 DO j=1,polh_new(1,i)
2473 IF (j/=jmax) THEN
2474 np=np+1
2475 IF (np>nphmax) info=1
2476 IF (info==0) THEN
2477 jj=polh_new(2+j,i)
2478 polh_new(2+np,imax)=jj
2479 ity=ipoly_f(1,jj)
2480 IF (ity==2) THEN
2481 IF (ipoly_f(5,jj)==i) THEN
2482 ipoly_f(5,jj)=imax
2483 ELSEIF (ipoly_f(6,jj)==i) THEN
2484 ipoly_f(6,jj)=imax
2485 ENDIF
2486 ENDIF
2487 ELSE
2488 nphmax=max(nphmax,np)
2489 ENDIF
2490 ENDIF
2491 ENDDO
2492 polh_new(2,i)=imax
2493 IF (info==0) polh_new(1,imax)=np
2494 ENDIF
2495 ENDDO
2496 IF (info==1) GOTO 400
2497C
2498 vmin=ep30
2499 DO i=1,npolh
2500 IF (volu(i)<=zero) cycle
2501 IF (volu(i)<vmin) THEN
2502 vmin=volu(i)
2503 imin=i
2504 ENDIF
2505 DO j=1,polh_new(1,i)
2506 jj=polh_new(2+j,i)
2507 IF (ipoly_f(1,jj)==1.OR.rpoly_f(1,jj)<zero) cycle
2508 IF (ipoly_f(5,jj)==ipoly_f(6,jj)) rpoly_f(1,jj)=-one
2509 ENDDO
2510 ENDDO
2511 DEALLOCATE(volu_old, rpoly_f_old, ipoly_f_old, imerged)
2512C---------------------------------------------------
2513C Verifications
2514C---------------------------------------------------
2515 volph=zero
2516 areap=zero
2517 areael=zero
2518 npolhf=0
2519 DO i=1,npolh
2520 IF (volu(i)>zero) THEN
2521 npolhf=npolhf+1
2522 volph=volph+volu(i)
2523 ENDIF
2524 ENDDO
2525 nspoly=0
2526 ncpoly=0
2527 DO i=1,npoly
2528 ity=ipoly_f(1,i)
2529 IF (ity==1.AND.rpoly_f(1,i)>zero) THEN
2530 nspoly=nspoly+1
2531 areap=areap+rpoly_f(1,i)
2532 ELSEIF (ity==2.AND.rpoly_f(1,i)>zero) THEN
2533 ncpoly=ncpoly+1
2534 ENDIF
2535 ENDDO
2536 DO iel=1,nel
2537 areael=areael+elarea(iel)
2538 ENDDO
2539C---------------------------------------------------
2540C Coordonnees elementaires des nouveaux noeuds
2541C---------------------------------------------------
2542 ALLOCATE(fvdata(ifv)%IFVNOD(3,nns+nns2),
2543 . fvdata(ifv)%RFVNOD(2,nns+nns2))
2544 DO i=1,nns+nns2
2545 fvdata(ifv)%IFVNOD(1,i)=0
2546 ENDDO
2547C
2548 nns_old=nns
2549 nns=0
2550 DO i=1,npoly
2551 nnp=ipoly_f(2,i)
2552 DO j=1,nnp
2553 IF (ipoly_f(6+j,i)>0) THEN
2554 nns=nns+1
2555 IF (ifvnod(nns)<0) THEN
2556 fvdata(ifv)%IFVNOD(1,nns)=2
2557 fvdata(ifv)%IFVNOD(2,nns)=-ifvnod(nns)
2558 cycle
2559 ENDIF
2560 fvdata(ifv)%IFVNOD(1,nns)=1
2561 iel=ifvnod(nns)
2562 fvdata(ifv)%IFVNOD(2,nns)=iel
2563 xx=rpoly_f(4+3*(j-1)+1,i)
2564 yy=rpoly_f(4+3*(j-1)+2,i)
2565 zz=rpoly_f(4+3*(j-1)+3,i)
2566C
2567 n1=elema(1,iel)
2568 n2=elema(2,iel)
2569 n3=elema(3,iel)
2570 IF (tagela(iel)>0) THEN
2571 x1=xxx(1,n1)
2572 y1=xxx(2,n1)
2573 z1=xxx(3,n1)
2574 x2=xxx(1,n2)
2575 y2=xxx(2,n2)
2576 z2=xxx(3,n2)
2577 x3=xxx(1,n3)
2578 y3=xxx(2,n3)
2579 z3=xxx(3,n3)
2580 ELSEIF (tagela(iel)<0) THEN
2581 x1=xxxa(1,n1)
2582 y1=xxxa(2,n1)
2583 z1=xxxa(3,n1)
2584 x2=xxxa(1,n2)
2585 y2=xxxa(2,n2)
2586 z2=xxxa(3,n2)
2587 x3=xxxa(1,n3)
2588 y3=xxxa(2,n3)
2589 z3=xxxa(3,n3)
2590 ENDIF
2591C
2592 vpx=xx-x1
2593 vpy=yy-y1
2594 vpz=zz-z1
2595 v1x=x2-x1
2596 v1y=y2-y1
2597 v1z=z2-z1
2598 v2x=x3-x1
2599 v2y=y3-y1
2600 v2z=z3-z1
2601 CALL coorloc(vpx, vpy, vpz, v1x, v1y,
2602 . v1z, v2x, v2y, v2z, ksi,
2603 . eta)
2604C
2605 fvdata(ifv)%RFVNOD(1,nns)=ksi
2606 fvdata(ifv)%RFVNOD(2,nns)=eta
2607 ELSE
2608 jj=-ipoly_f(6+j,i)
2609 fvdata(ifv)%IFVNOD(1,nns_old+jj)=3
2610 fvdata(ifv)%IFVNOD(2,nns_old+jj)=ifvnod2(1,jj)
2611 fvdata(ifv)%IFVNOD(3,nns_old+jj)=ifvnod2(2,jj)
2612 fvdata(ifv)%RFVNOD(1,nns_old+jj)=rfvnod2(1,jj)
2613 ENDIF
2614 ENDDO
2615 ENDDO
2616 DO i=1,nns2
2617 IF (fvdata(ifv)%IFVNOD(1,nns+i)==0) THEN
2618 fvdata(ifv)%IFVNOD(1,nns+i)=3
2619 fvdata(ifv)%IFVNOD(2,nns+i)=1
2620 fvdata(ifv)%IFVNOD(3,nns+i)=2
2621 fvdata(ifv)%RFVNOD(1,nns+i)=one
2622 ENDIF
2623 ENDDO
2624C---------------------------------------------------
2625C Triangulation des polygones
2626C---------------------------------------------------
2627 nntr=0
2628 DO i=1,npoly
2629 nn=ipoly_f(2,i)
2630 nntr=nntr+nn
2631 ENDDO
2632 ALLOCATE(fvdata(ifv)%IFVTRI(6,nntr),
2633 . fvdata(ifv)%IFVPOLY(nntr),
2634 . fvdata(ifv)%IFVTADR(npoly+1))
2635C
2636 nntr=0
2637 npoly_old=npoly
2638 npoly=0
2639 ALLOCATE(redir_poly(npoly_old))
2640 DO i=1,npoly_old
2641 redir_poly(i)=0
2642 ENDDO
2643C
2644 nns=0
2645 DO i=1,npoly_old
2646 IF (rpoly_f(1,i)<=zero) THEN
2647 DO j=1,ipoly_f(2,i)
2648 IF (ipoly_f(6+j,i)>0) nns=nns+1
2649 ENDDO
2650 cycle
2651 ENDIF
2652 npoly=npoly+1
2653 redir_poly(i)=npoly
2654 nnp=ipoly_f(2,i)
2655 nhol=0
2656 IF (ipoly_f(1,i)==2) nhol=ipoly_f(6+nnp+1,i)
2657 ALLOCATE(pnodes(2,nnp), pseg(2,nnp), pholes(2,nhol),
2658 . ptri(3,nnp), redir(nnp))
2659 ptri(1:3,1:nnp)=0
2660C
2661 DO j=1,nnp
2662 IF (ipoly_f(6+j,i)>0) THEN
2663 nns=nns+1
2664 redir(j)=nns
2665 ELSE
2666 jj=-ipoly_f(6+j,i)
2667 redir(j)=nns_old+jj
2668 ENDIF
2669 ENDDO
2670 IF (ipoly_f(1,i)==1) THEN
2671 ipsurf=ipoly_f(3,i)
2672 ic1=0
2673 ic2=0
2674 ELSEIF (ipoly_f(1,i)==2) THEN
2675 ipsurf=0
2676 ic1=ipoly_f(5,i)
2677 ic2=ipoly_f(6,i)
2678 ENDIF
2679C
2680C Coordonnees des sommets dans le plan du polygone
2681 nx=rpoly_f(2,i)
2682 ny=rpoly_f(3,i)
2683 nz=rpoly_f(4,i)
2684C
2685 x0=rpoly_f(5,i)
2686 y0=rpoly_f(6,i)
2687 z0=rpoly_f(7,i)
2688 x1=rpoly_f(8,i)
2689 y1=rpoly_f(9,i)
2690 z1=rpoly_f(10,i)
2691 nrm1=sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
2692 vx1=(x1-x0)/nrm1
2693 vy1=(y1-y0)/nrm1
2694 vz1=(z1-z0)/nrm1
2695 vx2=ny*vz1-nz*vy1
2696 vy2=nz*vx1-nx*vz1
2697 vz2=nx*vy1-ny*vx1
2698C
2699 pnodes(1,1)=zero
2700 pnodes(2,1)=zero
2701 IF (nhol==0) THEN
2702 DO j=2,nnp
2703 xx=rpoly_f(4+3*(j-1)+1,i)
2704 yy=rpoly_f(4+3*(j-1)+2,i)
2705 zz=rpoly_f(4+3*(j-1)+3,i)
2706 vx=xx-x0
2707 vy=yy-y0
2708 vz=zz-z0
2709 pnodes(1,j)=vx*vx1+vy*vy1+vz*vz1
2710 pnodes(2,j)=vx*vx2+vy*vy2+vz*vz2
2711 pseg(1,j-1)=j-1
2712 pseg(2,j-1)=j
2713 ENDDO
2714 pseg(1,nnp)=nnp
2715 pseg(2,nnp)=1
2716 nseg=nnp
2717 ELSE
2718 ALLOCATE(adr(nhol+1))
2719 DO j=1,nhol
2720 adr(j)=ipoly_f(6+nnp+1+j,i)
2721 ENDDO
2722 adr(nhol+1)=nnp
2723 nseg=0
2724 DO j=2,adr(1)
2725 xx=rpoly_f(4+3*(j-1)+1,i)
2726 yy=rpoly_f(4+3*(j-1)+2,i)
2727 zz=rpoly_f(4+3*(j-1)+3,i)
2728 vx=xx-x0
2729 vy=yy-y0
2730 vz=zz-z0
2731 pnodes(1,j)=vx*vx1+vy*vy1+vz*vz1
2732 pnodes(2,j)=vx*vx2+vy*vy2+vz*vz2
2733 nseg=nseg+1
2734 pseg(1,nseg)=j-1
2735 pseg(2,nseg)=j
2736 ENDDO
2737 nseg=nseg+1
2738 pseg(1,nseg)=adr(1)
2739 pseg(2,nseg)=1
2740 DO j=1,nhol
2741 xx=rpoly_f(4+3*adr(j)+1,i)
2742 yy=rpoly_f(4+3*adr(j)+2,i)
2743 zz=rpoly_f(4+3*adr(j)+3,i)
2744 vx=xx-x0
2745 vy=yy-y0
2746 vz=zz-z0
2747 pnodes(1,adr(j)+1)=vx*vx1+vy*vy1+vz*vz1
2748 pnodes(2,adr(j)+1)=vx*vx2+vy*vy2+vz*vz2
2749 DO k=adr(j)+2,adr(j+1)
2750 xx=rpoly_f(4+3*(k-1)+1,i)
2751 yy=rpoly_f(4+3*(k-1)+2,i)
2752 zz=rpoly_f(4+3*(k-1)+3,i)
2753 vx=xx-x0
2754 vy=yy-y0
2755 vz=zz-z0
2756 pnodes(1,k)=vx*vx1+vy*vy1+vz*vz1
2757 pnodes(2,k)=vx*vx2+vy*vy2+vz*vz2
2758 nseg=nseg+1
2759 pseg(1,nseg)=k-1
2760 pseg(2,nseg)=k
2761 ENDDO
2762 nseg=nseg+1
2763 pseg(1,nseg)=adr(j+1)
2764 pseg(2,nseg)=adr(j)+1
2765C
2766 xx=rpoly_f(4+3*nnp+3*(j-1)+1,i)
2767 yy=rpoly_f(4+3*nnp+3*(j-1)+2,i)
2768 zz=rpoly_f(4+3*nnp+3*(j-1)+3,i)
2769 vx=xx-x0
2770 vy=yy-y0
2771 vz=zz-z0
2772 pholes(1,j)=vx*vx1+vy*vy1+vz*vz1
2773 pholes(2,j)=vx*vx2+vy*vy2+vz*vz2
2774 ENDDO
2775 ENDIF
2776C
2777 nelp = 0
2778 CALL c_tricall(pnodes, pseg, pholes, ptri, nnp,
2779 . nseg, nhol, nelp )
2780C
2781 fvdata(ifv)%IFVTADR(npoly)=nntr+1
2782 DO j=1,nelp
2783 nntr=nntr+1
2784 n1=ptri(1,j)
2785 n2=ptri(2,j)
2786 n3=ptri(3,j)
2787C Verification de la normale
2788 x1=rpoly_f(4+3*(n1-1)+1,i)
2789 y1=rpoly_f(4+3*(n1-1)+2,i)
2790 z1=rpoly_f(4+3*(n1-1)+3,i)
2791 x2=rpoly_f(4+3*(n2-1)+1,i)
2792 y2=rpoly_f(4+3*(n2-1)+2,i)
2793 z2=rpoly_f(4+3*(n2-1)+3,i)
2794 x3=rpoly_f(4+3*(n3-1)+1,i)
2795 y3=rpoly_f(4+3*(n3-1)+2,i)
2796 z3=rpoly_f(4+3*(n3-1)+3,i)
2797 x12=x2-x1
2798 y12=y2-y1
2799 z12=z2-z1
2800 x13=x3-x1
2801 y13=y3-y1
2802 z13=z3-z1
2803 nrx=y12*z13-z12*y13
2804 nry=z12*x13-x12*z13
2805 nrz=x12*y13-y12*x13
2806 ss=nrx*nx+nry*ny+nrz*nz
2807C
2808 IF (ss>zero) THEN
2809 fvdata(ifv)%IFVTRI(1,nntr)=redir(n1)
2810 fvdata(ifv)%IFVTRI(2,nntr)=redir(n2)
2811 fvdata(ifv)%IFVTRI(3,nntr)=redir(n3)
2812 ELSE
2813 fvdata(ifv)%IFVTRI(1,nntr)=redir(n1)
2814 fvdata(ifv)%IFVTRI(2,nntr)=redir(n3)
2815 fvdata(ifv)%IFVTRI(3,nntr)=redir(n2)
2816 ENDIF
2817 fvdata(ifv)%IFVTRI(4,nntr)=ipsurf
2818 fvdata(ifv)%IFVTRI(5,nntr)=ic1
2819 fvdata(ifv)%IFVTRI(6,nntr)=ic2
2820 fvdata(ifv)%IFVPOLY(nntr)=nntr
2821 ENDDO
2822C
2823 DEALLOCATE(pnodes, pseg, pholes, ptri, redir)
2824 ENDDO
2825 fvdata(ifv)%IFVTADR(npoly+1)=nntr+1
2826C---------------------------------------------------
2827C Structure des polyhedres
2828C---------------------------------------------------
2829 ALLOCATE(fvdata(ifv)%IFVPOLH(2*npoly),
2830 . fvdata(ifv)%IFVPADR(npolh+1),
2831 . fvdata(ifv)%IDPOLH(npolh),
2832 . fvdata(ifv)%IBPOLH(npolh))
2833 nnp=0
2834 npolh_old=npolh
2835 npolh=0
2836 ALLOCATE(redir_polh(npolh_old))
2837 DO i=1,npolh_old
2838 redir_polh(i)=0
2839 ENDDO
2840C
2841 DO i=1,npolh_old
2842 IF (volu(i)<=zero) cycle
2843 npolh=npolh+1
2844 redir_polh(i)=npolh
2845 fvdata(ifv)%IFVPADR(npolh)=nnp+1
2846 DO j=1,polh_new(1,i)
2847 jj=redir_poly(polh_new(2+j,i))
2848 IF (jj==0) cycle
2849 nnp=nnp+1
2850 fvdata(ifv)%IFVPOLH(nnp)=redir_poly(polh_new(2+j,i))
2851 ENDDO
2852C
2853 fvdata(ifv)%IDPOLH(npolh)=npolh
2854C
2855 ii=polh_new(2,i)
2856 IF (ii>=0) ii=0
2857 IF (ii<0.AND.ibsa(-ii)==1) ii=-ii
2858 fvdata(ifv)%IBPOLH(npolh)=ii
2859 ENDDO
2860 fvdata(ifv)%IFVPADR(npolh+1)=nnp+1
2861C
2862 DO i=1,nntr
2863 IF (fvdata(ifv)%IFVTRI(4,i)==0) THEN
2864 ic1=fvdata(ifv)%IFVTRI(5,i)
2865 ic2=fvdata(ifv)%IFVTRI(6,i)
2866 fvdata(ifv)%IFVTRI(5,i)=redir_polh(ic1)
2867 fvdata(ifv)%IFVTRI(6,i)=redir_polh(ic2)
2868 ENDIF
2869 ENDDO
2870C
2871 ivolu(46)=nns+nns2
2872 ivolu(47)=nntr
2873 ivolu(48)=npoly
2874 ivolu(49)=npolh
2875C
2876 fvdata(ifv)%NNS=nns+nns2
2877 fvdata(ifv)%NNTR=nntr
2878 fvdata(ifv)%LENP=fvdata(ifv)%IFVTADR(npoly+1)-1
2879 fvdata(ifv)%NPOLY=npoly
2880 fvdata(ifv)%LENH=fvdata(ifv)%IFVPADR(npolh+1)-1
2881 fvdata(ifv)%NPOLH=npolh
2882C
2883 ALLOCATE(fvdata(ifv)%MPOLH(npolh), fvdata(ifv)%QPOLH(3,npolh),
2884 . fvdata(ifv)%EPOLH(npolh), fvdata(ifv)%PPOLH(npolh),
2885 . fvdata(ifv)%RPOLH(npolh), fvdata(ifv)%GPOLH(npolh),
2886 . fvdata(ifv)%CPAPOLH(npolh), fvdata(ifv)%CPBPOLH(npolh),
2887 . fvdata(ifv)%CPCPOLH(npolh), fvdata(ifv)%RMWPOLH(npolh),
2888 . fvdata(ifv)%VPOLH_INI(npolh), fvdata(ifv)%DTPOLH(npolh),
2889 . fvdata(ifv)%TPOLH(npolh), fvdata(ifv)%CPDPOLH(npolh),
2890 . fvdata(ifv)%CPEPOLH(npolh), fvdata(ifv)%CPFPOLH(npolh))
2891C
2892 DO i=1,npolh_old
2893 ii=redir_polh(i)
2894 IF (ii==0) cycle
2895 fvdata(ifv)%VPOLH_INI(ii)=volu(i)
2896 ENDDO
2897C
2898 DEALLOCATE(ipoly_f, rpoly_f, polh_f, volu, redir_poly, redir_polh,
2899 . polh_new, ifvnod, ifvnod2, rfvnod2, xns, xns2)
2900C---------------------------------------------------
2901C Impressions
2902C---------------------------------------------------
2903 nstr=0
2904 nctr=0
2905 DO i=1,nntr
2906 IF (fvdata(ifv)%IFVTRI(4,i)>0) THEN
2907 nstr=nstr+1
2908 ELSE
2909 nctr=nctr+1
2910 ENDIF
2911 ENDDO
2912C
2913 WRITE(istdo,*)
2914 WRITE(istdo,'(A25,I10,A24)')
2915 .' ** MONITORED VOLUME ID: ',ivolu(1),' - FINITE VOLUME MESH **'
2916 WRITE(istdo,'(A42,I8)')
2917 . ' NUMBER OF SURFACE POLYGONS : ',nspoly
2918 WRITE(istdo,'(A42,I8)')
2919 . ' NUMBER OF SURFACE TRIANGLES : ',nstr
2920 WRITE(istdo,'(A42,I8)')
2921 . ' NUMBER OF COMMUNICATION POLYGONS : ',ncpoly
2922 WRITE(istdo,'(A42,I8)')
2923 . ' NUMBER OF COMMUNICATION TRIANGLES : ',nctr
2924 WRITE(istdo,'(A42,I8)')
2925 . ' NUMBER OF POLYHEDRA (FINITE VOLUMES): ',npolhf
2926 IF (ilvout>=1) WRITE(istdo,'(A29,G16.9,A8,I8,A1)')
2927 . ' MIN. POLYHEDRON VOLUME : ',vmin,' (POLY. ',imin,')'
2928 IF (ilvout>=1) WRITE(istdo,'(A29,G16.9)')
2929 . ' INITIAL MERGING VOLUME : ',volumin
2930 WRITE(istdo,'(A29,G16.9,A17,G16.9,A1)')
2931 .' SUM VOLUME POLYHEDRA : ',volph,' (VOLUME AIRBAG: ',volg,')'
2932 WRITE(istdo,'(A29,G16.9,A17,G16.9,A1)')
2933 .' SUM AREA SURF. POLYGONS: ',areap,
2934 . ' (AREA AIRBAG : ',areael,')'
2935 WRITE(istdo,*)
2936 WRITE(iout,*)
2937 WRITE(iout,'(A25,I10,A24)')
2938 .' ** MONITORED VOLUME ID: ',ivolu(1),' - FINITE VOLUME MESH **'
2939 WRITE(iout,'(A42,I8)')
2940 . ' NUMBER OF SURFACE POLYGONS : ',nspoly
2941 WRITE(iout,'(A42,I8)')
2942 . ' NUMBER OF COMMUNICATION POLYGONS : ',ncpoly
2943 WRITE(iout,'(A42,I8)')
2944 . ' NUMBER OF POLYHEDRA (FINITE VOLUMES): ',npolhf
2945 IF (ilvout>=1) WRITE(iout,'(A29,G16.9,A8,I8,A1)')
2946 . ' MIN. POLYHEDRON VOLUME : ',vmin,' (POLY. ',imin,')'
2947 IF (ilvout>=1) WRITE(iout,'(A29,G16.9)')
2948 . ' INITIAL MERGING VOLUME : ',volumin
2949 WRITE(iout,'(A29,G16.9,A17,G16.9,A1)')
2950 .' SUM VOLUME POLYHEDRA : ',volph,' (VOLUME AIRBAG: ',volg,')'
2951 WRITE(iout,'(A29,G16.9,A17,G16.9,A1)')
2952 .' SUM AREA SURF. POLYGONS: ',areap,
2953 . ' (AREA AIRBAG : ',areael,')'
2954 WRITE(iout,*)
2955C---------------------------------------------------
2956C Anim des volumes finis
2957C---------------------------------------------------
2958 fvdata(ifv)%NPOLH_ANIM=npolh
2959 lenp=fvdata(ifv)%IFVTADR(npoly+1)
2960 lenh=fvdata(ifv)%IFVPADR(npolh+1)
2961 ALLOCATE(fvdata(ifv)%IFVPOLY_ANIM(lenp),
2962 . fvdata(ifv)%IFVTADR_ANIM(npoly+1),
2963 . fvdata(ifv)%IFVPOLH_ANIM(lenh),
2964 . fvdata(ifv)%IFVPADR_ANIM(npolh+1),
2965 . fvdata(ifv)%IFVTRI_ANIM(6,nntr),
2966 . fvdata(ifv)%REDIR_ANIM(nns+nns2),
2967 . fvdata(ifv)%NOD_ANIM(3,nns+nns2),
2968 . redir(nns+nns2), itagt(nntr))
2969 DO i=1,lenp
2970 fvdata(ifv)%IFVPOLY_ANIM(i)=fvdata(ifv)%IFVPOLY(i)
2971 ENDDO
2972 DO i=1,npoly+1
2973 fvdata(ifv)%IFVTADR_ANIM(i)=fvdata(ifv)%IFVTADR(i)
2974 ENDDO
2975 DO i=1,lenh
2976 fvdata(ifv)%IFVPOLH_ANIM(i)=fvdata(ifv)%IFVPOLH(i)
2977 ENDDO
2978 DO i=1,npolh+1
2979 fvdata(ifv)%IFVPADR_ANIM(i)=fvdata(ifv)%IFVPADR(i)
2980 ENDDO
2981C Elimination des noeuds doubles pour anim
2982 tole=em05*fac_length
2983 tole2=tole**2
2984 nns_anim=0
2985 IF (ilvout/=0) WRITE(istdo,'(A25,I10,A39)')
2986 . ' ** MONITORED VOLUME ID: ',ivolu(1),
2987 . ' - MERGING COINCIDENT NODES FOR ANIM **'
2988 ALLOCATE(pnodes(3,nns+nns2))
2989 DO i=1,nns
2990 IF (fvdata(ifv)%IFVNOD(1,i)==1) THEN
2991 iel=fvdata(ifv)%IFVNOD(2,i)
2992 ksi=fvdata(ifv)%RFVNOD(1,i)
2993 eta=fvdata(ifv)%RFVNOD(2,i)
2994 n1=elema(1,iel)
2995 n2=elema(2,iel)
2996 n3=elema(3,iel)
2997 IF (tagela(iel)>0) THEN
2998 x1=xxx(1,n1)
2999 y1=xxx(2,n1)
3000 z1=xxx(3,n1)
3001 x2=xxx(1,n2)
3002 y2=xxx(2,n2)
3003 z2=xxx(3,n2)
3004 x3=xxx(1,n3)
3005 y3=xxx(2,n3)
3006 z3=xxx(3,n3)
3007 ELSEIF (tagela(iel)<0) THEN
3008 x1=xxxa(1,n1)
3009 y1=xxxa(2,n1)
3010 z1=xxxa(3,n1)
3011 x2=xxxa(1,n2)
3012 y2=xxxa(2,n2)
3013 z2=xxxa(3,n2)
3014 x3=xxxa(1,n3)
3015 y3=xxxa(2,n3)
3016 z3=xxxa(3,n3)
3017 ENDIF
3018 pnodes(1,i)=(one-ksi-eta)*x1+ksi*x2+eta*x3
3019 pnodes(2,i)=(one-ksi-eta)*y1+ksi*y2+eta*y3
3020 pnodes(3,i)=(one-ksi-eta)*z1+ksi*z2+eta*z3
3021 ELSEIF (fvdata(ifv)%IFVNOD(1,i)==2) THEN
3022c IF (NSPMD == 1) THEN
3023c II=FVDATA(IFV)%IFVNOD(2,I)
3024c PNODES(1,I)=X(1,II)
3025c PNODES(2,I)=X(2,II)
3026c PNODES(3,I)=X(3,II)
3027c ELSE
3028 ii=fvdata(ifv)%IFVNOD(2,i)
3029 pnodes(1,i)=xxxsa(1,ii)
3030 pnodes(2,i)=xxxsa(2,ii)
3031 pnodes(3,i)=xxxsa(3,ii)
3032c ENDIF
3033 ENDIF
3034 ENDDO
3035 DO i=1,nns2
3036 ii=nns+i
3037 i1=fvdata(ifv)%IFVNOD(2,ii)
3038 i2=fvdata(ifv)%IFVNOD(3,ii)
3039 alpha=fvdata(ifv)%RFVNOD(1,ii)
3040 pnodes(1,ii)=alpha*pnodes(1,i1)+(one-alpha)*pnodes(1,i2)
3041 pnodes(2,ii)=alpha*pnodes(2,i1)+(one-alpha)*pnodes(2,i2)
3042 pnodes(3,ii)=alpha*pnodes(3,i1)+(one-alpha)*pnodes(3,i2)
3043 ENDDO
3044 DO i=1,nns+nns2
3045 IF (ilvout<0) CALL progbar_c(i,nns+nns2)
3046C
3047 xx=pnodes(1,i)
3048 yy=pnodes(2,i)
3049 zz=pnodes(3,i)
3050 ifound=0
3051 DO j=1,nns_anim
3052 xn(1)=fvdata(ifv)%NOD_ANIM(1,j)
3053 xn(2)=fvdata(ifv)%NOD_ANIM(2,j)
3054 xn(3)=fvdata(ifv)%NOD_ANIM(3,j)
3055 dd2=(xx-xn(1))**2+(yy-xn(2))**2+(zz-xn(3))**2
3056 IF (dd2<=tole2) THEN
3057 ifound=j
3058 EXIT
3059 ENDIF
3060 ENDDO
3061 IF (ifound==0) THEN
3062 nns_anim=nns_anim+1
3063 redir(i)=nns_anim
3064 fvdata(ifv)%REDIR_ANIM(nns_anim)=i
3065 fvdata(ifv)%NOD_ANIM(1,nns_anim)=xx
3066 fvdata(ifv)%NOD_ANIM(2,nns_anim)=yy
3067 fvdata(ifv)%NOD_ANIM(3,nns_anim)=zz
3068 ELSE
3069 redir(i)=ifound
3070 ENDIF
3071 ENDDO
3072 fvdata(ifv)%NNS_ANIM=nns_anim
3073 fvdata(ifv)%ID=ivolu(1)
3074C
3075 DO i=1,nntr
3076 n1=fvdata(ifv)%IFVTRI(1,i)
3077 n2=fvdata(ifv)%IFVTRI(2,i)
3078 n3=fvdata(ifv)%IFVTRI(3,i)
3079 fvdata(ifv)%IFVTRI_ANIM(1,i)=redir(n1)
3080 fvdata(ifv)%IFVTRI_ANIM(2,i)=redir(n2)
3081 fvdata(ifv)%IFVTRI_ANIM(3,i)=redir(n3)
3082 fvdata(ifv)%IFVTRI_ANIM(4,i)=
3083 . fvdata(ifv)%IFVTRI(4,i)
3084 fvdata(ifv)%IFVTRI_ANIM(5,i)=
3085 . fvdata(ifv)%IFVTRI(5,i)
3086 fvdata(ifv)%IFVTRI_ANIM(6,i)=
3087 . fvdata(ifv)%IFVTRI(6,i)
3088 ENDDO
3089 DEALLOCATE(pnodes)
3090C
3091 DEALLOCATE(redir, itagt)
3092 DEALLOCATE(xxxsa)
3093C
3094 RETURN
3095C
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
horipoly
subroutine horipoly(inedge, rnedge, ledge, nedge, ipoly, rpoly, iz, ielnod, npoly, nx, ny, nz, inz, ibric)
Definition fvmesh.F:4984
polyhedr
subroutine polyhedr(ipoly, rpoly, polb, npolb, polh, npolh, nrpmax, nphmax, ibric, lmin, info, npolhmax, nppmax, tole2)
Definition fvmesh.F:4119
horiedge
subroutine horiedge(ipoly, rpoly, nx, ny, nz, nbnedge, inedge, rnedge, x0, y0, z0, inz, nns3, nref, aref, nnsp)
Definition fvmesh.F:4906
facepoly
subroutine facepoly(quad, tri, ntri, ipoly, rpoly, n, normf, normt, nsmax, nnp, nrpmax, nb, nv, lmin, npolmax, nppmax, info, ielnod, xxx, elema, ibuf, nela, ibric, ifac, mvid, ilvout, ibufa, tagela, xxxa)
Definition fvmesh.F:3257
coorloc
subroutine coorloc(vpx, vpy, vpz, v1x, v1y, v1z, v2x, v2y, v2z, ksi, eta)
Definition fvmesh.F:4278
gpolcut
subroutine gpolcut(ipoly, rpoly, ipoly_old, rpoly_old, inedge, rnedge, nbnedge, nx, ny, nz, x0, y0, z0, ins, rns, nn, nhol, inz, iz, nns3, npoly, ns, ielnod, ielnod_old)
Definition fvmesh.F:4319
itribox
subroutine itribox(tri, box, norm, nverts, poly, nvmax)
Definition fvmesh.F:3107
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
min
#define min(a, b)
Definition macros.h:20
fvbag_mod
Definition fvbag_mod.F:44
fvbag_mod::fvspmd
type(fvbag_spmd), dimension(:), allocatable fvspmd
Definition fvbag_mod.F:129
fvbag_mod::fvdata
type(fvbag_data), dimension(:), allocatable fvdata
Definition fvbag_mod.F:128
progbar_c
void progbar_c(int *icur, int *imax)
Definition progbar_c.c:30
spmd_fvb_gath
subroutine spmd_fvb_gath(ifv, x, xxx, xxxa, xxxsa, ido)
Definition spmd_fvb_gath.F:40
facnor
subroutine facnor(x, d, ii, xnorm, cdg, invert)
Definition facnor.F:30
c_tricall
subroutine c_tricall(pnodes, pseg, pholes, ptri, nnp, nseg, nhol, nelp)
Definition fvmbags_stub.F:34
tritri3
subroutine tritri3(icut, fv0, fv1, fv2, fu0, fu1, fu2)
Definition fvmbags_stub.F:74
tribox3
subroutine tribox3(icut, bcenter, bhalfsize, tverts)
Definition fvmbags_stub.F:55

◆ gpolcut()

subroutine gpolcut ( integer, dimension(6+2*nn+1+nhol,*) ipoly,
rpoly,
integer, dimension(*) ipoly_old,
rpoly_old,
integer, dimension(6,*) inedge,
rnedge,
integer nbnedge,
nx,
ny,
nz,
x0,
y0,
z0,
integer, dimension(2,*) ins,
rns,
integer nn,
integer nhol,
integer inz,
integer, dimension(3,*) iz,
integer nns3,
integer npoly,
integer ns,
integer, dimension(2*nn,*) ielnod,
integer, dimension(*) ielnod_old )

Definition at line 4314 of file fvmesh.F.

4319C-----------------------------------------------
4320C I m p l i c i t T y p e s
4321C-----------------------------------------------
4322#include "implicit_f.inc"
4323C-----------------------------------------------
4324C D u m m y A r g u m e n t s
4325C-----------------------------------------------
4326 INTEGER NN, NHOL, IPOLY(6+2*NN+1+NHOL,*), IPOLY_OLD(*),
4327 . INEDGE(6,*), NBNEDGE, INS(2,*), INZ, IZ(3,*), NNS3,
4328 . NPOLY, NS, IELNOD(2*NN,*), IELNOD_OLD(*)
4329 my_real
4330 . rpoly(4+3*2*nn+3*nhol,*), rpoly_old(*), rnedge(6,*), nx,
4331 . ny, nz, x0, y0, z0, rns(4,*)
4332C-----------------------------------------------
4333C L o c a l V a r i a b l e s
4334C-----------------------------------------------
4335 INTEGER NN1, I, IADHOL(NHOL+1), NSEG, ITAG(2*NN+1), II,
4336 . TSEG(3,2*NN), NHOL_OLD, ICUT, J, JJ, NSEG_INI, NSEG1,
4337 . REDIR(NN), N1, N2, ISTOP, NNP, ICLOSE, ISEG,
4338 . POLY(2*NN,NN), IN, IN1, IN2, LENPOLY(NN), K, KK,
4339 . THOL(NHOL), JN, NHOLP, IADHOLP(NHOL), KN, I1, I2,
4340 . ITAG2(2*NN), ITAGN(NN), ITAGS(NN), ISEG_OLD
4341 my_real
4342 . x1, y1, z1, x2, y2, z2, zl1, zl2, alpha, npx, npy, npz,
4343 . xp0, yp0, zp0, vx, vy, vz, xx, yy, zz, xl(nn), xlmin,
4344 . xlc, xx1, yy1, zz1, xx2, yy2, zz2, vx1, vy1,
4345 . vz1, vx2, vy2, vz2, nr1, nr2, ss, vvx, vvy, vvz,
4346 . ss1, zl, xns(3,nn), tole, ll, zlm
4347C
4348C
4349 tole = zero
4350 tole = epsilon(tole)
4351C
4352C Liste des segments
4353 IF (nhol==0) THEN
4354 nn1=nn
4355 ELSE
4356 nn1=ipoly_old(6+nn+1+1)-1
4357 DO i=1,nhol
4358 iadhol(i)=ipoly_old(6+nn+1+i)
4359 ENDDO
4360 iadhol(nhol+1)=nn+1
4361 ENDIF
4362 ns=0
4363 nseg=0
4364C
4365 DO i=1,2*nn
4366 itag(i)=0
4367 itag2(i)=0
4368 ENDDO
4369 itag(2*nn+1)=0
4370 DO i=1,nn
4371 itagn(i)=0
4372 itags(i)=0
4373 ENDDO
4374C Segments du contour exterieur
4375 DO i=1,nn1
4376 ii=i+1
4377 IF (i==nn1) ii=1
4378 x1=rpoly_old(4+3*(i-1)+1)
4379 y1=rpoly_old(4+3*(i-1)+2)
4380 z1=rpoly_old(4+3*(i-1)+3)
4381 x2=rpoly_old(4+3*(ii-1)+1)
4382 y2=rpoly_old(4+3*(ii-1)+2)
4383 z2=rpoly_old(4+3*(ii-1)+3)
4384 zl1=(x1-x0)*nx+(y1-y0)*ny+(z1-z0)*nz
4385 zl2=(x2-x0)*nx+(y2-y0)*ny+(z2-z0)*nz
4386 IF (zl1-zl2/=zero) THEN
4387 alpha=zl2/(zl2-zl1)
4388 IF (alpha>zero.AND.alpha<one) THEN
4389 ns=ns+1
4390 ins(1,ns)=ipoly_old(6+i)
4391 ins(2,ns)=ipoly_old(6+ii)
4392 rns(1,ns)=alpha
4393 xns(1,ns)=alpha*x1+(one-alpha)*x2
4394 xns(2,ns)=alpha*y1+(one-alpha)*y2
4395 xns(3,ns)=alpha*z1+(one-alpha)*z2
4396 nseg=nseg+1
4397 tseg(1,nseg)=i
4398 tseg(2,nseg)=-ns
4399 tseg(3,nseg)=nseg+1
4400 nseg=nseg+1
4401 tseg(1,nseg)=-ns
4402 tseg(2,nseg)=ii
4403 tseg(3,nseg)=nseg+1
4404 ELSEIF (alpha==zero) THEN
4405 IF (itagn(ii)==0) THEN
4406 ns=ns+1
4407 itagn(ii)=ns
4408 ins(1,ns)=ipoly_old(6+i)
4409 ins(2,ns)=ipoly_old(6+ii)
4410 rns(1,ns)=zero
4411 xns(1,ns)=x2
4412 xns(2,ns)=y2
4413 xns(3,ns)=z2
4414 itags(ns)=1
4415 nseg=nseg+1
4416 tseg(1,nseg)=i
4417 tseg(2,nseg)=-ns
4418 tseg(3,nseg)=nseg+1
4419 nseg=nseg+1
4420 tseg(1,nseg)=-ns
4421 tseg(2,nseg)=ii
4422 tseg(3,nseg)=nseg+1
4423 ELSE
4424 nseg=nseg+1
4425 tseg(1,nseg)=i
4426 tseg(2,nseg)=ii
4427 tseg(3,nseg)=nseg+1
4428 ENDIF
4429 ELSEIF (alpha==one) THEN
4430 IF (itagn(i)==0) THEN
4431 ns=ns+1
4432 itagn(i)=ns
4433 ins(1,ns)=ipoly_old(6+i)
4434 ins(2,ns)=ipoly_old(6+ii)
4435 rns(1,ns)=one
4436 xns(1,ns)=x1
4437 xns(2,ns)=y1
4438 xns(3,ns)=z1
4439 itags(ns)=1
4440 nseg=nseg+1
4441 tseg(1,nseg)=i
4442 tseg(2,nseg)=-ns
4443 tseg(3,nseg)=nseg+1
4444 nseg=nseg+1
4445 tseg(1,nseg)=-ns
4446 tseg(2,nseg)=ii
4447 tseg(3,nseg)=nseg+1
4448 ELSE
4449 nseg=nseg+1
4450 tseg(1,nseg)=i
4451 tseg(2,nseg)=ii
4452 tseg(3,nseg)=nseg+1
4453 ENDIF
4454 ELSE
4455 nseg=nseg+1
4456 tseg(1,nseg)=i
4457 tseg(2,nseg)=ii
4458 tseg(3,nseg)=nseg+1
4459 ENDIF
4460 ELSE
4461 nseg=nseg+1
4462 tseg(1,nseg)=i
4463 tseg(2,nseg)=ii
4464 tseg(3,nseg)=nseg+1
4465 ENDIF
4466 ENDDO
4467 tseg(3,nseg)=1
4468C Segments des trous
4469 nhol_old=nhol
4470 nhol=0
4471 DO i=1,nhol_old
4472 icut=0
4473 nseg_ini=nseg
4474 DO j=iadhol(i),iadhol(i+1)-1
4475 jj=j+1
4476 IF (j==(iadhol(i+1)-1)) jj=iadhol(i)
4477 x1=rpoly_old(4+3*(j-1)+1)
4478 y1=rpoly_old(4+3*(j-1)+2)
4479 z1=rpoly_old(4+3*(j-1)+3)
4480 x2=rpoly_old(4+3*(jj-1)+1)
4481 y2=rpoly_old(4+3*(jj-1)+2)
4482 z2=rpoly_old(4+3*(jj-1)+3)
4483 zl1=(x1-x0)*nx+(y1-y0)*ny+(z1-z0)*nz
4484 zl2=(x2-x0)*nx+(y2-y0)*ny+(z2-z0)*nz
4485 IF (zl1-zl2/=zero) THEN
4486 alpha=zl2/(zl2-zl1)
4487 IF (alpha>zero.AND.alpha<one) THEN
4488 icut=1
4489 ns=ns+1
4490 ins(1,ns)=ipoly_old(6+j)
4491 ins(2,ns)=ipoly_old(6+jj)
4492 rns(1,ns)=alpha
4493 xns(1,ns)=alpha*x1+(one-alpha)*x2
4494 xns(2,ns)=alpha*y1+(one-alpha)*y2
4495 xns(3,ns)=alpha*z1+(one-alpha)*z2
4496 nseg=nseg+1
4497 tseg(1,nseg)=j
4498 tseg(2,nseg)=-ns
4499 tseg(3,nseg)=nseg+1
4500 nseg=nseg+1
4501 tseg(1,nseg)=-ns
4502 tseg(2,nseg)=jj
4503 tseg(3,nseg)=nseg+1
4504 ELSEIF (alpha==zero) THEN
4505 icut=1
4506 IF (itagn(ii)==0) THEN
4507 ns=ns+1
4508 itagn(jj)=ns
4509 ins(1,ns)=ipoly_old(6+j)
4510 ins(2,ns)=ipoly_old(6+jj)
4511 rns(1,ns)=zero
4512 xns(1,ns)=x2
4513 xns(2,ns)=y2
4514 xns(3,ns)=z2
4515 itags(ns)=1
4516 nseg=nseg+1
4517 tseg(1,nseg)=j
4518 tseg(2,nseg)=-ns
4519 tseg(3,nseg)=nseg+1
4520 nseg=nseg+1
4521 tseg(1,nseg)=-ns
4522 tseg(2,nseg)=jj
4523 tseg(3,nseg)=nseg+1
4524 ELSE
4525 nseg=nseg+1
4526 tseg(1,nseg)=j
4527 tseg(2,nseg)=jj
4528 tseg(3,nseg)=nseg+1
4529 ENDIF
4530 ELSEIF (alpha==one) THEN
4531 icut=1
4532 IF (itagn(i)==0) THEN
4533 ns=ns+1
4534 itagn(j)=ns
4535 ins(1,ns)=ipoly_old(6+j)
4536 ins(2,ns)=ipoly_old(6+jj)
4537 rns(1,ns)=one
4538 xns(1,ns)=x1
4539 xns(2,ns)=y1
4540 xns(3,ns)=z1
4541 itags(ns)=1
4542 nseg=nseg+1
4543 tseg(1,nseg)=j
4544 tseg(2,nseg)=-ns
4545 tseg(3,nseg)=nseg+1
4546 nseg=nseg+1
4547 tseg(1,nseg)=-ns
4548 tseg(2,nseg)=jj
4549 tseg(3,nseg)=nseg+1
4550 ELSE
4551 nseg=nseg+1
4552 tseg(1,nseg)=j
4553 tseg(2,nseg)=jj
4554 tseg(3,nseg)=nseg+1
4555 ENDIF
4556 ELSE
4557 nseg=nseg+1
4558 tseg(1,nseg)=j
4559 tseg(2,nseg)=jj
4560 tseg(3,nseg)=nseg+1
4561 ENDIF
4562 ELSE
4563 nseg=nseg+1
4564 tseg(1,nseg)=j
4565 tseg(2,nseg)=jj
4566 tseg(3,nseg)=nseg+1
4567 ENDIF
4568 ENDDO
4569 tseg(3,nseg)=nseg_ini+1
4570C
4571 IF (icut==0) THEN
4572C Le trou n'est pas coupe
4573 nhol=nhol+1
4574 DO j=nseg_ini+1,nseg
4575 itag(j)=-nhol
4576 ENDDO
4577 ENDIF
4578 ENDDO
4579C
4580 nseg1=nseg
4581C Creation des nouvelles aretes dans le plan horizontal
4582 npx=rpoly_old(2)
4583 npy=rpoly_old(3)
4584 npz=rpoly_old(4)
4585 xp0=rpoly_old(5)
4586 yp0=rpoly_old(6)
4587 zp0=rpoly_old(7)
4588 vx=ny*npz-nz*npy
4589 vy=nz*npx-nx*npz
4590 vz=nx*npy-ny*npx
4591 DO i=1,ns
4592 xx=xns(1,i)
4593 yy=xns(2,i)
4594 zz=xns(3,i)
4595 xl(i)=(xx-xp0)*vx+(yy-yp0)*vy+(zz-zp0)*vz
4596 ENDDO
4597 xlmin=ep30
4598 DO i=1,ns
4599 redir(i)=i
4600 ENDDO
4601 DO i=1,ns
4602 xlmin=xl(redir(i))
4603 DO j=i+1,ns
4604 xlc=xl(redir(j))
4605 IF (xlc<xlmin) THEN
4606 jj=redir(j)
4607 redir(j)=redir(i)
4608 redir(i)=jj
4609 xlmin=xlc
4610 ENDIF
4611 ENDDO
4612 ENDDO
4613 ii=0
4614 DO i=1,ns
4615 ii=ii+1
4616 IF (ii==1) THEN
4617 n1=redir(i)
4618 n2=redir(i+1)
4619 x1=xns(1,n1)
4620 y1=xns(2,n1)
4621 z1=xns(3,n1)
4622 x2=xns(1,n2)
4623 y2=xns(2,n2)
4624 z2=xns(3,n2)
4625 ll=sqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)
4626C
4627 IF (ll>tole) THEN
4628 nseg=nseg+1
4629 tseg(1,nseg)=-n1
4630 tseg(2,nseg)=-n2
4631 tseg(3,nseg)=-1
4632 ENDIF
4633 ELSE
4634 ii=0
4635 ENDIF
4636 ENDDO
4637C Nouveaux polygones
4638 istop=0
4639 npoly=1
4640 nnp=0
4641 DO WHILE (istop==0)
4642 i=1
4643 DO WHILE (itag(i)/=0.AND.i<=nseg)
4644 i=i+1
4645 ENDDO
4646 IF (i==nseg+1) THEN
4647 istop=1
4648 cycle
4649 ENDIF
4650C
4651 iclose=0
4652 iseg=i
4653 DO WHILE (iclose==0)
4654 nnp=nnp+1
4655 poly(nnp,npoly)=tseg(1,iseg)
4656 in=tseg(2,iseg)
4657 IF (tseg(3,iseg)>0) THEN
4658 itag(iseg)=1
4659 IF (in<0) THEN
4660 iseg_old=iseg
4661 iseg=0
4662 DO j=nseg1+1,nseg
4663 IF (itag(j)/=0) cycle
4664 in1=tseg(1,j)
4665 in2=tseg(2,j)
4666 IF (in1==in) THEN
4667 iseg=j
4668 ELSEIF (in2==in) THEN
4669 iseg=j
4670 tseg(1,j)=in2
4671 tseg(2,j)=in1
4672 ENDIF
4673 ENDDO
4674 IF (iseg==0) iseg=tseg(3,iseg_old)
4675 ELSE
4676 iseg=tseg(3,iseg)
4677 ENDIF
4678 ELSE
4679 IF (itag2(iseg)==1) itag(iseg)=1
4680 itag2(iseg)=1
4681 iseg=0
4682 DO j=1,nseg1
4683 in1=tseg(1,j)
4684 IF (in1==in) iseg=j
4685 ENDDO
4686 ENDIF
4687C
4688 IF (itag(iseg)/=0) THEN
4689 iclose=1
4690 lenpoly(npoly)=nnp
4691 npoly=npoly+1
4692 nnp=0
4693 ENDIF
4694 ENDDO
4695 ENDDO
4696 npoly=npoly-1
4697C Attribution des trous
4698 DO i=1,nhol
4699 j=1
4700 DO WHILE (itag(j)/=-i)
4701 j=j+1
4702 ENDDO
4703 n1=tseg(1,j)
4704 xx=rpoly_old(4+3*(n1-1)+1)
4705 yy=rpoly_old(4+3*(n1-1)+2)
4706 zz=rpoly_old(4+3*(n1-1)+3)
4707 DO j=1,npoly
4708 alpha=zero
4709 DO k=1,lenpoly(j)
4710 kk=k+1
4711 IF (k==lenpoly(j)) kk=1
4712 n1=poly(k,i)
4713 n2=poly(kk,i)
4714 IF (n1>0) THEN
4715 xx1=rpoly_old(4+3*(n1-1)+1)
4716 yy1=rpoly_old(4+3*(n1-1)+2)
4717 zz1=rpoly_old(4+3*(n1-1)+3)
4718 ELSE
4719 xx1=xns(1,-n1)
4720 yy1=xns(2,-n1)
4721 zz1=xns(3,-n1)
4722 ENDIF
4723 IF (n2>0) THEN
4724 xx2=rpoly_old(4+3*(n2-1)+1)
4725 yy2=rpoly_old(4+3*(n2-1)+2)
4726 zz2=rpoly_old(4+3*(n2-1)+3)
4727 ELSE
4728 xx2=xns(1,-n2)
4729 yy2=xns(2,-n2)
4730 zz2=xns(3,-n2)
4731 ENDIF
4732 vx1=xx1-xx
4733 vy1=yy1-yy
4734 vz1=zz1-zz
4735 vx2=xx2-xx
4736 vy2=yy2-yy
4737 vz2=zz2-zz
4738 nr1=sqrt(vx1**2+vy1**2+vz1**2)
4739 nr2=sqrt(vx2**2+vy2**2+vz2**2)
4740 vx1=vx1/nr1
4741 vy1=vy1/nr1
4742 vz1=vz1/nr1
4743 vx2=vx2/nr2
4744 vy2=vy2/nr2
4745 vz2=vz2/nr2
4746 ss=vx1*vx2+vy1*vy2+vz1*vz2
4747 vvx=vy1*vz2-vz1*vy2
4748 vvy=vz1*vx2-vx1*vz2
4749 vvz=vx1*vy2-vy1*vx2
4750 ss1=npx*vvx+npy*vvy+npz*vvz
4751 IF (ss1>=zero) THEN
4752 alpha=alpha+acos(ss)
4753 ELSE
4754 alpha=alpha-acos(ss)
4755 ENDIF
4756 ENDDO
4757C
4758 IF (abs(alpha)>=one) thol(i)=j
4759 ENDDO
4760 ENDDO
4761C Creation des tableaux de sortie
4762 DO i=1,ns
4763 rns(2,i)=xns(1,i)
4764 rns(3,i)=xns(2,i)
4765 rns(4,i)=xns(3,i)
4766 ENDDO
4767C
4768 DO i=1,npoly
4769 ipoly(1,i)=ipoly_old(1)
4770 ipoly(3,i)=ipoly_old(3)
4771 ipoly(4,i)=ipoly_old(4)
4772 ipoly(5,i)=ipoly_old(5)
4773 ipoly(6,i)=ipoly_old(6)
4774 rpoly(1,i)=zero
4775 rpoly(2,i)=npx
4776 rpoly(3,i)=npy
4777 rpoly(4,i)=npz
4778 nnp=0
4779 DO j=1,lenpoly(i)
4780 jn=poly(j,i)
4781 IF (jn>0) THEN
4782 nnp=nnp+1
4783 ipoly(6+nnp,i)=ipoly_old(6+jn)
4784 rpoly(4+3*(nnp-1)+1,i)=rpoly_old(4+3*(jn-1)+1)
4785 rpoly(4+3*(nnp-1)+2,i)=rpoly_old(4+3*(jn-1)+2)
4786 rpoly(4+3*(nnp-1)+3,i)=rpoly_old(4+3*(jn-1)+3)
4787 ielnod(nnp,i)=ielnod_old(jn)
4788 ELSE
4789 IF (itags(-jn)==0) THEN
4790 nnp=nnp+1
4791 ipoly(6+nnp,i)=-nns3+jn
4792 rpoly(4+3*(nnp-1)+1,i)=xns(1,-jn)
4793 rpoly(4+3*(nnp-1)+2,i)=xns(2,-jn)
4794 rpoly(4+3*(nnp-1)+3,i)=xns(3,-jn)
4795 ielnod(nnp,i)=-1
4796 ELSE
4797 jj=poly(j+1,i)
4798 IF (j==lenpoly(i)) jj=poly(1,i)
4799 x1=xns(1,-jn)
4800 y1=xns(2,-jn)
4801 z1=xns(3,-jn)
4802 IF (jj>0) THEN
4803 x2=rpoly_old(4+3*(jj-1)+1)
4804 y2=rpoly_old(4+3*(jj-1)+2)
4805 z2=rpoly_old(4+3*(jj-1)+3)
4806 ELSE
4807 x2=xns(1,-jj)
4808 y2=xns(2,-jj)
4809 z2=xns(3,-jj)
4810 ENDIF
4811 ll=sqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)
4812 IF (ll>tole) THEN
4813 nnp=nnp+1
4814 ipoly(6+nnp,i)=-nns3+jn
4815 rpoly(4+3*(nnp-1)+1,i)=xns(1,-jn)
4816 rpoly(4+3*(nnp-1)+2,i)=xns(2,-jn)
4817 rpoly(4+3*(nnp-1)+3,i)=xns(3,-jn)
4818 ielnod(nnp,i)=-1
4819 ENDIF
4820 ENDIF
4821 ENDIF
4822 ENDDO
4823C
4824 nholp=0
4825
4826 DO j=1,nhol
4827 IF (thol(j)==i) THEN
4828 nholp=nholp+1
4829 iadholp(nholp)=nnp+1
4830 DO k=1,nseg
4831 IF (tseg(3,k)==-j) THEN
4832 nnp=nnp+1
4833 kn=tseg(1,k)
4834 ipoly(6+nnp,i)=kn
4835 rpoly(4+3*(nnp-1)+1,i)=rpoly_old(4+3*(kn-1)+1)
4836 rpoly(4+3*(nnp-1)+2,i)=rpoly_old(4+3*(kn-1)+2)
4837 rpoly(4+3*(nnp-1)+3,i)=rpoly_old(4+3*(kn-1)+3)
4838 ENDIF
4839 ENDDO
4840 ENDIF
4841 ENDDO
4842 ipoly(2,i)=nnp
4843 ipoly(6+nnp+1,i)=nholp
4844 DO j=1,nholp
4845 ipoly(6+nnp+1+j,i)=iadholp(j)
4846 ENDDO
4847 ENDDO
4848C
4849 DO i=nseg1+1,nseg
4850 nbnedge=nbnedge+1
4851 n1=-tseg(1,i)
4852 n2=-tseg(2,i)
4853 inedge(1,nbnedge)=ipoly_old(1)
4854 inedge(2,nbnedge)=nns3+n1
4855 inedge(3,nbnedge)=nns3+n2
4856 inedge(4,nbnedge)=ipoly_old(3)
4857 inedge(5,nbnedge)=ipoly_old(4)
4858 inedge(6,nbnedge)=inz
4859C
4860 xx1=xns(1,n1)
4861 yy1=xns(2,n1)
4862 zz1=xns(3,n1)
4863C
4864 xx2=xns(1,n2)
4865 yy2=xns(2,n2)
4866 zz2=xns(3,n2)
4867C
4868 rnedge(1,nbnedge)=xx1
4869 rnedge(2,nbnedge)=yy1
4870 rnedge(3,nbnedge)=zz1
4871 rnedge(4,nbnedge)=xx2
4872 rnedge(5,nbnedge)=yy2
4873 rnedge(6,nbnedge)=zz2
4874 ENDDO
4875C
4876 DO i=1,npoly
4877 zl=zero
4878 zlm=zero
4879 DO j=1,ipoly(2,i)
4880 xx=rpoly(4+3*(j-1)+1,i)
4881 yy=rpoly(4+3*(j-1)+2,i)
4882 zz=rpoly(4+3*(j-1)+3,i)
4883 zl=(xx-x0)*nx+(yy-y0)*ny+(zz-z0)*nz
4884 IF (abs(zl)>abs(zlm)) zlm=zl
4885 ENDDO
4886 iz(1,i)=1
4887 IF (zlm>zero) THEN
4888 iz(2,i)=inz+1
4889 ELSEIF (zlm<zero) THEN
4890 iz(2,i)=inz
4891 ENDIF
4892 ENDDO
4893C
4894 RETURN

◆ horiedge()

subroutine horiedge ( integer, dimension(*) ipoly,
rpoly,
nx,
ny,
nz,
integer nbnedge,
integer, dimension(6,*) inedge,
rnedge,
x0,
y0,
z0,
integer inz,
integer nns3,
integer, dimension(2,*) nref,
aref,
integer nnsp )

Definition at line 4902 of file fvmesh.F.

4906C-----------------------------------------------
4907C I m p l i c i t T y p e s
4908C-----------------------------------------------
4909#include "implicit_f.inc"
4910C-----------------------------------------------
4911C D u m m y A r g u m e n t s
4912C-----------------------------------------------
4913 INTEGER IPOLY(*), NBNEDGE, INEDGE(6,*), INZ, NNS3, NREF(2,*), NNSP
4914 my_real
4915 . rpoly(*), nx, ny, nz, rnedge(6,*), x0, y0, z0, aref(4,*)
4916C-----------------------------------------------
4917C L o c a l V a r i a b l e s
4918C-----------------------------------------------
4919 INTEGER NN, I, II
4920 my_real
4921 . x1, y1, z1, x2, y2, z2, zl1, zl2, tole, dd
4922C
4923* TOLE=DLAMCH('Epsmach')
4924 tole=em10
4925 nn=ipoly(2)
4926 nnsp=0
4927 DO i=1,nn
4928 ii=i+1
4929 IF (i==nn) ii=1
4930 x1=rpoly(4+3*(i-1)+1)
4931 y1=rpoly(4+3*(i-1)+2)
4932 z1=rpoly(4+3*(i-1)+3)
4933 x2=rpoly(4+3*(ii-1)+1)
4934 y2=rpoly(4+3*(ii-1)+2)
4935 z2=rpoly(4+3*(ii-1)+3)
4936 dd=(x1-x2)**2+(y1-y2)**2+(z1-z2)**2
4937 zl1=(x1-x0)*nx+(y1-y0)*ny+(z1-z0)*nz
4938 zl2=(x2-x0)*nx+(y2-y0)*ny+(z2-z0)*nz
4939 IF (zl1==zero.AND.zl2==zero.AND.dd>=tole) THEN
4940 nbnedge=nbnedge+1
4941C
4942 nnsp=nnsp+1
4943 nref(1,nnsp)=ipoly(6+i)
4944 nref(2,nnsp)=ipoly(6+ii)
4945 aref(1,nnsp)=one
4946 aref(2,nnsp)=x1
4947 aref(3,nnsp)=y1
4948 aref(4,nnsp)=z1
4949 nnsp=nnsp+1
4950 nref(1,nnsp)=ipoly(6+i)
4951 nref(2,nnsp)=ipoly(6+ii)
4952 aref(1,nnsp)=zero
4953 aref(2,nnsp)=x2
4954 aref(3,nnsp)=y2
4955 aref(4,nnsp)=z2
4956C
4957 inedge(1,nbnedge)=ipoly(1)
4958 inedge(2,nbnedge)=nns3+nnsp-1
4959 inedge(3,nbnedge)=nns3+nnsp
4960 inedge(4,nbnedge)=ipoly(3)
4961 inedge(5,nbnedge)=ipoly(4)
4962 inedge(6,nbnedge)=inz
4963C
4964 rnedge(1,nbnedge)=rpoly(4+3*(i-1)+1)
4965 rnedge(2,nbnedge)=rpoly(4+3*(i-1)+2)
4966 rnedge(3,nbnedge)=rpoly(4+3*(i-1)+3)
4967 rnedge(4,nbnedge)=rpoly(4+3*(ii-1)+1)
4968 rnedge(5,nbnedge)=rpoly(4+3*(ii-1)+2)
4969 rnedge(6,nbnedge)=rpoly(4+3*(ii-1)+3)
4970 ENDIF
4971 ENDDO
4972C
4973 RETURN

◆ horipoly()

subroutine horipoly ( integer, dimension(6,*) inedge,
rnedge,
integer, dimension(*) ledge,
integer nedge,
integer, dimension(6+2*nedge+1+nedge,*) ipoly,
rpoly,
integer, dimension(3,*) iz,
integer, dimension(nedge,*) ielnod,
integer npoly,
nx,
ny,
nz,
integer inz,
integer ibric )

Definition at line 4981 of file fvmesh.F.

4984C-----------------------------------------------
4985C I m p l i c i t T y p e s
4986C-----------------------------------------------
4987#include "implicit_f.inc"
4988C-----------------------------------------------
4989C D u m m y A r g u m e n t s
4990C-----------------------------------------------
4991 INTEGER INEDGE(6,*), LEDGE(*), NEDGE, IPOLY(6+2*NEDGE+1+NEDGE,*),
4992 . IZ(3,*), IELNOD(NEDGE,*), NPOLY, INZ, IBRIC
4993 my_real
4994 . rnedge(6,*), rpoly(4+6*nedge+3*nedge,*), nx, ny, nz
4995C-----------------------------------------------
4996C L o c a l V a r i a b l e s
4997C-----------------------------------------------
4998 INTEGER NN, I, II, TNOD(2*NEDGE), TSEG(2,NEDGE), NN_OLD, JFOUND,
4999 . J, JJ, REDIR1(2*NEDGE), REDIR2(2*NEDGE), ITAG(2*NEDGE),
5000 . ITAGSEG(NEDGE+1), ISTOP, ICLOSE, N1, N2, IN, NNP,
5001 . POLY(2*NEDGE,NEDGE), ISEG, LENPOLY(NEDGE), IFOUND,
5002 . IADHOL(NEDGE), NHOL, K, KK, IPOUT, M, MM, NSEC, KSMIN,
5003 . NEDGE_OLD, TSEG_OLD(2,NEDGE)
5004 my_real
5005 . tole, xnod(3,2*nedge), xx1, yy1, zz1, xx2, yy2, zz2,
5006 . dd, xloc(2,nedge), xsec(nedge), phol(3,nedge), alpha,
5007 . x1, y1, z1, vx1, vy1, vz1, vx2, vy2, vz2, nr1, nr2,
5008 . ss, vvx, vvy, vvz, ss1, x2, y2, z2, ylmin, ylmax, xsmin1,
5009 . xsmin2, xx, yy, zz, ylsec, xs, ys, lmax, ll
5010C
5011 my_real
5012 . dlamch
5013 EXTERNAL dlamch
5014C
5015* TOLE=DLAMCH('Epsmach')
5016 tole=em10
5017C Creation de la liste des noeuds
5018 nn=0
5019 lmax=zero
5020 DO i=1,nedge
5021 ii=ledge(i)
5022 nn=nn+1
5023 tnod(nn)=inedge(2,ii)
5024 xnod(1,nn)=rnedge(1,ii)
5025 xnod(2,nn)=rnedge(2,ii)
5026 xnod(3,nn)=rnedge(3,ii)
5027 tseg(1,i)=nn
5028 nn=nn+1
5029 tnod(nn)=inedge(3,ii)
5030 xnod(1,nn)=rnedge(4,ii)
5031 xnod(2,nn)=rnedge(5,ii)
5032 xnod(3,nn)=rnedge(6,ii)
5033 tseg(2,i)=nn
5034 ll=(rnedge(1,ii)-rnedge(4,ii))**2+
5035 . (rnedge(2,ii)-rnedge(5,ii))**2+
5036 . (rnedge(3,ii)-rnedge(6,ii))**2
5037 lmax=max(lmax,ll)
5038 ENDDO
5039* TOLE=TOLE*LMAX
5040C Elimination des noeuds doubles
5041 nn_old=nn
5042 nn=0
5043 DO i=1,nn_old
5044 xx1=xnod(1,i)
5045 yy1=xnod(2,i)
5046 zz1=xnod(3,i)
5047 jfound=0
5048 DO j=1,nn
5049 jj=redir1(j)
5050 xx2=xnod(1,jj)
5051 yy2=xnod(2,jj)
5052 zz2=xnod(3,jj)
5053 dd=sqrt((xx1-xx2)**2+(yy1-yy2)**2+(zz1-zz2)**2)
5054 IF (dd<=tole) jfound=j
5055 ENDDO
5056 IF (jfound==0) THEN
5057 nn=nn+1
5058 redir1(nn)=i
5059 redir2(i)=nn
5060 ELSE
5061 redir2(i)=jfound
5062 ENDIF
5063 ENDDO
5064C
5065 DO i=1,nedge
5066 n1=tseg(1,i)
5067 n2=tseg(2,i)
5068 tseg(1,i)=redir2(n1)
5069 tseg(2,i)=redir2(n2)
5070 ENDDO
5071C Elimination des segments de longueur nulle
5072 nedge_old=nedge
5073 nedge=0
5074 DO i=1,nedge_old
5075 tseg_old(1,i)=tseg(1,i)
5076 tseg_old(2,i)=tseg(2,i)
5077 ENDDO
5078 DO i=1,nedge_old
5079 IF (tseg_old(1,i)/=tseg_old(2,i)) THEN
5080 nedge=nedge+1
5081 tseg(1,nedge)=tseg_old(1,i)
5082 tseg(2,nedge)=tseg_old(2,i)
5083 ENDIF
5084 ENDDO
5085C Construction des polygones
5086 DO i=1,nn
5087 itag(i)=0
5088 ENDDO
5089 DO i=1,nedge+1
5090 itagseg(i)=0
5091 ENDDO
5092 npoly=1
5093 istop=0
5094 DO WHILE (istop==0)
5095 i=1
5096 DO WHILE (itagseg(i)==1.AND.i<=nedge)
5097 i=i+1
5098 ENDDO
5099 IF (i==nedge+1) THEN
5100 istop=1
5101 cycle
5102 ENDIF
5103C
5104 iclose=0
5105 iseg=i
5106 itagseg(iseg)=1
5107 n1=tseg(1,iseg)
5108 n2=tseg(2,iseg)
5109 itag(n1)=1
5110 itag(n2)=1
5111 in=n2
5112 nnp=1
5113 poly(1,npoly)=redir1(n1)
5114 DO WHILE (iclose==0)
5115 ifound=0
5116 i=0
5117 DO WHILE (ifound==0)
5118 i=i+1
5119 IF (itagseg(i)==1) cycle
5120 n1=tseg(1,i)
5121 n2=tseg(2,i)
5122 IF (n1==in) THEN
5123 ifound=1
5124 IF (itag(n2)==1) iclose=1
5125 iseg=i
5126 in=n2
5127 nnp=nnp+1
5128 poly(nnp,npoly)=redir1(n1)
5129 itag(n2)=1
5130 ELSEIF (n2==in) THEN
5131 ifound=1
5132 IF (itag(n1)==1) iclose=1
5133 iseg=i
5134 in=n1
5135 nnp=nnp+1
5136 poly(nnp,npoly)=redir1(n2)
5137 itag(n1)=1
5138 ENDIF
5139 ENDDO
5140 itagseg(iseg)=1
5141 ENDDO
5142C
5143 IF (iclose==1) THEN
5144 lenpoly(npoly)=nnp
5145 npoly=npoly+1
5146 ENDIF
5147 ENDDO
5148 npoly=npoly-1
5149C Creation des tableaux de sortie
5150 DO i=1,npoly
5151 ipoly(1,i)=2
5152 ipoly(2,i)=lenpoly(i)
5153 ipoly(3,i)=ibric
5154 ipoly(4,i)=ibric
5155 ipoly(5,i)=0
5156 ipoly(6,i)=0
5157 rpoly(1,i)=zero
5158 rpoly(2,i)=nx
5159 rpoly(3,i)=ny
5160 rpoly(4,i)=nz
5161 DO j=1,lenpoly(i)
5162 jj=poly(j,i)
5163 ipoly(6+j,i)=-tnod(jj)
5164 rpoly(4+3*(j-1)+1,i)=xnod(1,jj)
5165 rpoly(4+3*(j-1)+2,i)=xnod(2,jj)
5166 rpoly(4+3*(j-1)+3,i)=xnod(3,jj)
5167 ielnod(j,i)=-1
5168 ENDDO
5169 ipoly(6+lenpoly(i)+1,i)=0
5170C
5171 iz(1,i)=2
5172 iz(2,i)=inz
5173 iz(3,i)=inz+1
5174 ENDDO
5175C
5176C Recherche des trous
5177 DO i=1,npoly
5178C
5179 nhol=0
5180 DO j=1,npoly
5181 iadhol(j)=0
5182 ENDDO
5183C
5184 DO j=1,npoly
5185 IF (i==j) cycle
5186 alpha=zero
5187 x1=rpoly(5,j)
5188 y1=rpoly(6,j)
5189 z1=rpoly(7,j)
5190 DO k=1,lenpoly(i)
5191 kk=k+1
5192 IF (k==lenpoly(i)) kk=1
5193 xx1=rpoly(4+3*(k-1)+1,i)
5194 yy1=rpoly(4+3*(k-1)+2,i)
5195 zz1=rpoly(4+3*(k-1)+3,i)
5196 xx2=rpoly(4+3*(kk-1)+1,i)
5197 yy2=rpoly(4+3*(kk-1)+2,i)
5198 zz2=rpoly(4+3*(kk-1)+3,i)
5199 vx1=xx1-x1
5200 vy1=yy1-y1
5201 vz1=zz1-z1
5202 vx2=xx2-x1
5203 vy2=yy2-y1
5204 vz2=zz2-z1
5205 nr1=sqrt(vx1**2+vy1**2+vz1**2)
5206 nr2=sqrt(vx2**2+vy2**2+vz2**2)
5207 vx1=vx1/nr1
5208 vy1=vy1/nr1
5209 vz1=vz1/nr1
5210 vx2=vx2/nr2
5211 vy2=vy2/nr2
5212 vz2=vz2/nr2
5213 ss=vx1*vx2+vy1*vy2+vz1*vz2
5214 vvx=vy1*vz2-vz1*vy2
5215 vvy=vz1*vx2-vx1*vz2
5216 vvz=vx1*vy2-vy1*vx2
5217 ss1=nx*vvx+ny*vvy+nz*vvz
5218 IF (ss1>=zero) THEN
5219 alpha=alpha+acos(ss)
5220 ELSE
5221 alpha=alpha-acos(ss)
5222 ENDIF
5223 ENDDO
5224C
5225 IF (abs(alpha)>=one) THEN
5226C Le premier point du polygone i est a l'interieur du polygone j
5227C On teste tous les autres
5228 ipout=0
5229 DO k=2,lenpoly(j)
5230 x2=rpoly(4+3*(k-1)+1,j)
5231 y2=rpoly(4+3*(k-1)+2,j)
5232 z2=rpoly(4+3*(k-1)+3,j)
5233 alpha=zero
5234 DO m=1,lenpoly(i)
5235 mm=m+1
5236 IF (m==lenpoly(i)) mm=1
5237 xx1=rpoly(4+3*(m-1)+1,i)
5238 yy1=rpoly(4+3*(m-1)+2,i)
5239 zz1=rpoly(4+3*(m-1)+3,i)
5240 xx2=rpoly(4+3*(mm-1)+1,i)
5241 yy2=rpoly(4+3*(mm-1)+2,i)
5242 zz2=rpoly(4+3*(mm-1)+3,i)
5243 vx1=xx1-x1
5244 vy1=yy1-y1
5245 vz1=zz1-z1
5246 vx2=xx2-x1
5247 vy2=yy2-y1
5248 vz2=zz2-z1
5249 nr1=sqrt(vx1**2+vy1**2+vz1**2)
5250 nr2=sqrt(vx2**2+vy2**2+vz2**2)
5251 vx1=vx1/nr1
5252 vy1=vy1/nr1
5253 vz1=vz1/nr1
5254 vx2=vx2/nr2
5255 vy2=vy2/nr2
5256 vz2=vz2/nr2
5257 ss=vx1*vx2+vy1*vy2+vz1*vz2
5258 vvx=vy1*vz2-vz1*vy2
5259 vvy=vz1*vx2-vx1*vz2
5260 vvz=vx1*vy2-vy1*vx2
5261 ss1=nx*vvx+ny*vvy+nz*vvz
5262 IF (ss1>=zero) THEN
5263 alpha=alpha+acos(ss)
5264 ELSE
5265 alpha=alpha-acos(ss)
5266 ENDIF
5267 ENDDO
5268 IF (abs(alpha)<one) ipout=1
5269 ENDDO
5270C
5271 IF (ipout==1) cycle
5272C Le polygone j est un trou dans le polygone i
5273 ipoly(1,j)=-1
5274 nhol=nhol+1
5275 iadhol(nhol)=lenpoly(i)
5276 DO k=1,lenpoly(j)
5277 ipoly(6+iadhol(nhol)+k,i)=ipoly(6+k,j)
5278 ielnod(iadhol(nhol)+k,i)=ielnod(k,j)
5279 rpoly(4+3*iadhol(nhol)+3*(k-1)+1,i)=
5280 . rpoly(4+3*(k-1)+1,j)
5281 rpoly(4+3*iadhol(nhol)+3*(k-1)+2,i)=
5282 . rpoly(4+3*(k-1)+2,j)
5283 rpoly(4+3*iadhol(nhol)+3*(k-1)+3,i)=
5284 . rpoly(4+3*(k-1)+3,j)
5285 ENDDO
5286 lenpoly(i)=lenpoly(i)+lenpoly(j)
5287C Point interieur polygone j
5288 vx1=rpoly(5,j)-rpoly(8,j)
5289 vy1=rpoly(6,j)-rpoly(9,j)
5290 vz1=rpoly(7,j)-rpoly(10,j)
5291 ss=sqrt(vx1**2+vy1**2+vz1**2)
5292 vx1=vx1/ss
5293 vy1=vy1/ss
5294 vz1=vz1/ss
5295 vx2=ny*vz1-nz*vy1
5296 vy2=nz*vx1-nx*vz1
5297 vz2=nx*vy1-ny*vx1
5298 x1=rpoly(5,j)
5299 y1=rpoly(6,j)
5300 z1=rpoly(7,j)
5301 xloc(1,1)=zero
5302 xloc(2,1)=zero
5303 ylmin=ep30
5304 ylmax=-ep30
5305 DO k=2,lenpoly(j)
5306 xx=rpoly(4+3*(k-1)+1,j)
5307 yy=rpoly(4+3*(k-1)+2,j)
5308 zz=rpoly(4+3*(k-1)+3,j)
5309 vvx=xx-x1
5310 vvy=yy-y1
5311 vvz=zz-z1
5312 xloc(1,k)=vvx*vx1+vvy*vy1+vvz*vz1
5313 xloc(2,k)=vvx*vx2+vvy*vy2+vvz*vz2
5314 IF (xloc(2,k)<ylmin) ylmin=xloc(2,k)
5315 IF (xloc(2,k)>ylmax) ylmax=xloc(2,k)
5316 ENDDO
5317 ylsec=half*(ylmin+ylmax)
5318C
5319 nsec=0
5320 DO k=1,lenpoly(j)
5321 kk=k+1
5322 IF (k==lenpoly(j)) kk=1
5323 x1=xloc(1,k)
5324 y1=xloc(2,k)
5325 x2=xloc(1,kk)
5326 y2=xloc(2,kk)
5327 IF (y1-y2/=zero) THEN
5328 alpha=(ylsec-y2)/(y1-y2)
5329 IF (alpha>=zero.AND.alpha<=one) THEN
5330 nsec=nsec+1
5331 xsec(nsec)=alpha*x1+(one-alpha)*x2
5332 ENDIF
5333 ELSE
5334 IF (y1==ylsec) THEN
5335 nsec=nsec+1
5336 xsec(nsec)=x1
5337 nsec=nsec+1
5338 xsec(nsec)=x2
5339 ENDIF
5340 ENDIF
5341 ENDDO
5342C
5343 xsmin1=ep30
5344 DO k=1,nsec
5345 IF (xsec(k)<xsmin1) THEN
5346 xsmin1=xsec(k)
5347 ksmin=k
5348 ENDIF
5349 ENDDO
5350 xsmin2=ep30
5351 DO k=1,nsec
5352 IF (k==ksmin) cycle
5353 IF (xsec(k)<xsmin2) xsmin2=xsec(k)
5354 ENDDO
5355C
5356 xs=half*(xsmin1+xsmin2)
5357 ys=ylsec
5358 phol(1,nhol)=rpoly(5,j)+xs*vx1+ys*vx2
5359 phol(2,nhol)=rpoly(6,j)+xs*vy1+ys*vy2
5360 phol(3,nhol)=rpoly(7,j)+xs*vz1+ys*vz2
5361 ENDIF
5362 ENDDO
5363C
5364 ipoly(2,i)=lenpoly(i)
5365 ipoly(6+lenpoly(i)+1,i)=nhol
5366 DO j=1,nhol
5367 ipoly(6+lenpoly(i)+1+j,i)=iadhol(j)
5368 rpoly(4+3*lenpoly(i)+3*(j-1)+1,i)=phol(1,j)
5369 rpoly(4+3*lenpoly(i)+3*(j-1)+2,i)=phol(2,j)
5370 rpoly(4+3*lenpoly(i)+3*(j-1)+3,i)=phol(3,j)
5371 ENDDO
5372 ENDDO
5373C
5374 RETURN
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69

◆ itribox()

subroutine itribox ( tri,
box,
norm,
integer nverts,
poly,
integer nvmax )

Definition at line 3105 of file fvmesh.F.

3107C-----------------------------------------------
3108C I m p l i c i t T y p e s
3109C-----------------------------------------------
3110#include "implicit_f.inc"
3111C-----------------------------------------------
3112C D u m m y A r g u m e n t s
3113C-----------------------------------------------
3114 INTEGER NVERTS, NVMAX
3115 my_real
3116 . tri(3,*), box(3,*), norm(3,*), poly(3,*)
3117C-----------------------------------------------
3118C L o c a l V a r i a b l e s
3119C-----------------------------------------------
3120 INTEGER I, NPOUT, J
3121 my_real
3122 . a(3), n(3), polyout(3,nvmax)
3123 INTEGER P_REF(6)
3124 DATA p_ref /1,5,1,2,3,4/
3125C
3126C Intersection triangle-box
3127 nverts=3
3128 DO i=1,nverts
3129 poly(1,i)=tri(1,i)
3130 poly(2,i)=tri(2,i)
3131 poly(3,i)=tri(3,i)
3132 ENDDO
3133C
3134 DO i=1,6
3135 a(1)=box(1,p_ref(i))
3136 a(2)=box(2,p_ref(i))
3137 a(3)=box(3,p_ref(i))
3138 n(1)=norm(1,i)
3139 n(2)=norm(2,i)
3140 n(3)=norm(3,i)
3141 CALL polclip(poly, nverts, a, n, polyout,
3142 . npout)
3143 nverts=npout
3144 DO j=1,nverts
3145 poly(1,j)=polyout(1,j)
3146 poly(2,j)=polyout(2,j)
3147 poly(3,j)=polyout(3,j)
3148 ENDDO
3149 ENDDO
3150C
3151 RETURN
polclip
subroutine polclip(polyin, npin, a, n, polyout, npout)
Definition fvmesh.F:3161

◆ polclip()

subroutine polclip ( polyin,
integer npin,
a,
n,
polyout,
integer npout )

Definition at line 3159 of file fvmesh.F.

3161C-----------------------------------------------
3162C I m p l i c i t T y p e s
3163C-----------------------------------------------
3164#include "implicit_f.inc"
3165C-----------------------------------------------
3166C D u m m y A r g u m e n t s
3167C-----------------------------------------------
3168 INTEGER NPIN, NPOUT
3169 my_real
3170 . polyin(3,*), a(*), n(*), polyout(3,*)
3171C-----------------------------------------------
3172C L o c a l V a r i a b l e s
3173C-----------------------------------------------
3174 INTEGER I, II, IC1, IC2
3175 my_real
3176 . x1, y1, z1, x2, y2, z2, ss1, ss2, x12, y12, z12,
3177 . xa1, ya1, za1, alpha, xm, ym, zm, ssn
3178C
3179 my_real
3180 . dlamch
3181 EXTERNAL dlamch
3182C
3183C Clipping d'un polygone par un plan
3184 npout=0
3185 DO i=1,npin
3186 IF (i/=npin) THEN
3187 ii=i+1
3188 ELSE
3189 ii=1
3190 ENDIF
3191C
3192 x1=polyin(1,i)
3193 y1=polyin(2,i)
3194 z1=polyin(3,i)
3195 x2=polyin(1,ii)
3196 y2=polyin(2,ii)
3197 z2=polyin(3,ii)
3198C
3199 ss1=(x1-a(1))*n(1)+(y1-a(2))*n(2)+(z1-a(3))*n(3)
3200 ss2=(x2-a(1))*n(1)+(y2-a(2))*n(2)+(z2-a(3))*n(3)
3201 IF (ss1<zero.AND.ss2<zero) cycle
3202 IF (ss1>=zero.AND.ss2>=zero) THEN
3203 npout=npout+1
3204 polyout(1,npout)=x1
3205 polyout(2,npout)=y1
3206 polyout(3,npout)=z1
3207 cycle
3208 ENDIF
3209C
3210 x12=x2-x1
3211 y12=y2-y1
3212 z12=z2-z1
3213 xa1=x1-a(1)
3214 ya1=y1-a(2)
3215 za1=z1-a(3)
3216 ssn=x12*n(1)+y12*n(2)+z12*n(3)
3217 alpha=-(xa1*n(1)+ya1*n(2)+za1*n(3))/ssn
3218 xm=x1+alpha*x12
3219 ym=y1+alpha*y12
3220 zm=z1+alpha*z12
3221 IF (ss1>=zero) THEN
3222 npout=npout+1
3223 polyout(1,npout)=x1
3224 polyout(2,npout)=y1
3225 polyout(3,npout)=z1
3226 npout=npout+1
3227 polyout(1,npout)=xm
3228 polyout(2,npout)=ym
3229 polyout(3,npout)=zm
3230 ELSEIF (ss2>=zero) THEN
3231 npout=npout+1
3232 polyout(1,npout)=xm
3233 polyout(2,npout)=ym
3234 polyout(3,npout)=zm
3235 ENDIF
3236 ENDDO
3237C
3238 RETURN

◆ polyhedr()

subroutine polyhedr ( integer, dimension(6+nppmax,*) ipoly,
rpoly,
integer, dimension(*) polb,
integer npolb,
integer, dimension(nphmax+2,*) polh,
integer npolh,
integer nrpmax,
integer nphmax,
integer ibric,
lmin,
integer info,
integer npolhmax,
integer nppmax,
tole2 )

Definition at line 4116 of file fvmesh.F.

4119C-----------------------------------------------
4120C I m p l i c i t T y p e s
4121C-----------------------------------------------
4122#include "implicit_f.inc"
4123C-----------------------------------------------
4124C D u m m y A r g u m e n t s
4125C-----------------------------------------------
4126 INTEGER NPPMAX, IPOLY(6+NPPMAX,*), POLB(*), NPOLB, NPHMAX,
4127 . POLH(NPHMAX+2,*),NPOLH, NRPMAX, IBRIC, INFO, NPOLHMAX
4128 my_real
4129 . rpoly(nrpmax,*), lmin, tole2
4130C-----------------------------------------------
4131C L o c a l V a r i a b l e s
4132C-----------------------------------------------
4133 INTEGER I, ITAG(NPOLB), ICMAX, ICUR, II, J, JJ, K, KK, ISTOP,
4134 . L, LL, ICUR_OLD, ITY, JMIN, PMIN, POLD
4135 my_real
4136 . x1, y1, z1, x2, y2, z2, xx1, yy1, zz1, xx2, yy2, zz2,
4137 . dd11, dd12, dd21, dd22, tole
4138C
4139 tole=tole2*lmin
4140C
4141 DO i=1,npolb
4142 itag(i)=0
4143 ENDDO
4144C
4145 icmax=0
4146 DO i=1,npolb
4147 IF (itag(i)==0) THEN
4148 icmax=icmax+1
4149 itag(i)=icmax
4150 icur=icmax
4151 ELSE
4152 icur=itag(i)
4153 ENDIF
4154 ii=polb(i)
4155 DO j=1,ipoly(2,ii)
4156 IF (j/=ipoly(2,ii)) THEN
4157 jj=j+1
4158 ELSE
4159 jj=1
4160 ENDIF
4161 x1=rpoly(4+3*(j-1)+1,ii)
4162 y1=rpoly(4+3*(j-1)+2,ii)
4163 z1=rpoly(4+3*(j-1)+3,ii)
4164 x2=rpoly(4+3*(jj-1)+1,ii)
4165 y2=rpoly(4+3*(jj-1)+2,ii)
4166 z2=rpoly(4+3*(jj-1)+3,ii)
4167 DO k=1,npolb
4168 IF (k==i) cycle
4169 kk=polb(k)
4170 istop=0
4171 l=0
4172 DO WHILE (istop==0.AND.l<ipoly(2,kk))
4173 l=l+1
4174 IF (l/=ipoly(2,kk)) THEN
4175 ll=l+1
4176 ELSE
4177 ll=1
4178 ENDIF
4179 xx1=rpoly(4+3*(l-1)+1,kk)
4180 yy1=rpoly(4+3*(l-1)+2,kk)
4181 zz1=rpoly(4+3*(l-1)+3,kk)
4182 xx2=rpoly(4+3*(ll-1)+1,kk)
4183 yy2=rpoly(4+3*(ll-1)+2,kk)
4184 zz2=rpoly(4+3*(ll-1)+3,kk)
4185 dd11=(xx1-x1)**2+(yy1-y1)**2+(zz1-z1)**2
4186 dd21=(xx2-x1)**2+(yy2-y1)**2+(zz2-z1)**2
4187 dd12=(xx1-x2)**2+(yy1-y2)**2+(zz1-z2)**2
4188 dd22=(xx2-x2)**2+(yy2-y2)**2+(zz2-z2)**2
4189 IF ((dd11<=tole.AND.dd22<=tole).OR.
4190 . (dd21<=tole.AND.dd12<=tole)) istop=l
4191 ENDDO
4192 IF (istop/=0) THEN
4193 IF (itag(k)==0) THEN
4194 itag(k)=icur
4195 ELSE
4196 icur_old=icur
4197 icur=itag(k)
4198 DO l=1,npolb
4199 IF (itag(l)==icur_old) itag(l)=icur
4200 ENDDO
4201 ENDIF
4202 ENDIF
4203 ENDDO
4204 ENDDO
4205 ENDDO
4206C
4207 npolh=0
4208 DO i=1,icmax
4209 ii=0
4210 DO j=1,npolb
4211 IF (itag(j)==i) ii=ii+1
4212 ENDDO
4213 IF (ii/=0) npolh=npolh+1
4214 ENDDO
4215 IF (npolh>npolhmax) THEN
4216 info=1
4217 npolhmax=npolh
4218 RETURN
4219 ENDIF
4220C
4221 npolh=0
4222 DO i=1,icmax
4223 ii=0
4224 DO j=1,npolb
4225 IF (itag(j)==i) ii=ii+1
4226 ENDDO
4227 IF (ii/=0) THEN
4228 npolh=npolh+1
4229 polh(1,npolh)=ii
4230 polh(2,npolh)=ibric
4231 ii=0
4232 DO j=1,npolb
4233 IF (itag(j)==i) THEN
4234 ii=ii+1
4235 jj=polb(j)
4236 polh(2+ii,npolh)=jj
4237 ity=ipoly(1,jj)
4238 IF (ity==1) THEN
4239 ipoly(5,jj)=npolh
4240 cycle
4241 ENDIF
4242 IF (ipoly(5,jj)==0) THEN
4243 ipoly(5,jj)=npolh
4244 ELSE
4245 ipoly(6,jj)=npolh
4246 ENDIF
4247 ENDIF
4248 ENDDO
4249C Tri par ordre croissant de polh
4250 DO j=1,polh(1,npolh)
4251 jmin=j
4252 pmin=polh(2+j,npolh)
4253 DO k=j+1,polh(1,npolh)
4254 IF (polh(2+k,npolh)<pmin) THEN
4255 jmin=k
4256 pmin=polh(2+k,npolh)
4257 ENDIF
4258 ENDDO
4259 pold=polh(2+j,npolh)
4260 polh(2+j,npolh)=pmin
4261 polh(2+jmin,npolh)=pold
4262 ENDDO
4263 ENDIF
4264 ENDDO
4265C
4266 info=0
4267 RETURN