OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i22subvol.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| i22subvol ../engine/source/interfaces/int22/i22subvol.F
25!||--- called by ------------------------------------------------------
26!|| inttri ../engine/source/interfaces/intsort/inttri.F
27!||--- calls -----------------------------------------------------
28!|| i22aera ../engine/source/interfaces/int22/i22subvol.F
29!|| my_barrier ../engine/source/system/machine.F
30!||--- uses -----------------------------------------------------
31!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
32!|| element_mod ../common_source/modules/elements/element_mod.F90
33!|| groupdef_mod ../common_source/modules/groupdef_mod.F
34!|| i22bufbric_mod ../common_source/modules/interfaces/cut-cell-search_mod.F
35!|| i22edge_mod ../common_source/modules/interfaces/cut-cell-buffer_mod.F
36!|| i22tri_mod ../common_source/modules/interfaces/cut-cell-search_mod.F
37!||====================================================================
38 SUBROUTINE i22subvol(X ,NIN ,ITASK ,IPARI ,ITAB ,
39 . IXS ,IXTG ,V ,IPARG ,ELBUF_TAB,
40 . W ,IGRSH3N )
41C============================================================================
42C-----------------------------------------------
43C D e s c r i p t i o n
44C-----------------------------------------------
45C Interface Type22 (/INTER/TYPE22) is an FSI coupling method based on cut cell method.
46C This experimental cut cell method is not completed, abandoned, and is not an official option.
47C
48C CThis subroutine computes coordinates of intersection points and faces area
49C-----------------------------------------------
50C M o d u l e s
51C-----------------------------------------------
52 USE i22bufbric_mod
53 USE i22tri_mod
54 USE i22edge_mod
55 USE elbufdef_mod
56 USE groupdef_mod
57 use element_mod , only : nixs,nixtg
58C-----------------------------------------------
59C I m p l i c i t T y p e s
60C-----------------------------------------------
61#include "implicit_f.inc"
62#include "comlock.inc"
63C-----------------------------------------------
64C C o m m o n B l o c k s
65C-----------------------------------------------
66#include "com01_c.inc"
67#include "com04_c.inc"
68#include "param_c.inc"
69#include "task_c.inc"
70#include "subvolumes.inc"
71#include "inter22.inc"
72C-----------------------------------------------
73 INTERFACE
74 FUNCTION i22aera(NPTS,P, C)
75 INTEGER, INTENT(IN) :: NPTS
76 my_real, INTENT(IN) :: p(3,npts)
77 my_real, INTENT(INOUT) :: c(3)
78 my_real :: i22aera(3)
79 END FUNCTION
80 END INTERFACE
81C-----------------------------------------------
82C D u m m y A r g u m e n t s
83C-----------------------------------------------
84 INTEGER :: NIN, ITASK,IbugANIM22, IXS(NIXS,*), IPARI(48:50),
85 . ITAB(*),IXTG(NIXTG,*), IPARG(NPARG,*)
86 my_real, intent(inout), TARGET ::
87 . x(3,*),v(3,*),w(3,*)
88 TYPE(elbuf_struct_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
89!
90 TYPE (GROUP_) , DIMENSION(NGRSH3N) :: IGRSH3N
91C-----------------------------------------------
92C L o c a l V a r i a b l e s
93C-----------------------------------------------
94 INTEGER I, J, K,N, IAD(7), NBCUT, ICUT(8,7), IB, ICELL, NTRUS,INOD,
95 . IEDG(6) , KK , II , IAD0, IDSPRI, SECID, Fa(0:6),
96 . id, icode, idble,m, ie, nbf, nbl,nnod, mnod, nintpoint, m_sum,numpt,igr
97 my_real
98 . cutpoint(3), cutcoor, vmin, vmax, vol, vol1, vol2
99 CHARACTER*14 SECTYPE, SECTYPE2
100 my_real :: vF(1:3,0:6)
101! my_real,DIMENSION(:) , POINTER :: pface
102 TYPE(polyface_entity),dimension(:),pointer :: pFACE
103 my_real , POINTER :: pFACE0
104! INTEGER,DIMENSION(:) , POINTER :: pNNodCell
105! INTEGER,DIMENSION(:) , POINTER :: pNPointCell
106! INTEGER,DIMENSION(:) , POINTER :: pWhichCellNod
107 my_real,DIMENSION(:,:), POINTER :: pcut !=> BRICK_ENTITY%CUTCOOR(3,2)
108! my_real,DIMENSION(:) , POINTER :: pSUBVOL !=> BRICK_ENTITY%Vnew
109 TYPE(poly_entity), DIMENSION(:), POINTER :: pSUBVOL,pNNodCell,pNPointCell
110 TYPE(NODE_ENTITY), DIMENSION(:), POINTER :: pWhichCellNod
111 my_real,DIMENSION(:) , POINTER :: pvol !=> BRICK_ENTITY%VOL(2)
112 my_real,DIMENSION(:) , POINTER :: a1,a2,a3,a4,a5,a6
113 my_real :: wa1(3),wa2(3),wa3(3),wa4(3),wa5(3),wa6(3)
114 my_real :: wn1(3),wn2(3),wn3(3),wn4(3),wn5(3),wn6(3)
115 my_real :: wf1(3),wf2(3),wf3(3),wf4(3),wf5(3),wf6(3) , wf0(3)
116 my_real,DIMENSION(:) , POINTER :: n1,n2,n3,n4 ,n5,n6,n7,n8, nn
117 my_real :: c(1:3,0:6), scut, ctmp(3), valtmp(3)
118 my_real :: diag1(3), diag2(3), uncutvol, z(3,6), volratiomin
119
120 integer :: pNNodCell_bak(2), EdgNod1, EdgNod2
121 INTEGER :: Cod1, Cod2, Cod1_bis, Cod2_bis
122 LOGICAL :: bAMBIGOUS, IsSH3N, bTWICE
123 INTEGER :: INodTAG(8), IEtab(8),IEnod(1:4), IPCUT, ICRIT,ISHELL,NSH3N, NTRIA(1:3),NCELL
124 my_real :: pcut_vect(1:3,4), ievect(1:3,8), criteria(2), pt(3,3), volratio,norm,dist,tmp(3),scut1,scut2,beta
125 LOGICAL :: bCOMPL
126 INTEGER :: nFACE, IPOS, iCUTe(2,12), iPOSe(12), tagEDG, iCOMPL
127 INTEGER :: EdgePos(1:16,2), EdgeList(1:16) !12<2**4
128 INTEGER :: ITAG(1:8),NFL,NEL,NEL_B,NG,IDLOC
129 LOGICAL :: lFOUND
130 CHARACTER*14 :: tmp_SECTYPE(8)
131 INTEGER :: tmp_SecID_CELL(8),tmp_NumNOD(8),tmp_NumPOINT(8),NTAG
132
133 TYPE(poly_entity)tmpPOLY(8)
134 TYPE(CUT_PLANE)tmpCUT(8)
135C----------------------------------------------------------------
136
137!A TO PERFORM : SI PAS DE BRIQUE OU PLUS DE BRIQUE, METTRE QUAND MEME NODES A 0.
138 !----------------------------------------------------------------
139 !debug:
140 if (ibug22_subvol==1.AND.itask==0)then
141 print *,""
142 print *, "================================="
143 print *, "======== SUBVOL ==========="
144 print *, "================================="
145 endif
146
147 VolRatioMin = EM10
148
149 IGR = IPARI(49) ! sh3n group identifier
150 NEL = IPARI(50) ! sh3n group elements
151 II = 1 ! attention, first sh3n element of the group
152
153 NBF = 1+ITASK*NB/NTHREAD
154 NBL = (ITASK+1)*NB/NTHREAD
155
156 DO IB=NBF,NBL !1,NB
157
158 DO K=1,6
159 BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(K)%Center(1) = ZERO
160 BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(K)%Center(2) = ZERO
161 BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(K)%Center(3) = ZERO
162 ENDDO
163
164 ICODE = BRICK_LIST(NIN,IB)%ICODE
165 IDBLE = BRICK_LIST(NIN,IB)%IDBLE
166 bTWICE = .FALSE.
167 IF(ICODE==IDBLE)bTWICE = .TRUE.
168
169 SECID = 0
170 ID = BRICK_LIST(NIN,IB)%ID
171 pSUBVOL(1:9) => BRICK_LIST(NIN,IB)%POLY(1:9)!%Vnew
172 pSUBVOL(1:9)%Vnew = ZERO
173 pNNodCell(1:9) => BRICK_LIST(NIN,IB)%POLY(1:9)!%NumNOD
174 pNNodCell(1:9)%NumNOD = 0
175 pNPointCell(1:9) => BRICK_LIST(NIN,IB)%POLY(1:9)!%NumPOINT
176 pNPointCell(1:9)%NumPOINT = 0
177 !pListNodID => BRICK_LIST(NIN,IB)%POLY(1:9)%ListNodID(1:8)
178
179 BRICK_LIST(NIN,IB)%POLY(1:9)%ListNodID(1) = 0
180 BRICK_LIST(NIN,IB)%POLY(1:9)%ListNodID(2) = 0
181 BRICK_LIST(NIN,IB)%POLY(1:9)%ListNodID(3) = 0
182 BRICK_LIST(NIN,IB)%POLY(1:9)%ListNodID(4) = 0
183 BRICK_LIST(NIN,IB)%POLY(1:9)%ListNodID(5) = 0
184 BRICK_LIST(NIN,IB)%POLY(1:9)%ListNodID(6) = 0
185 BRICK_LIST(NIN,IB)%POLY(1:9)%ListNodID(7) = 0
186 BRICK_LIST(NIN,IB)%POLY(1:9)%ListNodID(8) = 0
187
188 BRICK_LIST(NIN,IB)%POLY(1:9)%FACE0%Surf = ZERO
189 DO K=1,6
190 BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(K)%Surf = ZERO
191 ENDDO
192
193 NINTPOINT = 0
194 M_SUM = 0
195 pWhichCellNod(1:8) => BRICK_LIST(NIN,IB)%NODE(1:8)!%WhichCell
196 pWhichCellNod(1:8)%WhichCell = 0
197 EdgeList(1:16) = 0
198 EdgePos(1:16,1:2) = 0
199
200 DO K=1,6
201 BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(K)%Vel(1) = ZERO
202 BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(K)%Vel(2) = ZERO
203 BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(K)%Vel(3) = ZERO
204 ENDDO
205 BRICK_LIST(NIN,IB)%PCUT(1:9)%Vel(1) = ZERO
206 BRICK_LIST(NIN,IB)%PCUT(1:9)%Vel(2) = ZERO
207 BRICK_LIST(NIN,IB)%PCUT(1:9)%Vel(3) = ZERO
208
209 !--------------------------------------------------------------------------------!
210 ! 1. CHECK FIRST IF AMBIGUOUS COMBINATION !
211 !--------------------------------------------------------------------------------!
212 !bAMBIGOUS == .TRUE. => ambiguous combination detected
213 ! (code1_bis, code2_bis) is then the alternative combination
214 IF(BRICK_LIST(NIN,IB)%NBCUT == 2)THEN
215 SECTYPE = BRICK_LIST(NIN,IB)%SECTYPE(1)
216 Cod1 = IABS(BRICK_LIST(NIN,IB)%SecID_CELL(1))
217 bAMBIGOUS = .FALSE.
218 IF(SECTYPE(1:5)=='PENTA')THEN
219 SECTYPE = BRICK_LIST(NIN,IB)%SECTYPE(2)
220 Cod2 = IABS(BRICK_LIST(NIN,IB)%SecID_CELL(2))
221 IF(SECTYPE(1:5)=='PENTA')THEN
222.AND. IF (Cod1==10 Cod2==17 )THEN
223 bAMBIGOUS = .TRUE.
224 Cod1_bis = 12
225 Cod2_bis = 19
226.AND. ELSEIF (Cod1==11 Cod2==14 )THEN
227 bAMBIGOUS = .TRUE.
228 Cod1_bis = 15
229 Cod2_bis = 18
230.AND. ELSEIF (Cod1==12 Cod2==19 )THEN
231 bAMBIGOUS = .TRUE.
232 Cod1_bis = 10
233 Cod2_bis = 17
234.AND. ELSEIF (Cod1==13 Cod2==16 )THEN
235 bAMBIGOUS = .TRUE.
236 Cod1_bis = 09
237 Cod2_bis = 20
238.AND. ELSEIF (Cod1==15 Cod2==18 )THEN
239 bAMBIGOUS = .TRUE.
240 Cod1_bis = 11
241 Cod2_bis = 14
242.AND. ELSEIF (Cod1== 09 Cod2==20 )THEN
243 bAMBIGOUS = .TRUE.
244 Cod1_bis = 13
245 Cod2_bis = 16
246 ENDIF
247 ENDIF
248 ENDIF
249 !CHECK SURFACE NORMAL AND SOLVE THE AMBIGUITY
250.NOT. IF(bAMBIGOUS)THEN
251 !print *, "not ambiguous"
252 ELSE
253 !print *, "ambiguous=", Cod1,Cod2
254 !----------------------------------------------------------------
255 !RETRIEVING LAGRANGIAN FACES & CUT PLANE NORMAL
256 !----------------------------------------------------------------
257 SECID = Cod1
258 IAD(1:4) = (/((IB-1)*12+IABS(Gcorner(K,SECID)),K=1,4)/)
259 IAD(5) = IAD(1)
260 ICUT(1:2,1:4) = 1
261 DO K=1,4
262.AND. IF(EDGE_LIST(NIN,IAD(K))%NBCUT > 1 Gcorner(K,SECID) > 0)ICUT(1,K) = 2 !should be 1 but formulation might be updated later, keep general form using edge orientation (ICUT=2 possible).
263CC print *,"---1", IAD(K),ICUT(1,K), EDGE_LIST(NIN,IAD(K))%CUTSHELL(ICUT(1,K))
264 IEtab(K) = EDGE_LIST(NIN,IAD(K))%CUTSHELL(ICUT(1,K))
265 CUTCOOR = EDGE_LIST(NIN,IAD(K))%CUTCOOR(ICUT(1,K)) ! manage double intersections see comment above
266 CUTPOINT(1:3) = X(1:3, EDGE_LIST(NIN,IAD(K))%NODE(1) ) + CUTCOOR * (EDGE_LIST(NIN,IAD(K))%VECTOR(1:3))
267 EDGE_LIST(NIN,IAD(K))%CUTPOINT(1:3,ICUT(1,1)) = CUTPOINT(1:3)
268 END DO !NEXT K
269 a1 => EDGE_LIST(NIN,IAD(1))%CUTPOINT(1:3,ICUT(1,1))
270 a2 => EDGE_LIST(NIN,IAD(2))%CUTPOINT(1:3,ICUT(1,2))
271 a3 => EDGE_LIST(NIN,IAD(3))%CUTPOINT(1:3,ICUT(1,3))
272 a4 => EDGE_LIST(NIN,IAD(4))%CUTPOINT(1:3,ICUT(1,4))
273 vF(:,0) = I22AERA(4,(/a1,a2,a3,a4/),C(1:3, 0 )) !Cut surface normal vector (not unitary : 2S factor)
274 ! cannot be zero if intersections are limited to cutcoor = em06 and one-em06 in i22intersect.f (0 and 1.00 impossible => no degenerate penta with scut=0.000, and no division by 0)
275 PCut_vect(1:3,1) = vF(:,0) / SQRT(SUM(vF(:,0)*vF(:,0)))
276 !----------------------------------------------------------------
277 SECID = Cod2
278 IAD(1:4) = (/((IB-1)*12+IABS(Gcorner(K,SECID)),K=1,4)/)
279 IAD(5) = IAD(1)
280 DO K=1,4
281.AND. IF(EDGE_LIST(NIN,IAD(K))%NBCUT > 1 Gcorner(K,SECID) > 0)ICUT(2,K) = 2 !should be 1 but formulation might be updated later, keep general form using edge orientation (ICUT=2 possible).
282CC print *,"---2", IAD(K),ICUT(2,K), EDGE_LIST(NIN,IAD(K))%CUTSHELL(ICUT(2,K))
283 IEtab(4+K) = EDGE_LIST(NIN,IAD(K))%CUTSHELL(ICUT(2,K))
284 CUTCOOR = EDGE_LIST(NIN,IAD(K))%CUTCOOR(ICUT(2,K)) ! manage double intersections see comment above
285 CUTPOINT(1:3) = X(1:3, EDGE_LIST(NIN,IAD(K))%NODE(1) ) + CUTCOOR * (EDGE_LIST(NIN,IAD(K))%VECTOR(1:3))
286 EDGE_LIST(NIN,IAD(K))%CUTPOINT(1:3,ICUT(2,1)) = CUTPOINT(1:3)
287 END DO !NEXT K
288 a1 => EDGE_LIST(NIN,IAD(1))%CUTPOINT(1:3,ICUT(2,1))
289 a2 => EDGE_LIST(NIN,IAD(2))%CUTPOINT(1:3,ICUT(2,2))
290 a3 => EDGE_LIST(NIN,IAD(3))%CUTPOINT(1:3,ICUT(2,3))
291 a4 => EDGE_LIST(NIN,IAD(4))%CUTPOINT(1:3,ICUT(2,4))
292 vF(:,0) = I22AERA(4,(/a1,a2,a3,a4/),C(1:3, 0 )) !Cut surface normal vector (not unitary : 2S factor)
293 ! cannot be zero if intersections are limited to cutcoor = em06 and one-em06 in i22intersect.f (0 and 1.00 impossible => no degenerate penta with scut=0.000, and no division by 0)
294 PCut_vect(1:3,2) = vF(:,0) / SQRT(SUM(vF(:,0)*vF(:,0)))
295 !----------------------------------------------------------------
296 SECID = Cod1_bis
297 IAD(1:4) = (/((IB-1)*12+IABS(Gcorner(K,SECID)),K=1,4)/)
298 IAD(5) = IAD(1)
299 ICUT(1,1:4) = 1
300 DO K=1,4
301.AND. IF(EDGE_LIST(NIN,IAD(K))%NBCUT > 1 Gcorner(K,SECID) > 0)ICUT(1,K) = 2 !should be 1 but formulation might be updated later, keep general form using edge orientation (ICUT=2 possible).
302 !IEtab(4+K) = EDGE_LIST(NIN,IAD(K))%CUTSHELL(ICUT(1,K))
303 CUTCOOR = EDGE_LIST(NIN,IAD(K))%CUTCOOR(ICUT(1,K)) ! manage double intersections see comment above
304 CUTPOINT(1:3) = X(1:3, EDGE_LIST(NIN,IAD(K))%NODE(1) ) + CUTCOOR * (EDGE_LIST(NIN,IAD(K))%VECTOR(1:3))
305 EDGE_LIST(NIN,IAD(K))%CUTPOINT(1:3,ICUT(1,1)) = CUTPOINT(1:3)
306 END DO !NEXT K
307 a1 => EDGE_LIST(NIN,IAD(1))%CUTPOINT(1:3,ICUT(1,1))
308 a2 => EDGE_LIST(NIN,IAD(2))%CUTPOINT(1:3,ICUT(1,2))
309 a3 => EDGE_LIST(NIN,IAD(3))%CUTPOINT(1:3,ICUT(1,3))
310 a4 => EDGE_LIST(NIN,IAD(4))%CUTPOINT(1:3,ICUT(1,4))
311 vF(:,0) = I22AERA(4,(/a1,a2,a3,a4/),C(1:3, 0 )) !Cut surface normal vector (not unitary : 2S factor)
312 ! cannot be zero if intersections are limited to cutcoor = em06 and one-em06 in i22intersect.f (0 and 1.00 impossible => no degenerate penta with scut=0.000, and no division by 0)
313 PCut_vect(1:3,3) = vF(:,0) / SQRT(SUM(vF(:,0)*vF(:,0)))
314 !----------------------------------------------------------------
315 SECID = Cod2_bis
316 IAD(1:4) = (/((IB-1)*12+IABS(Gcorner(K,SECID)),K=1,4)/)
317 IAD(5) = IAD(1)
318 ICUT(2,1:4) = 1
319 DO K=1,4
320.AND. IF(EDGE_LIST(NIN,IAD(K))%NBCUT > 1 Gcorner(K,SECID) > 0)ICUT(2,K) = 2 !should be 1 but formulation might be updated later, keep general form using edge orientation (ICUT=2 possible).
321 !IEtab(4+K) = EDGE_LIST(NIN,IAD(K))%CUTSHELL(ICUT(2,K))
322 CUTCOOR = EDGE_LIST(NIN,IAD(K))%CUTCOOR(ICUT(2,K)) ! manage double intersections see comment above
323 CUTPOINT(1:3) = X(1:3, EDGE_LIST(NIN,IAD(K))%NODE(1) ) + CUTCOOR * (EDGE_LIST(NIN,IAD(K))%VECTOR(1:3))
324 EDGE_LIST(NIN,IAD(K))%CUTPOINT(1:3,ICUT(2,1)) = CUTPOINT(1:3)
325 END DO !NEXT K
326 a1 => EDGE_LIST(NIN,IAD(1))%CUTPOINT(1:3,ICUT(2,1))
327 a2 => EDGE_LIST(NIN,IAD(2))%CUTPOINT(1:3,ICUT(2,2))
328 a3 => EDGE_LIST(NIN,IAD(3))%CUTPOINT(1:3,ICUT(2,3))
329 a4 => EDGE_LIST(NIN,IAD(4))%CUTPOINT(1:3,ICUT(2,4))
330 vF(:,0) = I22AERA(4,(/a1,a2,a3,a4/),C(1:3, 0 )) !Cut surface normal vector (not unitary : 2S factor)
331 ! cannot be zero if intersections are limited to cutcoor = em06 and one-em06 in i22intersect.f (0 and 1.00 impossible => no degenerate penta with scut=0.000, and no division by 0)
332 PCut_vect(1:3,4) = vF(:,0) / SQRT(SUM(vF(:,0)*vF(:,0)))
333 !----------------------------------------------------------------
334 !COMPUTING NORMALS AT LAGRANGIAN FACES
335 !----------------------------------------------------------------
336 DO I=1,8
337 IE = IEtab(I)
338 IEnod(3:4)=IRECT_L(3:4,IE)
339 IF(IEnod(3)==IEnod(4))THEN
340 IsSH3N=.true.
341 Diag1(1:3) = IRECT_L((/5,9,13/),IE) - IRECT_L((/6,10,14/),IE)
342 Diag2(1:3) = IRECT_L((/5,9,13/),IE) - IRECT_L((/7,11,15/),IE)
343 IEvect(1,I) = + Diag1(2)*Diag2(3) - diag1(3)*Diag2(2)
344 IEvect(2,I) = + Diag1(3)*Diag2(1) - diag1(1)*Diag2(3)
345 IEvect(3,I) = + Diag1(1)*Diag2(2) - diag1(2)*Diag2(1)
346 IEvect(1:3,I) = HALF * IEvect(1:3,I) / SQRT(SUM( IEvect(1:3,I) * IEvect(1:3,I) ))
347 ELSE
348 IsSH3N=.false.
349 Diag1(1:3) = IRECT_L((/5, 9,13/),IE) - IRECT_L((/7,11,15/),IE)
350 Diag2(1:3) = IRECT_L((/6,10,14/),IE) - IRECT_L((/8,12,16/),IE)
351 IEvect(1,I) = + Diag1(2)*Diag2(3) - diag1(3)*Diag2(2)
352 IEvect(2,I) = + Diag1(3)*Diag2(1) - diag1(1)*Diag2(3)
353 IEvect(3,I) = + Diag1(1)*Diag2(2) - diag1(2)*Diag2(1)
354 IEvect(1:3,I) = HALF * IEvect(1:3,I) / SQRT(SUM( IEvect(1:3,I) * IEvect(1:3,I) ))
355 ENDIF
356 ENDDO
357 !----------------------------------------------------------------
358 !COMPUTING DETERMINATION CRITERIA
359 !----------------------------------------------------------------
360 CRITERIA(:) = ZERO
361 DO ICRIT = 1,2
362 DO IPCUT = 1,2
363 DO J = 1,4
364 CRITERIA(ICRIT) = CRITERIA(ICRIT) + ABS( SUM(IEvect(1:3,(IPCUT-1)*4+J)*PCut_vect( 1:3 , (ICRIT-1)*2 + IPCUT) ) )
365 ENDDO! J
366 ENDDO! PCUT
367 ENDDO! ICRIT
368 !----------------------------------------------------------------
369 !DETERMINING THE CORRECT COMBINATION
370 !----------------------------------------------------------------
371 IF (CRITERIA(1) > CRITERIA(2))THEN
372 !keep it
373 ELSEIF(CRITERIA(2) > CRITERIA(1))THEN
374 !change it
375 BRICK_LIST(NIN,IB)%SecID_CELL(1:2) = (/ Cod1_bis, Cod2_bis /)
376 BRICK_LIST(NIN,IB)%SECTYPE(1) = StrCode(Cod1_bis)
377 BRICK_LIST(NIN,IB)%SECTYPE(2) = StrCode(Cod2_bis)
378CC print *, "CHANGING RETAINED COMBINATION TO :",StrCode(cod1_bis),StrCode(cod2_bis)
379 ELSE
380CC print *, "inter22 : ambigus intersection not solved" ;
381 ENDIF
382 !----------------------------------------------------------------
383 ENDIF!(bAMBIGOUS)
384 ENDIF!(BRICK_LIST(NIN,IB)%NBCUT == 2)
385
386
387 !--------------------------------------------------------------------------------!
388 ! 2. COMPUTE POLYHEDRON VOLUME AND FACES !
389 !--------------------------------------------------------------------------------!
390 NCELL = BRICK_LIST(NIN,IB)%NBCUT
391 ICUTe(1:2,1:12) = 0
392
393 !normilized volume
394 UNCUTVOL = ONE
395 DO NG=1,NGROUP
396 NEL_B=IPARG(2,NG)
397 NFL =IPARG(3,NG)
398.AND. IF( ID>NFL ID<=NFL+NEL_B)THEN
399 IDLOC = ID - NFL
400 UNCUTVOL = ELBUF_TAB(NG)%GBUF%VOL(IDLOC)
401 EXIT
402 ENDIF
403 ENDDO!next ng
404
405 ITAG(1:8) = 0
406
407 DO ICELL=1, NCELL
408 M = 0
409 NTRUS = 0
410 NNOD = 0
411 MNOD = 0
412 SECTYPE = BRICK_LIST(NIN,IB)%SECTYPE(ICELL)
413 SECID = BRICK_LIST(NIN,IB)%SecID_Cell(ICELL)
414 bCOMPL = .FALSE.
415 IF(SECID<0)THEN
416 bCOMPL = .TRUE.
417 SECID = - SECID
418 ENDIF
419 INodTAG(1:8) = 0
420 IF(SECTYPE(1:5)=='TETRA')THEN
421 NTRUS= C_TETRA
422 M = M_TETRA
423 NNOD = N_TETRA+C_TETRA
424 MNOD = N_TETRA
425 ELSEIF(SECTYPE(1:5)=='PENTA')THEN
426 NTRUS= C_PENTA
427 M = M_PENTA
428 NNOD = N_PENTA+C_PENTA
429 MNOD = N_PENTA
430 ELSEIF(SECTYPE(1:5)=='POLY3')THEN
431 NTRUS= C_POLY3
432 M = M_POLY3
433 NNOD = N_POLY3+C_POLY3
434 MNOD = N_POLY3
435 ELSEIF(SECTYPE(1:5)=='HEXAE')THEN
436 NTRUS= C_HEXAE
437 M = M_HEXAE
438 NNOD = N_HEXAE+C_HEXAE
439 MNOD = N_HEXAE
440 ELSEIF(SECTYPE(1:6)=='POLY4 ')THEN
441 NTRUS= C_POLY4
442 M = M_POLY4
443 NNOD = N_POLY4+C_POLY4
444 MNOD = N_POLY4
445 ELSEIF(SECTYPE(1:6)=='POLY4A')THEN
446 NTRUS= C_POLY4A
447 M = M_POLY4A
448 NNOD = N_POLY4A+C_POLY4A
449 MNOD = N_POLY4A
450 ELSEIF(SECTYPE(1:6)=='POLY4B')THEN
451 NTRUS= C_POLY4B
452 M = M_POLY4B
453 NNOD = N_POLY4B+C_POLY4B
454 MNOD = N_POLY4B
455 ELSEIF(SECTYPE(1:5)=='POLYC')THEN
456 NTRUS= C_POLYC
457 M = M_POLYC
458 NNOD = N_POLYC+C_POLYC
459 MNOD = N_POLYc
460 END IF
461 M_SUM = M_SUM + M
462 NINTPOINT = NINTPOINT + NTRUS
463 IAD(1:NTRUS) = (/((IB-1)*12+IABS(Gcorner(K,SECID)),K=1,NTRUS)/)
464 IAD(NTRUS+1) = IAD(1)
465 DO K=1,NTRUS
466 EdgeList(IABS(Gcorner(K,SECID))) = EdgeList(IABS(Gcorner(K,SECID))) + 1
467 ENDDO
468
469 !----------------------------------------------------------------
470 !calculation of intersection point coordinates
471 !----------------------------------------------------------------
472 !if twice the same intersections, the 2nd one is the biggest one choose farthest intersection points
473 DO K=1,NTRUS
474 tagEDG = Gcorner(K,SECID)
475 IF(tagEDG<0)THEN
476 !Negative orientation: we take the closest
477 ipos = 1
478 tagEDG = - tagEDG
479 ELSE
480 !We take the most
481 IF(EDGE_LIST(NIN,IAD(K))%NBCUT==1)THEN
482 ipos = 1
483 ELSE
484 ipos = 2
485 ENDIF
486 ENDIF
487 EdgePos(IABS(tagEdg), ipos) = 1
488 IF(iCUTe(ipos,tagEDG)==1)ipos = MOD(ipos,2)+1 !already taken by Polyhedre precedent, we take the other
489 CUTCOOR = EDGE_LIST(NIN,IAD(K))%CUTCOOR(ipos)
490 CUTPOINT(1:3) = X(1:3, EDGE_LIST(NIN,IAD(K))%NODE(1) ) + CUTCOOR * (EDGE_LIST(NIN,IAD(K))%VECTOR(1:3))
491 IE = EDGE_LIST(NIN,IAD(K))%CUTSHELL(ipos)
492 EDGE_LIST(NIN,IAD(K))%CUTPOINT(1:3,ipos) = CUTPOINT(1:3)
493 EDGE_LIST(NIN,IAD(K))%CUTVEL(1:3,ipos) = IRECT_L(24:26,IE)
494 ICUTe(ipos, IABS(Gcorner(K,SECID))) = 1 !We will take the following intersection point.
495 iPOSe(K) = ipos
496 END DO !K
497! print *, "gcorner = ", gcorner(1:ntrus,secid)
498! print *, "iPOSe = ", iPOSe(1:NTRUS)
499
500 if( brick_list(nin,ib)%ICODE==0.and.icell>=1)then
501 print *, "error i22subvol."
502 stop 221
503 end if
504 !----------------------------------------------------------------
505 ! ANIMATION (positioning)
506 !----------------------------------------------------------------
507 nsh3n = getnumtria(getpolyhedratype(secid))
508 IF(ii<=nel-nsh3n.AND.nspmd<=1)THEN !ii=iad0=0 if no grshel, otherwise means there are still unused shells in the group, at least 4 to make a poly4 section
509
510C if (IBUG22_subvol==1)WRITE (*,FMT='(A20,I10)')
511C . "path intersection : ",IXS(11,ID)
512#include "lockon.inc"
513 !position the shells
514 !NSH3N = GetNumTria(GetPolyhedraType(SECID))
515c print *,"i22subvol : brickID=", IXS(11,BRICK_LIST(NIN,IB)%ID)
516c print *,"i22subvol : icode =" , BRICK_LIST(NIN,IB)%ICODE
517c print *,"i22subvol : idble =" , BRICK_LIST(NIN,IB)%IDBLE
518c print *,"i22subvol : bTWIC =" , bTWICE
519c print *,"i22subvol : icell =" , icell
520c print *,"i22subvol : typ =" , BRICK_LIST(NIN,IB)%SECTYPE(ICELL) , BRICK_LIST(NIN,IB)%SECID_Cell(ICELL)
521
522 DO k=1,nsh3n
523 ntria(1:3) = gtria(1:3,k,getpolyhedratype(secid))
524 x(1:3,ixtg(2,igrsh3n(igr)%ENTITY(ii))) = edge_list(nin,iad(ntria(1)) )%CUTPOINT(1:3,ipose(ntria(1)) )
525 x(1:3,ixtg(3,igrsh3n(igr)%ENTITY(ii))) = edge_list(nin,iad(ntria(2)) )%CUTPOINT(1:3,ipose(ntria(2)) )
526 x(1:3,ixtg(4,igrsh3n(igr)%ENTITY(ii))) = edge_list(nin,iad(ntria(3)) )%CUTPOINT(1:3,ipose(ntria(3)) )
527 v(1:3,ixtg(2,igrsh3n(igr)%ENTITY(ii))) = edge_list(nin,iad(ntria(1)) )%CUTVEL(1:3,ipose(ntria(1)) )
528 v(1:3,ixtg(3,igrsh3n(igr)%ENTITY(ii))) = edge_list(nin,iad(ntria(2)) )%CUTVEL(1:3,ipose(ntria(2)) )
529 v(1:3,ixtg(4,igrsh3n(igr)%ENTITY(ii))) = edge_list(nin,iad(ntria(3)) )%CUTVEL(1:3,ipose(ntria(3)) )
530C print *, " : sh3n_id=", IXTG(NIXTG,IBUFSSG(II))
531
532c if (IBUG22_subvol==1) then
533c pt(1,:) = EDGE_LIST(NIN,IAD(Ntria(1)) )%CUTPOINT(1:3,iPOSe(Ntria(1)))
534c pt(2,:) = EDGE_LIST(NIN,IAD(Ntria(2)) )%CUTPOINT(1:3,iPOSe(Ntria(2)))
535c pt(3,:) = EDGE_LIST(NIN,IAD(Ntria(3)) )%CUTPOINT(1:3,iPOSe(Ntria(3)))
536c WRITE(*,FMT='(I12,A,F16.9,A,F16.9,A,F16.9,A)') IXTG(NIXTG,IBUFSSG(II)),":(",pt(1,1),",",pt(1,2),",",pt(1,3),")"
537c WRITE(*,FMT='(I12,A,F16.9,A,F16.9,A,F16.9,A)') IXTG(NIXTG,IBUFSSG(II)),":(",pt(2,1),",",pt(2,2),",",pt(2,3),")"
538c WRITE(*,FMT='(I12,A,F16.9,A,F16.9,A,F16.9,A)') IXTG(NIXTG,IBUFSSG(II)),":(",pt(3,1),",",pt(3,2),",",pt(3,3),")"
539c endif
540C print *," ntria =", K
541C print *," II =", II
542
543 ii=ii+1 ! next sh3n
544 END DO !next K
545C if (IBUG22_subvol==1)WRITE (*,FMT='(A1)') " "
546#include "lockoff.inc"
547 END IF !IAD>0
548 !----------------------------------------------------------------
549
550 !----------------------------------------------------------------
551 ! SUBVOLUME AERAS & VOLUMES CALCULATION
552 ! It is first necessary to calculate vector areas
553 !----------------------------------------------------------------
554 fa(0) = 0
555 !OPTIM: Remove ABS if graph all oriented outgoing normal.
556
557 IF(sectype(1:5)=='TETRA')THEN
558 nface = f_tetra
559 !intersection points
560 a1 => edge_list(nin,iad(1))%CUTPOINT(1:3,ipose(1))
561 a2 => edge_list(nin,iad(2))%CUTPOINT(1:3,ipose(2))
562 a3 => edge_list(nin,iad(3))%CUTPOINT(1:3,ipose(3))
563 n1 => x(1:3,ixs(gnode(1,secid)+1,id))
564 fa(1:3) = gface(1:3,secid)
565 vf(:,1) = i22aera(3,(/n1,a2,a1/), c(1:3,fa(1)))
566 vf(:,2) = i22aera(3,(/n1,a3,a2/), c(1:3,fa(2)))
567 vf(:,3) = i22aera(3,(/n1,a1,a3/), c(1:3,fa(3)))
568 vf(:,0) = i22aera(3,(/a1,a2,a3/), c(1:3, 0 ))
569 vol = sum((/(vf(:,i)*c(:,fa(i)),i=0,3)/))
570 vol = third * abs(vol)
571 volratio = vol/uncutvol
572 IF(volratio <= volratiomin) THEN
573 itag(icell)=1
574 nintpoint = nintpoint - ntrus
575 !print *,"cleaning 1"
576 cycle
577 ENDIF
578 brick_list(nin,ib)%POLY(9)%FACE(fa(1))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(1))%Center(1:3) + a2 + a1 !put sum intersection point in icell 9 (will be used in sinit22_fvm.F to compute faceCenter for poly9
579 brick_list(nin,ib)%POLY(9)%FACE(fa(2))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(2))%Center(1:3) + a3 + a2
580 brick_list(nin,ib)%POLY(9)%FACE(fa(3))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(3))%Center(1:3) + a1 + a3
581 brick_list(nin,ib)%POLY(icell)%FACE(fa(1))%Center(1:3) = c(1:3,fa(1))
582 brick_list(nin,ib)%POLY(icell)%FACE(fa(2))%Center(1:3) = c(1:3,fa(2))
583 brick_list(nin,ib)%POLY(icell)%FACE(fa(3))%Center(1:3) = c(1:3,fa(3))
584 brick_list(nin,ib)%POLY(icell)%FACE0%Center(1:3) = c(1:3, 0 )
585 scut = sqrt(sum(vf(:,0)*vf(:,0)))
586 pface(1:6) => brick_list(nin,ib)%POLY(icell)%FACE(1:6)!%Surf
587 pface0=> brick_list(nin,ib)%POLY(icell)%FACE0%Surf
588 pface(fa(1))%Surf = sqrt(sum(vf(:,1)*vf(:,1))) !pointers indexes always start from 1
589 pface(fa(2))%Surf = sqrt(sum(vf(:,2)*vf(:,2)))
590 pface(fa(3))%Surf = sqrt(sum(vf(:,3)*vf(:,3)))
591 pface0 = scut
592 brick_list(nin,ib)%PCUT(icell)%SCUT = scut !area of the cutting surface
593 brick_list(nin,ib)%PCUT(icell)%B(1:3) = c(1:3,0) !centroid of the cutting surface (i.e. basis point)
594 brick_list(nin,ib)%PCUT(icell)%N(1:3) = vf(:,0) !Cut surface normal vector (not unitary : 2S factor)
595 brick_list(nin,ib)%PCUT(icell)%NP = 3 !polygon points (number)
596 brick_list(nin,ib)%PCUT(icell)%P(1:3,1)= a1 !polygon points (coordinates)
597 brick_list(nin,ib)%PCUT(icell)%P(1:3,2)= a2 !polygon points (coordinates)
598 brick_list(nin,ib)%PCUT(icell)%P(1:3,3)= a3 !polygon points (coordinates)
599 psubvol(icell)%Vnew = vol
600 pnnodcell(icell)%NumNOD = mnod !main(brick nodes)
601 pnpointcell(icell)%NumPOINT= nnod !main + intersection
602 brick_list(nin,ib)%POLY(icell)%CellCenter(1)= sum( (/a1(1),a2(1),a3(1),n1(1)/) ) / four
603 brick_list(nin,ib)%POLY(icell)%CellCenter(2)= sum( (/a1(2),a2(2),a3(2),n1(2)/) ) / four
604 brick_list(nin,ib)%POLY(icell)%CellCenter(3)= sum( (/a1(3),a2(3),a3(3),n1(3)/) ) / four
605 brick_list(nin,ib)%POLY(icell)%ListNodID(1:mnod) = gnode(1:mnod,secid) !IXS(Gnode(1,SECID)+1,ID)
606 inodtag(gnode(1:mnod,secid)) = 1 !tag local brick nodes [1 to 8] if used by polyhedron
607 pwhichcellnod(gnode(1:mnod,secid))%WhichCell = icell
608 !VELOCITY ON POLYHEDRA FACES
609 a1 => edge_list(nin,iad(1))%CUTVEL(1:3,ipose(1))
610 a2 => edge_list(nin,iad(2))%CUTVEL(1:3,ipose(2))
611 a3 => edge_list(nin,iad(3))%CUTVEL(1:3,ipose(3))
612 n1 => v(1:3,ixs(gnode(1,secid)+1,id))
613 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(1),ICELL) = THIRD* (n1(1:3) + a2(1:3) + a1(1:3))
614 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(2),ICELL) = THIRD* (n1(1:3) + a3(1:3) + a2(1:3))
615 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(3),ICELL) = THIRD* (n1(1:3) + a1(1:3) + a3(1:3))
616 !BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE0%Vel(1:3) = THIRD* (a1(1:3) + a2(1:3) + a3(1:3))
617 brick_list(nin,ib)%PCUT(icell)%Vel(1:3) = third* (a1(1:3) + a2(1:3) + a3(1:3))
618
619 ELSEIF(sectype(1:5)=='PENTA')THEN
620 nface = f_penta
621 a1 => edge_list(nin,iad(1))%CUTPOINT(1:3,ipose(1))
622 a2 => edge_list(nin,iad(2))%CUTPOINT(1:3,ipose(2))
623 a3 => edge_list(nin,iad(3))%CUTPOINT(1:3,ipose(3))
624 a4 => edge_list(nin,iad(4))%CUTPOINT(1:3,ipose(4))
625 n1 => x(1:3,ixs(gnode(1,secid)+1,id))
626 n2 => x(1:3,ixs(gnode(2,secid)+1,id))
627 fa(1:4) = gface(1:4,secid)
628 vf(:,1) = i22aera(3,(/n1,a2,a1/), c(1:3,fa(1)))
629 vf(:,2) = i22aera(3,(/n2,a4,a3/), c(1:3,fa(2)))
630 vf(:,3) = i22aera(4,(/n1,n2,a3,a2/),c(1:3,fa(3)))
631 vf(:,4) = i22aera(4,(/n2,n1,a1,a4/),c(1:3,fa(4)))
632 vf(:,0) = i22aera(4,(/a1,a2,a3,a4/),c(1:3, 0 ))
633 vol = sum((/(vf(:,i)*c(1:3,fa(i)),i=0,4)/))
634 vol = third * abs(vol)
635 volratio = vol/uncutvol
636 brick_list(nin,ib)%POLY(9)%FACE(fa(1))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(1))%Center(1:3)+ a2 + a1 !put sum intersection point in icell 9 (will be used in sinit22_fvm.F to compute faceCenter for poly9
637 brick_list(nin,ib)%POLY(9)%FACE(fa(2))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(2))%Center(1:3)+ a4 + a3
638 brick_list(nin,ib)%POLY(9)%FACE(fa(3))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(3))%Center(1:3)+ a3 + a2
639 brick_list(nin,ib)%POLY(9)%FACE(fa(4))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(4))%Center(1:3)+ a1 + a4
640 brick_list(nin,ib)%POLY(icell)%FACE(fa(1))%Center(1:3) = c(1:3,fa(1))
641 brick_list(nin,ib)%POLY(icell)%FACE(fa(2))%Center(1:3) = c(1:3,fa(2))
642 brick_list(nin,ib)%POLY(icell)%FACE(fa(3))%Center(1:3) = c(1:3,fa(3))
643 brick_list(nin,ib)%POLY(icell)%FACE(fa(4))%Center(1:3) = c(1:3,fa(4))
644 brick_list(nin,ib)%POLY(icell)%FACE0%Center(1:3) = c(1:3, 0 )
645 scut = sqrt(sum(vf(:,0)*vf(:,0)))
646 pface(1:6) => brick_list(nin,ib)%POLY(icell)%FACE(1:6)!%Surf
647 pface0=> brick_list(nin,ib)%POLY(icell)%FACE0%Surf
648 pface(fa(1))%Surf = sqrt(sum(vf(:,1)*vf(:,1)))
649 pface(fa(2))%Surf = sqrt(sum(vf(:,2)*vf(:,2)))
650 pface(fa(3))%Surf = sqrt(sum(vf(:,3)*vf(:,3)))
651 pface(fa(4))%Surf = sqrt(sum(vf(:,4)*vf(:,4)))
652 pface0 = scut
653 IF(volratio <= volratiomin) THEN
654 IF(scut / exp(0.666*log(uncutvol)) <= em03)THEN
655 itag(icell)=1
656 nintpoint = nintpoint - ntrus
657 !print *,"cleaning 2"
658 cycle
659 ENDIF
660 ENDIF
661 brick_list(nin,ib)%PCUT(icell)%SCUT = scut !area of the cutting surface
662 brick_list(nin,ib)%PCUT(icell)%B(1:3) = c(1:3,0) !centroid of the cutting surface (i.e. basis point)
663 brick_list(nin,ib)%PCUT(icell)%N(1:3) = vf(:,0) !Cut surface normal vector (not unitary : 2S factor)
664 brick_list(nin,ib)%PCUT(icell)%NP = 4 !polygon points (number)
665 brick_list(nin,ib)%PCUT(icell)%P(1:3,1)= a1 !polygon points (coordinates)
666 brick_list(nin,ib)%PCUT(icell)%P(1:3,2)= a2 !polygon points (coordinates)
667 brick_list(nin,ib)%PCUT(icell)%P(1:3,3)= a3 !polygon points (coordinates)
668 brick_list(nin,ib)%PCUT(icell)%P(1:3,4)= a4 !polygon points (coordinates)
669 psubvol(icell)%Vnew = vol
670 pnnodcell(icell)%NumNOD = mnod !main(brick nodes)
671 pnpointcell(icell)%NumPOINT= nnod
672 brick_list(nin,ib)%POLY(icell)%CellCenter(1)= sum( (/a1(1),a2(1),a3(1),a4(1),n1(1),n2(1)/) ) / six
673 brick_list(nin,ib)%POLY(icell)%CellCenter(2)= sum( (/a1(2),a2(2),a3(2),a4(2),n1(2),n2(2)/) ) / six
674 brick_list(nin,ib)%POLY(icell)%CellCenter(3)= sum( (/a1(3),a2(3),a3(3),a4(3),n1(3),n2(3)/) ) / six
675 brick_list(nin,ib)%POLY(icell)%ListNodID(1:mnod) = gnode(1:mnod,secid) !IXS(Gnode(1,SECID)+1,ID)
676 inodtag(gnode(1:mnod,secid)) = 1 !tag local brick nodes [1 to 8] if used by polyhedron
677 pwhichcellnod(gnode(1:mnod,secid))%WhichCell = icell
678 !VELOCITY ON POLYHEDRA FACES
679 a1 => edge_list(nin,iad(1))%CUTVEL(1:3,ipose(1))
680 a2 => edge_list(nin,iad(2))%CUTVEL(1:3,ipose(2))
681 a3 => edge_list(nin,iad(3))%CUTVEL(1:3,ipose(3))
682 a4 => edge_list(nin,iad(4))%CUTVEL(1:3,ipose(4))
683 n1 => v(1:3,ixs(gnode(1,secid)+1,id))
684 n2 => v(1:3,ixs(gnode(2,secid)+1,id))
685 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(1),ICELL) = THIRD* (n1(1:3) + a2(1:3) + a1(1:3))
686 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(2),ICELL) = THIRD* (n2(1:3) + a4(1:3) + a3(1:3))
687 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(3),ICELL) = FOURTH* (n1(1:3) + n2(1:3) + a3(1:3) + a2(1:3))
688 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(4),ICELL) = FOURTH* (n2(1:3) + n1(1:3) + a3(1:3) + a2(1:3))
689 !BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE0%Vel(1:3) = FOURTH* (a1(1:3) + a2(1:3) + a3(1:3) + a4(1:3))
690 brick_list(nin,ib)%PCUT(icell)%Vel(1:3) = fourth* (a1(1:3) + a2(1:3) + a3(1:3) + a4(1:3))
691
692 ELSEIF(sectype(1:5)=='POLY3')THEN
693 nface = f_poly3
694 a1 => edge_list(nin,iad(1))%CUTPOINT(1:3,ipose(1))
695 a2 => edge_list(nin,iad(2))%CUTPOINT(1:3,ipose(2))
696 a3 => edge_list(nin,iad(3))%CUTPOINT(1:3,ipose(3))
697 a4 => edge_list(nin,iad(4))%CUTPOINT(1:3,ipose(4))
698 a5 => edge_list(nin,iad(5))%CUTPOINT(1:3,ipose(5))
699 n1 => x(1:3,ixs(gnode(1,secid)+1,id))
700 n2 => x(1:3,ixs(gnode(2,secid)+1,id))
701 n3 => x(1:3,ixs(gnode(3,secid)+1,id))
702 fa(1:5) = gface(1:5,secid)
703 vf(:,1) = i22aera(5,(/n1,n2,n3,a3,a2/),c(1:3,fa(1)))
704 vf(:,2) = i22aera(3,(/n3,a4,a3/), c(1:3,fa(2)))
705 vf(:,3) = i22aera(4,(/n3,n2,a5,a4/), c(1:3,fa(3)))
706 vf(:,4) = i22aera(4,(/n2,n1,a1,a5/), c(1:3,fa(4)))
707 vf(:,5) = i22aera(3,(/n1,a2,a1/), c(1:3,fa(5)))
708 vf(:,0) = i22aera(5,(/a1,a2,a3,a4,a5/),c(1:3, 0 ))
709 brick_list(nin,ib)%POLY(9)%FACE(fa(1))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(1))%Center(1:3) + a3 + a2 !put sum intersection point in icell 9 (will be used in sinit22_fvm.F to compute faceCenter for poly9
710 brick_list(nin,ib)%POLY(9)%FACE(fa(2))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(2))%Center(1:3) + a4 + a3
711 brick_list(nin,ib)%POLY(9)%FACE(fa(3))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(3))%Center(1:3) + a5 + a4
712 brick_list(nin,ib)%POLY(9)%FACE(fa(4))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(4))%Center(1:3) + a1 + a5
713 brick_list(nin,ib)%POLY(9)%FACE(fa(5))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(5))%Center(1:3) + a2 + a1
714 brick_list(nin,ib)%POLY(icell)%FACE(fa(1))%Center(1:3) = c(1:3,fa(1))
715 brick_list(nin,ib)%POLY(icell)%FACE(fa(2))%Center(1:3) = c(1:3,fa(2))
716 brick_list(nin,ib)%POLY(icell)%FACE(fa(3))%Center(1:3) = c(1:3,fa(3))
717 brick_list(nin,ib)%POLY(icell)%FACE(fa(4))%Center(1:3) = c(1:3,fa(4))
718 brick_list(nin,ib)%POLY(icell)%FACE(fa(5))%Center(1:3) = c(1:3,fa(5))
719 brick_list(nin,ib)%POLY(icell)%FACE0%Center(1:3) = c(1:3, 0 )
720 scut = sqrt(sum(vf(:,0)*vf(:,0)))
721 pface(1:6) => brick_list(nin,ib)%POLY(icell)%FACE(1:6)!%Surf
722 pface0=> brick_list(nin,ib)%POLY(icell)%FACE0%Surf
723 pface(fa(1))%Surf = sqrt(sum(vf(:,1)*vf(:,1)))
724 pface(fa(2))%Surf = sqrt(sum(vf(:,2)*vf(:,2)))
725 pface(fa(3))%Surf = sqrt(sum(vf(:,3)*vf(:,3)))
726 pface(fa(4))%Surf = sqrt(sum(vf(:,4)*vf(:,4)))
727 pface(fa(5))%Surf = sqrt(sum(vf(:,5)*vf(:,5)))
728 pface0 = scut
729 brick_list(nin,ib)%PCUT(icell)%SCUT = scut !area of the cutting surface
730 brick_list(nin,ib)%PCUT(icell)%B(1:3) = c(1:3,0) !centroid of the cutting surface (i.e. basis point)
731 brick_list(nin,ib)%PCUT(icell)%N(1:3) = vf(:,0) !Cut surface normal vector (not unitary : 2S factor)
732 brick_list(nin,ib)%PCUT(icell)%NP = 5 !polygon points (number)
733 brick_list(nin,ib)%PCUT(icell)%P(1:3,1)= a1 !polygon points (coordinates)
734 brick_list(nin,ib)%PCUT(icell)%P(1:3,2)= a2 !polygon points (coordinates)
735 brick_list(nin,ib)%PCUT(icell)%P(1:3,3)= a3 !polygon points (coordinates)
736 brick_list(nin,ib)%PCUT(icell)%P(1:3,4)= a4 !polygon points (coordinates)
737 brick_list(nin,ib)%PCUT(icell)%P(1:3,5)= a5 !polygon points (coordinates)
738 vol = sum((/(vf(:,i)*c(1:3,fa(i)),i=0,5)/))
739 vol = third * abs(vol)
740 psubvol(icell)%Vnew = vol
741 pnnodcell(icell)%NumNOD = mnod !main(brick nodes)
742 pnpointcell(icell)%NumPOINT= nnod
743 brick_list(nin,ib)%POLY(icell)%CellCenter(1)= sum( (/a1(1),a2(1),a3(1),a4(1),a5(1),n1(1),n2(1),n3(1)/) ) / eight
744 brick_list(nin,ib)%POLY(icell)%CellCenter(2)= sum( (/a1(2),a2(2),a3(2),a4(2),a5(2),n1(2),n2(2),n3(2)/) ) / eight
745 brick_list(nin,ib)%POLY(icell)%CellCenter(3)= sum( (/a1(3),a2(3),a3(3),a4(3),a5(3),n1(3),n2(3),n3(3)/) ) / eight
746 brick_list(nin,ib)%POLY(icell)%ListNodID(1:mnod) = gnode(1:mnod,secid) !IXS(Gnode(1,SECID)+1,ID)
747 inodtag(gnode(1:mnod,secid)) = 1 !tag local brick nodes [1 to 8] if used by polyhedron
748 pwhichcellnod(gnode(1:mnod,secid))%WhichCell = icell
749 !VELOCITY ON POLYHEDRA FACES
750 a1 => edge_list(nin,iad(1))%CUTVEL(1:3,ipose(1))
751 a2 => edge_list(nin,iad(2))%CUTVEL(1:3,ipose(2))
752 a3 => edge_list(nin,iad(3))%CUTVEL(1:3,ipose(3))
753 a4 => edge_list(nin,iad(4))%CUTVEL(1:3,ipose(4))
754 a5 => edge_list(nin,iad(5))%CUTVEL(1:3,ipose(5))
755 n1 => v(1:3,ixs(gnode(1,secid)+1,id))
756 n2 => v(1:3,ixs(gnode(2,secid)+1,id))
757 n3 => v(1:3,ixs(gnode(3,secid)+1,id))
758 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(1),ICELL) = ONE_FIFTH * (n1+n2+n3+a3+a2)
759 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(2),ICELL) = THIRD * (n3+a4+a3)
760 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(3),ICELL) = FOURTH * (n3+n2+a5+a4)
761 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(4),ICELL) = FOURTH * (n2+n1+a1+a5)
762 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(5),ICELL) = THIRD * (n1+a2+a1)
763 !BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE0%Vel(1:3) = ONE_FIFTH * (a1+a2+a3+a4+a5)
764 brick_list(nin,ib)%PCUT(icell)%Vel(1:3) = one_fifth * (a1+a2+a3+a4+a5)
765
766 ELSEIF(sectype(1:5)=='HEXAE')THEN
767 nface = f_hexae
768 a1 => edge_list(nin,iad(1))%CUTPOINT(1:3,ipose(1))
769 a2 => edge_list(nin,iad(2))%CUTPOINT(1:3,ipose(2))
770 a3 => edge_list(nin,iad(3))%CUTPOINT(1:3,ipose(3))
771 a4 => edge_list(nin,iad(4))%CUTPOINT(1:3,ipose(4))
772 n1 => x(1:3,ixs(gnode(1,secid)+1,id))
773 n2 => x(1:3,ixs(gnode(2,secid)+1,id))
774 n3 => x(1:3,ixs(gnode(3,secid)+1,id))
775 n4 => x(1:3,ixs(gnode(4,secid)+1,id))
776 fa(1:5) = gface(1:5,secid)
777 vf(:,1) = i22aera(4,(/n4,n3,n2,n1/), c(1:3,fa(1)))
778 vf(:,2) = i22aera(4,(/n1,n2,a2,a1/), c(1:3,fa(2)))
779 vf(:,3) = i22aera(4,(/n2,n3,a3,a2/), c(1:3,fa(3)))
780 vf(:,4) = i22aera(4,(/n3,n4,a4,a3/), c(1:3,fa(4)))
781 vf(:,5) = i22aera(4,(/n4,n1,a1,a4/), c(1:3,fa(5)))
782 vf(:,0) = i22aera(4,(/a1,a2,a3,a4/), c(1:3, 0 ))
783 brick_list(nin,ib)%POLY(9)%FACE(fa(1))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(1))%Center(1:3) + zero !put sum intersection point in icell 9 (will be used in sinit22_fvm.F to compute faceCenter for poly9
784 brick_list(nin,ib)%POLY(9)%FACE(fa(2))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(2))%Center(1:3) + a2 + a1
785 brick_list(nin,ib)%POLY(9)%FACE(fa(3))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(3))%Center(1:3) + a3 + a2
786 brick_list(nin,ib)%POLY(9)%FACE(fa(4))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(4))%Center(1:3) + a4 + a3
787 brick_list(nin,ib)%POLY(9)%FACE(fa(5))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(5))%Center(1:3) + a1 + a4
788 brick_list(nin,ib)%POLY(icell)%FACE(fa(1))%Center(1:3) = c(1:3,fa(1))
789 brick_list(nin,ib)%POLY(icell)%FACE(fa(2))%Center(1:3) = c(1:3,fa(2))
790 brick_list(nin,ib)%POLY(icell)%FACE(fa(3))%Center(1:3) = c(1:3,fa(3))
791 brick_list(nin,ib)%POLY(icell)%FACE(fa(4))%Center(1:3) = c(1:3,fa(4))
792 brick_list(nin,ib)%POLY(icell)%FACE(fa(5))%Center(1:3) = c(1:3,fa(5))
793 brick_list(nin,ib)%POLY(icell)%FACE0%Center(1:3) = c(1:3, 0 )
794 scut = sqrt(sum(vf(:,0)*vf(:,0)))
795 pface(1:6) => brick_list(nin,ib)%POLY(icell)%FACE(1:6)!%Surf
796 pface0=> brick_list(nin,ib)%POLY(icell)%FACE0%Surf
797 pface(fa(1))%Surf = sqrt(sum(vf(:,1)*vf(:,1)))
798 pface(fa(2))%Surf = sqrt(sum(vf(:,2)*vf(:,2)))
799 pface(fa(3))%Surf = sqrt(sum(vf(:,3)*vf(:,3)))
800 pface(fa(4))%Surf = sqrt(sum(vf(:,4)*vf(:,4)))
801 pface(fa(5))%Surf = sqrt(sum(vf(:,5)*vf(:,5)))
802 pface0 = scut
803 brick_list(nin,ib)%PCUT(icell)%SCUT = scut !area of the cutting surface
804 brick_list(nin,ib)%PCUT(icell)%B(1:3) = c(1:3,0) !centroid of the cutting surface (i.e. basis point)
805 brick_list(nin,ib)%PCUT(icell)%N(1:3) = vf(:,0) !Cut surface normal vector (not unitary : 2S factor)
806 brick_list(nin,ib)%PCUT(icell)%NP = 4 !polygon points (number)
807 brick_list(nin,ib)%PCUT(icell)%P(1:3,1)= a1 !polygon points (coordinates)
808 brick_list(nin,ib)%PCUT(icell)%P(1:3,2)= a2 !polygon points (coordinates)
809 brick_list(nin,ib)%PCUT(icell)%P(1:3,3)= a3 !polygon points (coordinates)
810 brick_list(nin,ib)%PCUT(icell)%P(1:3,4)= a4 !polygon points (coordinates)
811 vol = sum((/(vf(:,i)*c(1:3,fa(i)),i=0,5)/))
812 vol = third * abs(vol)
813 psubvol(icell)%Vnew = vol
814 pnnodcell(icell)%NumNOD = mnod !main(brick nodes)
815 pnpointcell(icell)%NumPOINT= nnod
816 brick_list(nin,ib)%POLY(icell)%CellCenter(1)= sum( (/a1(1),a2(1),a3(1),a4(1),n1(1),n2(1),n3(1),n4(1)/) ) / eight
817 brick_list(nin,ib)%POLY(icell)%CellCenter(2)= sum( (/a1(2),a2(2),a3(2),a4(2),n1(2),n2(2),n3(2),n4(2)/) ) / eight
818 brick_list(nin,ib)%POLY(icell)%CellCenter(3)= sum( (/a1(3),a2(3),a3(3),a4(3),n1(3),n2(3),n3(3),n4(3)/) ) / eight
819 brick_list(nin,ib)%POLY(icell)%ListNodID(1:mnod) = gnode(1:mnod,secid) !IXS(Gnode(1,SECID)+1,ID)
820 inodtag(gnode(1:mnod,secid)) = 1 !tag local brick nodes [1 to 8] if used by polyhedron
821 pwhichcellnod(gnode(1:mnod,secid))%WhichCell = icell
822 !VELOCITY ON POLYHEDRA FACES
823 a1 => edge_list(nin,iad(1))%CUTVEL(1:3,ipose(1))
824 a2 => edge_list(nin,iad(2))%CUTVEL(1:3,ipose(2))
825 a3 => edge_list(nin,iad(3))%CUTVEL(1:3,ipose(3))
826 a4 => edge_list(nin,iad(4))%CUTVEL(1:3,ipose(4))
827 n1 => v(1:3,ixs(gnode(1,secid)+1,id))
828 n2 => v(1:3,ixs(gnode(2,secid)+1,id))
829 n3 => v(1:3,ixs(gnode(3,secid)+1,id))
830 n4 => v(1:3,ixs(gnode(4,secid)+1,id))
831 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(1),ICELL) = FOURTH * (n4+n3+n2+n1)
832 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(2),ICELL) = FOURTH * (n1+n2+a2+a1)
833 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(3),ICELL) = FOURTH * (n2+n3+a3+a2)
834 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(4),ICELL) = FOURTH * (n3+n4+a4+a3)
835 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(5),ICELL) = FOURTH * (n4+n1+a1+a4)
836 !BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE0%Vel(1:3) = FOURTH * (a1+a2+a3+a4)
837 brick_list(nin,ib)%PCUT(icell)%Vel(1:3) = fourth * (a1+a2+a3+a4)
838
839 ELSEIF(sectype(1:6)=='POLY4 ')THEN
840 nface = f_poly4
841 a1 => edge_list(nin,iad(1))%CUTPOINT(1:3,ipose(1))
842 a2 => edge_list(nin,iad(2))%CUTPOINT(1:3,ipose(2))
843 a3 => edge_list(nin,iad(3))%CUTPOINT(1:3,ipose(3))
844 a4 => edge_list(nin,iad(4))%CUTPOINT(1:3,ipose(4))
845 a5 => edge_list(nin,iad(5))%CUTPOINT(1:3,ipose(5))
846 a6 => edge_list(nin,iad(6))%CUTPOINT(1:3,ipose(6))
847 n1 => x(1:3,ixs(gnode(1,secid)+1,id))
848 n2 => x(1:3,ixs(gnode(2,secid)+1,id))
849 n3 => x(1:3,ixs(gnode(3,secid)+1,id))
850 n4 => x(1:3,ixs(gnode(4,secid)+1,id))
851 fa(1:6) = gface(1:6,secid)
852 vf(:,1) =i22aera(3,(/n1,a2,a1/) , c(1:3,fa(1)))
853 vf(:,2) =i22aera(5,(/n2,n4,a5,a4,n3/), c(1:3,fa(2)))
854 vf(:,3) =i22aera(5,(/a1,a6,n4,n2,n1/), c(1:3,fa(3)))
855 vf(:,4) =i22aera(5,(/n1,n2,n3,a3,a2/), c(1:3,fa(4)))
856 vf(:,5) =i22aera(3,(/n4,a6,a5/), c(1:3,fa(5)))
857 vf(:,6) =i22aera(3,(/n3,a4,a3/), c(1:3,fa(6)))
858 vf(:,0) =i22aera(6,(/a1,a2,a3,a4,a5,a6/),c(1:3, 0 ))
859 vol = sum((/(vf(:,i)*c(1:3,fa(i)),i=0,6)/))
860 vol = third * abs(vol)
861 if(icell==2)then
862 vol1 = brick_list(nin,ib)%poly(1)%vnew
863 vol2 = vol
864 if((vol1+vol2)/uncutvol>=one-em04)then
865 itag(icell)=1
866 nintpoint = nintpoint - ntrus
867 !print *,"cleaning 3"
868 cycle
869 elseif( abs((vol1-vol2)/uncutvol)<=em04 )then
870 scut1 = brick_list(nin,ib)%Pcut(1)%Scut(1)
871 scut2 = brick_list(nin,ib)%Pcut(2)%Scut(1)
872 tmp(1) = c(1, 0 ) - brick_list(nin,ib)%POLY(1)%FACE0%Center(1)
873 tmp(2) = c(2, 0 ) - brick_list(nin,ib)%POLY(1)%FACE0%Center(2)
874 tmp(3) = c(3, 0 ) - brick_list(nin,ib)%POLY(1)%FACE0%Center(3)
875 dist = tmp(1)*tmp(1)+tmp(2)*tmp(2)+tmp(3)*tmp(3)
876 dist = sqrt(dist)
877 if (dist / exp(0.333*log(uncutvol)) <= em04)THEN
878 itag(icell)=1
879 nintpoint = nintpoint - ntrus
880 !print *,"cleaning 4"
881 cycle
882 ENDIF
883 endif
884 endif
885 brick_list(nin,ib)%POLY(9)%FACE(fa(1))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(1))%Center(1:3) + a2 + a1 !put sum intersection point in icell 9 (will be used in sinit22_fvm.F to compute faceCenter for poly9
886 brick_list(nin,ib)%POLY(9)%FACE(fa(2))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(2))%Center(1:3) + a5 + a4
887 brick_list(nin,ib)%POLY(9)%FACE(fa(3))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(3))%Center(1:3) + a1 + a6
888 brick_list(nin,ib)%POLY(9)%FACE(fa(4))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(4))%Center(1:3) + a3 + a2
889 brick_list(nin,ib)%POLY(9)%FACE(fa(5))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(5))%Center(1:3) + a6 + a5
890 brick_list(nin,ib)%POLY(9)%FACE(fa(6))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(6))%Center(1:3) + a4 + a3
891 brick_list(nin,ib)%POLY(icell)%FACE(fa(1))%Center(1:3) = c(1:3,fa(1))
892 brick_list(nin,ib)%POLY(icell)%FACE(fa(2))%Center(1:3) = c(1:3,fa(2))
893 brick_list(nin,ib)%POLY(icell)%FACE(fa(3))%Center(1:3) = c(1:3,fa(3))
894 brick_list(nin,ib)%POLY(icell)%FACE(fa(4))%Center(1:3) = c(1:3,fa(4))
895 brick_list(nin,ib)%POLY(icell)%FACE(fa(5))%Center(1:3) = c(1:3,fa(5))
896 brick_list(nin,ib)%POLY(icell)%FACE(fa(6))%Center(1:3) = c(1:3,fa(6))
897 brick_list(nin,ib)%POLY(icell)%FACE0%Center(1:3) = c(1:3, 0 )
898 scut = sqrt(sum(vf(:,0)*vf(:,0)))
899 pface(1:6) => brick_list(nin,ib)%POLY(icell)%FACE(1:6)!%Surf
900 pface0=> brick_list(nin,ib)%POLY(icell)%FACE0%Surf
901 pface(fa(1))%Surf = sqrt(sum(vf(:,1)*vf(:,1)))
902 pface(fa(2))%Surf = sqrt(sum(vf(:,2)*vf(:,2)))
903 pface(fa(3))%Surf = sqrt(sum(vf(:,3)*vf(:,3)))
904 pface(fa(4))%Surf = sqrt(sum(vf(:,4)*vf(:,4)))
905 pface(fa(5))%Surf = sqrt(sum(vf(:,5)*vf(:,5)))
906 pface(fa(6))%Surf = sqrt(sum(vf(:,6)*vf(:,6)))
907 pface0 = scut
908 brick_list(nin,ib)%PCUT(icell)%SCUT = scut !area of the cutting surface
909 brick_list(nin,ib)%PCUT(icell)%B(1:3) = c(1:3,0) !centroid of the cutting surface (i.e. basis point)
910 brick_list(nin,ib)%PCUT(icell)%N(1:3) = vf(:,0) !Cut surface normal vector (not unitary : 2S factor)
911 brick_list(nin,ib)%PCUT(icell)%NP = 6 !polygon points (number)
912 brick_list(nin,ib)%PCUT(icell)%P(1:3,1)= a1 !polygon points (coordinates)
913 brick_list(nin,ib)%PCUT(icell)%P(1:3,2)= a2 !polygon points (coordinates)
914 brick_list(nin,ib)%PCUT(icell)%P(1:3,3)= a3 !polygon points (coordinates)
915 brick_list(nin,ib)%PCUT(icell)%P(1:3,4)= a4 !polygon points (coordinates)
916 brick_list(nin,ib)%PCUT(icell)%P(1:3,5)= a5 !polygon points (coordinates)
917 brick_list(nin,ib)%PCUT(icell)%P(1:3,6)= a6 !polygon points (coordinates)
918 psubvol(icell)%Vnew = vol
919 pnnodcell(icell)%NumNOD = mnod !main(brick nodes)
920 pnpointcell(icell)%NumPOINT= nnod
921 brick_list(nin,ib)%POLY(icell)%CellCenter(1)=sum((/a1(1),a2(1),a3(1),a4(1),a5(1),a6(1),n1(1),n2(1),n3(1),n4(1)/))/ten
922 brick_list(nin,ib)%POLY(icell)%CellCenter(2)=sum((/a1(2),a2(2),a3(2),a4(2),a5(2),a6(2),n1(2),n2(2),n3(2),n4(2)/))/ten
923 brick_list(nin,ib)%POLY(icell)%CellCenter(3)=sum((/a1(3),a2(3),a3(3),a4(3),a5(3),a6(3),n1(3),n2(3),n3(3),n4(3)/))/ten
924 brick_list(nin,ib)%POLY(icell)%ListNodID(1:mnod) = gnode(1:mnod,secid) !IXS(Gnode(1,SECID)+1,ID)
925 inodtag(gnode(1:mnod,secid)) = 1 !tag local brick nodes [1 to 8] if used by polyhedron
926 pwhichcellnod(gnode(1:mnod,secid))%WhichCell = icell
927 !VELOCITY ON POLYHEDRA FACES
928 a1 => edge_list(nin,iad(1))%CUTVEL(1:3,ipose(1))
929 a2 => edge_list(nin,iad(2))%CUTVEL(1:3,ipose(2))
930 a3 => edge_list(nin,iad(3))%CUTVEL(1:3,ipose(3))
931 a4 => edge_list(nin,iad(4))%CUTVEL(1:3,ipose(4))
932 a5 => edge_list(nin,iad(5))%CUTVEL(1:3,ipose(5))
933 a6 => edge_list(nin,iad(6))%CUTVEL(1:3,ipose(6))
934 n1 => v(1:3,ixs(gnode(1,secid)+1,id))
935 n2 => v(1:3,ixs(gnode(2,secid)+1,id))
936 n3 => v(1:3,ixs(gnode(3,secid)+1,id))
937 n4 => v(1:3,ixs(gnode(4,secid)+1,id))
938 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(1),ICELL) = THIRD * (n1+a2+a1)
939 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(2),ICELL) = ONE_FIFTH * (n2+n4+a5+a4+n3)
940 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(3),ICELL) = ONE_FIFTH * (a1+a6+n4+n2+n1)
941 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(4),ICELL) = ONE_FIFTH * (n1+n2+n3+a3+a2)
942 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(5),ICELL) = THIRD * (n4+a6+a5)
943 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(6),ICELL) = THIRD * (n3+a4+a3)
944 !BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE0%Vel(1:3) = ONE_OVER_6 * (a1+a2+a3+a4+a5+a6 )
945 brick_list(nin,ib)%PCUT(icell)%Vel(1:3) = one_over_6 * (a1+a2+a3+a4+a5+a6 )
946
947 ELSEIF(sectype(1:6)=='POLY4A')THEN
948 nface = f_poly4a
949 a1 => edge_list(nin,iad(1))%CUTPOINT(1:3,ipose(1))
950 a2 => edge_list(nin,iad(2))%CUTPOINT(1:3,ipose(2))
951 a3 => edge_list(nin,iad(3))%CUTPOINT(1:3,ipose(3))
952 a4 => edge_list(nin,iad(4))%CUTPOINT(1:3,ipose(4))
953 a5 => edge_list(nin,iad(5))%CUTPOINT(1:3,ipose(5))
954 a6 => edge_list(nin,iad(6))%CUTPOINT(1:3,ipose(6))
955 n1 => x(1:3,ixs(gnode(1,secid)+1,id))
956 n2 => x(1:3,ixs(gnode(2,secid)+1,id))
957 n3 => x(1:3,ixs(gnode(3,secid)+1,id))
958 n4 => x(1:3,ixs(gnode(4,secid)+1,id))
959 fa(1:6) = gface(1:6,secid)
960 vf(:,1) =i22aera(3,(/n1,a1,a6/) , c(1:3,fa(1)))
961 vf(:,2) =i22aera(4,(/n1,n2,a2,a1/) , c(1:3,fa(2)))
962 vf(:,3) =i22aera(5,(/n2,n3,n4,a3,a2/), c(1:3,fa(3)))
963 vf(:,4) =i22aera(3,(/n4,a4,a3/) , c(1:3,fa(4)))
964 vf(:,5) =i22aera(4,(/n4,n3,a5,a4/) , c(1:3,fa(5)))
965 vf(:,6) =i22aera(5,(/n3,n2,n1,a6,a5/), c(1:3,fa(6)))
966 vf(:,0) =i22aera(6,(/a1,a2,a3,a4,a5,a6/),c(1:3, 0 ))
967 brick_list(nin,ib)%POLY(9)%FACE(fa(1))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(1))%Center(1:3) + a1 + a6 !put sum intersection point in icell 9 (will be used in sinit22_fvm.F to compute faceCenter for poly9
968 brick_list(nin,ib)%POLY(9)%FACE(fa(2))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(2))%Center(1:3) + a2 + a1
969 brick_list(nin,ib)%POLY(9)%FACE(fa(3))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(3))%Center(1:3) + a3 + a2
970 brick_list(nin,ib)%POLY(9)%FACE(fa(4))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(4))%Center(1:3) + a4 + a3
971 brick_list(nin,ib)%POLY(9)%FACE(fa(5))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(5))%Center(1:3) + a5 + a4
972 brick_list(nin,ib)%POLY(9)%FACE(fa(6))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(6))%Center(1:3) + a6 + a5
973 brick_list(nin,ib)%POLY(icell)%FACE(fa(1))%Center(1:3) = c(1:3,fa(1))
974 brick_list(nin,ib)%POLY(icell)%FACE(fa(2))%Center(1:3) = c(1:3,fa(2))
975 brick_list(nin,ib)%POLY(icell)%FACE(fa(3))%Center(1:3) = c(1:3,fa(3))
976 brick_list(nin,ib)%POLY(icell)%FACE(fa(4))%Center(1:3) = c(1:3,fa(4))
977 brick_list(nin,ib)%POLY(icell)%FACE(fa(5))%Center(1:3) = c(1:3,fa(5))
978 brick_list(nin,ib)%POLY(icell)%FACE(fa(6))%Center(1:3) = c(1:3,fa(6))
979 brick_list(nin,ib)%POLY(icell)%FACE0%Center(1:3) = c(1:3, 0 )
980 scut = sqrt(sum(vf(:,0)*vf(:,0)))
981 pface(1:6) => brick_list(nin,ib)%POLY(icell)%FACE(1:6)!%Surf
982 pface0=> brick_list(nin,ib)%POLY(icell)%FACE0%Surf
983 pface(fa(1))%Surf = sqrt(sum(vf(:,1)*vf(:,1)))
984 pface(fa(2))%Surf = sqrt(sum(vf(:,2)*vf(:,2)))
985 pface(fa(3))%Surf = sqrt(sum(vf(:,3)*vf(:,3)))
986 pface(fa(4))%Surf = sqrt(sum(vf(:,4)*vf(:,4)))
987 pface(fa(5))%Surf = sqrt(sum(vf(:,5)*vf(:,5)))
988 pface(fa(6))%Surf = sqrt(sum(vf(:,6)*vf(:,6)))
989 pface0 = scut
990 brick_list(nin,ib)%PCUT(icell)%SCUT = scut !area of the cutting surface
991 brick_list(nin,ib)%PCUT(icell)%B(1:3) = c(1:3,0) !centroid of the cutting surface (i.e. basis point)
992 brick_list(nin,ib)%PCUT(icell)%N(1:3) = vf(:,0) !Cut surface normal vector (not unitary : 2S factor)
993 brick_list(nin,ib)%PCUT(icell)%NP = 6 !polygon points (number)
994 brick_list(nin,ib)%PCUT(icell)%P(1:3,1)= a1 !polygon points (coordinates)
995 brick_list(nin,ib)%PCUT(icell)%P(1:3,2)= a2 !polygon points (coordinates)
996 brick_list(nin,ib)%PCUT(icell)%P(1:3,3)= a3 !polygon points (coordinates)
997 brick_list(nin,ib)%PCUT(icell)%P(1:3,4)= a4 !polygon points (coordinates)
998 brick_list(nin,ib)%PCUT(icell)%P(1:3,5)= a5 !polygon points (coordinates)
999 brick_list(nin,ib)%PCUT(icell)%P(1:3,6)= a6 !polygon points (coordinates)
1000 vol = sum((/(vf(:,i)*c(1:3,fa(i)),i=0,6)/))
1001 vol = third * abs(vol)
1002 psubvol(icell)%Vnew = vol
1003 pnnodcell(icell)%NumNOD = mnod !main(brick nodes)
1004 pnpointcell(icell)%NumPOINT= nnod
1005 brick_list(nin,ib)%POLY(icell)%CellCenter(1)= sum((/a1(1),a2(1),a3(1),a4(1),a5(1),a6(1),n1(1),n2(1),n3(1),n4(1)/))/ten
1006 brick_list(nin,ib)%POLY(icell)%CellCenter(2)= sum((/a1(2),a2(2),a3(2),a4(2),a5(2),a6(2),n1(2),n2(2),n3(2),n4(2)/))/ten
1007 brick_list(nin,ib)%POLY(icell)%CellCenter(3)= sum((/a1(3),a2(3),a3(3),a4(3),a5(3),a6(3),n1(3),n2(3),n3(3),n4(3)/))/ten
1008 brick_list(nin,ib)%POLY(icell)%ListNodID(1:mnod) = gnode(1:mnod,secid) !IXS(Gnode(1,SECID)+1,ID)
1009 inodtag(gnode(1:mnod,secid)) = 1 !tag local brick nodes [1 to 8] if used by polyhedron
1010 pwhichcellnod(gnode(1:mnod,secid))%WhichCell = icell
1011 !VELOCITY ON POLYHEDRA FACES
1012 a1 => edge_list(nin,iad(1))%CUTVEL(1:3,ipose(1))
1013 a2 => edge_list(nin,iad(2))%CUTVEL(1:3,ipose(2))
1014 a3 => edge_list(nin,iad(3))%CUTVEL(1:3,ipose(3))
1015 a4 => edge_list(nin,iad(4))%CUTVEL(1:3,ipose(4))
1016 a5 => edge_list(nin,iad(5))%CUTVEL(1:3,ipose(5))
1017 a6 => edge_list(nin,iad(6))%CUTVEL(1:3,ipose(6))
1018 n1 => v(1:3,ixs(gnode(1,secid)+1,id))
1019 n2 => v(1:3,ixs(gnode(2,secid)+1,id))
1020 n3 => v(1:3,ixs(gnode(3,secid)+1,id))
1021 n4 => v(1:3,ixs(gnode(4,secid)+1,id))
1022 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(1),ICELL) = THIRD * (n1+a1+a6)
1023 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(2),ICELL) = FOURTH * (n1+n2+a2+a1)
1024 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(3),ICELL) = ONE_FIFTH * (n2+n3+n4+a3+a2)
1025 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(4),ICELL) = THIRD * (n4+a4+a3)
1026 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(5),ICELL) = FOURTH * (n4+n3+a5+a4)
1027 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(6),ICELL) = ONE_FIFTH * (n3+n2+n1+a6+a5)
1028 !BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE0%Vel(1:3) = ONE_OVER_6 * (a1+a2+a3+a4+a5+a6 )
1029 brick_list(nin,ib)%PCUT(icell)%Vel(1:3) = one_over_6 * (a1+a2+a3+a4+a5+a6 )
1030
1031 ELSEIF(sectype(1:6)=='POLY4B')THEN
1032 nface = f_poly4b
1033 a1 => edge_list(nin,iad(1))%CUTPOINT(1:3,ipose(1))
1034 a2 => edge_list(nin,iad(2))%CUTPOINT(1:3,ipose(2))
1035 a3 => edge_list(nin,iad(3))%CUTPOINT(1:3,ipose(3))
1036 a4 => edge_list(nin,iad(4))%CUTPOINT(1:3,ipose(4))
1037 a5 => edge_list(nin,iad(5))%CUTPOINT(1:3,ipose(5))
1038 a6 => edge_list(nin,iad(6))%CUTPOINT(1:3,ipose(6))
1039 n1 => x(1:3,ixs(gnode(1,secid)+1,id))
1040 n2 => x(1:3,ixs(gnode(2,secid)+1,id))
1041 n3 => x(1:3,ixs(gnode(3,secid)+1,id))
1042 n4 => x(1:3,ixs(gnode(4,secid)+1,id))
1043 fa(1:6) = gface(1:6,secid)
1044 vf(:,1) = - i22aera(3,(/n1,a1,a6/) , c(1:3,fa(1)))
1045 vf(:,2) = - i22aera(4,(/n1,n2,a2,a1/) , c(1:3,fa(2)))
1046 vf(:,3) = - i22aera(5,(/n2,n3,n4,a3,a2/), c(1:3,fa(3)))
1047 vf(:,4) = - i22aera(3,(/n4,a4,a3/) , c(1:3,fa(4)))
1048 vf(:,5) = - i22aera(4,(/n4,n3,a5,a4/) , c(1:3,fa(5)))
1049 vf(:,6) = - i22aera(5,(/n3,n2,n1,a6,a5/), c(1:3,fa(6)))
1050 vf(:,0) = - i22aera(6,(/a1,a2,a3,a4,a5,a6/),c(1:3, 0 ))
1051 brick_list(nin,ib)%POLY(9)%FACE(fa(1))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(1))%Center(1:3) + a1 + a6 !put sum intersection point in icell 9 (will be used in sinit22_fvm.F to compute faceCenter for poly9
1052 brick_list(nin,ib)%POLY(9)%FACE(fa(2))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(2))%Center(1:3) + a2 + a1
1053 brick_list(nin,ib)%POLY(9)%FACE(fa(3))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(3))%Center(1:3) + a3 + a2
1054 brick_list(nin,ib)%POLY(9)%FACE(fa(4))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(4))%Center(1:3) + a4 + a3
1055 brick_list(nin,ib)%POLY(9)%FACE(fa(5))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(5))%Center(1:3) + a5 + a4
1056 brick_list(nin,ib)%POLY(9)%FACE(fa(6))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(6))%Center(1:3) + a6 + a5
1057 brick_list(nin,ib)%POLY(icell)%FACE(fa(1))%Center(1:3) = c(1:3,fa(1))
1058 brick_list(nin,ib)%POLY(icell)%FACE(fa(2))%Center(1:3) = c(1:3,fa(2))
1059 brick_list(nin,ib)%POLY(icell)%FACE(fa(3))%Center(1:3) = c(1:3,fa(3))
1060 brick_list(nin,ib)%POLY(icell)%FACE(fa(4))%Center(1:3) = c(1:3,fa(4))
1061 brick_list(nin,ib)%POLY(icell)%FACE(fa(5))%Center(1:3) = c(1:3,fa(5))
1062 brick_list(nin,ib)%POLY(icell)%FACE(fa(6))%Center(1:3) = c(1:3,fa(6))
1063 brick_list(nin,ib)%POLY(icell)%FACE0%Center(1:3) = c(1:3, 0 )
1064 scut = sqrt(sum(vf(:,0)*vf(:,0)))
1065 pface(1:6) => brick_list(nin,ib)%POLY(icell)%FACE(1:6)!%Surf
1066 pface0=> brick_list(nin,ib)%POLY(icell)%FACE0%Surf
1067 pface(fa(1))%Surf = sqrt(sum(vf(:,1)*vf(:,1)))
1068 pface(fa(2))%Surf = sqrt(sum(vf(:,2)*vf(:,2)))
1069 pface(fa(3))%Surf = sqrt(sum(vf(:,3)*vf(:,3)))
1070 pface(fa(4))%Surf = sqrt(sum(vf(:,4)*vf(:,4)))
1071 pface(fa(5))%Surf = sqrt(sum(vf(:,5)*vf(:,5)))
1072 pface(fa(6))%Surf = sqrt(sum(vf(:,6)*vf(:,6)))
1073 pface0 = scut
1074 brick_list(nin,ib)%PCUT(icell)%SCUT = scut !area of the cutting surface
1075 brick_list(nin,ib)%PCUT(icell)%B(1:3) = c(1:3,0) !centroid of the cutting surface (i.e. basis point)
1076 brick_list(nin,ib)%PCUT(icell)%N(1:3) = vf(:,0) !Cut surface normal vector (not unitary : 2S factor)
1077 brick_list(nin,ib)%PCUT(icell)%NP = 6 !polygon points (number)
1078 brick_list(nin,ib)%PCUT(icell)%P(1:3,1)= a1 !polygon points (coordinates)
1079 brick_list(nin,ib)%PCUT(icell)%P(1:3,2)= a2 !polygon points (coordinates)
1080 brick_list(nin,ib)%PCUT(icell)%P(1:3,3)= a3 !polygon points (coordinates)
1081 brick_list(nin,ib)%PCUT(icell)%P(1:3,4)= a4 !polygon points (coordinates)
1082 brick_list(nin,ib)%PCUT(icell)%P(1:3,5)= a5 !polygon points (coordinates)
1083 brick_list(nin,ib)%PCUT(icell)%P(1:3,6)= a6 !polygon points (coordinates)
1084 vol = sum((/(vf(:,i)*c(1:3,fa(i)),i=0,6)/))
1085 vol = third * abs(vol)
1086 psubvol(icell)%Vnew = vol
1087 pnnodcell(icell)%NumNOD = mnod !main(brick nodes)
1088 pnpointcell(icell)%NumPOINT= nnod
1089 brick_list(nin,ib)%POLY(icell)%CellCenter(1)=sum((/a1(1),a2(1),a3(1),a4(1),a5(1),a6(1),n1(1),n2(1),n3(1),n4(1)/))/ten
1090 brick_list(nin,ib)%POLY(icell)%CellCenter(2)=sum((/a1(2),a2(2),a3(2),a4(2),a5(2),a6(2),n1(2),n2(2),n3(2),n4(2)/))/ten
1091 brick_list(nin,ib)%POLY(icell)%CellCenter(3)=sum((/a1(3),a2(3),a3(3),a4(3),a5(3),a6(3),n1(3),n2(3),n3(3),n4(3)/))/ten
1092 brick_list(nin,ib)%POLY(icell)%ListNodID(1:mnod) = gnode(1:mnod,secid) !IXS(Gnode(1,SECID)+1,ID)
1093 inodtag(gnode(1:mnod,secid)) = 1 !tag local brick nodes [1 to 8] if used by polyhedron
1094 pwhichcellnod(gnode(1:mnod,secid))%WhichCell = icell
1095 !VELOCITY ON POLYHEDRA FACES
1096 a1 => edge_list(nin,iad(1))%CUTVEL(1:3,ipose(1))
1097 a2 => edge_list(nin,iad(2))%CUTVEL(1:3,ipose(2))
1098 a3 => edge_list(nin,iad(3))%CUTVEL(1:3,ipose(3))
1099 a4 => edge_list(nin,iad(4))%CUTVEL(1:3,ipose(4))
1100 a5 => edge_list(nin,iad(5))%CUTVEL(1:3,ipose(5))
1101 a6 => edge_list(nin,iad(6))%CUTVEL(1:3,ipose(6))
1102 n1 => v(1:3,ixs(gnode(1,secid)+1,id))
1103 n2 => v(1:3,ixs(gnode(2,secid)+1,id))
1104 n3 => v(1:3,ixs(gnode(3,secid)+1,id))
1105 n4 => v(1:3,ixs(gnode(4,secid)+1,id))
1106 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(1),ICELL) = THIRD * (n1+a1+a6)
1107 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(2),ICELL) = FOURTH * (n1+n2+a2+a1)
1108 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(3),ICELL) = ONE_FIFTH * (n2+n3+n4+a3+a2)
1109 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(4),ICELL) = THIRD * (n4+a4+a3)
1110 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(5),ICELL) = FOURTH * (n4+n3+a5+a4)
1111 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(6),ICELL) = ONE_FIFTH * (n3+n2+n1+a6+a5)
1112 !BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE0%Vel(1:3) = ONE_OVER_6 * (a1+a2+a3+a4+a5+a6 )
1113 brick_list(nin,ib)%PCUT(icell)%Vel(1:3) = one_over_6 * (a1+a2+a3+a4+a5+a6 )
1114
1115 ELSEIF(sectype(1:5)=='POLYC')THEN
1116 nface = f_polyc
1117 a1 => edge_list(nin,iad(1))%CUTPOINT(1:3,ipose(1))
1118 a2 => edge_list(nin,iad(2))%CUTPOINT(1:3,ipose(2))
1119 a3 => edge_list(nin,iad(3))%CUTPOINT(1:3,ipose(3))
1120 a4 => edge_list(nin,iad(4))%CUTPOINT(1:3,ipose(4))
1121 a5 => edge_list(nin,iad(5))%CUTPOINT(1:3,ipose(5)) !points 5,6 on same edge : IAD(5)=IAD(6) BUT ICUT=1
1122 a6 => edge_list(nin,iad(6))%CUTPOINT(1:3,ipose(6)) !points 5,6 on same edge : IAD(5)=IAD(6) BUT ICUT=2
1123 n1 => x(1:3,ixs(gnode(1,secid)+1,id))
1124 n2 => x(1:3,ixs(gnode(2,secid)+1,id))
1125 n3 => x(1:3,ixs(gnode(3,secid)+1,id))
1126 n4 => x(1:3,ixs(gnode(4,secid)+1,id))
1127 fa(1:6) = gface(1:6,secid)
1128 vf(:,1) =i22aera(4,(/n1,n2,n3,n4/) , c(1:3,fa(1)))
1129 vf(:,2) =i22aera(4,(/n1,a1,a2,n2/) , c(1:3,fa(2)))
1130 vf(:,3) =i22aera(4,(/n2,a2,a3,n3/) , c(1:3,fa(3)))
1131 vf(:,4) =i22aera(4,(/n1,n4,a4,a1/) , c(1:3,fa(4)))
1132 vf(:,5) =i22aera(7,(/n4,n3,a3,a6,a5,a4,n4/) , c(1:3,fa(5)))
1133 vf(:,0) =i22aera(6,(/a1,a4,a5,a6,a3,a2/) , c(1:3, 0 )) !TO BE UPDATED LATER IF NEEDED : this is an approximation
1134 brick_list(nin,ib)%POLY(9)%FACE(fa(1))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(1))%Center(1:3) + zero !put sum intersection point in icell 9 (will be used in sinit22_fvm.F to compute faceCenter for poly9
1135 brick_list(nin,ib)%POLY(9)%FACE(fa(2))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(2))%Center(1:3) +a1+a2
1136 brick_list(nin,ib)%POLY(9)%FACE(fa(3))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(3))%Center(1:3) +a2+a3
1137 brick_list(nin,ib)%POLY(9)%FACE(fa(4))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(4))%Center(1:3) +a4+a1
1138 brick_list(nin,ib)%POLY(9)%FACE(fa(5))%Center(1:3) = brick_list(nin,ib)%POLY(9)%FACE(fa(5))%Center(1:3) +a3+a6+a5+a4
1139 brick_list(nin,ib)%POLY(icell)%FACE(fa(1))%Center(1:3) = c(1:3,fa(1))
1140 brick_list(nin,ib)%POLY(icell)%FACE(fa(2))%Center(1:3) = c(1:3,fa(2))
1141 brick_list(nin,ib)%POLY(icell)%FACE(fa(3))%Center(1:3) = c(1:3,fa(3))
1142 brick_list(nin,ib)%POLY(icell)%FACE(fa(4))%Center(1:3) = c(1:3,fa(4))
1143 brick_list(nin,ib)%POLY(icell)%FACE(fa(5))%Center(1:3) = c(1:3,fa(5))
1144 brick_list(nin,ib)%POLY(icell)%FACE0%Center(1:3) = c(1:3, 0 )
1145 scut = sqrt(sum(vf(:,0)*vf(:,0)))
1146 pface(1:6) => brick_list(nin,ib)%POLY(icell)%FACE(1:6)!%Surf
1147 pface0=> brick_list(nin,ib)%POLY(icell)%FACE0%Surf
1148 pface(fa(1))%Surf = sqrt(sum(vf(:,1)*vf(:,1)))
1149 pface(fa(2))%Surf = sqrt(sum(vf(:,2)*vf(:,2)))
1150 pface(fa(3))%Surf = sqrt(sum(vf(:,3)*vf(:,3)))
1151 pface(fa(4))%Surf = sqrt(sum(vf(:,4)*vf(:,4)))
1152 pface(fa(5))%Surf = sqrt(sum(vf(:,5)*vf(:,5)))
1153 pface0 = scut
1154 brick_list(nin,ib)%PCUT(icell)%SCUT = scut !area of the cutting surface
1155 brick_list(nin,ib)%PCUT(icell)%B(1:3) = c(1:3,0) !centroid of the cutting surface (i.e. basis point)
1156 brick_list(nin,ib)%PCUT(icell)%N(1:3) = vf(:,0) !Cut surface normal vector (not unitary : 2S factor)
1157 brick_list(nin,ib)%PCUT(icell)%NP = 6 !polygon points (number)
1158 brick_list(nin,ib)%PCUT(icell)%P(1:3,1)= a1 !polygon points (coordinates)
1159 brick_list(nin,ib)%PCUT(icell)%P(1:3,2)= a2 !polygon points (coordinates)
1160 brick_list(nin,ib)%PCUT(icell)%P(1:3,3)= a3 !polygon points (coordinates)
1161 brick_list(nin,ib)%PCUT(icell)%P(1:3,4)= a4 !polygon points (coordinates)
1162 brick_list(nin,ib)%PCUT(icell)%P(1:3,5)= a5 !polygon points (coordinates)
1163 brick_list(nin,ib)%PCUT(icell)%P(1:3,6)= a6 !polygon points (coordinates)
1164 vol = sum((/(vf(:,i)*c(1:3,fa(i)),i=0,5)/))
1165 vol = third * abs(vol)
1166 psubvol(icell)%Vnew = vol
1167 pnnodcell(icell)%NumNOD = mnod !main(brick nodes)
1168 pnpointcell(icell)%NumPOINT= nnod
1169 brick_list(nin,ib)%POLY(icell)%CellCenter(1)=sum((/a1(1),a2(1),a3(1),a4(1),a5(1),a6(1),n1(1),n2(1),n3(1),n4(1)/))/ten
1170 brick_list(nin,ib)%POLY(icell)%CellCenter(2)=sum((/a1(2),a2(2),a3(2),a4(2),a5(2),a6(2),n1(2),n2(2),n3(2),n4(2)/))/ten
1171 brick_list(nin,ib)%POLY(icell)%CellCenter(3)=sum((/a1(3),a2(3),a3(3),a4(3),a5(3),a6(3),n1(3),n2(3),n3(3),n4(3)/))/ten
1172 brick_list(nin,ib)%POLY(icell)%ListNodID(1:mnod) = gnode(1:mnod,secid) !IXS(Gnode(1,SECID)+1,ID)
1173 inodtag(gnode(1:mnod,secid)) = 1 !tag local brick nodes [1 to 8] if used by polyhedron
1174 pwhichcellnod(gnode(1:mnod,secid))%WhichCell = icell
1175 !VELOCITY ON POLYHEDRA FACES
1176 a1 => edge_list(nin,iad(1))%CUTVEL(1:3,ipose(1))
1177 a2 => edge_list(nin,iad(2))%CUTVEL(1:3,ipose(2))
1178 a3 => edge_list(nin,iad(3))%CUTVEL(1:3,ipose(3))
1179 a4 => edge_list(nin,iad(4))%CUTVEL(1:3,ipose(4))
1180 a5 => edge_list(nin,iad(5))%CUTVEL(1:3,ipose(5))
1181 a6 => edge_list(nin,iad(6))%CUTVEL(1:3,ipose(6))
1182 n1 => v(1:3,ixs(gnode(1,secid)+1,id))
1183 n2 => v(1:3,ixs(gnode(2,secid)+1,id))
1184 n3 => v(1:3,ixs(gnode(3,secid)+1,id))
1185 n4 => v(1:3,ixs(gnode(4,secid)+1,id))
1186 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(1),ICELL) = FOURTH * (n1+n2+n3+n4)
1187 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(2),ICELL) = FOURTH * (n1+a1+a2+n2)
1188 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(3),ICELL) = FOURTH * (n2+a2+a3+n3)
1189 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(4),ICELL) = FOURTH * (n1+n4+a4+a1)
1190 !!BRICK_LIST(NIN,IB)%Vel_Face(1:3,Fa(5),ICELL) = (n4+n3+a3+a6+a5+a4+n4) / SEVEN
1191 !BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE0%Vel(1:3) = ONE_OVER_6 * (a1+a4+a5+a6+a3+a2)
1192 brick_list(nin,ib)%PCUT(icell)%Vel(1:3) = one_over_6 * (a1+a4+a5+a6+a3+a2)
1193 END IF
1194
1195 !Treating complementary polyhedron (example : icode=2842;idble=2312)
1196 IF(bcompl)THEN
1197 pface(1:6) => brick_list(nin,ib)%POLY(icell)%FACE(1:6)!%Surf
1198 pface0=> brick_list(nin,ib)%POLY(icell)%FACE0%Surf
1199 DO j=1,nface
1200 pface(fa(j))%Surf =- pface(fa(j))%Surf
1201 ENDDO
1202 psubvol(icell)%Vnew =- psubvol(icell)%Vnew
1203 pnnodcell(icell)%NumNOD = 8-mnod
1204 ENDIF
1205
1206 !------------------------------------------------------------!
1207 ! GRID VELOCITIES ON POLYEDRA FACES !
1208 !------------------------------------------------------------!
1209 brick_list(nin,ib)%POLY(icell)%FACE0%W(1:3) = zero
1210 brick_list(nin,ib)%POLY(icell)%FACE(fa(1))%W(1:3) = zero
1211 IF(i22_aleul == 1)THEN
1212 IF(sectype(1:5)=='TETRA')THEN
1213 nface = f_tetra
1214 fa(1:3) = gface(1:3,secid)
1215 !intersection points
1216 a1 => edge_list(nin,iad(1))%CUTPOINT(1:3,ipose(1))
1217 a2 => edge_list(nin,iad(2))%CUTPOINT(1:3,ipose(2))
1218 a3 => edge_list(nin,iad(3))%CUTPOINT(1:3,ipose(3))
1219 !interpolate w_a1
1220 beta = edge_list(nin,iad(1))%CUTCOOR(ipose(1));
1221 j = iad(1)-(ib-1)*12;
1222 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1223 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1224 wa1(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1225 wa1(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1226 wa1(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1227 !interpolate w_a2
1228 beta = edge_list(nin,iad(2))%CUTCOOR(ipose(2));
1229 j = iad(2)-(ib-1)*12;
1230 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1231 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1232 wa2(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1233 wa2(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1234 wa2(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1235 !interpolate w_a3
1236 beta = edge_list(nin,iad(3))%CUTCOOR(ipose(3));
1237 j = iad(3)-(ib-1)*12;
1238 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1239 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1240 wa3(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1241 wa3(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1242 wa3(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1243 !brick nodes
1244 n1 => x(1:3,ixs(gnode(1,secid)+1,id))
1245 wn1(1)= w(1,ixs(gnode(1,secid)+1,id))
1246 wn1(2)= w(2,ixs(gnode(1,secid)+1,id))
1247 wn1(3)= w(3,ixs(gnode(1,secid)+1,id))
1248 !Grid velocity on face1 :
1249 wf1(1) = third*(wn1(1)+wa2(1)+wa1(1))
1250 wf1(2) = third*(wn1(2)+wa2(2)+wa1(2))
1251 wf1(3) = third*(wn1(3)+wa2(3)+wa1(3))
1252 !Grid velocity on face2 :
1253 wf2(1) = third*(wn1(1)+wa3(1)+wa2(1))
1254 wf2(2) = third*(wn1(2)+wa3(2)+wa2(2))
1255 wf2(3) = third*(wn1(3)+wa3(3)+wa2(3))
1256 !Grid velocity on face3 :
1257 wf3(1) = third*(wn1(1)+wa1(1)+wa3(1))
1258 wf3(2) = third*(wn1(2)+wa1(2)+wa3(2))
1259 wf3(3) = third*(wn1(3)+wa1(3)+wa3(3))
1260 !Grid velocity on face0 :
1261 wf0(1) = third*(wa1(1)+wa2(1)+wa3(1))
1262 wf0(2) = third*(wa1(2)+wa2(2)+wa3(2))
1263 wf0(3) = third*(wa1(3)+wa2(3)+wa3(3))
1264 !storage
1265 brick_list(nin,ib)%POLY(icell)%FACE0%W(1:3) = wf0(1:3)
1266 brick_list(nin,ib)%POLY(icell)%FACE(fa(1))%W(1:3) = wf1(1:3)
1267 brick_list(nin,ib)%POLY(icell)%FACE(fa(2))%W(1:3) = wf2(1:3)
1268 brick_list(nin,ib)%POLY(icell)%FACE(fa(3))%W(1:3) = wf3(1:3)
1269
1270
1271 ELSEIF(sectype(1:5)=='PENTA')THEN
1272 nface = f_penta
1273 fa(1:4) = gface(1:4,secid)
1274 a1 => edge_list(nin,iad(1))%CUTPOINT(1:3,ipose(1))
1275 a2 => edge_list(nin,iad(2))%CUTPOINT(1:3,ipose(2))
1276 a3 => edge_list(nin,iad(3))%CUTPOINT(1:3,ipose(3))
1277 a4 => edge_list(nin,iad(4))%CUTPOINT(1:3,ipose(4))
1278 n1 => x(1:3,ixs(gnode(1,secid)+1,id))
1279 n2 => x(1:3,ixs(gnode(2,secid)+1,id))
1280 !interpolate w_a1
1281 beta = edge_list(nin,iad(1))%CUTCOOR(ipose(1));
1282 j = iad(1)-(ib-1)*12;
1283 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1284 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1285 wa1(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1286 wa1(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1287 wa1(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1288 !interpolate w_a2
1289 beta = edge_list(nin,iad(2))%CUTCOOR(ipose(2));
1290 j = iad(2)-(ib-1)*12;
1291 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1292 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1293 wa2(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1294 wa2(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1295 wa2(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1296 !interpolate w_a3
1297 beta = edge_list(nin,iad(3))%CUTCOOR(ipose(3));
1298 j = iad(3)-(ib-1)*12;
1299 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1300 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1301 wa3(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1302 wa3(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1303 wa3(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1304 !interpolate w_a4
1305 beta = edge_list(nin,iad(4))%CUTCOOR(ipose(4));
1306 j = iad(4)-(ib-1)*12;
1307 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1308 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1309 wa4(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1310 wa4(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1311 wa4(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1312 !brick nodes
1313 wn1(1)= w(1,ixs(gnode(1,secid)+1,id))
1314 wn1(2)= w(2,ixs(gnode(1,secid)+1,id))
1315 wn1(3)= w(3,ixs(gnode(1,secid)+1,id))
1316 !brick nodes
1317 wn2(1)= w(1,ixs(gnode(2,secid)+1,id))
1318 wn2(2)= w(2,ixs(gnode(2,secid)+1,id))
1319 wn2(3)= w(3,ixs(gnode(2,secid)+1,id))
1320 !Grid velocity on face1 :
1321 wf1(1) = third*(wn1(1)+wa2(1)+wa1(1))
1322 wf1(2) = third*(wn1(2)+wa2(2)+wa1(2))
1323 wf1(3) = third*(wn1(3)+wa2(3)+wa1(3))
1324 !Grid velocity on face2 :
1325 wf2(1) = third*(wn2(1)+wa4(1)+wa3(1))
1326 wf2(2) = third*(wn2(2)+wa4(2)+wa3(2))
1327 wf2(3) = third*(wn2(3)+wa4(3)+wa3(3))
1328 !Grid velocity on face3 :
1329 wf3(1) = fourth*(wn1(1)+wn2(1)+wa3(1)+wa2(1))
1330 wf3(2) = fourth*(wn1(2)+wn2(2)+wa3(2)+wa2(2))
1331 wf3(3) = fourth*(wn1(3)+wn2(3)+wa3(3)+wa2(3))
1332 !Grid velocity on face4 :
1333 wf4(1) = fourth*(wn2(1)+wn1(1)+wa1(1)+wa4(1))
1334 wf4(2) = fourth*(wn2(2)+wn1(2)+wa1(2)+wa4(2))
1335 wf4(3) = fourth*(wn2(3)+wn1(3)+wa1(3)+wa4(3))
1336 !Grid velocity on face0 :
1337 wf0(1) = fourth*(wa1(1)+wa2(1)+wa3(1)+wa4(1))
1338 wf0(2) = fourth*(wa1(2)+wa2(2)+wa3(2)+wa4(2))
1339 wf0(3) = fourth*(wa1(3)+wa2(3)+wa3(3)+wa4(3))
1340 !storage
1341 brick_list(nin,ib)%POLY(icell)%FACE0%W(1:3) = wf0(1:3)
1342 brick_list(nin,ib)%POLY(icell)%FACE(fa(1))%W(1:3) = wf1(1:3)
1343 brick_list(nin,ib)%POLY(icell)%FACE(fa(2))%W(1:3) = wf2(1:3)
1344 brick_list(nin,ib)%POLY(icell)%FACE(fa(3))%W(1:3) = wf3(1:3)
1345 brick_list(nin,ib)%POLY(icell)%FACE(fa(4))%W(1:3) = wf4(1:3)
1346
1347 ELSEIF(sectype(1:5)=='POLY3')THEN
1348 nface = f_poly3
1349 fa(1:5) = gface(1:5,secid)
1350 a1 => edge_list(nin,iad(1))%CUTPOINT(1:3,ipose(1))
1351 a2 => edge_list(nin,iad(2))%CUTPOINT(1:3,ipose(2))
1352 a3 => edge_list(nin,iad(3))%CUTPOINT(1:3,ipose(3))
1353 a4 => edge_list(nin,iad(4))%CUTPOINT(1:3,ipose(4))
1354 a5 => edge_list(nin,iad(5))%CUTPOINT(1:3,ipose(5))
1355 n1 => x(1:3,ixs(gnode(1,secid)+1,id))
1356 n2 => x(1:3,ixs(gnode(2,secid)+1,id))
1357 n3 => x(1:3,ixs(gnode(3,secid)+1,id))
1358 !interpolate w_a1
1359 beta = edge_list(nin,iad(1))%CUTCOOR(ipose(1));
1360 j = iad(1)-(ib-1)*12;
1361 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1362 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1363 wa1(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1364 wa1(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1365 wa1(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1366 !interpolate w_a2
1367 beta = edge_list(nin,iad(2))%CUTCOOR(ipose(2));
1368 j = iad(2)-(ib-1)*12;
1369 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1370 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1371 wa2(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1372 wa2(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1373 wa2(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1374 !interpolate w_a3
1375 beta = edge_list(nin,iad(3))%CUTCOOR(ipose(3));
1376 j = iad(3)-(ib-1)*12;
1377 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1378 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1379 wa3(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1380 wa3(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1381 wa3(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1382 !interpolate w_a4
1383 beta = edge_list(nin,iad(4))%CUTCOOR(ipose(4));
1384 j = iad(4)-(ib-1)*12;
1385 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1386 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1387 wa4(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1388 wa4(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1389 wa4(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1390 !interpolate w_a5
1391 beta = edge_list(nin,iad(4))%CUTCOOR(ipose(5));
1392 j = iad(5)-(ib-1)*12;
1393 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1394 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1395 wa5(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1396 wa5(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1397 wa5(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1398 !brick nodes
1399 wn1(1)= w(1,ixs(gnode(1,secid)+1,id))
1400 wn1(2)= w(2,ixs(gnode(1,secid)+1,id))
1401 wn1(3)= w(3,ixs(gnode(1,secid)+1,id))
1402 !brick nodes
1403 wn2(1)= w(1,ixs(gnode(2,secid)+1,id))
1404 wn2(2)= w(2,ixs(gnode(2,secid)+1,id))
1405 wn2(3)= w(3,ixs(gnode(2,secid)+1,id))
1406 !brick nodes
1407 wn3(1)= w(1,ixs(gnode(3,secid)+1,id))
1408 wn3(2)= w(2,ixs(gnode(3,secid)+1,id))
1409 wn3(3)= w(3,ixs(gnode(3,secid)+1,id))
1410 !Grid velocity on face1 :
1411 wf1(1) = one_fifth*(wn1(1)+wn2(1)+wn3(1)+wa3(1)+wa2(1))
1412 wf1(2) = one_fifth*(wn1(2)+wn2(2)+wn3(2)+wa3(2)+wa2(2))
1413 wf1(3) = one_fifth*(wn1(3)+wn2(3)+wn3(3)+wa3(3)+wa2(3))
1414 !Grid velocity on face2 :
1415 wf2(1) = third*(wn3(1)+wa4(1)+wa3(1))
1416 wf2(2) = third*(wn3(2)+wa4(2)+wa3(2))
1417 wf2(3) = third*(wn3(3)+wa4(3)+wa3(3))
1418 !Grid velocity on face3 :
1419 wf3(1) = fourth*(wn3(1)+wn2(1)+wa5(1)+wa4(1))
1420 wf3(2) = fourth*(wn3(2)+wn2(2)+wa5(2)+wa4(2))
1421 wf3(3) = fourth*(wn3(3)+wn2(3)+wa5(3)+wa4(3))
1422 !Grid velocity on face4 :
1423 wf4(1) = fourth*(wn2(1)+wn1(1)+wa1(1)+wa5(1))
1424 wf4(2) = fourth*(wn2(2)+wn1(2)+wa1(2)+wa5(2))
1425 wf4(3) = fourth*(wn2(3)+wn1(3)+wa1(3)+wa5(3))
1426 !Grid velocity on face5 :
1427 wf5(1) = third*(wn1(1)+wa2(1)+wa1(1))
1428 wf5(2) = third*(wn1(2)+wa2(2)+wa1(2))
1429 wf5(3) = third*(wn1(3)+wa2(3)+wa1(3))
1430 !Grid velocity on face0 :
1431 wf0(1) = one_fifth*(wa1(1)+wa2(1)+wa3(1)+wa4(1)+wa5(1))
1432 wf0(2) = one_fifth*(wa1(2)+wa2(2)+wa3(2)+wa4(2)+wa5(2))
1433 wf0(3) = one_fifth*(wa1(3)+wa2(3)+wa3(3)+wa4(3)+wa5(3))
1434 !storage
1435 brick_list(nin,ib)%POLY(icell)%FACE0%W(1:3) = wf0(1:3)
1436 brick_list(nin,ib)%POLY(icell)%FACE(fa(1))%W(1:3) = wf1(1:3)
1437 brick_list(nin,ib)%POLY(icell)%FACE(fa(2))%W(1:3) = wf2(1:3)
1438 brick_list(nin,ib)%POLY(icell)%FACE(fa(3))%W(1:3) = wf3(1:3)
1439 brick_list(nin,ib)%POLY(icell)%FACE(fa(4))%W(1:3) = wf4(1:3)
1440 brick_list(nin,ib)%POLY(icell)%FACE(fa(5))%W(1:3) = wf5(1:3)
1441
1442 ELSEIF(sectype(1:5)=='HEXAE')THEN
1443 nface = f_hexae
1444 fa(1:5) = gface(1:5,secid)
1445 a1 => edge_list(nin,iad(1))%CUTPOINT(1:3,ipose(1))
1446 a2 => edge_list(nin,iad(2))%CUTPOINT(1:3,ipose(2))
1447 a3 => edge_list(nin,iad(3))%CUTPOINT(1:3,ipose(3))
1448 a4 => edge_list(nin,iad(4))%CUTPOINT(1:3,ipose(4))
1449 n1 => x(1:3,ixs(gnode(1,secid)+1,id))
1450 n2 => x(1:3,ixs(gnode(2,secid)+1,id))
1451 n3 => x(1:3,ixs(gnode(3,secid)+1,id))
1452 n4 => x(1:3,ixs(gnode(4,secid)+1,id))
1453 !interpolate w_a1
1454 beta = edge_list(nin,iad(1))%CUTCOOR(ipose(1));
1455 j = iad(1)-(ib-1)*12;
1456 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1457 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1458 wa1(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1459 wa1(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1460 wa1(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1461 !interpolate w_a2
1462 beta = edge_list(nin,iad(2))%CUTCOOR(ipose(2));
1463 j = iad(2)-(ib-1)*12;
1464 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1465 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1466 wa2(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1467 wa2(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1468 wa2(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1469 !interpolate w_a3
1470 beta = edge_list(nin,iad(3))%CUTCOOR(ipose(3));
1471 j = iad(3)-(ib-1)*12;
1472 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1473 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1474 wa3(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1475 wa3(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1476 wa3(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1477 !interpolate w_a4
1478 beta = edge_list(nin,iad(4))%CUTCOOR(ipose(4));
1479 j = iad(4)-(ib-1)*12;
1480 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1481 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1482 wa4(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1483 wa4(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1484 wa4(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1485 !brick nodes
1486 wn1(1)= w(1,ixs(gnode(1,secid)+1,id))
1487 wn1(2)= w(2,ixs(gnode(1,secid)+1,id))
1488 wn1(3)= w(3,ixs(gnode(1,secid)+1,id))
1489 !brick nodes
1490 wn2(1)= w(1,ixs(gnode(2,secid)+1,id))
1491 wn2(2)= w(2,ixs(gnode(2,secid)+1,id))
1492 wn2(3)= w(3,ixs(gnode(2,secid)+1,id))
1493 !brick nodes
1494 wn3(1)= w(1,ixs(gnode(3,secid)+1,id))
1495 wn3(2)= w(2,ixs(gnode(3,secid)+1,id))
1496 wn3(3)= w(3,ixs(gnode(3,secid)+1,id))
1497 !brick nodes
1498 wn4(1)= w(1,ixs(gnode(4,secid)+1,id))
1499 wn4(2)= w(2,ixs(gnode(4,secid)+1,id))
1500 wn4(3)= w(3,ixs(gnode(4,secid)+1,id))
1501 !Grid velocity on face1 :
1502 wf1(1) = fourth*(wn4(1)+wn3(1)+wn2(1)+wn1(1))
1503 wf1(2) = fourth*(wn4(2)+wn3(2)+wn2(2)+wn1(2))
1504 wf1(3) = fourth*(wn4(3)+wn3(3)+wn2(3)+wn1(3))
1505 !Grid velocity on face2 :
1506 wf2(1) = fourth*(wn1(1)+wn2(1)+wa2(1)+wa1(1))
1507 wf2(2) = fourth*(wn1(2)+wn2(2)+wa2(2)+wa1(2))
1508 wf2(3) = fourth*(wn1(3)+wn2(3)+wa2(3)+wa1(3))
1509 !Grid velocity on face3 :
1510 wf3(1) = fourth*(wn2(1)+wn3(1)+wa3(1)+wa2(1))
1511 wf3(2) = fourth*(wn2(2)+wn3(2)+wa3(2)+wa2(2))
1512 wf3(3) = fourth*(wn2(3)+wn3(3)+wa3(3)+wa2(3))
1513 !Grid velocity on face4 :
1514 wf4(1) = fourth*(wn3(1)+wn4(1)+wa4(1)+wa3(1))
1515 wf4(2) = fourth*(wn3(2)+wn4(2)+wa4(2)+wa3(2))
1516 wf4(3) = fourth*(wn3(3)+wn4(3)+wa4(3)+wa3(3))
1517 !Grid velocity on face5 :
1518 wf5(1) = fourth*(wn4(1)+wn1(1)+wa1(1)+wa4(1))
1519 wf5(2) = fourth*(wn4(2)+wn1(2)+wa1(2)+wa4(2))
1520 wf5(3) = fourth*(wn4(3)+wn1(3)+wa1(3)+wa4(3))
1521 !Grid velocity on face0 :
1522 wf0(1) = one_fifth*(wa1(1)+wa2(1)+wa3(1)+wa4(1))
1523 wf0(2) = one_fifth*(wa1(2)+wa2(2)+wa3(2)+wa4(2))
1524 wf0(3) = one_fifth*(wa1(3)+wa2(3)+wa3(3)+wa4(3))
1525 !storage
1526 brick_list(nin,ib)%POLY(icell)%FACE0%W(1:3) = wf0(1:3)
1527 brick_list(nin,ib)%POLY(icell)%FACE(fa(1))%W(1:3) = wf1(1:3)
1528 brick_list(nin,ib)%POLY(icell)%FACE(fa(2))%W(1:3) = wf2(1:3)
1529 brick_list(nin,ib)%POLY(icell)%FACE(fa(3))%W(1:3) = wf3(1:3)
1530 brick_list(nin,ib)%POLY(icell)%FACE(fa(4))%W(1:3) = wf4(1:3)
1531 brick_list(nin,ib)%POLY(icell)%FACE(fa(5))%W(1:3) = wf5(1:3)
1532
1533
1534 ELSEIF(sectype(1:6)=='POLY4 ')THEN
1535 nface = f_poly4
1536 fa(1:6) = gface(1:6,secid)
1537 a1 => edge_list(nin,iad(1))%CUTPOINT(1:3,ipose(1))
1538 a2 => edge_list(nin,iad(2))%CUTPOINT(1:3,ipose(2))
1539 a3 => edge_list(nin,iad(3))%CUTPOINT(1:3,ipose(3))
1540 a4 => edge_list(nin,iad(4))%CUTPOINT(1:3,ipose(4))
1541 a5 => edge_list(nin,iad(5))%CUTPOINT(1:3,ipose(5))
1542 a6 => edge_list(nin,iad(6))%CUTPOINT(1:3,ipose(6))
1543 n1 => x(1:3,ixs(gnode(1,secid)+1,id))
1544 n2 => x(1:3,ixs(gnode(2,secid)+1,id))
1545 n3 => x(1:3,ixs(gnode(3,secid)+1,id))
1546 n4 => x(1:3,ixs(gnode(4,secid)+1,id))
1547 !interpolate w_a1
1548 beta = edge_list(nin,iad(1))%CUTCOOR(ipose(1));
1549 j = iad(1)-(ib-1)*12;
1550 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1551 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1552 wa1(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1553 wa1(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1554 wa1(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1555 !interpolate w_a2
1556 beta = edge_list(nin,iad(2))%CUTCOOR(ipose(2));
1557 j = iad(2)-(ib-1)*12;
1558 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1559 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1560 wa2(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1561 wa2(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1562 wa2(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1563 !interpolate w_a3
1564 beta = edge_list(nin,iad(3))%CUTCOOR(ipose(3));
1565 j = iad(3)-(ib-1)*12;
1566 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1567 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1568 wa3(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1569 wa3(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1570 wa3(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1571 !interpolate w_a4
1572 beta = edge_list(nin,iad(4))%CUTCOOR(ipose(4));
1573 j = iad(4)-(ib-1)*12;
1574 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1575 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1576 wa4(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1577 wa4(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1578 wa4(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1579 !interpolate w_a5
1580 beta = edge_list(nin,iad(5))%CUTCOOR(ipose(5));
1581 j = iad(5)-(ib-1)*12;
1582 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1583 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1584 wa5(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1585 wa5(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1586 wa5(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1587 !interpolate w_a6
1588 beta = edge_list(nin,iad(6))%CUTCOOR(ipose(6));
1589 j = iad(6)-(ib-1)*12;
1590 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1591 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1592 wa6(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1593 wa6(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1594 wa6(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1595 !brick nodes
1596 wn1(1)= w(1,ixs(gnode(1,secid)+1,id))
1597 wn1(2)= w(2,ixs(gnode(1,secid)+1,id))
1598 wn1(3)= w(3,ixs(gnode(1,secid)+1,id))
1599 !brick nodes
1600 wn2(1)= w(1,ixs(gnode(2,secid)+1,id))
1601 wn2(2)= w(2,ixs(gnode(2,secid)+1,id))
1602 wn2(3)= w(3,ixs(gnode(2,secid)+1,id))
1603 !brick nodes
1604 wn3(1)= w(1,ixs(gnode(3,secid)+1,id))
1605 wn3(2)= w(2,ixs(gnode(3,secid)+1,id))
1606 wn3(3)= w(3,ixs(gnode(3,secid)+1,id))
1607 !brick nodes
1608 wn4(1)= w(1,ixs(gnode(4,secid)+1,id))
1609 wn4(2)= w(2,ixs(gnode(4,secid)+1,id))
1610 wn4(3)= w(3,ixs(gnode(4,secid)+1,id))
1611 !Grid velocity on face1 :
1612 wf1(1) = third*(wn1(1)+wa1(1)+wa6(1))
1613 wf1(2) = third*(wn1(2)+wa1(2)+wa6(2))
1614 wf1(3) = third*(wn1(3)+wa1(3)+wa6(3))
1615 !Grid velocity on face2 :
1616 wf2(1) = fourth*(wn1(1)+wn2(1)+wa2(1)+wa1(1))
1617 wf2(2) = fourth*(wn1(2)+wn2(2)+wa2(2)+wa1(2))
1618 wf2(3) = fourth*(wn1(3)+wn2(3)+wa2(3)+wa1(3))
1619 !Grid velocity on face3 :
1620 wf3(1) = one_fifth*(wn2(1)+wn3(1)+wn4(1)+wa3(1)+wa2(1))
1621 wf3(2) = one_fifth*(wn2(2)+wn3(2)+wn4(2)+wa3(2)+wa2(2))
1622 wf3(3) = one_fifth*(wn2(3)+wn3(3)+wn4(3)+wa3(3)+wa2(3))
1623 !Grid velocity on face4 :
1624 wf4(1) = third*(wn4(1)+wa4(1)+wa3(1))
1625 wf4(2) = third*(wn4(2)+wa4(2)+wa3(2))
1626 wf4(3) = third*(wn4(3)+wa4(3)+wa3(3))
1627 !Grid velocity on face5 :
1628 wf5(1) = fourth*(wn4(1)+wn3(1)+wa5(1)+wa4(1))
1629 wf5(2) = fourth*(wn4(2)+wn3(2)+wa5(2)+wa4(2))
1630 wf5(3) = fourth*(wn4(3)+wn3(3)+wa5(3)+wa4(3))
1631 !Grid velocity on face6 :
1632 wf6(1) = fourth*(wn3(1)+wn2(1)+wn1(1)+wa6(1)+wa5(1))
1633 wf6(2) = fourth*(wn3(2)+wn2(2)+wn1(2)+wa6(2)+wa5(2))
1634 wf6(3) = fourth*(wn3(3)+wn2(3)+wn1(3)+wa6(3)+wa5(3))
1635 !Grid velocity on face0 :
1636 wf0(1) = one_over_6*(wa1(1)+wa2(1)+wa3(1)+wa4(1)+wa5(1)+wa6(1))
1637 wf0(2) = one_over_6*(wa1(2)+wa2(2)+wa3(2)+wa4(2)+wa5(2)+wa6(2))
1638 wf0(3) = one_over_6*(wa1(3)+wa2(3)+wa3(3)+wa4(3)+wa5(3)+wa6(3))
1639 !storage
1640 brick_list(nin,ib)%POLY(icell)%FACE0%W(1:3) = wf0(1:3)
1641 brick_list(nin,ib)%POLY(icell)%FACE(fa(1))%W(1:3) = wf1(1:3)
1642 brick_list(nin,ib)%POLY(icell)%FACE(fa(2))%W(1:3) = wf2(1:3)
1643 brick_list(nin,ib)%POLY(icell)%FACE(fa(3))%W(1:3) = wf3(1:3)
1644 brick_list(nin,ib)%POLY(icell)%FACE(fa(4))%W(1:3) = wf4(1:3)
1645 brick_list(nin,ib)%POLY(icell)%FACE(fa(5))%W(1:3) = wf5(1:3)
1646 brick_list(nin,ib)%POLY(icell)%FACE(fa(6))%W(1:3) = wf6(1:3)
1647
1648
1649 ELSEIF(sectype(1:6)=='POLY4A')THEN
1650 nface = f_poly4a
1651 fa(1:6) = gface(1:6,secid)
1652 a1 => edge_list(nin,iad(1))%CUTPOINT(1:3,ipose(1))
1653 a2 => edge_list(nin,iad(2))%CUTPOINT(1:3,ipose(2))
1654 a3 => edge_list(nin,iad(3))%CUTPOINT(1:3,ipose(3))
1655 a4 => edge_list(nin,iad(4))%CUTPOINT(1:3,ipose(4))
1656 a5 => edge_list(nin,iad(5))%CUTPOINT(1:3,ipose(5))
1657 a6 => edge_list(nin,iad(6))%CUTPOINT(1:3,ipose(6))
1658 n1 => x(1:3,ixs(gnode(1,secid)+1,id))
1659 n2 => x(1:3,ixs(gnode(2,secid)+1,id))
1660 n3 => x(1:3,ixs(gnode(3,secid)+1,id))
1661 n4 => x(1:3,ixs(gnode(4,secid)+1,id))
1662 !interpolate w_a1
1663 beta = edge_list(nin,iad(1))%CUTCOOR(ipose(1));
1664 j = iad(1)-(ib-1)*12;
1665 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1666 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1667 wa1(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1668 wa1(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1669 wa1(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1670 !interpolate w_a2
1671 beta = edge_list(nin,iad(2))%CUTCOOR(ipose(2));
1672 j = iad(2)-(ib-1)*12;
1673 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1674 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1675 wa2(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1676 wa2(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1677 wa2(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1678 !interpolate w_a3
1679 beta = edge_list(nin,iad(3))%CUTCOOR(ipose(3));
1680 j = iad(3)-(ib-1)*12;
1681 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1682 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1683 wa3(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1684 wa3(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1685 wa3(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1686 !interpolate w_a4
1687 beta = edge_list(nin,iad(4))%CUTCOOR(ipose(4));
1688 j = iad(4)-(ib-1)*12;
1689 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1690 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1691 wa4(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1692 wa4(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1693 wa4(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1694 !interpolate w_a5
1695 beta = edge_list(nin,iad(5))%CUTCOOR(ipose(5));
1696 j = iad(5)-(ib-1)*12;
1697 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1698 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1699 wa5(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1700 wa5(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1701 wa5(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1702 !interpolate w_a6
1703 beta = edge_list(nin,iad(6))%CUTCOOR(ipose(6));
1704 j = iad(6)-(ib-1)*12;
1705 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1706 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1707 wa6(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1708 wa6(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1709 wa6(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1710 !brick nodes
1711 wn1(1)= w(1,ixs(gnode(1,secid)+1,id))
1712 wn1(2)= w(2,ixs(gnode(1,secid)+1,id))
1713 wn1(3)= w(3,ixs(gnode(1,secid)+1,id))
1714 !brick nodes
1715 wn2(1)= w(1,ixs(gnode(2,secid)+1,id))
1716 wn2(2)= w(2,ixs(gnode(2,secid)+1,id))
1717 wn2(3)= w(3,ixs(gnode(2,secid)+1,id))
1718 !brick nodes
1719 wn3(1)= w(1,ixs(gnode(3,secid)+1,id))
1720 wn3(2)= w(2,ixs(gnode(3,secid)+1,id))
1721 wn3(3)= w(3,ixs(gnode(3,secid)+1,id))
1722 !brick nodes
1723 wn4(1)= w(1,ixs(gnode(4,secid)+1,id))
1724 wn4(2)= w(2,ixs(gnode(4,secid)+1,id))
1725 wn4(3)= w(3,ixs(gnode(4,secid)+1,id))
1726 !Grid velocity on face1 :
1727 wf1(1) = third*(wn1(1)+wa1(1)+wa6(1))
1728 wf1(2) = third*(wn1(2)+wa1(2)+wa6(2))
1729 wf1(3) = third*(wn1(3)+wa1(3)+wa6(3))
1730 !Grid velocity on face2 :
1731 wf2(1) = fourth*(wn1(1)+wn2(1)+wa2(1)+wa1(1))
1732 wf2(2) = fourth*(wn1(2)+wn2(2)+wa2(2)+wa1(2))
1733 wf2(3) = fourth*(wn1(3)+wn2(3)+wa2(3)+wa1(3))
1734 !Grid velocity on face3 :
1735 wf3(1) = one_fifth*(wn2(1)+wn3(1)+wn4(1)+wa3(1)+wa2(1))
1736 wf3(2) = one_fifth*(wn2(2)+wn3(2)+wn4(2)+wa3(2)+wa2(2))
1737 wf3(3) = one_fifth*(wn2(3)+wn3(3)+wn4(3)+wa3(3)+wa2(3))
1738 !Grid velocity on face4 :
1739 wf4(1) = third*(wn4(1)+wa4(1)+wa3(1))
1740 wf4(2) = third*(wn4(2)+wa4(2)+wa3(2))
1741 wf4(3) = third*(wn4(3)+wa4(3)+wa3(3))
1742 !Grid velocity on face5 :
1743 wf5(1) = fourth*(wn4(1)+wn3(1)+wa5(1)+wa4(1))
1744 wf5(2) = fourth*(wn4(2)+wn3(2)+wa5(2)+wa4(2))
1745 wf5(3) = fourth*(wn4(3)+wn3(3)+wa5(3)+wa4(3))
1746 !Grid velocity on face6 :
1747 wf6(1) = fourth*(wn3(1)+wn2(1)+wn1(1)+wa6(1)+wa5(1))
1748 wf6(2) = fourth*(wn3(2)+wn2(2)+wn1(2)+wa6(2)+wa5(2))
1749 wf6(3) = fourth*(wn3(3)+wn2(3)+wn1(3)+wa6(3)+wa5(3))
1750 !Grid velocity on face0 :
1751 wf0(1) = one_over_6*(wa1(1)+wa2(1)+wa3(1)+wa4(1)+wa5(1)+wa6(1))
1752 wf0(2) = one_over_6*(wa1(2)+wa2(2)+wa3(2)+wa4(2)+wa5(2)+wa6(2))
1753 wf0(3) = one_over_6*(wa1(3)+wa2(3)+wa3(3)+wa4(3)+wa5(3)+wa6(3))
1754 !storage
1755 brick_list(nin,ib)%POLY(icell)%FACE0%W(1:3) = wf0(1:3)
1756 brick_list(nin,ib)%POLY(icell)%FACE(fa(1))%W(1:3) = wf1(1:3)
1757 brick_list(nin,ib)%POLY(icell)%FACE(fa(2))%W(1:3) = wf2(1:3)
1758 brick_list(nin,ib)%POLY(icell)%FACE(fa(3))%W(1:3) = wf3(1:3)
1759 brick_list(nin,ib)%POLY(icell)%FACE(fa(4))%W(1:3) = wf4(1:3)
1760 brick_list(nin,ib)%POLY(icell)%FACE(fa(5))%W(1:3) = wf5(1:3)
1761 brick_list(nin,ib)%POLY(icell)%FACE(fa(6))%W(1:3) = wf6(1:3)
1762
1763 ELSEIF(sectype(1:6)=='POLY4B')THEN
1764 nface = f_poly4b
1765 fa(1:6) = gface(1:6,secid)
1766 a1 => edge_list(nin,iad(1))%CUTPOINT(1:3,ipose(1))
1767 a2 => edge_list(nin,iad(2))%CUTPOINT(1:3,ipose(2))
1768 a3 => edge_list(nin,iad(3))%CUTPOINT(1:3,ipose(3))
1769 a4 => edge_list(nin,iad(4))%CUTPOINT(1:3,ipose(4))
1770 a5 => edge_list(nin,iad(5))%CUTPOINT(1:3,ipose(5))
1771 a6 => edge_list(nin,iad(6))%CUTPOINT(1:3,ipose(6))
1772 n1 => x(1:3,ixs(gnode(1,secid)+1,id))
1773 n2 => x(1:3,ixs(gnode(2,secid)+1,id))
1774 n3 => x(1:3,ixs(gnode(3,secid)+1,id))
1775 n4 => x(1:3,ixs(gnode(4,secid)+1,id))
1776 !interpolate w_a1
1777 beta = edge_list(nin,iad(1))%CUTCOOR(ipose(1));
1778 j = iad(1)-(ib-1)*12;
1779 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1780 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1781 wa1(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1782 wa1(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1783 wa1(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1784 !interpolate w_a2
1785 beta = edge_list(nin,iad(2))%CUTCOOR(ipose(2));
1786 j = iad(2)-(ib-1)*12;
1787 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1788 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1789 wa2(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1790 wa2(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1791 wa2(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1792 !interpolate w_a3
1793 beta = edge_list(nin,iad(3))%CUTCOOR(ipose(3));
1794 j = iad(3)-(ib-1)*12;
1795 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1796 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1797 wa3(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1798 wa3(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1799 wa3(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1800 !interpolate w_a4
1801 beta = edge_list(nin,iad(4))%CUTCOOR(ipose(4));
1802 j = iad(4)-(ib-1)*12;
1803 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1804 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1805 wa4(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1806 wa4(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1807 wa4(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1808 !interpolate w_a5
1809 beta = edge_list(nin,iad(5))%CUTCOOR(ipose(5));
1810 j = iad(5)-(ib-1)*12;
1811 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1812 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1813 wa5(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1814 wa5(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1815 wa5(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1816 !interpolate w_a6
1817 beta = edge_list(nin,iad(6))%CUTCOOR(ipose(6));
1818 j = iad(6)-(ib-1)*12;
1819 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1820 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1821 wa6(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1822 wa6(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1823 wa6(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1824 !brick nodes
1825 wn1(1)= w(1,ixs(gnode(1,secid)+1,id))
1826 wn1(2)= w(2,ixs(gnode(1,secid)+1,id))
1827 wn1(3)= w(3,ixs(gnode(1,secid)+1,id))
1828 !brick nodes
1829 wn2(1)= w(1,ixs(gnode(2,secid)+1,id))
1830 wn2(2)= w(2,ixs(gnode(2,secid)+1,id))
1831 wn2(3)= w(3,ixs(gnode(2,secid)+1,id))
1832 !brick nodes
1833 wn3(1)= w(1,ixs(gnode(3,secid)+1,id))
1834 wn3(2)= w(2,ixs(gnode(3,secid)+1,id))
1835 wn3(3)= w(3,ixs(gnode(3,secid)+1,id))
1836 !brick nodes
1837 wn4(1)= w(1,ixs(gnode(4,secid)+1,id))
1838 wn4(2)= w(2,ixs(gnode(4,secid)+1,id))
1839 wn4(3)= w(3,ixs(gnode(4,secid)+1,id))
1840 !Grid velocity on face1 :
1841 wf1(1) = third*(wn1(1)+wa1(1)+wa6(1))
1842 wf1(2) = third*(wn1(2)+wa1(2)+wa6(2))
1843 wf1(3) = third*(wn1(3)+wa1(3)+wa6(3))
1844 !Grid velocity on face2 :
1845 wf2(1) = fourth*(wn1(1)+wn2(1)+wa2(1)+wa1(1))
1846 wf2(2) = fourth*(wn1(2)+wn2(2)+wa2(2)+wa1(2))
1847 wf2(3) = fourth*(wn1(3)+wn2(3)+wa2(3)+wa1(3))
1848 !Grid velocity on face3 :
1849 wf3(1) = one_fifth*(wn2(1)+wn3(1)+wn4(1)+wa3(1)+wa2(1))
1850 wf3(2) = one_fifth*(wn2(2)+wn3(2)+wn4(2)+wa3(2)+wa2(2))
1851 wf3(3) = one_fifth*(wn2(3)+wn3(3)+wn4(3)+wa3(3)+wa2(3))
1852 !Grid velocity on face4 :
1853 wf4(1) = third*(wn4(1)+wa4(1)+wa3(1))
1854 wf4(2) = third*(wn4(2)+wa4(2)+wa3(2))
1855 wf4(3) = third*(wn4(3)+wa4(3)+wa3(3))
1856 !Grid velocity on face5 :
1857 wf5(1) = fourth*(wn4(1)+wn3(1)+wa5(1)+wa4(1))
1858 wf5(2) = fourth*(wn4(2)+wn3(2)+wa5(2)+wa4(2))
1859 wf5(3) = fourth*(wn4(3)+wn3(3)+wa5(3)+wa4(3))
1860 !Grid velocity on face6 :
1861 wf6(1) = fourth*(wn3(1)+wn2(1)+wn1(1)+wa6(1)+wa5(1))
1862 wf6(2) = fourth*(wn3(2)+wn2(2)+wn1(2)+wa6(2)+wa5(2))
1863 wf6(3) = fourth*(wn3(3)+wn2(3)+wn1(3)+wa6(3)+wa5(3))
1864 !Grid velocity on face0 :
1865 wf0(1) = one_over_6*(wa1(1)+wa2(1)+wa3(1)+wa4(1)+wa5(1)+wa6(1))
1866 wf0(2) = one_over_6*(wa1(2)+wa2(2)+wa3(2)+wa4(2)+wa5(2)+wa6(2))
1867 wf0(3) = one_over_6*(wa1(3)+wa2(3)+wa3(3)+wa4(3)+wa5(3)+wa6(3))
1868 !storage
1869 brick_list(nin,ib)%POLY(icell)%FACE0%W(1:3) = wf0(1:3)
1870 brick_list(nin,ib)%POLY(icell)%FACE(fa(1))%W(1:3) = wf1(1:3)
1871 brick_list(nin,ib)%POLY(icell)%FACE(fa(2))%W(1:3) = wf2(1:3)
1872 brick_list(nin,ib)%POLY(icell)%FACE(fa(3))%W(1:3) = wf3(1:3)
1873 brick_list(nin,ib)%POLY(icell)%FACE(fa(4))%W(1:3) = wf4(1:3)
1874 brick_list(nin,ib)%POLY(icell)%FACE(fa(5))%W(1:3) = wf5(1:3)
1875 brick_list(nin,ib)%POLY(icell)%FACE(fa(6))%W(1:3) = wf6(1:3)
1876
1877 ELSEIF(sectype(1:5)=='POLYC')THEN
1878 nface = f_polyc
1879 fa(1:6) = gface(1:6,secid)
1880 a1 => edge_list(nin,iad(1))%CUTPOINT(1:3,ipose(1))
1881 a2 => edge_list(nin,iad(2))%CUTPOINT(1:3,ipose(2))
1882 a3 => edge_list(nin,iad(3))%CUTPOINT(1:3,ipose(3))
1883 a4 => edge_list(nin,iad(4))%CUTPOINT(1:3,ipose(4))
1884 a5 => edge_list(nin,iad(5))%CUTPOINT(1:3,ipose(5)) !points 5,6 on same edge : IAD(5)=IAD(6) BUT ICUT=1
1885 a6 => edge_list(nin,iad(6))%CUTPOINT(1:3,ipose(6)) !points 5,6 on same edge : IAD(5)=IAD(6) BUT ICUT=2
1886 n1 => x(1:3,ixs(gnode(1,secid)+1,id))
1887 n2 => x(1:3,ixs(gnode(2,secid)+1,id))
1888 n3 => x(1:3,ixs(gnode(3,secid)+1,id))
1889 n4 => x(1:3,ixs(gnode(4,secid)+1,id))
1890 !interpolate w_a1
1891 beta = edge_list(nin,iad(1))%CUTCOOR(ipose(1));
1892 j = iad(1)-(ib-1)*12;
1893 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1894 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1895 wa1(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1896 wa1(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1897 wa1(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1898 !interpolate w_a2
1899 beta = edge_list(nin,iad(2))%CUTCOOR(ipose(2));
1900 j = iad(2)-(ib-1)*12;
1901 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1902 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1903 wa2(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1904 wa2(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1905 wa2(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1906 !interpolate w_a3
1907 beta = edge_list(nin,iad(3))%CUTCOOR(ipose(3));
1908 j = iad(3)-(ib-1)*12;
1909 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1910 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1911 wa3(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1912 wa3(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1913 wa3(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1914 !interpolate w_a4
1915 beta = edge_list(nin,iad(4))%CUTCOOR(ipose(4));
1916 j = iad(4)-(ib-1)*12;
1917 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1918 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1919 wa4(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1920 wa4(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1921 wa4(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1922 !interpolate w_a5
1923 beta = edge_list(nin,iad(5))%CUTCOOR(ipose(5));
1924 j = iad(5)-(ib-1)*12;
1925 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1926 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1927 wa5(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1928 wa5(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1929 wa5(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1930 !interpolate w_a6
1931 beta = edge_list(nin,iad(6))%CUTCOOR(ipose(6));
1932 j = iad(6)-(ib-1)*12;
1933 edgnod1 = ixs(1+int22_buf%iEDGE(1,j),id);
1934 edgnod2 = ixs(1+int22_buf%iEDGE(2,j),id)
1935 wa6(1) = beta * w(1,edgnod1) + (one-beta)*w(1,edgnod2)
1936 wa6(2) = beta * w(2,edgnod1) + (one-beta)*w(2,edgnod2)
1937 wa6(3) = beta * w(3,edgnod1) + (one-beta)*w(3,edgnod2)
1938 !brick nodes
1939 wn1(1)= w(1,ixs(gnode(1,secid)+1,id))
1940 wn1(2)= w(2,ixs(gnode(1,secid)+1,id))
1941 wn1(3)= w(3,ixs(gnode(1,secid)+1,id))
1942 !brick nodes
1943 wn2(1)= w(1,ixs(gnode(2,secid)+1,id))
1944 wn2(2)= w(2,ixs(gnode(2,secid)+1,id))
1945 wn2(3)= w(3,ixs(gnode(2,secid)+1,id))
1946 !brick nodes
1947 wn3(1)= w(1,ixs(gnode(3,secid)+1,id))
1948 wn3(2)= w(2,ixs(gnode(3,secid)+1,id))
1949 wn3(3)= w(3,ixs(gnode(3,secid)+1,id))
1950 !brick nodes
1951 wn4(1)= w(1,ixs(gnode(4,secid)+1,id))
1952 wn4(2)= w(2,ixs(gnode(4,secid)+1,id))
1953 wn4(3)= w(3,ixs(gnode(4,secid)+1,id))
1954 !Grid velocity on face1 :
1955 wf1(1) = fourth*(wn1(1)+wn2(1)+wn3(1)+wn4(1))
1956 wf1(2) = fourth*(wn1(2)+wn2(2)+wn3(2)+wn4(2))
1957 wf1(3) = fourth*(wn1(3)+wn2(3)+wn3(3)+wn4(3))
1958 !Grid velocity on face2 :
1959 wf2(1) = fourth*(wn1(1)+wa1(1)+wa2(1)+wn2(1))
1960 wf2(2) = fourth*(wn1(2)+wa1(2)+wa2(2)+wn2(2))
1961 wf2(3) = fourth*(wn1(3)+wa1(3)+wa2(3)+wn2(3))
1962 !Grid velocity on face3 :
1963 wf3(1) = fourth*(wn2(1)+wa2(1)+wa3(1)+wn3(1))
1964 wf3(2) = fourth*(wn2(2)+wa2(2)+wa3(2)+wn3(2))
1965 wf3(3) = fourth*(wn2(3)+wa2(3)+wa3(3)+wn3(3))
1966 !Grid velocity on face4 :
1967 wf4(1) = fourth*(wn1(1)+wn4(1)+wa4(1)+wa1(1))
1968 wf4(2) = fourth*(wn1(2)+wn4(2)+wa4(2)+wa1(2))
1969 wf4(3) = fourth*(wn1(3)+wn4(3)+wa4(3)+wa1(3))
1970 !Grid velocity on face5 :
1971 wf5(1) = (wn4(1)+wn3(1)+wa3(1)+wa6(1)+wa5(1)+wa4(1)+wn4(1))/seven
1972 wf5(2) = (wn4(2)+wn3(2)+wa3(2)+wa6(2)+wa5(2)+wa4(2)+wn4(2))/seven
1973 wf5(3) = (wn4(3)+wn3(3)+wa3(3)+wa6(3)+wa5(3)+wa4(3)+wn4(3))/seven
1974 !Grid velocity on face0 :
1975 wf0(1) = one_over_6*(wa1(1)+wa2(1)+wa3(1)+wa4(1)+wa5(1)+wa6(1))
1976 wf0(2) = one_over_6*(wa1(2)+wa2(2)+wa3(2)+wa4(2)+wa5(2)+wa6(2))
1977 wf0(3) = one_over_6*(wa1(3)+wa2(3)+wa3(3)+wa4(3)+wa5(3)+wa6(3))
1978 !storage
1979 brick_list(nin,ib)%POLY(icell)%FACE0%W(1:3) = wf0(1:3)
1980 brick_list(nin,ib)%POLY(icell)%FACE(fa(1))%W(1:3) = wf1(1:3)
1981 brick_list(nin,ib)%POLY(icell)%FACE(fa(2))%W(1:3) = wf2(1:3)
1982 brick_list(nin,ib)%POLY(icell)%FACE(fa(3))%W(1:3) = wf3(1:3)
1983 brick_list(nin,ib)%POLY(icell)%FACE(fa(4))%W(1:3) = wf4(1:3)
1984 brick_list(nin,ib)%POLY(icell)%FACE(fa(5))%W(1:3) = wf5(1:3)
1985
1986 END IF
1987 ENDIF ! ILAE/=0
1988
1989
1990
1991
1992
1993
1994 END DO !ICELL
1995
1996
1997 ntag = sum(itag(1:8))
1998 IF(ntag>0)THEN
1999 !CLEAN DEGENERATED TETRA & HEXA (USELESS AND WOULD NEED TO BE ATTACHED TO A main CELL WHICH COULD NOT EXIST)
2000 numpt = 0
2001 DO k=1,8
2002 tmp_sectype(k) = brick_list(nin,ib)%SECTYPE(k)
2003 tmp_secid_cell(k) = brick_list(nin,ib)%SecID_CELL(k)
2004 !tmp_NumNOD(K) = BRICK_LIST(NIN,IB)%POLY(K)%NumNOD
2005 !tmp_NumPOINT(K) = BRICK_LIST(NIN,IB)%POLY(K)%NumPOINT
2006 tmppoly(k) = brick_list(nin,ib)%POLY(k)
2007 tmpcut(k) = brick_list(nin,ib)%PCUT(k)
2008 ENDDO
2009 !polyhedron K deletion, updating identifiers (resorting) for polyhedrons ids K+1 (if K==8, nothing to do, last one deleted, nothing to resort)
2010 DO k=1,8
2011 IF(itag(k)==1)THEN
2012 DO n=k,7
2013 brick_list(nin,ib)%SECTYPE(n) = tmp_sectype(n+1)
2014 brick_list(nin,ib)%SecID_CELL(n) = tmp_secid_cell(n+1)
2015 !BRICK_LIST(NIN,IB)%POLY(N)%NumNOD = tmp_NumNOD(N+1)
2016 !BRICK_LIST(NIN,IB)%POLY(N)%NumPOINT = tmp_NumPOINT(N+1)
2017 brick_list(nin,ib)%POLY(n) = tmppoly(n+1)
2018 brick_list(nin,ib)%PCUT(n) = tmpcut(n+1)
2019 ENDDO
2020 numpt=numpt+1
2021 !updating buffer whichcellnode. find a cell_id from a local node_id. This changed due to polyhedra k deletion and resorting
2022 DO n=1,8
2023 inod = brick_list(nin,ib)%NODE(n)%WhichCell
2024 IF(inod >= k)THEN
2025 brick_list(nin,ib)%NODE(n)%WhichCell = inod - ntag
2026 ENDIF
2027 ENDDO
2028 ENDIF
2029 ENDDO
2030 brick_list(nin,ib)%NBCUT = brick_list(nin,ib)%NBCUT - numpt
2031 DO k=brick_list(nin,ib)%NBCUT+1,8
2032 brick_list(nin,ib)%PCUT(k)%NP = 0
2033 brick_list(nin,ib)%POLY(k)%NumNOD = 0
2034 brick_list(nin,ib)%POLY(k)%NumPOINT = 0
2035 brick_list(nin,ib)%SecID_CELL(k) = 0
2036 brick_list(nin,ib)%SECTYPE(k) = '______________'
2037 ENDDO
2038 ENDIF !(NTAG>0)
2039
2040 !----------------------------------------------------------------
2041 ! Removing Unused Intersection Point.
2042 ! (Otherwise troube with complementary cell #9 in sinit22.F)
2043 !----------------------------------------------------------------
2044 DO j=1,12
2045 k = (ib-1)*12+j
2046 !EdgeList(J) : real number of points used on edge J
2047 !EdgePos(J,1:2) : their positiono n edge 1,or 2, or 1&2
2048 IF(edgelist(j) /= edge_list(nin,k)%NBCUT)THEN
2049 IF(edgelist(j)==1)THEN
2050 edge_list(nin,k)%NBCUT = 1
2051 !keep only relevant node
2052 IF(edgepos(j,1)==1)THEN
2053 !do nothing, keep first one
2054 ELSEIF(edgepos(j,2)==1)THEN
2055 !switch
2056 edge_list(nin,k)%CUTSHELL(1) = edge_list(nin,k)%CUTSHELL(2)
2057 edge_list(nin,k)%CUTCOOR(1) = edge_list(nin,k)%CUTCOOR(2)
2058 edge_list(nin,k)%CUTPOINT(1:3,1) = edge_list(nin,k)%CUTPOINT(1:3,2)
2059 edge_list(nin,k)%CUTVEL(1:3,1) = edge_list(nin,k)%CUTVEL(1:3,2)
2060 ENDIF
2061 ELSEIF(edgelist(j)==0)THEN
2062 edge_list(nin,k)%NBCUT = 0
2063 ENDIF
2064 ENDIF
2065 ENDDO
2066
2067
2068 !----------------------------------------------------------------
2069 ! Build Cell 9 which is complementary one.
2070 !----------------------------------------------------------------
2071 ! moved to sinit22.F
2072
2073
2074 IF(brick_list(nin,ib)%NBCUT == 0)THEN
2075 icell = 1
2076 pnnodcell(icell)%NumNOD = 8 !main(brick nodes)
2077 pnpointcell(icell)%NumPOINT = 8 !main + intersection
2078 n1 => x(1:3, ixs(1+1,brick_list(nin,ib)%ID ))
2079 n2 => x(1:3, ixs(1+2,brick_list(nin,ib)%ID ))
2080 n3 => x(1:3, ixs(1+3,brick_list(nin,ib)%ID ))
2081 n4 => x(1:3, ixs(1+4,brick_list(nin,ib)%ID ))
2082 n5 => x(1:3, ixs(1+5,brick_list(nin,ib)%ID ))
2083 n6 => x(1:3, ixs(1+6,brick_list(nin,ib)%ID ))
2084 n7 => x(1:3, ixs(1+7,brick_list(nin,ib)%ID ))
2085 n8 => x(1:3, ixs(1+8,brick_list(nin,ib)%ID ))
2086 brick_list(nin,ib)%POLY(icell)%CellCenter(1)= sum( (/n1(1),n2(1),n3(1),n4(1),n5(1),n6(1),n7(1),n8(1)/) ) / eight
2087 brick_list(nin,ib)%POLY(icell)%CellCenter(2)= sum( (/n1(2),n2(2),n3(2),n4(2),n5(2),n6(2),n7(2),n8(2)/) ) / eight
2088 brick_list(nin,ib)%POLY(icell)%CellCenter(3)= sum( (/n1(3),n2(3),n3(3),n4(3),n5(3),n6(3),n7(3),n8(3)/) ) / eight
2089 brick_list(nin,ib)%POLY(icell)%ListNodID(1:8) = (/ 1,2,3,4,5,6,7,8/)
2090 brick_list(nin,ib)%NODE(1:8)%WhichCell = icell
2091 ENDIF
2092
2093 END DO !next IB (cut cell)
2094
2095
2096 !----------------------------------------------------------------
2097 ! double Multiplicity Intersection.
2098 !----------------------------------------------------------------
2099 DO ib=nbf,nbl !1,NB
2100 IF(brick_list(nin,ib)%NBCUT == 2)THEN
2101 sectype = brick_list(nin,ib)%SECTYPE(1)
2102 sectype2 = brick_list(nin,ib)%SECTYPE(2)
2103 cod1 = iabs(brick_list(nin,ib)%SecID_CELL(1))
2104 cod2 = iabs(brick_list(nin,ib)%SecID_CELL(2))
2105 IF(cod1==cod2)THEN
2106 vol1 = brick_list(nin,ib)%POLY(1)%Vnew
2107 vol2 = brick_list(nin,ib)%POLY(2)%Vnew
2108 numpt = brick_list(nin,ib)%POLY(1)%NumPOINT
2109
2110 IF(vol1>vol2)THEN
2111
2112 valtmp(1:3) = brick_list(nin,ib)%POLY(1)%CellCenter(1:3)
2113
2114 !backup poly1
2115 brick_list(nin,ib)%POLY(3)%FACE(1:6)%Surf = brick_list(nin,ib)%POLY(1)%FACE(1:6)%Surf
2116 brick_list(nin,ib)%PCUT(3)%SCUT = brick_list(nin,ib)%PCUT(1)%SCUT
2117 brick_list(nin,ib)%PCUT(3)%B(1:3) = brick_list(nin,ib)%PCUT(1)%B(1:3)
2118 brick_list(nin,ib)%PCUT(3)%N(1:3) = brick_list(nin,ib)%PCUT(1)%N(1:3)
2119 brick_list(nin,ib)%PCUT(3)%NP = brick_list(nin,ib)%PCUT(1)%NP
2120 brick_list(nin,ib)%PCUT(3)%P(1:3,1) = brick_list(nin,ib)%PCUT(1)%P(1:3,1)
2121 brick_list(nin,ib)%PCUT(3)%P(1:3,2) = brick_list(nin,ib)%PCUT(1)%P(1:3,2)
2122 brick_list(nin,ib)%PCUT(3)%P(1:3,3) = brick_list(nin,ib)%PCUT(1)%P(1:3,3)
2123 brick_list(nin,ib)%PCUT(3)%P(1:3,4) = brick_list(nin,ib)%PCUT(1)%P(1:3,4)
2124 brick_list(nin,ib)%PCUT(3)%P(1:3,5) = brick_list(nin,ib)%PCUT(1)%P(1:3,5)
2125 brick_list(nin,ib)%PCUT(3)%P(1:3,6) = brick_list(nin,ib)%PCUT(1)%P(1:3,6)
2126 brick_list(nin,ib)%POLY(3)%Vnew = brick_list(nin,ib)%POLY(1)%Vnew
2127 brick_list(nin,ib)%POLY(3)%NumNOD = brick_list(nin,ib)%POLY(1)%NumNOD
2128 brick_list(nin,ib)%POLY(3)%NumPOINT = brick_list(nin,ib)%POLY(1)%NumPOINT
2129 brick_list(nin,ib)%POLY(3)%ListNodID(1:8) = brick_list(nin,ib)%POLY(1)%ListNodID(1:8)
2130 brick_list(nin,ib)%PCUT(3)%Vel(1:3) = brick_list(nin,ib)%PCUT(1)%Vel(1:3)
2131
2132 !switch poly2 into poly1
2133 !1. poly1 <- poly2
2134 brick_list(nin,ib)%POLY(1)%FACE(1:6)%Surf = brick_list(nin,ib)%POLY(2)%FACE(1:6)%Surf
2135 brick_list(nin,ib)%PCUT(1)%SCUT = brick_list(nin,ib)%PCUT(2)%SCUT
2136 brick_list(nin,ib)%PCUT(1)%B(1:3) = brick_list(nin,ib)%PCUT(2)%B(1:3)
2137 brick_list(nin,ib)%PCUT(1)%N(1:3) = brick_list(nin,ib)%PCUT(2)%N(1:3)
2138 brick_list(nin,ib)%PCUT(1)%NP = brick_list(nin,ib)%PCUT(2)%NP
2139 brick_list(nin,ib)%PCUT(1)%P(1:3,1) = brick_list(nin,ib)%PCUT(2)%P(1:3,1)
2140 brick_list(nin,ib)%PCUT(1)%P(1:3,2) = brick_list(nin,ib)%PCUT(2)%P(1:3,2)
2141 brick_list(nin,ib)%PCUT(1)%P(1:3,3) = brick_list(nin,ib)%PCUT(2)%P(1:3,3)
2142 brick_list(nin,ib)%PCUT(1)%P(1:3,4) = brick_list(nin,ib)%PCUT(2)%P(1:3,4)
2143 brick_list(nin,ib)%PCUT(1)%P(1:3,5) = brick_list(nin,ib)%PCUT(2)%P(1:3,5)
2144 brick_list(nin,ib)%PCUT(1)%P(1:3,6) = brick_list(nin,ib)%PCUT(2)%P(1:3,6)
2145 brick_list(nin,ib)%POLY(1)%Vnew = brick_list(nin,ib)%POLY(2)%Vnew
2146 brick_list(nin,ib)%POLY(1)%NumNOD = brick_list(nin,ib)%POLY(2)%NumNOD
2147 brick_list(nin,ib)%POLY(1)%NumPOINT = brick_list(nin,ib)%POLY(2)%NumPOINT
2148 brick_list(nin,ib)%POLY(1)%ListNodID(1:8) = brick_list(nin,ib)%POLY(2)%ListNodID(1:8)
2149 brick_list(nin,ib)%PCUT(1)%Vel(1:3) = brick_list(nin,ib)%PCUT(2)%Vel(1:3)
2150
2151 !2. poly2 <- poly3
2152 brick_list(nin,ib)%POLY(2)%FACE(1:6)%Surf = -brick_list(nin,ib)%POLY(3)%FACE(1:6)%Surf
2153 brick_list(nin,ib)%PCUT(2)%SCUT = brick_list(nin,ib)%PCUT(3)%SCUT
2154 brick_list(nin,ib)%PCUT(2)%B(1:3) = brick_list(nin,ib)%PCUT(3)%B(1:3)
2155 brick_list(nin,ib)%PCUT(2)%N(1:3) = brick_list(nin,ib)%PCUT(3)%N(1:3)
2156 brick_list(nin,ib)%PCUT(2)%NP = brick_list(nin,ib)%PCUT(3)%NP
2157 brick_list(nin,ib)%PCUT(2)%P(1:3,1) = brick_list(nin,ib)%PCUT(3)%P(1:3,1)
2158 brick_list(nin,ib)%PCUT(2)%P(1:3,2) = brick_list(nin,ib)%PCUT(3)%P(1:3,2)
2159 brick_list(nin,ib)%PCUT(2)%P(1:3,3) = brick_list(nin,ib)%PCUT(3)%P(1:3,3)
2160 brick_list(nin,ib)%PCUT(2)%P(1:3,4) = brick_list(nin,ib)%PCUT(3)%P(1:3,4)
2161 brick_list(nin,ib)%PCUT(2)%P(1:3,5) = brick_list(nin,ib)%PCUT(3)%P(1:3,5)
2162 brick_list(nin,ib)%PCUT(2)%P(1:3,6) = brick_list(nin,ib)%PCUT(3)%P(1:3,6)
2163 brick_list(nin,ib)%POLY(2)%Vnew = -brick_list(nin,ib)%POLY(3)%Vnew
2164 brick_list(nin,ib)%POLY(2)%NumNOD = 8-brick_list(nin,ib)%POLY(3)%NumNOD
2165 brick_list(nin,ib)%POLY(2)%NumPOINT = brick_list(nin,ib)%POLY(1)%NumPOINT-brick_list(nin,ib)%POLY(1)%NumNOD
2166 . + 8-brick_list(nin,ib)%POLY(3)%NumNOD
2167 brick_list(nin,ib)%POLY(2)%ListNodID(1:8) = brick_list(nin,ib)%POLY(3)%ListNodID(1:8)
2168 brick_list(nin,ib)%PCUT(2)%Vel(1:3) = brick_list(nin,ib)%PCUT(3)%Vel(1:3)
2169
2170 !3.reset poly3.
2171 brick_list(nin,ib)%POLY(3)%FACE(1:6)%Surf = zero
2172 brick_list(nin,ib)%POLY(3)%Vnew = zero
2173 brick_list(nin,ib)%POLY(3)%NumNOD = 0
2174 brick_list(nin,ib)%POLY(3)%NumPOINT = 0
2175 brick_list(nin,ib)%POLY(3)%CellCenter(1) = zero
2176 brick_list(nin,ib)%POLY(3)%CellCenter(2) = zero
2177 brick_list(nin,ib)%POLY(3)%CellCenter(3) = zero
2178 brick_list(nin,ib)%POLY(3)%ListNodID(1:8) = 0
2179 DO k=0,6
2180 brick_list(nin,ib)%POLY(3)%FACE(k)%Vel(1) = zero
2181 brick_list(nin,ib)%POLY(3)%FACE(k)%Vel(2) = zero
2182 brick_list(nin,ib)%POLY(3)%FACE(k)%Vel(3) = zero
2183 ENDDO
2184
2185 DO i=1,8
2186 IF(brick_list(nin,ib)%NODE(i)%WhichCell == 2)THEN
2187
2188 ELSE
2189 brick_list(nin,ib)%NODE(i)%WhichCell = 1
2190 ENDIF
2191 ENDDO
2192
2193 !cell center update
2194 !from previous calculation based on intersection points + nodes, remove nodes and add new ones.
2195 brick_list(nin,ib)%POLY(1)%CellCenter(1:3) = brick_list(nin,ib)%POLY(1)%CellCenter(1:3) * numpt !(direct sum, and no longer an average)
2196 DO i=1,8
2197 nn => x(1:3, ixs(1+i,brick_list(nin,ib)%ID ))
2198 IF(brick_list(nin,ib)%NODE(i)%WhichCell == 2)THEN
2199 brick_list(nin,ib)%POLY(1)%CellCenter(1) = brick_list(nin,ib)%POLY(1)%CellCenter(1) - nn(1)
2200 brick_list(nin,ib)%POLY(1)%CellCenter(2) = brick_list(nin,ib)%POLY(1)%CellCenter(2) - nn(2)
2201 brick_list(nin,ib)%POLY(1)%CellCenter(3) = brick_list(nin,ib)%POLY(1)%CellCenter(3) - nn(3)
2202 numpt = numpt - 1
2203 ELSE
2204 brick_list(nin,ib)%POLY(1)%CellCenter(1) = brick_list(nin,ib)%POLY(1)%CellCenter(1) + nn(1)
2205 brick_list(nin,ib)%POLY(1)%CellCenter(2) = brick_list(nin,ib)%POLY(1)%CellCenter(2) + nn(2)
2206 brick_list(nin,ib)%POLY(1)%CellCenter(3) = brick_list(nin,ib)%POLY(1)%CellCenter(3) + nn(3)
2207 numpt = numpt + 1
2208 ENDIF
2209 ENDDO
2210 brick_list(nin,ib)%POLY(1)%CellCenter(1:3) = brick_list(nin,ib)%POLY(1)%CellCenter(1:3) / numpt !we get an average
2211
2212
2213
2214 !IF(Vol1>Vol2)THEN
2215 ELSE
2216 valtmp(1:3) = brick_list(nin,ib)%POLY(2)%CellCenter(1:3)
2217 brick_list(nin,ib)%POLY(2)%FACE(1:6)%Surf = -brick_list(nin,ib)%POLY(2)%FACE(1:6)%Surf
2218 brick_list(nin,ib)%POLY(2)%Vnew = -brick_list(nin,ib)%POLY(2)%Vnew
2219 brick_list(nin,ib)%POLY(2)%NumNOD = 8-brick_list(nin,ib)%POLY(2)%NumNOD
2220 brick_list(nin,ib)%POLY(2)%NumPOINT = brick_list(nin,ib)%POLY(2)%NumPOINT-brick_list(nin,ib)%POLY(1)%NumNOD
2221 . + 8-brick_list(nin,ib)%POLY(1)%NumNOD
2222 DO i=1,8
2223 IF(brick_list(nin,ib)%NODE(i)%WhichCell == 2)THEN
2224 brick_list(nin,ib)%NODE(i)%WhichCell = 1
2225 ELSE
2226 brick_list(nin,ib)%NODE(i)%WhichCell = 2
2227 ENDIF
2228 ENDDO
2229
2230 !cell center update
2231 !from previous calculation based on intersection points + nodes, remove nodes and add new ones.
2232 numpt = brick_list(nin,ib)%POLY(2)%NumPOINT
2233 brick_list(nin,ib)%POLY(2)%CellCenter(1:3) = brick_list(nin,ib)%POLY(2)%CellCenter(1:3) * numpt !(direct sum, and no longer an average)
2234 DO i=1,8
2235 IF(brick_list(nin,ib)%NODE(i)%WhichCell == 1)THEN
2236 nn => x(1:3, ixs(1+i,brick_list(nin,ib)%ID ))
2237 brick_list(nin,ib)%POLY(2)%CellCenter(1) = brick_list(nin,ib)%POLY(2)%CellCenter(1) - nn(1)
2238 brick_list(nin,ib)%POLY(2)%CellCenter(2) = brick_list(nin,ib)%POLY(2)%CellCenter(2) - nn(2)
2239 brick_list(nin,ib)%POLY(2)%CellCenter(3) = brick_list(nin,ib)%POLY(2)%CellCenter(3) - nn(3)
2240 numpt = numpt - 1
2241 ELSE
2242 nn => x(1:3, ixs(1+i,brick_list(nin,ib)%ID ))
2243 brick_list(nin,ib)%POLY(2)%CellCenter(1) = brick_list(nin,ib)%POLY(2)%CellCenter(1) + nn(1)
2244 brick_list(nin,ib)%POLY(2)%CellCenter(2) = brick_list(nin,ib)%POLY(2)%CellCenter(2) + nn(2)
2245 brick_list(nin,ib)%POLY(2)%CellCenter(3) = brick_list(nin,ib)%POLY(2)%CellCenter(3) + nn(3)
2246 numpt = numpt + 1
2247 ENDIF
2248 ENDDO
2249 brick_list(nin,ib)%POLY(2)%CellCenter(1:3) = brick_list(nin,ib)%POLY(2)%CellCenter(1:3) / numpt !we get an average
2250
2251 ENDIF !Vol1>Vol2
2252
2253
2254 !set new list nod (complementary from poly 1)
2255 inodtag(1:8) = 0
2256 DO i=1,brick_list(nin,ib)%POLY(1)%NumNOD
2257 inodtag( brick_list(nin,ib)%POLY(1)%ListNodID(i) ) = 1
2258 ENDDO
2259 j = 1
2260 DO i=1,8
2261 IF(inodtag(i) == 0)THEN
2262 brick_list(nin,ib)%POLY(2)%ListNodID(j) = i
2263 j = j + 1
2264 ENDIF
2265 ENDDO
2266
2267 brick_list(nin,ib)%POLY(9)%CellCenter(1:3) = valtmp(1:3)
2268 brick_list(nin,ib)%POLY(9)%NumNOD = 0
2269 brick_list(nin,ib)%POLY(9)%NumPOINT = 0
2270
2271 ENDIF !(Cod1==Cod2)
2272 ENDIF !IF(BRICK_LIST(NIN,IB)%NBCUT == 2)
2273 ENDDO !next IB
2274
2275
2276 !----------------------------------------------------------------
2277 ! save edge data in cut cell buffer for velocity interpolation.
2278 !----------------------------------------------------------------
2279 DO ib=nbf,nbl
2280 nbcut = brick_list(nin,ib)%NBCUT
2281 !IF(NBCUT==0)CYCLE !need EDGE_LIST(NIN,K)%NBCUT for velocity interpolation
2282 DO j=1,12
2283 k = (ib-1)*12+j
2284 brick_list(nin,ib)%Edge(j)%NODE(1:2) = edge_list(nin,k)%NODE(1:2)
2285 brick_list(nin,ib)%Edge(j)%NBCUT = edge_list(nin,k)%NBCUT
2286 brick_list(nin,ib)%Edge(j)%CUTSHELL(1:2) = edge_list(nin,k)%CUTSHELL(1:2)
2287 brick_list(nin,ib)%Edge(j)%CUTCOOR(1:2) = edge_list(nin,k)%CUTCOOR(1:2)
2288 brick_list(nin,ib)%Edge(j)%CUTPOINT(1:3,1) = edge_list(nin,k)%CUTPOINT(1:3,1)
2289 brick_list(nin,ib)%Edge(j)%CUTPOINT(1:3,2) = edge_list(nin,k)%CUTPOINT(1:3,2)
2290 brick_list(nin,ib)%Edge(j)%CUTVEL(1:3,1) = edge_list(nin,k)%CUTVEL(1:3,1)
2291 brick_list(nin,ib)%Edge(j)%CUTVEL(1:3,2) = edge_list(nin,k)%CUTVEL(1:3,2)
2292 brick_list(nin,ib)%Edge(j)%VECTOR(1:3) = edge_list(nin,k)%VECTOR(1:3)
2293 brick_list(nin,ib)%Edge(j)%LEN = edge_list(nin,k)%LEN
2294 enddo!next J
2295 enddo!next IB
2296
2297
2298 !----------------------------------------------------------------
2299 ! 3. ANIMATION (reset others)
2300 !----------------------------------------------------------------
2301 CALL my_barrier
2302
2303 IF (ii<=nel.AND.nspmd<=1.AND.itask==0) THEN
2304#include "lockon.inc"
2305 !initialization of grsh3n nodes to zero
2306 DO j=1,nel
2307 !EXIT
2308c print *, "idsh3n=", IXTG(NIXTG,IBUFSSG(J)) !IDSPRI=
2309 x(1:3,ixtg(2,igrsh3n(igr)%ENTITY(j)))=(/zero,zero,zero/)
2310 x(1:3,ixtg(3,igrsh3n(igr)%ENTITY(j)))=(/zero,zero,zero/)
2311 x(1:3,ixtg(4,igrsh3n(igr)%ENTITY(j)))=(/zero,zero,zero/)
2312 v(1:3,ixtg(2,igrsh3n(igr)%ENTITY(j)))=(/zero,zero,zero/)
2313 v(1:3,ixtg(3,igrsh3n(igr)%ENTITY(j)))=(/zero,zero,zero/)
2314 v(1:3,ixtg(4,igrsh3n(igr)%ENTITY(j)))=(/zero,zero,zero/)
2315 END DO
2316#include "lockoff.inc"
2317 END IF
2318 !----------------------------------------------------------------
2319
2320
2321
2322 !----------------------------------------------------------------
2323 !debug:
2324 if(itask==0)then
2325
2326 if (ibug22_subvol==1)then
2327 do ib=1,nb
2328 id = brick_list(nin,ib)%ID
2329 ncell = brick_list(nin,ib)%NBCUT
2330 write (*,fmt='(A,2I12)') "=== BRICK ID===", ixs(11,id) , id
2331 icell = 1
2332 do while (icell <= ncell .OR. icell == 9 ) !loop on ICELL
2333 mnod = brick_list(nin,ib)%POLY(icell)%NumNOD
2334 write (*,fmt='(A )') "|"
2335 if(icell==9)then
2336 write (*,fmt='(A,I1,A,A ,A1,I2,A1 )') "+== ICELL= ", icell , ", SecType=", 'COMPLEMENT' ,
2337 . "(", 0 , ")"
2338 else
2339 write (*,fmt='(A,I1,A,A14,A1,I2,A1 )') "+== ICELL= ", icell , ", SecType=", brick_list(nin,ib)%SECTYPE(icell) ,
2340 . "(", brick_list(nin,ib)%SecID_Cell(icell) , ")"
2341 endif
2342 pface(1:6) => brick_list(nin,ib)%POLY(icell)%FACE(1:6)!%Surf
2343 pface0=> brick_list(nin,ib)%POLY(icell)%FACE0%Surf
2344 write (*,fmt='(A )') "| |"
2345 write (*,fmt='(A,F20.14)') "| +======SUVOLUMES=", brick_list(nin,ib)%POLY(icell)%Vnew
2346 write (*,fmt='(A,6F20.14)') "| +=======SUBFACES=", brick_list(nin,ib)%POLY(icell)%FACE(1:6)%Surf
2347 write (*,fmt='(A,F20.14)') "| +=======CUT AERA=", brick_list(nin,ib)%POLY(icell)%FACE0%Surf
2348 write (*,fmt='(A,A,I2)') "| +======NUM POINT="," ",brick_list(nin,ib)%POLY(icell)%NumPOINT
2349 write (*,fmt='(A,A,I1,A,8I12)')"| +======NODE LIST="," (",mnod,") ", brick_list(nin,ib)%POLY(icell)%ListNodID(1:mnod)
2350 write (*,fmt='(A,A ,8I12)') "| radIDs=", " ",
2351 . ixs(1+brick_list(nin,ib)%POLY(icell)%ListNodID(1:mnod),id)
2352 write (*,fmt='(A,A ,8I12)') "| userIDs="," ",
2353 . itab(ixs(1+brick_list(nin,ib)%POLY(icell)%ListNodID(1:mnod),id))
2354! do i=0,6
2355! print *, "Sn=", vF(:,I)
2356! print *, " C=", C(1:3, Fa(I))
2357! end do
2358 icell = icell + 1
2359 if(icell>ncell .and. icell/=10) icell = 9 ! list is also : {1:ncell} U {9}
2360 enddo!next icell
2361 write (*,fmt='(A )') " "
2362
2363 end do!next I
2364
2365 endif!(IBUG22_subvol==1)
2366 end if
2367 !----------------------------------------------------------------
2368
2369
2370 RETURN
2371 END
2372
2373
2374!||====================================================================
2375!|| i22aera ../engine/source/interfaces/int22/i22subvol.F
2376!||--- called by ------------------------------------------------------
2377!|| i22for3 ../engine/source/interfaces/int22/i22for3.F
2378!|| i22subvol ../engine/source/interfaces/int22/i22subvol.f
2379!|| sinit22_fvm ../engine/source/interfaces/int22/sinit22_fvm.F
2380!||====================================================================
2381 FUNCTION i22aera(NPTS, P, C)
2382C=======================================================================
2383C-----------------------------------------------
2384C D e s c r i p t i o n
2385C-----------------------------------------------
2386C compute oriented areas of polygons p0, p1, .., pn with n=npts
2387C Precondition : NPTS>=3, & P0, P1, P3 non colineaires
2388C OPTIMISATION FOR LA SUITE : ne pas calculer les centroides ici car on duplique certaines coperations.
2389C-----------------------------------------------
2390C I m p l i c i t T y p e s
2391C-----------------------------------------------
2392#include "implicit_f.inc"
2393C-----------------------------------------------
2394C D u m m y A r g u m e n t s
2395C-----------------------------------------------
2396 my_real :: i22aera(3)
2397C-----------------------------------------------
2398 INTEGER, INTENT(IN) :: NPTS
2399 my_real, INTENT(IN) :: p(3,npts)
2400 my_real, INTENT(INOUT) :: c(3)
2401C-----------------------------------------------
2402C L o c a l V a r i a b l e s
2403C-----------------------------------------------
2404 INTEGER h, k,l,r,i,imax
2405 my_real
2406 . hh, diag13(3), diag24(3), s2n(3), norml2, s2,
2407 . v1(1:3), v2(1:3), v3(1:3), v4(1:3)
2408C-----------------------------------------------
2409
2410 IF(npts<=2)THEN
2411 !print *, "**error in inter22 : contact force algorithm"
2412 i22aera = zero
2413 stop 222
2414 return
2415 ELSEIF (npts<=4)THEN
2416 !produit vectoriel des diagonales diag(13) x diag(24)
2417 imax = max(npts,3)
2418 diag13(1:3) = p(1:3,3) - p(1:3,1)
2419 diag24(1:3) = p(1:3,imax) - p(1:3,2)
2420 s2n(1) = diag13(2)*diag24(3) - diag13(3)*diag24(2)
2421 s2n(2) = diag13(3)*diag24(1) - diag13(1)*diag24(3)
2422 s2n(3) = diag13(1)*diag24(2) - diag13(2)*diag24(1)
2423 s2n(1:3) = half * s2n(1:3)
2424 i22aera(1:3) = s2n(1:3)
2425 !I22AERA = SQRT(SUM(S2n(1:3)*S2n(1:3)))
2426 ELSEIF (npts>=5) THEN
2427 k=npts
2428 hh=half*(k-1)
2429 h=int(hh)
2430 IF( mod(k,2) == 1 )THEN
2431 l=0
2432 ELSE
2433 l=k-1
2434 END IF
2435 DO i=1, h-1
2436 v1(1:3) = p(1:3, 2*i+1) - p(1:3, 1)
2437 v2(1:3) = p(1:3,2*(i+1)) - p(1:3, 2*i)
2438 v3(1:3) = p(1:3, 2*h+1) - p(1:3, 1)
2439 v4(1:3) = p(1:3, l+1) - p(1:3, 2*h)
2440 s2n(1) = v1(2)*v2(3) - v1(3)*v2(2)
2441 s2n(2) = v1(3)*v2(1) - v1(1)*v2(3)
2442 s2n(3) = v1(1)*v2(2) - v1(2)*v2(1)
2443 s2n(1) = v3(2)*v4(3) - v3(3)*v4(2) + s2n(1)
2444 s2n(2) = v3(3)*v4(1) - v3(1)*v4(3) + s2n(2)
2445 s2n(3) = v3(1)*v4(2) - v3(2)*v4(1) + s2n(3)
2446 s2n(1:3) = half*s2n(1:3)
2447 i22aera = s2n(1:3)
2448 !I22AERA = SQRT(SUM(S2n(1:3)*S2n(1:3)))
2449 END DO
2450 END IF
2451
2452 c(1:3)=p(1:3,1)
2453 DO i=2,npts
2454 c(1:3)=c(1:3)+p(1:3,i)
2455 END DO
2456 c(1:3)=c(1:3)/npts
2457
2458 RETURN
2459 END
2460
my_real
#define my_real
Definition cppsort.cpp:32
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
i22subvol
subroutine i22subvol(x, nin, itask, ipari, itab, ixs, ixtg, v, iparg, elbuf_tab, w, igrsh3n)
Definition i22subvol.F:41
i22aera
function i22aera(npts, p, c)
Definition i22subvol.F:2382
max
#define max(a, b)
Definition macros.h:21
id
initmumps id
Definition multiplerhs_example.m:4
groupdef_mod
Definition groupdef_mod.F:668
i22bufbric_mod
Definition cut-cell-search_mod.F:95
i22bufbric_mod::edge_list
type(edge_entity), dimension(:,:), allocatable, target edge_list
Definition cut-cell-search_mod.F:106
i22bufbric_mod::brick_list
type(brick_entity), dimension(:,:), allocatable, target brick_list
Definition cut-cell-search_mod.F:105
i22bufbric_mod::int22_buf
type(int22_buf_t) int22_buf
Definition cut-cell-search_mod.F:130
i22edge_mod
Definition cut-cell-buffer_mod.F:48
i22tri_mod
Definition cut-cell-search_mod.F:212
i22tri_mod::nb
integer nb
Definition cut-cell-search_mod.F:487
i22tri_mod::ibug22_subvol
integer ibug22_subvol
Definition cut-cell-search_mod.F:246
my_barrier
subroutine my_barrier
Definition machine.F:31
i22edge_mod::poly_entity
Definition cut-cell-buffer_mod.F:152
i22edge_mod::polyface_entity
Definition cut-cell-buffer_mod.F:114