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

Go to the source code of this file.

Functions/Subroutines

subroutine fvmesh1 (ibuf, elem, x, ivolu, bric, xb, rvolu, nel, neli, nbric, tbric, sfac, dxm, nba, nela, tba, tfaca, tagels, ibufa, elema, tagela, ixs, id, titr, nb_node, ityp)
subroutine coorloc (vpx, vpy, vpz, v1x, v1y, v1z, v2x, v2y, v2z, ksi, eta)
subroutine fvdim (t_monvol)
subroutine fvdeal ()

Function/Subroutine Documentation

◆ coorloc()

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

Definition at line 3416 of file fvmesh.F.

3419C-----------------------------------------------
3420C I m p l i c i t T y p e s
3421C-----------------------------------------------
3422#include "implicit_f.inc"
3423C-----------------------------------------------
3424C D u m m y A r g u m e n t s
3425C-----------------------------------------------
3426 my_real
3427 . vpx, vpy, vpz, v1x, v1y,
3428 . v1z, v2x, v2y, v2z, ksi,
3429 . eta
3430C-----------------------------------------------
3431C L o c a l V a r i a b l e s
3432C-----------------------------------------------
3433 my_real
3434 . sp1, sp2, s11, s22, s12, d
3435C
3436 sp1=vpx*v1x+vpy*v1y+vpz*v1z
3437 sp2=vpx*v2x+vpy*v2y+vpz*v2z
3438 s11=v1x*v1x+v1y*v1y+v1z*v1z
3439 s22=v2x*v2x+v2y*v2y+v2z*v2z
3440 s12=v1x*v2x+v1y*v2y+v1z*v2z
3441C
3442 d=s11*s22-s12*s12
3443C
3444 ksi=(sp1*s22-sp2*s12)/d
3445 eta=(sp2*s11-sp1*s12)/d
3446C
3447 RETURN
my_real
#define my_real
Definition cppsort.cpp:32

◆ fvdeal()

subroutine fvdeal

Definition at line 3500 of file fvmesh.F.

3501C-----------------------------------------------
3502C M o d u l e s
3503C-----------------------------------------------
3504 USE fvbag_mod
3505C-----------------------------------------------
3506C I m p l i c i t T y p e s
3507C-----------------------------------------------
3508#include "implicit_f.inc"
3509C-----------------------------------------------
3510C C o m m o n B l o c k s
3511C-----------------------------------------------
3512C-----------------------------------------------
3513C L o c a l V a r i a b l e s
3514C-----------------------------------------------
3515 INTEGER I
3516C
3517 IF (ALLOCATED(fvdata)) THEN
3518 DO i=1,nfvbag
3519 IF (ASSOCIATED(fvdata(i)%IFVNOD)) DEALLOCATE(fvdata(i)%IFVNOD)
3520 IF (ASSOCIATED(fvdata(i)%RFVNOD)) DEALLOCATE(fvdata(i)%RFVNOD)
3521 IF (ASSOCIATED(fvdata(i)%IFVTRI)) DEALLOCATE(fvdata(i)%IFVTRI)
3522 IF (ASSOCIATED(fvdata(i)%IFVPOLY)) DEALLOCATE(fvdata(i)%IFVPOLY)
3523 IF (ASSOCIATED(fvdata(i)%IFVTADR)) DEALLOCATE(fvdata(i)%IFVTADR)
3524 IF (ASSOCIATED(fvdata(i)%IFVPOLH)) DEALLOCATE(fvdata(i)%IFVPOLH)
3525 IF (ASSOCIATED(fvdata(i)%IFVPADR)) DEALLOCATE(fvdata(i)%IFVPADR)
3526 IF (ASSOCIATED(fvdata(i)%IDPOLH)) DEALLOCATE(fvdata(i)%IDPOLH)
3527 IF (ASSOCIATED(fvdata(i)%IBPOLH)) DEALLOCATE(fvdata(i)%IBPOLH)
3528 IF (ASSOCIATED(fvdata(i)%IFVPOLY_ANIM)) DEALLOCATE(fvdata(i)%IFVPOLY_ANIM)
3529 IF (ASSOCIATED(fvdata(i)%IFVTADR_ANIM)) DEALLOCATE(fvdata(i)%IFVTADR_ANIM)
3530 IF (ASSOCIATED(fvdata(i)%IFVPOLH_ANIM)) DEALLOCATE(fvdata(i)%IFVPOLH_ANIM)
3531 IF (ASSOCIATED(fvdata(i)%IFVPADR_ANIM)) DEALLOCATE(fvdata(i)%IFVPADR_ANIM)
3532 IF (ASSOCIATED(fvdata(i)%REDIR_ANIM)) DEALLOCATE(fvdata(i)%REDIR_ANIM)
3533 IF (ASSOCIATED(fvdata(i)%NOD_ANIM)) DEALLOCATE(fvdata(i)%NOD_ANIM)
3534 IF (ASSOCIATED(fvdata(i)%MPOLH)) DEALLOCATE(fvdata(i)%MPOLH)
3535 IF (ASSOCIATED(fvdata(i)%QPOLH)) DEALLOCATE(fvdata(i)%QPOLH)
3536 IF (ASSOCIATED(fvdata(i)%EPOLH)) DEALLOCATE(fvdata(i)%EPOLH)
3537 IF (ASSOCIATED(fvdata(i)%PPOLH)) DEALLOCATE(fvdata(i)%PPOLH)
3538 IF (ASSOCIATED(fvdata(i)%RPOLH)) DEALLOCATE(fvdata(i)%RPOLH)
3539 IF (ASSOCIATED(fvdata(i)%GPOLH)) DEALLOCATE(fvdata(i)%GPOLH)
3540 IF (ASSOCIATED(fvdata(i)%CPAPOLH)) DEALLOCATE(fvdata(i)%CPAPOLH)
3541 IF (ASSOCIATED(fvdata(i)%CPBPOLH)) DEALLOCATE(fvdata(i)%CPBPOLH)
3542 IF (ASSOCIATED(fvdata(i)%CPCPOLH)) DEALLOCATE(fvdata(i)%CPCPOLH)
3543 IF (ASSOCIATED(fvdata(i)%RMWPOLH)) DEALLOCATE(fvdata(i)%RMWPOLH)
3544 IF (ASSOCIATED(fvdata(i)%VPOLH_INI)) DEALLOCATE(fvdata(i)%VPOLH_INI)
3545 IF (ASSOCIATED(fvdata(i)%DTPOLH)) DEALLOCATE(fvdata(i)%DTPOLH)
3546C
3547 IF (ASSOCIATED(fvdata(i)%BRIC)) DEALLOCATE(fvdata(i)%BRIC)
3548 IF (ASSOCIATED(fvdata(i)%TBRIC)) DEALLOCATE(fvdata(i)%TBRIC)
3549 IF (ASSOCIATED(fvdata(i)%XB)) DEALLOCATE(fvdata(i)%XB)
3550 IF (ASSOCIATED(fvdata(i)%SFAC)) DEALLOCATE(fvdata(i)%SFAC)
3551C
3552 IF (ASSOCIATED(fvdata(i)%TPOLH)) DEALLOCATE(fvdata(i)%TPOLH)
3553 IF (ASSOCIATED(fvdata(i)%CPDPOLH)) DEALLOCATE(fvdata(i)%CPDPOLH)
3554 IF (ASSOCIATED(fvdata(i)%CPEPOLH)) DEALLOCATE(fvdata(i)%CPEPOLH)
3555 IF (ASSOCIATED(fvdata(i)%CPFPOLH)) DEALLOCATE(fvdata(i)%CPFPOLH)
3556 ENDDO
3557 DEALLOCATE(fvdata)
3558 ENDIF
3559C
3560 IF (ALLOCATED(fvspmd)) THEN
3561 DO i=1,nfvbag
3562 IF (ASSOCIATED(fvspmd(i)%IBUF_L)) DEALLOCATE(fvspmd(i)%IBUF_L)
3563 IF (ASSOCIATED(fvspmd(i)%IBUFA_L)) DEALLOCATE(fvspmd(i)%IBUFA_L)
3564 IF (ASSOCIATED(fvspmd(i)%IBUFSA_L)) DEALLOCATE(fvspmd(i)%IBUFSA_L)
3565 IF (ASSOCIATED(fvspmd(i)%IXSA)) DEALLOCATE(fvspmd(i)%IXSA)
3566 IF (ASSOCIATED(fvspmd(i)%ELEMSA)) DEALLOCATE(fvspmd(i)%ELEMSA)
3567 ENDDO
3568 DEALLOCATE(fvspmd)
3569 ENDIF
3570C
3571 RETURN
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
fvbag_mod::nfvbag
integer nfvbag
Definition fvbag_mod.F:127

