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