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 4277 of file fvmesh.F.

4280C-----------------------------------------------
4281C I m p l i c i t T y p e s
4282C-----------------------------------------------
4283#include "implicit_f.inc"
4284C-----------------------------------------------
4285C D u m m y A r g u m e n t s
4286C-----------------------------------------------
4287 my_real
4288 . vpx, vpy, vpz, v1x, v1y,
4289 . v1z, v2x, v2y, v2z, ksi,
4290 . eta
4291C-----------------------------------------------
4292C L o c a l V a r i a b l e s
4293C-----------------------------------------------
4294 my_real
4295 . sp1, sp2, s11, s22, s12, d
4296C
4297 sp1=vpx*v1x+vpy*v1y+vpz*v1z
4298 sp2=vpx*v2x+vpy*v2y+vpz*v2z
4299 s11=v1x*v1x+v1y*v1y+v1z*v1z
4300 s22=v2x*v2x+v2y*v2y+v2z*v2z
4301 s12=v1x*v2x+v1y*v2y+v1z*v2z
4302C
4303 d=s11*s22-s12*s12
4304C
4305 ksi=(sp1*s22-sp2*s12)/d
4306 eta=(sp2*s11-sp1*s12)/d
4307C
4308 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 3253 of file fvmesh.F.

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

◆ 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 45 of file fvmesh.F.

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

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

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

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

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

◆ polclip()

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

Definition at line 3161 of file fvmesh.F.

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

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