◆ fvdim()

subroutine fvdim ( type(monvol_struct_), dimension(nvolu), intent(inout) t_monvol)

Definition at line 3457 of file fvmesh.F.

3458C-----------------------------------------------
3459C M o d u l e s
3460C-----------------------------------------------
3461 USE fvbag_mod
3462 USE monvol_struct_mod
3463C-----------------------------------------------
3464C I m p l i c i t T y p e s
3465C-----------------------------------------------
3466#include "implicit_f.inc"
3467C-----------------------------------------------
3468C C o m m o n B l o c k s
3469C-----------------------------------------------
3470#include "com04_c.inc"
3471C-----------------------------------------------
3472C D u m m y A r g u m e n t s
3473C-----------------------------------------------
3474 TYPE(MONVOL_STRUCT_), DIMENSION(NVOLU), INTENT(INOUT) :: T_MONVOL
3475C-----------------------------------------------
3476C L o c a l V a r i a b l e s
3477C-----------------------------------------------
3478 INTEGER N, ITYP
3479C
3480 ALLOCATE(fvid(nvolu))
3481 fvid(:) = -1
3482 nfvbag=0
3483 DO n = 1, nvolu
3484 ityp = t_monvol(n)%TYPE
3485 IF (ityp == 6 .OR. ityp == 8) THEN
3486 nfvbag = nfvbag+1
3487 t_monvol(n)%IVOLU(45) = nfvbag
3488 t_monvol(n)%IVOLU(57) = 0
3489 fvid(nfvbag) = n
3490 ENDIF
3491 ENDDO
3492C
3493 RETURN
fvbag_mod::fvid
integer, dimension(:), allocatable fvid
Definition fvbag_mod.F:130
monvol_struct_mod
Definition monvol_struct_mod.F:66

◆ fvmesh1()

subroutine fvmesh1 ( integer, dimension(*) ibuf,
integer, dimension(3,*) elem,
x,
integer, dimension(*) ivolu,
integer, dimension(8,*) bric,
xb,
rvolu,
integer nel,
integer neli,
integer nbric,
integer, dimension(13,*) tbric,
sfac,
dxm,
integer nba,
integer nela,
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 id,
character(len=nchartitle) titr,
integer nb_node,
integer ityp )

Definition at line 48 of file fvmesh.F.

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