OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sinit22_fvm.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!|| sinit22_fvm ../engine/source/interfaces/int22/sinit22_fvm.F
25!||--- called by ------------------------------------------------------
26!|| alemain ../engine/source/ale/alemain.F
27!||--- calls -----------------------------------------------------
28!|| destroy_cell ../engine/source/interfaces/int22/destroy_cell.F
29!|| get_group_id ../engine/source/interfaces/int22/get_group_id.F
30!|| get_unique_main_cell ../engine/source/interfaces/int22/get_unique_master_cell.F
31!|| i22aera ../engine/source/interfaces/int22/i22subvol.F
32!|| initbuf ../engine/share/resol/initbuf.F
33!|| link_with_unique_main_cell ../engine/source/interfaces/int22/link_with_unique_master_cell.F
34!|| my_barrier ../engine/source/system/machine.F
35!|| sigeps37_single_cell ../engine/source/interfaces/int22/sigeps37_single_cell.f
36!|| weighting_cell_nodes ../engine/source/interfaces/int22/weighting_cell_nodes.F
37!||--- uses -----------------------------------------------------
38!|| ale_connectivity_mod ../common_source/modules/ale/ale_connectivity_mod.F
39!|| alefvm_mod ../common_source/modules/ale/alefvm_mod.f
40!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.f90
41!|| groupdef_mod ../common_source/modules/groupdef_mod.F
42!|| i22bufbric_mod ../common_source/modules/interfaces/cut-cell-search_mod.F
43!|| i22tri_mod ../common_source/modules/interfaces/cut-cell-search_mod.F
44!|| initbuf_mod ../engine/share/resol/initbuf.F
45!|| multimat_param_mod ../common_source/modules/multimat_param_mod.F90
46!||====================================================================
47 SUBROUTINE sinit22_fvm (
48 . IXS , ELBUF_TAB ,IPARG ,ITAB ,ITASK ,
49 . BUFBRIC, NBRIC_L ,X ,ALE_CONNECTIVITY,V ,
50 . NV46 , VEUL ,IGRNOD,IPARI ,IGRTRUSS,
51 . IXT , BUFMAT ,IPM)
52C-----------------------------------------------
53C D e s c r i p t i o n
54C-----------------------------------------------
55C Interface Type22 (/INTER/TYPE22) is an FSI coupling method based on cut cell method.
56C This experimental cut cell method is not completed, abandoned, and is not an official option.
57C
58
59! this subroutine is called for fsi interface 22
60! It handles intersected bricks when lagrangian
61! surface moved in this new cycle.
62!S T E P S
63 !------------------------------------------------------------!
64 ! @00. OUTPUT CUT CELL BUFFER LIST !
65 ! @01. CINEATIC TIME STEP !
66 ! @02. INIT/PARAMETERS !
67 ! @03. INIT/MULTITHREADING !
68 ! @04. INIT/ALLOCATION !
69 ! @05. RESET TAG22 !
70 ! @06. COMPUTE UNEXTENDED main CELL VOLUMES !
71 ! @07. BUILDING BIJECTION APPLICATION for cut cells !
72 ! @08. ADJACENT BRICKS DATA IN CUT CELL BUFFER !
73 ! @09. AREA and VOLUME RATIO UPDATE for cut cells !
74 ! @10. BUILD CELL 9 !
75 ! @11. BUILD UNCUT CELL IN CUT CELL BUFFER !
76 ! @12. STORING BRICK REFERENCE VOLUME (UNCUT) !
77 ! @13. POLY 9 : MARK IF NO NODE ON A FACE !
78 ! @14. FILL MATERIAL BUFFER : %bakMAT%... !
79 ! @15. FINDING & STORING ADJACENT CELLS !
80 ! @16. BUILD POLYGONAL TRACE FOR POLY9 !
81 ! @17. PRE-TREATMENT FOR DBLE SECONDARY WITH BOTH V=0 !
82 ! @18. TAGGING SECONDARY CUT CELL !
83 ! @19. LINKING SECONDARY CELL TO A main ONE !
84 ! @20. LIST ALL SECONDARY CELLS LINKED TO A GIVEN main CELL !
85 ! @21. INDIRECT LINKING !
86 ! @22. INDIRECT LINKING !
87 ! @23. TAG SECONDARY CELLS WHICH ARE BASED ON FORMER main NODE !
88 ! @24. VOLUME MERGING BETWEEN SECONDARY CELL & ITS main ONE !
89 ! @25. MATERIAL MERGING (IF main HAS TO BE DEPORTED) !
90 ! @26. DEMERGING CRITERIA !
91 ! @27. MATERIAL DEMERGING !
92 ! @28. COMPUTE VOLD_Vnew_Cell for newly cell in buffer !
93 ! @29. COMPUTE VOLUME CHANGE FOR EACH CELL !
94 ! @30. LINK SWITCHING !
95 ! @31. LEVEL2 NEIGHBORHOOD (FOR FLUXES & CONVECTION) !
96 ! @32. COUNT POINTS ON A GIVEN POLYHEDRON FACE !
97 ! @33. LOCATE WHERE WAS main FOR NEXT CYCLE !
98 ! @34. SET IIAD22(1:NUMELS) = %ELBUG_TAB(NG)%TAG22(I) !
99 ! @35. SET MLW !
100 ! @36. RESET OLD SECONDARY LIST CONNECTIVITY !
101 ! @37. SET main STRONG NODE !
102 ! @38. DRAW SUPERCELL CONNEXION WITH 1D TRUSSES !
103 ! @39. COMPUTE DVOL PREDICTION FOR EACH CELL !
104 ! @40. BUILD SUPERCELL DVOL (PREDICTION) !
105 ! @41. FIX SUPERCELL DVOL (CORRECTION) !
106 ! @42. TIME ZERO VOLUME INIT (NOT DONE IN STARTER) !
107 ! @43. DEBUG OUTPUT !
108 ! @44. DEBUG - WRITE CUT CELL BUFFER !
109 ! @45. SUPERCELL CENTERS !
110 ! @46. MARK SUPER-CELL CENTERS WITH ORPHAN NODES !
111 ! @47. UNCUT CELLS + POLY 9 : CENTERS !
112 ! @48. MARK CELL CENTERS WITH ORPHAN NODES !
113 ! @49. MARK FACE CENTERS WITH ORPHAN NODES !
114 ! @50. DEALLOCATE !
115 !------------------------------------------------------------!
116
117!
118! BUFFERS :
119! BRICK_LIST : BUT CELL BUFFER(NIN,IB) : NIN interface id, IB:cut cell id
120C-----------------------------------------------
121C N o t e s F o r L a t e r
122C-----------------------------------------------
123C -1- Reset tag22 only for brick in previous buffer (check performances)
124C -2- remove pointers
125C -3- mettre les allocatables dans un module
126C-----------------------------------------------
127C M o d u l e s
128C-----------------------------------------------
129 USE initbuf_mod
130 USE i22bufbric_mod
131 USE i22tri_mod
132 USE elbufdef_mod
133 USE alefvm_mod
134 USE groupdef_mod
135 USE ale_connectivity_mod
136 USE multimat_param_mod , ONLY : m51_n0phas, m51_nvphas
137C-----------------------------------------------
138C I m p l i c i t T y p e s
139C-----------------------------------------------
140#include "implicit_f.inc"
141C-----------------------------------------------
142C C o m m o n B l o c k s
143C-----------------------------------------------
144#include "param_c.inc"
145#include "com01_c.inc"
146#include "com04_c.inc"
147#include "com08_c.inc"
148#include "task_c.inc"
149#include "vect01_c.inc"
150#include "inter22.inc"
151#include "mvsiz_p.inc"
152#include "comlock.inc"
153#include "subvolumes.inc"
154C-----------------------------------------------
155C D u m m y A r g u m e n t s
156C-----------------------------------------------
157 INTEGER,INTENT(IN) :: IXS(NIXS,*) ,IPARG(NPARG,*),ITAB(*) ,NV46 ,BUFBRIC(*),
158 . IPARI(*) ,IXT(NIXT,*) ,ITASK ,IPM(NPROPMI,*)
159 TYPE(ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
160 my_real,INTENT(IN) :: V(3,*), VEUL(LVEUL,*)
161 my_real,INTENT(IN),TARGET :: bufmat(*)
162 my_real,INTENT(INOUT) :: x(3,*)
163!
164 TYPE (GROUP_) , DIMENSION(NGRNOD) :: IGRNOD
165 TYPE (group_) , DIMENSION(NGRTRUS) :: igrtruss
166 TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECTIVITY
167C-----------------------------------------------
168C L o c a l V a r i a b l e s
169C-----------------------------------------------
170 TYPE(l_bufel_) ,POINTER :: LBUF1,LBUF2
171 INTEGER :: I,J,K0,K1,K2,JV,IDV, NBCUT, NBCUTv, NBCUTprev, NEL,NG,NFL,NBF,NBL,NBL1,NBRIC_L,NIN
172 INTEGER :: brickID,tNB,NTAG,IV,NGv,IAD22,NCELLv,ICV,IGR
173 INTEGER :: IB,IBV,IBv2,IBv_i,IBo,ICELL,ICELL2,MCELL,NCELL,MNOD,ID,ITAG(66)
174 INTEGER :: IPOS, LLT_, LLT_o,LLT_v, IDLOCv, IPOSf, IPOSiadj,ICELLv,ICELLv2
175 INTEGER :: INODES(8),INOD,INOD2,INODE,NNODES, NNODES2, ADD, ADD0, ITRIMAT, K, KV, Ko
176 LOGICAL :: lDONE,lStillTruss,lFOUND,lCYCLE,lTARGET,lStillNode, lCOND1, lCOND2
177 my_real :: volcell, volcell51(trimat), volsecnd51(trimat),volsecnd,volcell51_old(trimat)
178 my_real :: adjmain_vol(6), adjmain_face(6), adjmain_centroid(3,6),n(3,6), fac
179 INTEGER :: IDadj_MCELLv(6), IVadj_MCELLv(6), IBadj_MCELLv(6)
180 my_real,DIMENSION(:,:),ALLOCATABLE :: f
181 my_real, POINTER :: pvar, pvarv, pvaro
182 my_real, DIMENSION(:), POINTER :: pvar3
183 TYPE(g_bufel_) ,POINTER :: GBUF, GBUFv, GBUFo
184 TYPE(l_bufel_) ,POINTER :: LBUF
185 TYPE(poly_entity),DIMENSION(:), POINTER :: pIsMain ,pIsMainV
186 my_real,DIMENSION(:) , POINTER :: pfullface
187 TYPE(node_entity),DIMENSION(:) , POINTER :: pNodWasMain
188 !INTEGER,DIMENSION(:,:), POINTER :: pNAdjCELL
189 INTEGER,DIMENSION(:,:), POINTER :: pAdjBRICK, pAdjBRICKv
190 TYPE(node_entity),DIMENSION(:) , POINTER :: pWhereWasMain
191 TYPE(node_entity), DIMENSION(:), POINTER :: pWhichCellNod,pWhichCellNodv
192 INTEGER , POINTER :: pMainID
193 my_real,DIMENSION(:) , POINTER :: uparam !=> BRICK_ENTITY%POLY()%Vnew
194 TYPE(POLY_ENTITY), DIMENSION(:), POINTER :: pSUBVOLv
195 my_real , POINTER :: psubvold
196 TYPE(poly_entity), DIMENSION(:), POINTER :: pSUBVOL
197 CHARACTER*6 :: char1
198 LOGICAL :: bool1, bool2
199 CHARACTER*10 :: debugMAINSECND
200 CHARACTER*10 ,ALLOCATABLE :: debugMAINSECNDv(:,:,:)
201 INTEGER,ALLOCATABLE :: IsMainV(:,:,:)
202 TYPE(buf_mat_) ,POINTER :: MBUF,MBUFv, MBUFo
203 INTEGER :: ICRIT_MAT_DEPORTATION, ICRIT_DEMERGE
204! INTEGER,ALLOCATABLE,DIMENSION(:,:,:):: TARGET_DATA !index1:facet - index2:adjacent_cell_number - index3:ib(cut_cell_id)
205 INTEGER,ALLOCATABLE,DIMENSION(:,:,:):: ORIGIN_DATA !index1:facet - index2:adjacent_cell_number - index3:ib(cut_cell_id),J
206 INTEGER,ALLOCATABLE,DIMENSION(:) :: Ntarget,Norigin
207 INTEGER :: MTN_,ITAR, IORIG, MAIN_TAG(6,9) !index1:Face !index2:icv=cell number
208 INTEGER :: IPLA, ISOLNOD, FM, IBm,IBmCur,IBMo,IBMold, IDM, IE, IN, MT
209 INTEGER :: ICELLTAG(9), ICELLTAG_OLD(9),IC, NAdjCELL,NAdjCELLv, IFV, FV, FV2
210 INTEGER :: IVv,MCELLv, MCELLvv, NGvv,IDLOCvv, IBvv,IFVv
211 INTEGER :: GET_UNIQUE_MAIN_CELL, LINK_WITH_UNIQUE_MAIN_CELL,IADJ,NADJ,NADJv,IE_M
212 INTEGER :: NsecndNOD, NP, NP_NODE, JJ, NINTP(6,9), NNOD(6,9), NN, SECid, IDLOC, NewMnod(8), MLW
213 INTEGER :: ITASKv, NTRUS, NPOINT, IRES(2), NewInBuffer
214 INTEGER,ALLOCATABLE,DIMENSION(:,:) :: DESTROY_DATA
215 INTEGER :: I22LOOP,IT,IUNLINK,ITASKB, J1,J2, NTAR, INod_W, IDEMERGE, INod_W_old, ISECND, IAD0,II
216 INTEGER :: FV_old, NumSECND, NumSECND2, IC1, IC2, ITARGET, NumTarget, WasCut, LID, IFAC, IEDG
217 INTEGER :: Cod1, Cod2, Cod3, Icompl, Poly9woNodes, Poly9woNodesV, ISGN(2), ICODE, OLD_ICODE
218 INTEGER :: mainID,NC(8),MOLD,MNEW, IB2,NUM,IADBUF,NUPARAM,NGm,IDLOCm,ICELLm, NP_(9),NODE_ID
219 LOGICAL :: debug_outp,debug_outp2, lSKIP
220
221 my_real,ALLOCATABLE, DIMENSION(:,:) :: vol51, vol51v !index1:IB index2:ITRIMAT
222 my_real,ALLOCATABLE, DIMENSION(:) :: uvar,uvar_adv !(M51_N0PHAS+TRIMAT*M51_NVPHAS+I22LAW37) !NUVAR=M51_N0PHAS+TRIMAT*M51_NVPHAS
223 my_real :: vfrac(4),som, som_(4), somi, adjface(nv46,nadj_f), delta, ppoint(3)
224 my_real :: sumface, vectmp(3), volorig(24), pointtmp(3), point0(3), cut_point(3)
225 my_real :: cut_vel(3)
226 my_real :: eint, eintv, rho,rho_u(3),mom(3),rhov,sigv(6), sig(6),volv,vol,vol_m,vol_s
227 my_real :: var, var_,var__, var_6(6), var_vf(4), vold_phase(4)
228 my_real :: tmp(6),dxmin,vmax,ratiof,ratiofv,ratio,ratiov,uncutvol,uncutvolv
229 my_real :: vuncut, vi,vj, m(9,9) !index1:polyhedra in current conformation, index2:polyhedra in previous confromation
230 my_real :: vsum(3) , n_(3), vnew, vold, dvol_numeric, dvol_predic
231 my_real :: sgn, dvi, dvii, face, face9, m_tot, m_liq, m_toto,m_liqo, rho10, rho20, mfrac
232 my_real :: center(3,8)
233
234 my_real :: det, df11, df12, df21, df22, f1, f2, p1, p2, drho1, drho2, error, tol, r1, psh, pmin
235 my_real :: c1, gam, p0, p, rho1, rho2, mas, mas1, mas2, ssp, ssp1, ssp2, rhoc2_1, rhoc2_2, rhoc2
236 my_real :: alp, alpo, beta, volcellold, vel_m(3), surf_s, norm_s(3), adv, norm, vm, vs
237 my_real :: rho_adv, eint_adv, sig_adv(6), mom_adv(3), zm(3),zs(3), volratio
238 INTEGER :: ITER, NITER, LLTo, LLTm, NGo, IADV
239 INTEGER, ALLOCATABLE, DIMENSION(:) :: ORDER, VALUE
240
241 INTEGER :: NC1,NC2,NC3,NC4,NC5,NC6,NC7,NC8
242 INTEGER :: IAD2, IAD3, LGTH2, LGTH3
243
244 INTERFACE
245 FUNCTION i22aera(NPTS, P, C)
246 my_real :: i22aera(3)
247 INTEGER, INTENT(IN) :: NPTS
248 my_real, INTENT(IN) :: p(3,npts)
249 my_real, INTENT(INOUT) :: c(3)
250 END FUNCTION
251 END INTERFACE
252
253 TYPE pointer_array_r
254 my_real,DIMENSION(:), POINTER :: p
255 END TYPE
256
257 TYPE pointer_array_i
258 INTEGER,DIMENSION(:), POINTER :: p
259 END TYPE
260
261 TYPE(poly_entity), DIMENSION(:), POINTER :: pFACE ! WARNING : all pointers begin with 1 !!! index 0 impossible pFACE(1:7,:)=>%Face_Cell(0:6,:) <
262 TYPE(pointer_array_r) :: pFACEv(9)
263 TYPE(pointer_array_i) :: pListNodID(9)
264 TYPE(POINTER_ARRAY_I) :: pListNodIDv(9)
265 !TYPE(POINTER_ARRAY_I) :: pAdjCELL
266
267 INTEGER ICF(4,6)
268 DATA icf/1,4,3,2,3,4,8,7,5,6,7,8,1,2,6,5,2,3,7,6,1,5,8,4/
269
270 my_real ::
271 . aj7(mvsiz), aj8(mvsiz) , aj9(mvsiz),
272 . aji1, aji2, aji3,
273 . aji4, aji5, aji6,
274 . aji7, aji8, aji9,
275 . x17 , x28 , x35 , x46,
276 . y17 , y28 , y35 , y46,
277 . z17 , z28 , z35 , z46,
278 . jac_59_68(mvsiz), jac_67_49(mvsiz), jac_48_57(mvsiz),
279 . aj12, aj45, aj78,
280 . a17 , a28 ,
281 . b17 , b28 ,
282 . c17 , c28 ,
283 . x1(mvsiz),x2(mvsiz),x3(mvsiz),x4(mvsiz),x5(mvsiz),x6(mvsiz),x7(mvsiz),x8(mvsiz),
284 . y1(mvsiz),y2(mvsiz),y3(mvsiz),y4(mvsiz),y5(mvsiz),y6(mvsiz),y7(mvsiz),y8(mvsiz),
285 . z1(mvsiz),z2(mvsiz),z3(mvsiz),z4(mvsiz),z5(mvsiz),z6(mvsiz),z7(mvsiz),z8(mvsiz),
286 . dett(mvsiz),
287 . aj1(mvsiz),aj2(mvsiz),aj3(mvsiz),
288 . aj4(mvsiz),aj5(mvsiz),aj6(mvsiz),hxp,hyp,hzp,
289 . x1_,x2_,x3_,x4_,x5_,x6_,x7_,x8_,
290 . y1_,y2_,y3_,y4_,y5_,y6_,y7_,y8_,
291 . z1_,z2_,z3_,z4_,z5_,z6_,z7_,z8_
292
293C-----------------------------------------------
294C P r e c o n d i t i o n s
295C-----------------------------------------------
296 !IF(INT22 == 0) RETURN !already checked if this subroutine is called
297C-----------------------------------------------
298C S o u r c e L i n e s
299C-----------------------------------------------
300 !------------------------------------------------------------!
301 ! @00. OUTPUT CUT CELL BUFFER LIST !
302 !------------------------------------------------------------!
303 IF(ibug22_ccbuflist /= 0)THEN
304 nin = 1
305 IF(itask==0)THEN
306 write (*,fmt='(A, 1000I7)')"CUT CELL BUFFER : ", ixs(11,brick_list(nin,1:nb)%id)
307 ENDIF
308 ENDIF
309
310 !------------------------------------------------------------!
311 ! @01. CINEATIC TIME STEP !
312 !------------------------------------------------------------!
313 idt_int22 = 1 !kinematic time step by default. Will be replaced by 0 if exist i such as DTIX(I) < DT22_MIN (mqviscb.F)
314 dxmin = minval(dx22min_l(0:nthread-1))
315 vmax = maxval(v22max_l(0:nthread-1))
316 IF(vmax == zero)THEN
317 dt22_min = ep30
318 ELSE
319 dt22_min = dxmin/ncross22 / vmax
320 ENDIF
321 if(ibug22_dtmin==1)print *, "dt22_min = ", dt22_min
322.AND. !!sIF(DT1==ZERO ITASK==0) print *, "inter22, tolerance ratio ",RATIO22
323
324 !------------------------------------------------------------!
325 ! @02. INIT/PARAMETERS !
326 !------------------------------------------------------------!
327 IADV = 0 ! 0: advection step off 1: advection step on !When demerging a new main, advection is done from old main cell (possibly relocated due to a merging process) to new one.
328 NIN = 1
329 dbVOL = ZERO
330 dbMASS = ZERO
331 I22LOOP = 0 !infinite loop detected
332 IMERGEL(ITASK) = 0
333 IF(TRIMAT<=0)THEN
334 ALLOCATE(UVAR (I22LAW37))
335 ALLOCATE(UVAR_ADV (I22LAW37))
336 ELSE
337 ALLOCATE(UVAR (M51_N0PHAS+TRIMAT*M51_NVPHAS))
338 ALLOCATE(UVAR_ADV (M51_N0PHAS+TRIMAT*M51_NVPHAS))
339 ENDIF
340
341 !------------------------------------------------------------!
342 ! @03. INIT/MULTITHREADING !
343 !------------------------------------------------------------!
344 NBF = 1+ITASK*NB/NTHREAD
345 NBL = (ITASK+1)*NB/NTHREAD
346 NBL = MIN(NBL,NB)
347 tNB = NBL-NBF+1
348 I22_DEGENERATED = 0
349
350 !!!!!!!!!!!!!!!!!!!!!!DEBUG OUTPUT!!!!!!!!!!!!!!!!!!!!!!!
351 debug_outp = .false.
352 debug_outp2= .false.
353
354 do ib=nbf,nbl
355 ie = brick_list(nin,ib)%id
356.or. if (ie == ibug22_sinit ibug22_sinit == -1)then
357 debug_outp = .true.
358 exit
359 endif
360 enddo
361
362.AND. if(itask==0debug_outp)then
363 print *, ""
364 print *, " |----------i22sinit.f-----------|"
365 print *, " | initialization SUBROUTINE |"
366 print *, " |-------------------------------|"
367 Char1 = ' '
368 end if
369 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
370
371 !------------------------------------------------------------!
372 ! @04. INIT/ALLOCATION !
373 !------------------------------------------------------------!
374 ALLOCATE (debugMAINSECNDv(NBL-NBF+1,6,9))
375! ALLOCATE (TARGET_DATA(tNB,9,1:2)) !index1:IB, index2:1:Ntarget<=9, index3:1=localface_2=adjacent_cell
376 ALLOCATE (ORIGIN_DATA(NBF:NBL,9,1:3)) !index1:IB, index2:1:Ntarget<=9, index3:1=localface_2=adjacent_cell
377 ALLOCATE (IsMainV(NBF:NBL,6,9))
378 ALLOCATE (F(6,NBF:NBL))
379 ALLOCATE (VOL51(NBF:NBL,TRIMAT),VOL51v(NBF:NBL,TRIMAT))
380c ALLOCATE (Ntarget(NBF:NBL))
381 ALLOCATE (Norigin(NBF:NBL))
382 !-------------------!
383 CALL MY_BARRIER()
384 !-------------------!
385 !------------------------------------------------------------!
386 ! @05. RESET TAG22 !
387 !------------------------------------------------------------!
388 !TAG22(I) is adress in cut cell buffer "brick_list(NIN,TAG22(I))"
389 !for interface number NIN, and element I from current group.
390 !TAG22==0 means elem I is not in cut cell buffer.
391 DO NG=ITASK+1,NGROUP,NTHREAD
392 IF(IPARG(8,NG) /= 1) THEN
393 CALL INITBUF(
394 1 IPARG ,NG ,
395 2 MTN ,NEL ,NFT ,IAD ,ITY ,
396 3 NPT ,JALE ,ISMSTR ,JEUL ,JTUR ,
397 4 JTHE ,JLAG ,JMULT ,JHBE ,JIVF ,
398 5 NVAUX ,JPOR ,JCVT ,JCLOSE ,IPLA ,
399 6 IREP ,IINT ,IGTYP ,ISRAT ,ISROT ,
400 7 ICSEN ,ISORTH ,ISORTHG ,IFAILURE,JSMS )
401.AND. IF(JLAG /= 1 ITY<=2) THEN
402.AND. IF (MTN /= 0 IPARG(64,NG)==0) THEN
403 LFT = 1
404 LLT = NEL
405 ISOLNOD = IPARG(28,NG)
406.AND. IF (ITY == 1 ISOLNOD /= 4) THEN
407 GBUF => ELBUF_TAB(NG)%GBUF
408 GBUF%TAG22(LFT:LLT) = 0
409 ENDIF
410 ENDIF
411 ENDIF
412 ENDIF
413 ENDDO
414 !------------------------------------------------------------!
415 ! @06. COMPUTE UNEXTENDED main CELL VOLUMES !
416 !------------------------------------------------------------!
417 DO NG=ITASK+1,NGROUP,NTHREAD
418 IF(IPARG(8,NG) /= 1) THEN
419 CALL INITBUF(
420 1 IPARG ,NG ,
421 2 MTN ,NEL ,NFT ,IAD ,ITY ,
422 3 NPT ,JALE ,ISMSTR ,JEUL ,JTUR ,
423 4 JTHE ,JLAG ,JMULT ,JHBE ,JIVF ,
424 5 NVAUX ,JPOR ,JCVT ,JCLOSE ,IPLA ,
425 6 IREP ,IINT ,IGTYP ,ISRAT ,ISROT ,
426 7 ICSEN ,ISORTH ,ISORTHG ,IFAILURE,JSMS )
427.AND. IF(JLAG /= 1 ITY<=2) THEN
428.AND. IF (MTN /= 0 IPARG(64,NG)==0) THEN
429 LFT = 1
430 LLT = NEL
431 ISOLNOD = IPARG(28,NG)
432.AND. IF (ITY == 1 ISOLNOD /= 4) THEN
433 GBUF => ELBUF_TAB(NG)%GBUF
434 IF(JEUL/=0)THEN
435 IF(INTEG8==0)THEN
436 GBUF%VOL(LFT:LLT) = VEUL(32,NFT+LFT:NFT+LLT)
437 ELSEIF(INTEG8==1)THEN
438 GBUF%VOL(LFT:LLT) = VEUL(52,NFT+LFT:NFT+LLT)
439 ENDIF!INTEG8
440 ELSE
441 DO I=LFT,LLT
442 II = I+NFT
443 X1(I) = X(1,IXS(2,II)) ; Y1(I) = X(2,IXS(2,II)) ; Z1(I) = X(3,IXS(2,II)) ;
444 X2(I) = X(1,IXS(3,II)) ; Y2(I) = X(2,IXS(3,II)) ; Z2(I) = X(3,IXS(3,II)) ;
445 X3(I) = X(1,IXS(4,II)) ; Y3(I) = X(2,IXS(4,II)) ; Z3(I) = X(3,IXS(4,II)) ;
446 X4(I) = X(1,IXS(5,II)) ; Y4(I) = X(2,IXS(5,II)) ; Z4(I) = X(3,IXS(5,II)) ;
447 X5(I) = X(1,IXS(6,II)) ; Y5(I) = X(2,IXS(6,II)) ; Z5(I) = X(3,IXS(6,II)) ;
448 X6(I) = X(1,IXS(7,II)) ; Y6(I) = X(2,IXS(7,II)) ; Z6(I) = X(3,IXS(7,II)) ;
449 X7(I) = X(1,IXS(8,II)) ; Y7(I) = X(2,IXS(8,II)) ; Z7(I) = X(3,IXS(8,II)) ;
450 X8(I) = X(1,IXS(9,II)) ; Y8(I) = X(2,IXS(9,II)) ; Z8(I) = X(3,IXS(9,II)) ;
451 ENDDO
452 DO I=LFT,LLT
453 X17=X7(I)-X1(I)
454 X28=X8(I)-X2(I)
455 X35=X5(I)-X3(I)
456 X46=X6(I)-X4(I)
457 Y17=Y7(I)-Y1(I)
458 Y28=Y8(I)-Y2(I)
459 Y35=Y5(I)-Y3(I)
460 Y46=Y6(I)-Y4(I)
461 Z17=Z7(I)-Z1(I)
462 Z28=Z8(I)-Z2(I)
463 Z35=Z5(I)-Z3(I)
464 Z46=Z6(I)-Z4(I)
465 AJ1(I)=X17+X28-X35-X46
466 AJ2(I)=Y17+Y28-Y35-Y46
467 AJ3(I)=Z17+Z28-Z35-Z46
468 A17=X17+X46
469 A28=X28+X35
470 B17=Y17+Y46
471 B28=Y28+Y35
472 C17=Z17+Z46
473 C28=Z28+Z35
474 AJ4(I)=A17+A28
475 AJ5(I)=B17+B28
476 AJ6(I)=C17+C28
477 AJ7(I)=A17-A28
478 AJ8(I)=B17-B28
479 AJ9(I)=C17-C28
480 ENDDO
481 DO I=LFT,LLT
482 JAC_59_68(I)=AJ5(I)*AJ9(I)-AJ6(I)*AJ8(I)
483 JAC_67_49(I)=AJ6(I)*AJ7(I)-AJ4(I)*AJ9(I)
484 JAC_48_57(I)=AJ4(I)*AJ8(I)-AJ5(I)*AJ7(I)
485 ENDDO
486 DO I=LFT,LLT
487 DETT(I)=ONE_OVER_64*(AJ1(I)*JAC_59_68(I)+AJ2(I)*JAC_67_49(I)+AJ3(I)*JAC_48_57(I))
488 ENDDO
489 GBUF%VOL(LFT:LLT) = DETT(LFT:LLT)
490 ENDIF!JEUL/=0
491 ENDIF
492 ENDIF
493 ENDIF
494 ENDIF
495 ENDDO
496 !------------------------------------------------------------!
497 ! @07. BUILDING BIJECTION APPLICATION for cut cells !
498 ! (local array) <--> (global buffer) !
499 ! NBRIC <--> (I,NG) !
500 !------------------------------------------------------------!
501 !-------------------!
502 CALL MY_BARRIER() !need %TAG22 to be init to 0. (initialisation is multithreading segmenting 1:NG, not 1:NB
503 !-------------------!
504 DO IB=NBF,NBL !1,NB
505 IDB(IB) = BRICK_LIST(NIN,IB)%ID
506 BRICK_LIST(NIN,IB)%ITASK = ITASK
507 ENDDO
508 DO IB=NBF,NBL !1,NB
509 !recuperer NG,I
510 !From IB get I and NG
511 lDONE = .FALSE.
512 DO NG=1,NGROUP
513 NEL=IPARG(2,NG)
514 NFL=IPARG(3,NG)
515.AND. IF((IDB(IB)>NFL)(IDB(IB)<=NFL+NEL))THEN
516 IF(IPARG(11,NG)==1) JEUL = 1
517 IF(IPARG(7,NG)==1) JALE = 1
518 IDLOCB(IB) = IDB(IB) - NFL
519 NGB(IB) = NG
520 BRICK_LIST(NIN,IB)%NG = NGB(IB)
521 BRICK_LIST(NIN,IB)%IDLOC = IDLOCB(IB)
522 NELB(IB) = NEL
523 lDONE =.TRUE.
524 GBUF => ELBUF_TAB(NGB(IB))%GBUF
525 NBCUT = BRICK_LIST(NIN,IB)%NBCUT
526 GBUF%TAG22(IDLOCB(IB)) = IB ! les tags ont deja ete remis a zero
527 EXIT !next I
528 ELSE
529 CYCLE !next NG
530 ENDIF
531.NOT. if ((lDONE))then
532 write( *,*) "int 22 : error in group sorting"
533 stop 2201
534 end if
535 ENDDO !next NG
536 ENDDO !next i
537 !-------------------!
538 CALL MY_BARRIER()
539 !-------------------!
540 !------------------------------------------------------------!
541 ! @08. ADJACENT BRICKS DATA IN CUT CELL BUFFER !
542 !------------------------------------------------------------!
543 DO IB=NBF,NBL
544 brickID = IDB(IB)
545 IAD2 = ALE_CONNECTIVITY%ee_connect%iad_connect(brickID)
546 LGTH2 = ALE_CONNECTIVITY%ee_connect%iad_connect(brickID+1) -
547 . ALE_CONNECTIVITY%ee_connect%iad_connect(brickID)
548 DO J=1,LGTH2
549 IDV = ALE_CONNECTIVITY%ee_connect%connected(IAD2 + J - 1)
550 BRICK_LIST(NIN,IB)%Adjacent_Brick(J,1) = IDV !IV
551 BRICK_LIST(NIN,IB)%Adjacent_Brick(J,2) = 0 !NGV
552 BRICK_LIST(NIN,IB)%Adjacent_Brick(J,3) = 0 !IDLOCV
553 BRICK_LIST(NIN,IB)%Adjacent_Brick(J,4) = 0 !IBV
554 BRICK_LIST(NIN,IB)%Adjacent_Brick(J,5) = 0 !IFV
555 IF(IDV>0)THEN
556 IAD3 = ALE_CONNECTIVITY%ee_connect%iad_connect(IDV)
557 LGTH3 = ALE_CONNECTIVITY%ee_connect%iad_connect(IDV+1) -
558 . ALE_CONNECTIVITY%ee_connect%iad_connect(IDV)
559 DO JV=1,LGTH3
560 IF(ALE_CONNECTIVITY%ee_connect%connected(IAD3 + JV - 1)==brickID)THEN
561 BRICK_LIST(NIN,IB)%Adjacent_Brick(J,5) = JV
562 EXIT
563 ENDIF
564 ENDDO!next JV
565 CALL GET_GROUP_ID(IDV,NGv,IDLOCv,IPARG)
566 GBUFv => ELBUF_TAB(NGv)%GBUF
567 IAD22 = GBUFv%TAG22(IDLOCv)
568 BRICK_LIST(NIN,IB)%Adjacent_Brick(J,2) = NGv
569 BRICK_LIST(NIN,IB)%Adjacent_Brick(J,3) = IDLOCv
570 BRICK_LIST(NIN,IB)%Adjacent_Brick(J,4) = IAD22
571.AND. IF (IAD22==0 BRICK_LIST(NIN,IB)%NBCUT>0)THEN
572 print *, "**error : inter22 : lagrangian surface seems to"
573 print *, " reach eulerian boundary domain. "
574 print *, " check surface location and grbrick definition"
575 print *, " near related brick_id =", IXS(11,brick_list(nin,ib)%id)
576 !print *, "itask=", ITASK
577
578 stop 2202 !means that a real cut cell has an adjacent brick which is not in cut cell buffer : NOT EXPECTED unless if /INTER/TYPE22 GrBrick is not correctly defined.
579 ENDIF
580 ELSEIF(IDV==0)THEN
581 BRICK_LIST(NIN,IB)%Adjacent_Brick(J,2) = 0
582 BRICK_LIST(NIN,IB)%Adjacent_Brick(J,3) = 0
583 BRICK_LIST(NIN,IB)%Adjacent_Brick(J,4) = 0
584 ELSEIF(IDV<0)THEN
585 write (*,*) "unavailable in spmd"
586 stop 2203
587 ENDIF
588 ENDDO
589 ENDDO
590
591 !------------------------------------------------------------!
592 ! @09. AREA and VOLUME RATIO UPDATE for cut cells !
593 !------------------------------------------------------------!
594 DO IB=NBF,NBL
595 brickID = IDB(IB)
596 IF(I22_ALEUL==2)THEN
597 N(1:3,1) = (/ VEUL(14,brickID) , VEUL(20,brickID) , VEUL(26,brickID) /)
598 N(1:3,2) = (/ VEUL(15,brickID) , VEUL(21,brickID) , VEUL(27,brickID) /)
599 N(1:3,3) = (/ VEUL(16,brickID) , VEUL(22,brickID) , VEUL(28,brickID) /)
600 N(1:3,4) = (/ VEUL(17,brickID) , VEUL(23,brickID) , VEUL(29,brickID) /)
601 N(1:3,5) = (/ VEUL(18,brickID) , VEUL(24,brickID) , VEUL(30,brickID) /)
602 N(1:3,6) = (/ VEUL(19,brickID) , VEUL(25,brickID) , VEUL(31,brickID) /)
603 BRICK_LIST(NIN,IB)%N(1,1:3) = N(1:3,1) / SQRT(SUM(N(1:3,1)*N(1:3,1)))
604 BRICK_LIST(NIN,IB)%N(2,1:3) = N(1:3,2) / SQRT(SUM(N(1:3,2)*N(1:3,2)))
605 BRICK_LIST(NIN,IB)%N(3,1:3) = N(1:3,3) / SQRT(SUM(N(1:3,3)*N(1:3,3)))
606 BRICK_LIST(NIN,IB)%N(4,1:3) = N(1:3,4) / SQRT(SUM(N(1:3,4)*N(1:3,4)))
607 BRICK_LIST(NIN,IB)%N(5,1:3) = N(1:3,5) / SQRT(SUM(N(1:3,5)*N(1:3,5)))
608 BRICK_LIST(NIN,IB)%N(6,1:3) = N(1:3,6) / SQRT(SUM(N(1:3,6)*N(1:3,6)))
609 ELSEIF(I22_ALEUL==1)THEN
610 II = brickID
611 !---8 local node numbers NC1 TO NC8 for solid element I ---!
612 NC1=IXS(2,II)
613 NC2=IXS(3,II)
614 NC3=IXS(4,II)
615 NC4=IXS(5,II)
616 NC5=IXS(6,II)
617 NC6=IXS(7,II)
618 NC7=IXS(8,II)
619 NC8=IXS(9,II)
620 !
621 !---Coordinates of the 8 nodes
622 X1_=X(1,NC1)
623 Y1_=X(2,NC1)
624 Z1_=X(3,NC1)
625 !
626 X2_=X(1,NC2)
627 Y2_=X(2,NC2)
628 Z2_=X(3,NC2)
629 !
630 X3_=X(1,NC3)
631 Y3_=X(2,NC3)
632 Z3_=X(3,NC3)
633 !
634 X4_=X(1,NC4)
635 Y4_=X(2,NC4)
636 Z4_=X(3,NC4)
637 !
638 X5_=X(1,NC5)
639 Y5_=X(2,NC5)
640 Z5_=X(3,NC5)
641 !
642 X6_=X(1,NC6)
643 Y6_=X(2,NC6)
644 Z6_=X(3,NC6)
645 !
646 X7_=X(1,NC7)
647 Y7_=X(2,NC7)
648 Z7_=X(3,NC7)
649 !
650 X8_=X(1,NC8)
651 Y8_=X(2,NC8)
652 Z8_=X(3,NC8)
653 !
654 ! Face-1
655 N(1,1)=(Y3_-Y1_)*(Z2_-Z4_) - (Z3_-Z1_)*(Y2_-Y4_)
656 N(2,1)=(Z3_-Z1_)*(X2_-X4_) - (X3_-X1_)*(Z2_-Z4_)
657 N(3,1)=(X3_-X1_)*(Y2_-Y4_) - (Y3_-Y1_)*(X2_-X4_)
658 ! Face-2
659 N(1,2)=(Y7_-Y4_)*(Z3_-Z8_) - (Z7_-Z4_)*(Y3_-Y8_)
660 N(2,2)=(Z7_-Z4_)*(X3_-X8_) - (X7_-X4_)*(Z3_-Z8_)
661 N(3,2)=(X7_-X4_)*(Y3_-Y8_) - (Y7_-Y4_)*(X3_-X8_)
662 ! Face-3
663 N(1,3)=(Y6_-Y8_)*(Z7_-Z5_) - (Z6_-Z8_)*(Y7_-Y5_)
664 N(2,3)=(Z6_-Z8_)*(X7_-X5_) - (X6_-X8_)*(Z7_-Z5_)
665 N(3,3)=(X6_-X8_)*(Y7_-Y5_) - (Y6_-Y8_)*(X7_-X5_)
666 ! Face-4
667 N(1,4)=(Y2_-Y5_)*(Z6_-Z1_) - (Z2_-Z5_)*(Y6_-Y1_)
668 N(2,4)=(Z2_-Z5_)*(X6_-X1_) - (X2_-X5_)*(Z6_-Z1_)
669 N(3,4)=(X2_-X5_)*(Y6_-Y1_) - (Y2_-Y5_)*(X6_-X1_)
670 ! Face-5
671 N(1,5)=(Y7_-Y2_)*(Z6_-Z3_) - (Z7_-Z2_)*(Y6_-Y3_)
672 N(2,5)=(Z7_-Z2_)*(X6_-X3_) - (X7_-X2_)*(Z6_-Z3_)
673 N(3,5)=(X7_-X2_)*(Y6_-Y3_) - (Y7_-Y2_)*(X6_-X3_)
674 ! Face-6
675 N(1,6)=(Y8_-Y1_)*(Z4_-Z5_) - (Z8_-Z1_)*(Y4_-Y5_)
676 N(2,6)=(Z8_-Z1_)*(X4_-X5_) - (X8_-X1_)*(Z4_-Z5_)
677 N(3,6)=(X8_-X1_)*(Y4_-Y5_) - (Y8_-Y1_)*(X4_-X5_)
678 !
679 BRICK_LIST(NIN,IB)%N(1,1:3) = N(1:3,1) / SQRT(SUM(N(1:3,1)*N(1:3,1)))
680 BRICK_LIST(NIN,IB)%N(2,1:3) = N(1:3,2) / SQRT(SUM(N(1:3,2)*N(1:3,2)))
681 BRICK_LIST(NIN,IB)%N(3,1:3) = N(1:3,3) / SQRT(SUM(N(1:3,3)*N(1:3,3)))
682 BRICK_LIST(NIN,IB)%N(4,1:3) = N(1:3,4) / SQRT(SUM(N(1:3,4)*N(1:3,4)))
683 BRICK_LIST(NIN,IB)%N(5,1:3) = N(1:3,5) / SQRT(SUM(N(1:3,5)*N(1:3,5)))
684 BRICK_LIST(NIN,IB)%N(6,1:3) = N(1:3,6) / SQRT(SUM(N(1:3,6)*N(1:3,6)))
685 ENDIF
686
687 DO J=1,NV46
688 ! compute Brick Face New Area !suplicated eflux3/aflux3, to optimize
689 F(J,IB) = HALF * SQRT( SUM( N(:,J) * N(:,J) ) ) ! N is 2Sn
690 ENDDO
691 pFACE(1:9) => BRICK_LIST(NIN,IB)%POLY(1:9)
692! pFACE(2)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(2)%FACE(1:6)%Surf
693! pFACE(3)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(3)%FACE(1:6)%Surf
694! pFACE(4)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(4)%FACE(1:6)%Surf
695! pFACE(5)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(5)%FACE(1:6)%Surf
696! pFACE(6)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(6)%FACE(1:6)%Surf
697! pFACE(7)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(7)%FACE(1:6)%Surf
698! pFACE(8)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(8)%FACE(1:6)%Surf
699! pFACE(9)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(9)%FACE(1:6)%Surf
700 pFullFACE => BRICK_LIST(NIN,IB)%Face_BRICK(1:6)
701 pFullFACE(1:6) = F(1:6,IB)
702 IF(BRICK_LIST(NIN,IB)%NBCUT==0)THEN
703 pFACE(1)%FACE(1:6)%Surf = F(1:6,IB)
704 CYCLE
705 ENDIF
706 NCELL = BRICK_LIST(NIN,IB)%NBCUT !cell 9 not yet taken into account.
707 DO ICELL=1,NCELL
708 !FACES UPDATE
709 IF(BRICK_LIST(NIN,IB)%POLY(ICELL)%Vnew<ZERO)THEN
710 DO J=1,NV46
711 pFACE(ICELL)%FACE(J)%Surf = F(J,IB) + pFACE(ICELL)%FACE(J)%Surf
712 if(pFACE(ICELL)%FACE(J)%Surf<ZERO)then
713 write (*,*) "**error : inter22 negative cell face"
714 endif
715 ENDDO! next J
716 ENDIF
717 !VOLUME UPDATE
718 VOL = BRICK_LIST(NIN,IB)%POLY(ICELL)%Vnew
719 IF(VOL<ZERO)THEN
720 GBUF => ELBUF_TAB(NGB(IB))%GBUF
721 BRICK_LIST(NIN,IB)%POLY(ICELL)%Vnew = GBUF%VOL(IDLOCB(IB)) + VOL
722 ENDIF
723 VOLratio = BRICK_LIST(NIN,IB)%POLY(ICELL)%Vnew / ELBUF_TAB(NGB(IB))%GBUF%VOL(IDLOCB(IB))
724 ! THRESHOLD VOLUME : Smax = O(V**(2/3)), dh= EM03 * O(V**(1/3)) => Veps = EM04*O(V) => VOLratio = Veps/V = EM04 (EM04 => dh = 0.01 % * average_edge)
725 IF(ABS(VOLratio) <= EM04)THEN
726 I22_DEGENERATED = 1
727 ENDIF
728 ENDDO!next ICELL
729 ENDDO !next IB (cut cell)
730
731
732
733
734 !------------------------------------------------------------!
735 ! @10. BUILD CELL 9 !
736 ! which is complementary one : !
737 ! (faces,volumes,brick nodes and num points) !
738 !------------------------------------------------------------!
739 DO IB=NBF,NBL
740 pListNodID(1)%p(1:8) => BRICK_LIST(NIN,IB)%POLY(1)%ListNodID(1:8)
741 pListNodID(2)%p(1:8) => BRICK_LIST(NIN,IB)%POLY(2)%ListNodID(1:8)
742 pListNodID(3)%p(1:8) => BRICK_LIST(NIN,IB)%POLY(3)%ListNodID(1:8)
743 pListNodID(4)%p(1:8) => BRICK_LIST(NIN,IB)%POLY(4)%ListNodID(1:8)
744 pListNodID(5)%p(1:8) => BRICK_LIST(NIN,IB)%POLY(5)%ListNodID(1:8)
745 pListNodID(6)%p(1:8) => BRICK_LIST(NIN,IB)%POLY(6)%ListNodID(1:8)
746 pListNodID(7)%p(1:8) => BRICK_LIST(NIN,IB)%POLY(7)%ListNodID(1:8)
747 pListNodID(8)%p(1:8) => BRICK_LIST(NIN,IB)%POLY(8)%ListNodID(1:8)
748 pListNodID(9)%p(1:8) => BRICK_LIST(NIN,IB)%POLY(9)%ListNodID(1:8)
749 pSUBVOL(1:9) => BRICK_LIST(NIN,IB)%POLY(1:9)
750 pWhichCellNod(1:8) => BRICK_LIST(NIN,IB)%NODE(1:8)
751 pFACE(1:9) => BRICK_LIST(NIN,IB)%POLY(1:9)!%FACE(1:6)%Surf
752! pFACE(2)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(2)%FACE(1:6)%Surf
753! pFACE(3)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(3)%FACE(1:6)%Surf
754! pFACE(4)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(4)%FACE(1:6)%Surf
755! pFACE(5)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(5)%FACE(1:6)%Surf
756! pFACE(6)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(6)%FACE(1:6)%Surf
757! pFACE(7)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(7)%FACE(1:6)%Surf
758! pFACE(8)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(8)%FACE(1:6)%Surf
759! pFACE(9)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(9)%FACE(1:6)%Surf
760 ID = BRICK_LISt(NIN,IB)%ID
761 NCELL = BRICK_LIST(NIN,IB)%NBCUT
762 Cod1 = IABS(BRICK_LIST(NIN,IB)%SecID_CELL(1))
763 Cod2 = IABS(BRICK_LIST(NIN,IB)%SecID_CELL(2))
764 IF(NCELL/=0)THEN !si decoupage de brique.
765 !---VOLUME AND CONNECTIVITY
766 BRICK_LIST(NIN,IB)%POLY(9)%FACE(1)%Surf = F(1,IB)-SUM( (/ (BRICK_LIST(NIN,IB)%POLY(K)%FACE(1)%Surf,K=1,NCELL)/) )
767 BRICK_LIST(NIN,IB)%POLY(9)%FACE(2)%Surf = F(2,IB)-SUM( (/ (BRICK_LIST(NIN,IB)%POLY(K)%FACE(2)%Surf,K=1,NCELL)/) )
768 BRICK_LIST(NIN,IB)%POLY(9)%FACE(3)%Surf = F(3,IB)-SUM( (/ (BRICK_LIST(NIN,IB)%POLY(K)%FACE(3)%Surf,K=1,NCELL)/) )
769 BRICK_LIST(NIN,IB)%POLY(9)%FACE(4)%Surf = F(4,IB)-SUM( (/ (BRICK_LIST(NIN,IB)%POLY(K)%FACE(4)%Surf,K=1,NCELL)/) )
770 BRICK_LIST(NIN,IB)%POLY(9)%FACE(5)%Surf = F(5,IB)-SUM( (/ (BRICK_LIST(NIN,IB)%POLY(K)%FACE(5)%Surf,K=1,NCELL)/) )
771 BRICK_LIST(NIN,IB)%POLY(9)%FACE(6)%Surf = F(6,IB)-SUM( (/ (BRICK_LIST(NIN,IB)%POLY(K)%FACE(6)%Surf,K=1,NCELL)/) )
772 BRICK_LIST(NIN,IB)%POLY(9)%FACE0%Surf = SUM( (/ (BRICK_LIST(NIN,IB)%POLY(K)%FACE0%Surf,K=1,NCELL) /) ) !sum of cut surfaces
773 BRICK_LIST(NIN,IB)%POLY(9)%Vnew = -SUM(BRICK_LIST(NIN,IB)%POLY(1:NCELL)%Vnew)
774 !print *, "psubvol(9)=", pSUBVOL(9)
775 BRICK_LIST(NIN,IB)%POLY(9)%NumNOD = 0 !should noty be needed
776 BRICK_LIST(NIN,IB)%POLY(9)%NumPOINT = 0 !should noty be needed
777 MNOD = 0
778 pPOINT(1:3) = ZERO
779 !saving brick nodes used by remaining polyhedron (complementary one)
780 DO J=1,8
781 IF(pWhichCellNod(J)%WhichCell/=0) CYCLE !a tagged brick node is already used by a polyhedron
782 pWhichCellNod(J)%WhichCell = 9
783 MNOD = MNOD + 1
784 pListNodID(9)%p(MNOD) = J !IXS(J+1,ID)
785 pPOINT(1) = pPOINT(1) + X(1,IXS(1+J,ID))
786 pPOINT(2) = pPOINT(2) + X(2,IXS(1+J,ID))
787 pPOINT(3) = pPOINT(3) + X(3,IXS(1+J,ID))
788 ENDDO
789 BRICK_LIST(NIN,IB)%POLY(9)%NumNOD = MNOD
790 BRICK_LIST(NIN,IB)%POLY(9)%NumPOINT =
791 . SUM(BRICK_LIST(NIN,IB)%POLY(1:NCELL)%NumPOINT ) - SUM(BRICK_LIST(NIN,IB)%POLY(1:NCELL)%NumNOD) + MNOD
792 GBUF => ELBUF_TAB(NGB(IB))%GBUF
793 BRICK_LIST(NIN,IB)%POLY(9)%Vnew = GBUF%VOL(IDLOCB(IB)) + BRICK_LIST(NIN,IB)%POLY(9)%Vnew
794 UNCUTvol = GBUF%VOL(IDLOCB(IB))
795 IF(ABS(BRICK_LIST(NIN,IB)%POLY(9)%Vnew /UNCUTVOL) <EM04
796.OR. . ABS(BRICK_LIST(NIN,IB)%POLY(1)%Vnew /UNCUTVOL) <EM04)THEN
797 I22_DEGENERATED = 1
798 ENDIF
799 !---AND ALSO CELL CENTER
800 !MNOD brick nodes above + K1 intersection points below.
801 NPOINT = BRICK_LIST(NIN,IB)%POLY(9)%NumPOINT
802 K1 = 0
803 DO J=1,12
804 NBCUT = BRICK_LIST(NIN,IB)%EDGE(J)%NBCUT
805 DO I=1,NBCUT
806 PointTMP(1:3) = BRICK_LIST(NIN,IB)%EDGE(J)%CUTPOINT(1:3,I)
807 pPOINT(1) = pPOINT(1) + PointTMP(1)
808 pPOINT(2) = pPOINT(2) + PointTMP(2)
809 pPOINT(3) = pPOINT(3) + PointTMP(3)
810 K1 = K1 + 1
811 ENDDO
812 ENDDO
813 IF(NPOINT/=0)THEN
814 pPOINT(1) = pPOINT(1) / (K1+MNOD) !NPOINT
815 pPOINT(2) = pPOINT(2) / (K1+MNOD) !NPOINT
816 pPOINT(3) = pPOINT(3) / (K1+MNOD) !NPOINT
817 BRICK_LIST(NIN,IB)%POLY(9)%CellCENTER(1:3) = pPOINT(1:3)
818 ENDIF
819 ENDIF!NCELL/=0
820 ENDDO!next IB
821
822 !------------------------------------------------------------!
823 ! @11. BUILD UNCUT CELL IN CUT CELL BUFFER !
824 !------------------------------------------------------------!
825 DO IB=NBF,NBL
826 IF(BRICK_LIST(NIN,IB)%NBCUT == 0)THEN
827 pSUBVOL(1:9) => BRICK_LIST(NIN,IB)%POLY(1:9)
828 pSUBVOL(1:9)%Vnew = 0
829 GBUF => ELBUF_TAB(NGB(IB))%GBUF
830 pSUBVOL(1)%Vnew = GBUF%VOL(IDLOCB(IB))
831 BRICK_LIST(NIN,IB)%Vnew_SCell = GBUF%VOL(IDLOCB(IB))
832 ENDIF
833 ENDDO
834
835 !------------------------------------------------------------!
836 ! @12. STORING BRICK REFERENCE VOLUME (UNCUT) !
837 !------------------------------------------------------------!
838 DO IB=NBF,NBL
839 BRICK_LIST(NIN,IB)%UncutVol = ELBUF_TAB(NGB(IB))%GBUF%VOL(IDLOCB(IB))
840 ENDDO
841
842 9 CONTINUE
843
844
845 !------------------------------------------------------------!
846 ! @13. POLY 9 : MARK IF NO NODE ON A FACE !
847 !------------------------------------------------------------!
848 !Poly9woNodes
849 DO IB=NBF,NBL
850 NCELL = BRICK_LIST(NIN,IB)%NBCUT
851 BRICK_LIST(NIN,IB)%Poly9woNodes(1:6,1:2) = 0
852 IF(NCELL<= 1) CYCLE
853 DO J=1,6
854 FACE = BRICK_LIST(NIN,IB)%POLY(9)%FACE(J)%Surf
855 IF(FACE ==ZERO) CYCLE
856 !loop on nodes (1,2,3,4,5,6,7,8) composing face J
857 lFOUND = .FALSE.
858 DO K=1,4
859 INOD = INT22_BUF%nodFACE(J,K)
860 ICELL = BRICK_LIST(NIN,IB)%NODE(INOD)%WhichCell
861 IF(ICELL == 9) lFOUND = .TRUE.
862 ENDDO
863 IF(lFOUND)CYCLE ! CELL 9 IS SHARING A NODE SO ITS ADJACENT CELL WAS FOUND ABOVE (See 14. FINDING & STORING ADJACENT CELLS)
864 !Here Cell 9 has a trace on face J but has no nodes on this face. An adjacent celle must be find (otherwise internal forces issue, missing area)
865 BRICK_LIST(NIN,IB)%Poly9woNodes(J,1) = 1
866 ENDDO!next J
867 !UPDATE FLUX FOR THIS CASE : NO NODES => NO FLUXES (CLOSED SURFACE HYPOTHESIS)
868 ENDDO! next IB
869
870
871
872 DO IB=NBF,NBL
873 BRICK_LIST(NIN,IB)%IsMergeTarget = 0
874 !print *, "0, ismergetarget", ixs(11,brick_list(nin,ib)%id)
875 ENDDO
876
877
878 !------------------------------------------------------------!
879 ! @14. FILL MATERIAL BUFFER : %bakMAT%... !
880 !------------------------------------------------------------!
881 DO IB=NBF,NBL
882 NCELL = BRICK_LIST(NIN,IB)%NBCUT
883 MCELL = BRICK_LIST(NIN,IB)%MainID
884 ID = BRICK_LIST(NIN,IB)%ID
885 NG = BRICK_LIST(NIN,IB)%NG
886 IDLOC = BRICK_LIST(NIN,IB)%IDLOC
887 MLW = BRIcK_LIST(NIN,IB)%MLW
888 GBUF => ELBUF_TAB(NG)%GBUF
889 LBUF => ELBUF_TAB(NG)%BUFLY(1)%LBUF(1,1,1)
890 MBUF => ELBUF_TAB(NG)%BUFLY(1)%MAT(1,1,1)
891 LLT_ = IPARG(2,NG)
892 ICELL = 0
893 IF(MCELL==0)THEN
894 BRICK_LIST(NIN,IB)%bakMAT%RHO = ZERO
895 BRICK_LIST(NIN,IB)%bakMAT%rhoE = ZERO
896 BRICK_LIST(NIN,IB)%bakMAT%rhoU(1) = ZERO
897 BRICK_LIST(NIN,IB)%bakMAT%rhoU(2) = ZERO
898 BRICK_LIST(NIN,IB)%bakMAT%rhoU(3) = ZERO
899 BRICK_LIST(NIN,IB)%bakMAT%ssp = ZERO
900 BRICK_LIST(NIN,IB)%bakMAT%SIG(1) = ZERO
901 BRICK_LIST(NIN,IB)%bakMAT%SIG(2) = ZERO
902 BRICK_LIST(NIN,IB)%bakMAT%SIG(3) = ZERO
903 BRICK_LIST(NIN,IB)%bakMAT%SIG(4) = ZERO
904 BRICK_LIST(NIN,IB)%bakMAT%SIG(5) = ZERO
905 BRICK_LIST(NIN,IB)%bakMAT%SIG(6) = ZERO
906 !------MULTI-MATERIAL-----------!
907 IF(I22LAW37+I22LAW51 == 0) CYCLE
908 BRICK_LIST(NIN,IB)%bakMAT%UVAR(1) = ZERO
909 BRICK_LIST(NIN,IB)%bakMAT%UVAR(2) = ZERO
910 BRICK_LIST(NIN,IB)%bakMAT%UVAR(3) = ZERO
911 BRICK_LIST(NIN,IB)%bakMAT%UVAR(4) = ZERO
912 BRICK_LIST(NIN,IB)%bakMAT%UVAR(5) = ZERO
913 IF(I22LAW51 == 0) CYCLE
914 DO K=6, I22LAW51
915 BRICK_LIST(NIN,IB)%bakMAT%UVAR(K) = ZERO
916 ENDDO!next K
917 ELSE
918 BRICK_LIST(NIN,IB)%bakMAT%RHO = GBUF%RHO(IDLOC)
919 BRICK_LIST(NIN,IB)%bakMAT%rhoE = GBUF%EINT(IDLOC)
920 BRICK_LIST(NIN,IB)%bakMAT%rhoU(1) = GBUF%MOM(LLT_*(1-1) +IDLOC)
921 BRICK_LIST(NIN,IB)%bakMAT%rhoU(2) = GBUF%MOM(LLT_*(2-1) +IDLOC)
922 BRICK_LIST(NIN,IB)%bakMAT%rhoU(3) = GBUF%MOM(LLT_*(3-1) +IDLOC)
923 BRICK_LIST(NIN,IB)%bakMAT%ssp = LBUF%SSP(IDLOC)
924 BRICK_LIST(NIN,IB)%bakMAT%SIG(1) = GBUF%SIG(LLT_*(1-1) +IDLOC)
925 BRICK_LIST(NIN,IB)%bakMAT%SIG(2) = GBUF%SIG(LLT_*(2-1) +IDLOC)
926 BRICK_LIST(NIN,IB)%bakMAT%SIG(3) = GBUF%SIG(LLT_*(3-1) +IDLOC)
927 BRICK_LIST(NIN,IB)%bakMAT%SIG(4) = GBUF%SIG(LLT_*(4-1) +IDLOC)
928 BRICK_LIST(NIN,IB)%bakMAT%SIG(5) = GBUF%SIG(LLT_*(5-1) +IDLOC)
929 BRICK_LIST(NIN,IB)%bakMAT%SIG(6) = GBUF%SIG(LLT_*(6-1) +IDLOC)
930 !------MULTI-MATERIAL-----------!
931 !IF(I22LAW37+I22LAW51 == 0) CYCLE
932.AND. IF(MLW/=37 MLW/=51)CYCLE
933 BRICK_LIST(NIN,IB)%bakMAT%UVAR(1) = MBUF%VAR((1-1)*LLT_+IDLOC)
934 BRICK_LIST(NIN,IB)%bakMAT%UVAR(2) = MBUF%VAR((2-1)*LLT_+IDLOC)
935 BRICK_LIST(NIN,IB)%bakMAT%UVAR(3) = MBUF%VAR((3-1)*LLT_+IDLOC)
936 BRICK_LIST(NIN,IB)%bakMAT%UVAR(4) = MBUF%VAR((4-1)*LLT_+IDLOC)
937 BRICK_LIST(NIN,IB)%bakMAT%UVAR(5) = MBUF%VAR((5-1)*LLT_+IDLOC)
938 IF(I22LAW51 == 0) CYCLE
939 DO K=6, I22LAW51
940 BRICK_LIST(NIN,IB)%bakMAT%UVAR(K) = MBUF%VAR((K-1)*LLT_+IDLOC)
941 ENDDO!next K
942 ENDIF
943 ENDDO!next IB
944
945 !-------------------!
946 CALL MY_BARRIER()
947 !-------------------!
948
949 !------------------------------------------------------------!
950 ! @15. FINDING & STORING ADJACENT CELLS !
951 !------------------------------------------------------------!
952 DO IB=NBF,NBL
953 pFACE(1:9) => BRICK_LIST(NIN,IB)%POLY(1:9)
954! pFACE(2)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(2)%FACE(1:6)%Surf
955! pFACE(3)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(3)%FACE(1:6)%Surf
956! pFACE(4)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(4)%FACE(1:6)%Surf
957! pFACE(5)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(5)%FACE(1:6)%Surf
958! pFACE(6)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(6)%FACE(1:6)%Surf
959! pFACE(7)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(7)%FACE(1:6)%Surf
960! pFACE(8)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(8)%FACE(1:6)%Surf
961! pFACE(9)%p(1:6) => BRICK_LIST(NIN,IB)%POLY(9)%FACE(1:6)%Surf
962 !pNAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%NAdjCell
963 pAdjBRICK => BRICK_LIST(NIN,IB)%Adjacent_Brick(1:6,1:5)
964 pListNodID(1)%p(1:8) => BRICK_LIST(NIN,IB)%POLY(1)%ListNodID(1:8)
965 pListNodID(2)%p(1:8) => BRICK_LIST(NIN,IB)%POLY(2)%ListNodID(1:8)
966 pListNodID(3)%p(1:8) => BRICK_LIST(NIN,IB)%POLY(3)%ListNodID(1:8)
967 pListNodID(4)%p(1:8) => BRICK_LIST(NIN,IB)%POLY(4)%ListNodID(1:8)
968 pListNodID(5)%p(1:8) => BRICK_LIST(NIN,IB)%POLY(5)%ListNodID(1:8)
969 pListNodID(6)%p(1:8) => BRICK_LIST(NIN,IB)%POLY(6)%ListNodID(1:8)
970 pListNodID(7)%p(1:8) => BRICK_LIST(NIN,IB)%POLY(7)%ListNodID(1:8)
971 pListNodID(8)%p(1:8) => BRICK_LIST(NIN,IB)%POLY(8)%ListNodID(1:8)
972 pListNodID(9)%p(1:8) => BRICK_LIST(NIN,IB)%POLY(9)%ListNodID(1:8)
973 ID = BRICK_LIST(NIN,IB)%ID
974 NCELL = BRICK_LIST(NIN,IB)%NBCUT
975 !pNAdjCELL(1:9,1:6) = 0
976 DO K=1,6
977 BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(K)%Adjacent_Cell(1) = 0
978 BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(K)%Adjacent_Cell(2) = 0
979 BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(K)%Adjacent_Cell(3) = 0
980 BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(K)%Adjacent_Cell(4) = 0
981 BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(K)%Adjacent_Cell(5) = 0
982 BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(K)%NAdjCell = 0
983 ENDDO
984 ICELL = 0
985 DO WHILE (ICELL<=NCELL) ! loop on polyhedron {1:NCELL} U {9}
986 ICELL = ICELL +1
987.AND. IF (ICELL>NCELL NCELL/=0)ICELL=9
988 DO J=1,NV46 ! loop on brick face
989 IF(pFACE(ICELL)%FACE(J)%Surf>ZERO)THEN ! if polyhedron has a non zero subface on this brick face
990 IV = pAdjBRICK(J,1)
991 IF(IV==0)CYCLE !if there is no adjacent brick on this face
992 !if there is an adjacent brick on this face
993 IAD22 = pAdjBRICK(J,4)
994 IF(IAD22 == 0)THEN
995 !adjacent brick is in cut cell buffer : it also uncut
996 BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(J)%NAdjCell = 1
997 BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(J)%Adjacent_Cell(1) = 1 !necessarily ICELLv=1
998 ELSEIF(BRICK_LIST(NIN,IAD22)%NBCUT==0)THEN
999 !adjacent brick is in cut cell buffer but is uncut
1000 BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(J)%NAdjCell = 1
1001 BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(J)%Adjacent_Cell(1) = 1 !necessarily ICELLv=1
1002 ELSE
1003
1004 pWhichCellNodv(1:8) => BRICK_LIST(NIN,IAD22)%NODE(1:8)
1005 !---brick nodes used by reference cell, must lay on current face J, and are identified by internal id : IXS(1+.,ID)
1006 NNODES = 0
1007 DO K1=1,BRICK_LIST(NIN,IB)%POLY(ICELL)%NumNOD !loop on polyhedron nodes !TODO : might be optimized by looping on nodFACE and cycle with pNNODCELL
1008 DO K2=1,4 !loop on face nodes
1009 IF(pListNodID(ICELL)%p(K1) == INT22_BUF%nodFACE(J,K2))THEN !polyhedron node on the current face
1010 NNODES = NNODES +1
1011 INODES(NNODES) = IXS(1+INT22_BUF%nodFACE(J,K2), ID)
1012 ENDIF
1013 ENDDO
1014 ENDDO
1015 !---search for local nodes id in adjacent brick.
1016 ICELLTAG(1:9) = 0
1017 JV = pAdjBRICK(J,5)
1018 DO K0=1,4 !loop on adjacent 'local' nodes (II,J)<->(IV,JV)
1019 K1 = INT22_BUF%nodFACE(JV,K0)
1020 DO IN=1,NNODES
1021 IF(IXS(1+K1,IV)==INODES(IN) )THEN
1022 ICELLv = pWhichCellNodv(K1)%WhichCell
1023 if (ICELLv==0)then
1024 print *, "itask,icellv,pwhichcellnodv(K1)",ITASK,ICELLv,pWhichCellNodv(K1)%WhichCell
1025 stop 2245
1026 endif
1027 IF(ICELLTAG(ICELLv)==0)THEN !check if cell is already stored as an adjacent cell, if not store it
1028 ICELLTAG( ICELLv ) = 1
1029 BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(J)%NAdjCell = BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(J)%NAdjCell + 1
1030 NAdjCell = BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(J)%NAdjCell
1031 BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(J)%Adjacent_Cell(NadjCell) = ICELLv
1032 ENDIF
1033 ENDIF
1034 ENDDO!next IN (polyhedron node on face)
1035 ENDDO!next K0 (adjacent 'local' node)
1036 ENDIF!testing IV & IBV
1037 ENDIF!(pFACE(ICELL)%p(J)>ZERO)
1038
1039 !COMPLEMENTARY POLYHEDRON WITH NO NODES ON A GIVEN FACE
1040 Poly9woNodes = BRICK_LIST(NIN,IB)%Poly9woNodes(J,1)
1041 ! forcemment IBv /= 0
1042 IF(Poly9woNodes == 1)THEN
1043 IV = BRICK_LIST(NIN,IB)%Adjacent_Brick(J,1)
1044 IBv = BRICK_LIST(NIN,IB)%Adjacent_Brick(J,4)
1045 IF(IBv==0)THEN
1046 BRICK_LIST(NIN,IB)%Poly9woNodes(J,2) = 0
1047 CYCLE
1048 ENDIF
1049 JV = BRICK_LIST(NIN,IB)%Adjacent_Brick(J,5)
1050 Poly9woNodesV = BRICK_LIST(NIN,IBv)%Poly9woNodes(Jv,1)
1051 IF(Poly9woNodesV==0)THEN
1052 NBCUTv = BRICK_LIST(NIN,IBv)%NBCUT
1053 IF(NBCUTv == 0)THEN
1054 !PAS DE VOISIN HYPOTHESE SURFACE FERMEE
1055 !le polyedre voisin est le polyedre 9
1056 ! +--------+--+-----------+
1057 ! | \1| !
1058 ! | | |
1059 ! | 9 | 1 |
1060 ! | | |
1061 ! | / | |
1062 ! | / 2 | |
1063 ! +------+----+-----------+
1064 !BRICK_LIST(NIN,IB)%POLY(9)%FACE(J)%NAdjCell = BRICK_LIST(NIN,IB)%PoLY(9)%FACE(J)%NAdjCell + 1
1065 !NAdjCell = BRICK_LIST(NIN,IB)%POLY(9)%FACE(J)%NAdjCell
1066 !BRICK_LIST(NIN,IB)%Adjacent_Brick(J,9,NadjCell) = 1
1067 BRICK_LIST(NIN,IB)%Poly9woNodes(J,2) = 0
1068C !DESTRUCTION POLYEDRE VOISIN DETECTE AU DEBUT DU PARAGRAPHE @14.
1069C pNAdjCell(ICELL,J) = 0
1070C pAdjCELL(J,ICELL,1) = 0 !necessarily ICELLv=1
1071C non sinon cell trapped error et aussi on ne peut pas trouver de pression voisin poru le contact
1072 ELSE
1073 !le polyedre voisin utilise le point d'intersection le plus proche du noeud sommet => le polyedre voisin deveint systematiquement le polyedre 9. (faire dessins)
1074 ! +--------+--+-----------+ +--------+--+-----------+
1075 ! | \1| ! | \1| 1 !
1076 ! | + ! | +-----------+
1077 ! | 9 | 9 | | 9 | 9 |
1078 ! | +-----------+ | | |
1079 ! | / | 1 | | / | |
1080 ! | / 2 | | | / 2 | |
1081 ! +------+----+-----------+ +------+----+-----------+
1082 !autre hypothese => ... = 1
1083 BRICK_LIST(NIN,IB)%POLY(9)%FACE(J)%NAdjCell = BRICK_LIST(NIN,IB)%POLY(9)%FACE(J)%NAdjCell + 1
1084 NAdjCell = BRICK_LIST(NIN,IB)%POLY(9)%FACE(J)%NAdjCell
1085 BRICK_LIST(NIN,IB)%POLY(9)%FACE(J)%Adjacent_Cell(NadjCell) = 9
1086 BRICK_LIST(NIN,IB)%Poly9woNodes(J,2) = 9
1087 ENDIF
1088 ELSEIF(Poly9woNodesV/=0)THEN
1089 !le polyedre voisin est le polyedre 9
1090 ! +--------+--+-----------+
1091 ! | \1| 2 !
1092 ! | +-----------+
1093 ! | 9 | 9 |
1094 ! | +-----------+
1095 ! | / | 1 |
1096 ! | / 2 | |
1097 ! +------+----+-----------+
1098
1099 ! SIMPLY CONNECTING poly_9 <-> polyV_9
1100 BRICK_LIST(NIN,IB)%POLY(9)%FACE(J)%NAdjCell = BRICK_LIST(NIN,IB)%POLY(9)%FACE(J)%NAdjCell + 1
1101 NAdjCell = BRICK_LIST(NIN,IB)%POLY(9)%FACE(J)%NAdjCell
1102 BRICK_LIST(NIN,IB)%POLY(9)%FACE(J)%Adjacent_Cell(NadjCell) = 9
1103 BRICK_LIST(NIN,IB)%Poly9woNodes(J,2) = 9
1104 ENDIF!(Poly9woNodesV==0)
1105 ENDIF
1106 ENDDO!next J (face)
1107 ENDDO!next ICELL
1108 ENDDO!next IB
1109
1110
1111 !------------------------------------------------------------!
1112 ! @16. BUILD POLYGONAL TRACE FOR POLY9 !
1113 !------------------------------------------------------------!
1114 !CUT CELL BUFFER REMINDER (check common source for complete & update source)
1115 !my_real :: N(3) !normal vector
1116 !my_real :: B(3) !basis
1117 !my_real :: SCUT(1) !Cut Area
1118 !my_real :: SEFF(1) !Efficace Section = Wetted Surface : Cut Aera - Holes (computed by suming lagrangian projection)
1119 !my_real :: P(3,8) !Cut Surface Polygon (coordinates)
1120 !INTEGER :: NP !Cut Surface Polygon (number of points) ! NP<=6 | NP<=8 for additional Scut (closed surface hypothesis)
1121 !BRICK_LIST(NIN,IB)%Edge(J)%NODE(1:2) = EDGE_LIST(NIN,K)%NODE(1:2)
1122 !BRICK_LIST(NIN,IB)%Edge(J)%NBCUT = EDGE_LIST(NIN,K)%NBCUT
1123 !BRICK_LIST(NIN,IB)%Edge(J)%CUTSHELL(1:2) = EDGE_LIST(NIN,K)%CUTSHELL(1:2)
1124 !BRICK_LIST(NIN,IB)%Edge(J)%CUTCOOR(1:2) = EDGE_LIST(NIN,K)%CUTCOOR(1:2)
1125 !BRICK_LIST(NIN,IB)%Edge(J)%CUTPOINT(1:3,1) = EDGE_LIST(NIN,K)%CUTPOINT(1:3,1)
1126 !BRICK_LIST(NIN,IB)%Edge(J)%CUTPOINT(1:3,2) = EDGE_LIST(NIN,K)%CUTPOINT(1:3,2)
1127 !BRICK_LIST(NIN,IB)%Edge(J)%CUTVEL(1:3,1) = EDGE_LIST(NIN,K)%CUTVEL(1:3,1)
1128 !BRICK_LIST(NIN,IB)%Edge(J)%CUTVEL(1:3,2) = EDGE_LIST(NIN,K)%CUTVEL(1:3,2)
1129 !BRICK_LIST(NIN,IB)%Edge(J)%VECTOR(1:3) = EDGE_LIST(NIN,K)%VECTOR(1:3)
1130 !BRICK_LIST(NIN,IB)%Edge(J)%LEN = EDGE_LIST(NIN,K)%LEN
1131 !
1132 !
1133 ! +--------+--+ +-------+---+
1134 ! | \1| | | \ |
1135 ! | + Check WhichCellNode for each node of the face | | +
1136 ! | 9 | if digits 1,2, and 9 tagged then a polygon has => | | P |
1137 ! | + to be computed | | +
1138 ! | / | | | /|
1139 ! | / 2 | | |/ |
1140 ! +------+----+ +------+----+
1141 !
1142 !
1143 !must be counterclockwise from outward ( clockwise from point inside the brick )
1144 DO IB=NBF,NBL
1145 IE = BRICK_LIST(NIN,IB)%ID
1146 BRICK_LIST(NIN,IB)%ClosedSurf(1:6) = 0
1147 DO J=1,6
1148 ICODE = 0
1149 IV = BRICK_LIST(NIN,IB)%Adjacent_Brick(J,1)
1150 IF(IV==0)CYCLE
1151 DO K=1,4
1152 INOD = INT22_BUF%nodFACE(J,K)
1153 ICELL = BRICK_LIST(NIN,IB)%NODE(INOD)%WhichCell
1154 ICODE = IBSET(ICODE,ICELL)
1155 ENDDO
1156
1157 !ON EN PROFITE POUR MARQUER LES CAS DE FERMETURE (CLOSURE HYPOTHESIS)
1158.OR. IF(ICODE == 518 ICODE == 6) THEN
1159 FACE9 = BRICK_LIST(NIN,IB)%POLY(9)%FACE(J)%Surf
1160 IBV = BRICK_LIST(NIN,IB)%Adjacent_Brick(J,4)
1161 !necessairement IBV>0 sinon on ne peut pas avoir ICODE=518 (autrement dit : la brique IB est intersectee donc la voisine est dans le buffer)
1162 NBCUTv = BRICK_LIST(NIN,IBv)%NBCUT
1163.AND. IF(NBCUTv == 0 FACE9>ZERO) THEN
1164 !DISTINCITON IMPORTANTE CAR ON A BESOIN DU VOISIN POUR LES FORCES DE COTNACT MAIS SINON IL DOIT ETRE IGNORE POUR FINT ET CONVECTION
1165 BRICK_LIST(NIN,IB)%ClosedSurf(J) = 1
1166C print *, "brick_id=",ixs(11,iv)," face=",J, " has a closed surface"
1167C print *, "sinit22_fvm marked closure hypothesis."
1168 ENDIF
1169 ENDIF
1170
1171 Poly9woNodes = BRICK_LIST(NIN,IB)%Poly9woNodes(J,1)
1172 ICELLv = BRICK_LIST(NIN,IB)%Poly9woNodes(J,2)
1173.AND. IF(Poly9woNodes /=0 ICELLv /= 0)CYCLE
1174.AND. IF(Poly9woNodes==0 ICODE /= 518)CYCLE
1175 ! poly 1,2,9 on this current face => ICODE = b1000000110 = 518 (bits 1,2,9)
1176 NP = 0
1177 pPOINT(1:3) = ZERO
1178 CUT_Vel(1:3) = ZERO
1179 NBCUTprev = 0
1180 DO INOD=4,1,-1
1181 !base point
1182 pointTMP(1:3) = X(1:3,IXS(1+INT22_BUF%NodFace(J,INOD),IE))
1183 pPOINT(1) = pPOINT(1) + pointTMP(1)
1184 pPOINT(2) = pPOINT(2) + pointTMP(2)
1185 pPOINT(3) = pPOINT(3) + pointTMP(3)
1186 !polygon point
1187 iEDG = INT22_BUF%iLeftEdge(J,INOD)
1188 ISGN(1:2) = (/1,2/)
1189 IF(iEDG < 0) ISGN(1:2) = (/2,1/)
1190 iEDG = IABS(iEDG)
1191 NBCUT = BRICK_LIST(NIN,IB)%Edge(iEDG)%NBCUT
1192 IF(NBCUT == 0)THEN
1193 CYCLE
1194 ELSEIF(NBCUT == 1)THEN
1195 NP = NP + 1
1196 CUT_Point(1:3) = BRICK_LIST(NIN,IB)%Edge(iEDG)%CUTPOINT(1:3,1)
1197 BRICK_LIST(NIN,IB)%PCUT(8+J)%P(1:3,NP) = CUT_Point(1:3)
1198 CUT_Vel(1:3) = CUT_Vel(1:3) + BRICK_LIST(NIN,IB)%Edge(iEDG)%CUTVEL(1:3,1)
1199! IF(NBCUTprev==1)THEN
1200! !rajouter le point sommet au polygone 9
1201! NP = NP + 1
1202! NP _NODE = NP_NODE + 1
1203! CUT_Point(1:3) = X(1:3,IXS(1+iEDGE(ISGN(1),iEDG) ,IE))
1204! BRICK_LIST(NIN,IB)%PCUT(8+J)%P(1:3,NP) = CUT_Point(1:3)
1205! !do not increment CUTVEL
1206! print *, "pause-22- "
1207! pause
1208! ENDIF
1209 ELSE
1210 NP = NP + 1
1211 CUT_Point(1:3) = BRICK_LIST(NIN,IB)%Edge(iEDG)%CUTPOINT(1:3,ISGN(1))
1212 BRICK_LIST(NIN,IB)%PCUT(8+J)%P(1:3,NP) = CUT_Point(1:3)
1213 CUT_Vel(1:3) = CUT_Vel(1:3) + BRICK_LIST(NIN,IB)%Edge(iEDG)%CUTVEL(1:3,ISGN(1))
1214 NP = NP + 1
1215 CUT_Point(1:3) = BRICK_LIST(NIN,IB)%Edge(iEDG)%CUTPOINT(1:3,ISGN(2))
1216 CUT_Vel(1:3) = CUT_Vel(1:3) + BRICK_LIST(NIN,IB)%Edge(iEDG)%CUTVEL(1:3,ISGN(2))
1217 BRICK_LIST(NIN,IB)%PCUT(8+J)%P(1:3,NP) = CUT_Point(1:3)
1218 ENDIF
1219 NBCUTprev = NBCUT
1220 ENDDO!next INOD
1221 BRICK_LIST(NIN,IB)%PCUT(8+J)%NP = NP
1222 BRICK_LIST(NIN,IB)%PCUT(8+J)%B(1) = FOURTH * pPOINT(1)
1223 BRICK_LIST(NIN,IB)%PCUT(8+J)%B(2) = FOURTH * pPOINT(2)
1224 BRICK_LIST(NIN,IB)%PCUT(8+J)%B(3) = FOURTH * pPOINT(3)
1225 BRICK_LIST(NIN,IB)%PCUT(8+J)%N(1:3) = BRICK_LIST(NIN,IB)%N(J,1:3) * BRICK_LIST(NIN,IB)%POLY(9)%FACE(J)%Surf
1226 BRICK_LIST(NIN,IB)%PCUT(8+J)%SEFF = ZERO
1227 !BRICK_LIST(NIN,IB)%POLY(9)%FACE0%Vel(1:3) = CUT_Vel(1:3) / (NP*ONE)
1228 BRICK_LIST(NIN,IB)%PCUT(8+J)%Vel(1:3) = CUT_Vel(1:3) / (NP*ONE)
1229 pPOINT(1:3) = BRICK_LIST(NIN,IB)%PCUT(8+J)%B(1:3)
1230 VecTMP(1:3) = I22AERA(NP, BRICK_LIST(NIN,IB)%PCUT(8+J)%P(1:3,1:NP), pPOINT(1:3) )
1231 FACE = SQRT(SUM(VecTMP(1:3)*VecTMP(1:3)))
1232 BRICK_LIST(NIN,IB)%PCUT(8+J)%SCUT = FACE
1233 ENDDO!next J
1234 ENDDO!next IB
1235
1236
1237 !-------------------!
1238 CALL MY_BARRIER()
1239 !-------------------!
1240
1241 !------------------------------------------------------------!
1242 ! @17. PRE-TREATMENT FOR DBLE SECND WITH BOTH V=0 !
1243 !------------------------------------------------------------!
1244 NBL1 = NBL
1245 IF(I22_DEGENERATED == 1)THEN
1246 ALLOCATE (DESTROY_DATA(7,2*NB)) !on each thread
1247 ELSE
1248 NBL1 = 0
1249 ENDIF
1250 !-------------------!
1251 CALL MY_BARRIER()
1252 !-------------------!
1253 IPOS = 0
1254 DO IB=NBF,NBL1
1255 NCELL = BRICK_LIST(NIN,IB)%NBCUT
1256 ICELL = 0
1257 DO WHILE (ICELL<=NCELL) ! loop on polyhedron {1:NCELL} U {9}
1258 ICELL = ICELL +1
1259.AND. IF (ICELL>NCELL NCELL/=0)ICELL=9
1260 VOL = BRICK_LIST(NIN,IB)%POLY(ICELL)%Vnew
1261 UNCUTVOL = BRICK_LIST(NIN,IB)%UncutVol
1262 RATIO = VOL / UNCUTVOL
1263 IF(ABS(RATIO)>EM04)CYCLE !Vol=1.01e-05, UNCUTVOL=10 => ratio > EM06. So EM04 is retained. THRESHOLD VOLUME : Smax = O(V**(2/3)), dh= EM03 * O(V**(1/3)) => Veps = EM04*O(V) => ratio=Veps/V = EM04
1264 IDLOC = MAXLOC(BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(1:6)%SURF,1)
1265 RatioF = BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(IDLOC)%SURF / BRICK_LIST(NIN,IB)%Face_Brick(IDLOC)
1266 !volume degenere
1267 !on cherche une correspondance voisine
1268 J = IDLOC
1269 NAdjCELL = BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(J)%NAdjCell
1270 DO IADJ = 1,NAdjCELL
1271 ICELLv = BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(J)%Adjacent_Cell(IADJ)
1272 IBV = BRICK_LIST(NIN,IB)%Adjacent_Brick(J,4)
1273 JV = BRICK_LIST(NIN,IB)%Adjacent_Brick(J,5)
1274 ITASKv = BRICK_LIST(NIN,IBv)%ITASK
1275 VOLv = BRICK_LIST(NIN,IBv)%POLY(ICELLv)%Vnew
1276 UNCUTVOLv = BRICK_LIST(NIN,IBv)%UncutVOL
1277 RATIOv = VOLv/UncutVOLv
1278 !Vol=1.01e-05, UNCUTVOL=10 => ratio > EM06. So EM04 is retained. THRESHOLD VOLUME : Smax = O(V**(2/3)), dh= EM03 * O(V**(1/3)) => Veps = EM04*O(V) => ratio=Veps/V = EM04
1279.OR. IF(ABS(RATIOv)>EM04 IBv>IB )CYCLE
1280.AND..OR. !IF( (ITASKv==ITASKIB>IBv) (ITASKv<ITASK) )THEN
1281 IPOS = IPOS + 1
1282 DESTROY_DATA(1,IPOS) = NIN
1283 DESTROY_DATA(2,IPOS) = IB
1284 DESTROY_DATA(3,IPOS) = ICELL
1285 DESTROY_DATA(4,IPOS) = ICELLv
1286 DESTROY_DATA(5,IPOS) = IBv
1287 DESTROY_DATA(6,IPOS) = J
1288 DESTROY_DATA(7,IPOS) = Jv
1289 ! ENDIF
1290 ENDDO!next ICV
1291 !J was IDLOC
1292 ENDDO!next I
1293 ENDDO!next IB
1294
1295 !detach from the previous loop and launch after barrier to avoid multi-threading issue
1296.AND..and. if(itask==0 ibug22_destroy==1 IPOS>0)then
1297 print *, ""
1298 print *, " |------i22_destroy_cell.f-------|"
1299 print *, " | initialization subroutine |"
1300 print *, " |-------------------------------|"
1301 print *, ""
1302 endif
1303
1304 !-------------------!
1305 CALL MY_BARRIER()
1306 !-------------------!
1307
1308 DO I=1,IPOS
1309 IB = DESTROY_DATA(2,I)
1310 ICELL = DESTROY_DATA(3,I)
1311 ICELLv = DESTROY_DATA(4,I)
1312 IBv = DESTROY_DATA(5,I)
1313 J = DESTROY_DATA(6,I)
1314 Jv = DESTROY_DATA(7,I)
1315 !print *, "destroy", ixs(11,brick_list(nin,ib)%id),ICELL
1316 CALL DESTROY_CELL( NIN, IB,ICELL,ICELLv,IBv,J,Jv, IXS, ITASK)
1317 ENDDO
1318
1319 IF(I22_DEGENERATED == 1)THEN ! this is activated if SECONDARY/SECONDARY adjacency with both 0-volumes
1320 !-------------------!
1321 CALL MY_BARRIER()
1322 !-------------------!
1323
1324 I22LOOP = I22LOOP + 1
1325 if(I22LOOP >= 2)then
1326 print *, "**error : inter22, unexpected situation."
1327 endif
1328 I22_DEGENERATED = 0
1329 GOTO 9 !can be optimized later by updating adjacency in destroy_cell.F or looping only on bricks related to destroyed cell. (need to reach finish line)
1330 ENDIF
1331
1332
1333 !------------------------------------------------------------!
1334 ! @18. TAGGING SECONDARY CUT CELL !
1335 !------------------------------------------------------------!
1336 DO IB=NBF,NBL
1337 NG = BRICK_LIST(NIN,IB)%NG
1338 pMainID => BRICK_LIST(NIN,IB)%MainID
1339 pIsMain(1:9) => BRICK_LIST(NIN,IB)%POLY(1:9)
1340 pIsMain(1:9)%IsMain = 0
1341 pMainID = 0
1342 IF(BRICK_LIST(NIN,IB)%NBCUT==0)THEN
1343 pIsMain(1)%IsMain = 1
1344 pMainID = 1
1345 CYCLE
1346 ENDIF
1347 GBUF => ELBUF_TAB(NGB(IB))%GBUF
1348 VOL = BRICK_LIST(NIN,IB)%UncutVol
1349 NCELL = BRICK_LIST(NIN,IB)%NBCUT
1350 ICELL = 0
1351 !tag les potentiel mains
1352 DO WHILE (ICELL<=NCELL) ! loop on polyhedron {1:NCELL} U {9}
1353 ICELL = ICELL +1
1354.AND. IF (ICELL>NCELL NCELL/=0)ICELL=9
1355 !if lower than 50% then it is a SECONDARY.
1356 VolCELL = BRICK_LIST(NIN,IB)%POLY(ICELL)%Vnew
1357 FAC = VolCELL / VOL
1358 IF(FAC>CritMerge22)THEN
1359 !print *, "critmerge22=", CritMerge22
1360 pIsMain(ICELL)%IsMain = 1
1361 ENDIF
1362 ENDDO!next ICELL
1363 K = SUM(pIsMain(1:9)%IsMain)
1364 !In case of several potential main :Unicity with Reproductibility (indepandant from numbering)
1365 IF(K==1)THEN
1366 pMainID = MAXLOC(pIsMain(1:9)%IsMain,1)
1367 ELSEIF(K>=2)THEN
1368 MCELL = GET_UNIQUE_MAIN_CELL(NIN,IB,K)
1369 pIsMain(1:9)%IsMain = 0
1370 pIsMain(MCELL)%IsMain = 1
1371 pMainID = MCELL
1372 !In case of no cell satisfy criteria, choose cell 9 to avoid problem with partial cut
1373 ELSEIF(K==0)THEN
1374 pIsMain(9)%IsMain = 1
1375 pMainID = 9
1376 ENDIF
1377 ENDDO!next IB
1378
1379 !-------------------!
1380 CALL MY_BARRIER()
1381 !-------------------!
1382 !------------------------------------------------------------!
1383 ! @19. LINKING SECONDARY CELL TO A main ONE !
1384 !------------------------------------------------------------!
1385 N_UNLINKED_L(ITASK) = 0
1386 !if no main cell is found, it becomes itself a main cell (example : INTERF/INT_22/1BRICK-EOS/CEL_only)
1387 DO IB=NBF,NBL
1388 IE = BRICK_LIST(NIN,IB)%ID
1389 NCELL = BRICK_LIST(NIN,IB)%NBCUT
1390 pAdjBRICK => BRICK_LIST(NIN,IB)%Adjacent_Brick(1:6,1:5)
1391 !pAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%Adjacent_Cell(1:5)
1392 !pNAdjCell => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%NAdjCell
1393 pIsMain(1:9) => BRICK_LIST(NIN,IB)%POLY(1:9)
1394 pMainID => BRICK_LIST(NIN,IB)%MainID
1395 MCELL = pMainID
1396 ICELL = 0
1397 BRICK_LIST(NIN,IB)%POLY(1:9)%WhereIsMain(1) = 0
1398 BRICK_LIST(NIN,IB)%POLY(1:9)%WhereIsMain(2) = 0
1399 BRICK_LIST(NIN,IB)%POLY(1:9)%WhereIsMain(3) = 0
1400 BRICK_LIST(NIN,IB)%POLY(1:9)%WhereIsMain(4) = 0
1401 IF(NCELL == 0) THEN
1402 BRICK_LIST(NIN,IB)%POLY(MCELL)%WhereIsMain(3) = IE
1403 BRICK_LIST(NIN,IB)%POLY(MCELL)%WhereIsMain(4) = IB
1404 CYCLE
1405 ENDIF
1406 DO WHILE (ICELL<=NCELL) ! loop on polyhedron {1:NCELL} U {9}
1407 ICELL = ICELL +1
1408.AND. IF (ICELL>NCELL NCELL/=0)ICELL=9
1409 !--------------------------------------!
1410 !if cell needs to be linked to a main one
1411 IF(ICELL == MCELL)THEN
1412 BRICK_LIST(NIN,IB)%POLY(ICELL)%WhereIsMain(1:2) = 0
1413 BRICK_LIST(NIN,IB)%POLY(ICELL)%WhereIsMain(3) = IE
1414 BRICK_LIST(NIN,IB)%POLY(ICELL)%WhereIsMain(4) = IB
1415 ELSE
1416 adjMAIN_face(1:6) = ZERO
1417 adjMAIN_centroid(1:3,1:6) = ZERO
1418 IDadj_MCELLv(1:6) = 0
1419 IVadj_MCELLv(1:6) = 0
1420 IBadj_MCELLv(1:6) = 0
1421 !Loop on face to list adjacent main cells
1422 IF(SUM(BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(1:6)%NAdjCell) == 0)THEN
1423 print *, "**error : inter22 - cell trapped. check lagrangian surface surrounding brick id:",
1424 . IXS(11,BRICK_LIST(NIN,IB)%ID) ; stop 2206
1425 CYCLE
1426 ENDIF
1427 DO J=1,NV46
1428 NAdjCell = BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(J)%NAdjCell
1429 IF(NAdjCell == 0)CYCLE
1430 IV = pAdjBRICK(J,1)
1431 IBV = pAdjBRICK(J,4)
1432 IF(IBV==0)THEN
1433 adjMAIN_face(J) = BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(J)%Surf
1434 IDadj_MCELLv(J) = 1
1435 IVadj_MCELLv(J) = IV
1436 IBadj_MCELLv(J) = 0
1437 adjMAIN_centroid(1,J) = SUM(X(1,IXS(2:9,IV))) / EIGHT
1438 adjMAIN_centroid(2,J) = SUM(X(2,IXS(2:9,IV))) / EIGHT
1439 adjMAIN_centroid(3,J) = SUM(X(3,IXS(2:9,IV))) / EIGHT
1440 ELSE
1441 DO K=1,NAdjCell
1442 ICELLv = BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(J)%Adjacent_Cell(K)
1443 JV = pAdjBRICK(J,5)
1444 IF(BRICK_LIST(NIN,IBV)%POLY(ICELLv)%IsMain == 1)THEN
1445 adjMAIN_face(J) = MIN (BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(J)%Surf,
1446 . BRICK_LIST(NIN,IBv)%POLY(ICELLv)%FACE(JV)%Surf)
1447 adjMAIN_centroid(1:3,J) = BRICK_LIST(NIN,IBV)%POLY(ICELLv)%CellCenter(1:3)
1448 IDadj_MCELLv(J) = ICELLv
1449 IVadj_MCELLv(J) = IV
1450 IBadj_MCELLv(J) = IBV
1451 EXIT !no more main cell out this face (only one)
1452 ENDIF
1453 ENDDO!next K
1454 ENDIF
1455 ENDDO
1456 sumFACE = SUM(adjMAIN_face(1:6))
1457 IF(sumFACE==ZERO)THEN
1458c print *, "**Warning inter22 : no main cell to link with cell from BRICK ID:",IXS(11,BRICK_LIST(NIN,IB)%ID)
1459 N_UNLINKED_L(ITASK) = N_UNLINKED_L(ITASK) + 1
1460 UNLINKED_CELLS_L(ITASK,1,N_UNLINKED_L(ITASK)) = IB
1461 UNLINKED_CELLS_L(ITASK,2,N_UNLINKED_L(ITASK)) = ICELL
1462 CYCLE !number of adjacent cell on given face J : no adjacent cell on face J means no need to search out there
1463 ENDIF
1464! adjMAIN_face(1:6) = adjMAIN_face(1:6) / sumFACE
1465 !IPOS = LINK_WITH_UNIQUE_MAIN_CELL(adjMAIN_face, adjMAIN_centroid)
1466 IPOS = MAXLOC(adjMAIN_face(1:6),1)
1467 BRICK_LIST(NIN,IB)%POLY(ICELL)%WhereIsMain(1) = IPOS
1468 BRICK_LIST(NIN,IB)%POLY(ICELL)%WhereIsMain(2) = IDadj_MCELLv(IPOS)
1469 BRICK_LIST(NIN,IB)%POLY(ICELL)%WhereIsMain(3) = IVadj_MCELLv(IPOS)
1470 BRICK_LIST(NIN,IB)%POLY(ICELL)%WhereIsMain(4) = IBadj_MCELLv(IPOS)
1471 ENDIF
1472 ENDDO!next ICELL
1473 ENDDO
1474 !-----------------------------------------------------------!
1475 CALL MY_BARRIER() !
1476 !-----------------------------------------------------------!
1477.OR..AND. if((ibug22_UndirectLink==1ibug22_UndirectLink==-1) N_UNLINKED_L(ITASK)>0)
1478 . print *, " **sinit** : unlinked cell synthesis "
1479 DO i=1,n_unlinked_l(itask)
1480 ib = unlinked_cells_l(itask,1,i)
1481 icell = unlinked_cells_l(itask,2,i)
1482 if(ibug22_undirectlink==1.OR.ibug22_undirectlink==-1)write(*,fmt='(I10,A1,I1)') ixs(11,brick_list(nin,ib)%ID),".",icell
1483 ENDDO
1484 if(ibug22_undirectlink==1 .AND. n_unlinked_l(itask)>0)print *, ""
1485
1486
1487 !------------------------------------------------------------!
1488 ! @20. LIST ALL SECONDARY CELLS LINKED TO A GIVEN main CELL !
1489 !------------------------------------------------------------!
1490 DO ib=nbf,nbl
1491 nbcut = brick_list(nin,ib)%NBCUT
1492 mcell = brick_list(nin,ib)%MainID
1493 IF(mcell == 0) cycle
1494 !pNAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%NAdjCell
1495 ie = brick_list(nin,ib)%ID
1496 nsecndnod = 0
1497 numsecnd = 0
1498 brick_list(nin,ib)%SecndList%Num = 0
1499 brick_list(nin,ib)%SecndList%VOL_Unmerged = brick_list(nin,ib)%POLY(mcell)%Vnew
1500 DO j=1,nv46
1501 iv = brick_list(nin,ib)%Adjacent_Brick(j,1)
1502 ngv = brick_list(nin,ib)%Adjacent_Brick(j,2)
1503 idlocv = brick_list(nin,ib)%Adjacent_Brick(j,3)
1504 ibv = brick_list(nin,ib)%Adjacent_Brick(j,4) !necessarily different from zero because a cut cell(nbcut>0) has all its neighbor in cut cell buffer(tzinf increased)
1505 ifv = brick_list(nin,ib)%Adjacent_Brick(j,5)
1506 nadj = brick_list(nin,ib)%POLY(mcell)%FACE(j)%NAdjCell
1507 IF(ibv==0)cycle
1508 DO iadj = 1,nadj
1509 icellv = brick_list(nin,ib)%POLY(mcell)%FACE(j)%Adjacent_Cell(iadj)
1510 ie_m = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(3)
1511 IF(ie_m/=ie)cycle
1512 !secnd cell IBV.ICELLv is added into the list since it is linked to IB.MCELL
1513 numsecnd = numsecnd + 1
1514 nsecndnod = nsecndnod + brick_list(nin,ibv)%POLY(icellv)%NumNOD
1515 brick_list(nin,ib)%SecndList%FM(numsecnd) = j
1516 brick_list(nin,ib)%SecndList%FV(numsecnd) = ifv
1517 brick_list(nin,ib)%SecndList%IV(numsecnd) = iv
1518 brick_list(nin,ib)%SecndList%IBV(numsecnd) = ibv
1519 brick_list(nin,ib)%SecndList%ICELLv(numsecnd) = icellv
1520 brick_list(nin,ib)%SecndList%VOL(numsecnd) = brick_list(nin,ibv)%POLY(icellv)%vnew
1521 brick_list(nin,ib)%SecndList%NumNOD_cell(numsecnd) = brick_list(nin,ibv)%POLY(icellv)%NumNOD
1522 brick_list(nin,ib)%SecndList%ListNodID(numsecnd,1:8) = brick_list(nin,ibv)%POLY(icellv)%ListNodID(1:8)
1523 brick_list(nin,ib)%SecndList%SURF_v(numsecnd) = brick_list(nin,ibv)%POLY(icellv)%FACE(ifv)%SURF
1524 enddo!next IADJ
1525 brick_list(nin,ib)%SecndList%Num = numsecnd
1526 enddo!next J
1527 brick_list(nin,ib)%SecndList%NumSecndNodes = nsecndnod
1528 ENDDO
1529
1530 !------------------------------------------------------------!
1531 ! @21. INDIRECT LINKING !
1532 ! TAG RELATED CELL AND STORE ITS INDIRECT LINK !
1533 !------------------------------------------------------------!
1534 !CREER LES LIAISON DANS LE BUFFER CUT CELL
1535 DO iunlink=1,n_unlinked_l(itask)
1536 ib = unlinked_cells_l(itask,1,iunlink)
1537 icell = unlinked_cells_l(itask,2,iunlink)
1538 adjface(:,:) = zero
1539 ipos = 0
1540 DO j=1,nv46
1541 nadj = brick_list(nin,ib)%POLY(icell)%FACE(j)%NAdjCell
1542 ibv = brick_list(nin,ib)%Adjacent_Brick(j,4)
1543
1544 DO iadj=1,nadj
1545 jv = brick_list(nin,ib)%Adjacent_Brick(j,5)
1546 icellv = brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_Cell(iadj)
1547 ipos = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(1)
1548 IF(ipos==0)cycle
1549 IF(ibv/=0)THEN
1550 adjface(j,iadj) = min(brick_list(nin,ib)%POLY(icell)%FACE(j)%Surf,
1551 . brick_list(nin,ibv)%POLY(icellv)%FACE(jv)%Surf)
1552 ELSE
1553 adjface(j,iadj) = brick_list(nin,ib)%POLY(icell)%FACE(j)%Surf
1554 ENDIF
1555 ENDDO
1556 ENDDO
1557 ires(1:2) = maxloc(adjface)
1558 iposf = ires(1) !face number
1559 iposiadj = ires(2) !adjacent cell number
1560 ibv = brick_list(nin,ib)%Adjacent_Brick(iposf,4)
1561 icellv = brick_list(nin,ib)%POLY(icell)%FACE(iposf)%Adjacent_Cell(iposiadj)
1562 ipos = 0
1563 IF(icellv/=0)
1564 . ipos = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(1)
1565 IF(ipos==0 .OR. adjface(ires(1),ires(2))==zero)THEN
1566 print *, "***error : inter22 unable to treat cell ",ixs(11,brick_list(nin,ib)%ID)
1567 print *, " Cell seems trapped with no adjacent cells"
1568 stop 2207
1569 ENDIF
1570 IF(ipos>=10)THEN
1571 print *, "***error : inter22 unable to treat cell ",ixs(11,brick_list(nin,ib)%ID)
1572 print *, " Cell seems trapped with no adjacent cells"
1573 print *, " Fluid mesh in interaction area should be refined"
1574 stop 2208
1575 ENDIF
1576 j = iposf*10 + ipos
1577 brick_list(nin,ib)%POLY(icell)%WhereIsMain(1) = j
1578 brick_list(nin,ib)%POLY(icell)%WhereIsMain(2) = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(2)
1579 brick_list(nin,ib)%POLY(icell)%WhereIsMain(3) = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(3)
1580 brick_list(nin,ib)%POLY(icell)%WhereIsMain(4) = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(4)
1581 IF((ibug22_undirectlink==1.OR.ibug22_undirectlink==-1) .AND. n_unlinked_l(itask)>0)THEN
1582 write(*,fmt='(A,I10,A1,I1,A,I2,I2,A1,I2,A2)') "unlinked cell:", ixs(11,brick_list(nin,ib)%ID),
1583 . ".",icell, " is now linked to faces ", iposf, ipos,"(",j," )"
1584 ENDIF
1585 ENDDO
1586
1587 !------------------------------------------------------------!
1588 ! @22. INDIRECT LINKING !
1589 ! UPDATE main CELL BUFFER !
1590 !------------------------------------------------------------!
1591 !need to wait all thread finishing to link before updating their supercell volume and their buffer
1592 !-------------------!
1593 CALL my_barrier()
1594 !-------------------!
1595
1596 !MERGE DES VOLUMES ET MISE A JOUR BUFFER main (main secnd peut etre sur un autre thread : traitement OpenMP)
1597 DO it=0,nthread-1
1598 DO iunlink=1,n_unlinked_l(it)
1599 ib = unlinked_cells_l(it,1,iunlink)
1600 icell = unlinked_cells_l(it,2,iunlink)
1601 ibm = brick_list(nin,ib)%POLY(icell)%WhereIsMain(4)
1602 itaskb = brick_list(nin,ibm)%ITASK
1603
1604! write (*,FMT='(A,I5,A,I5,I5,A,I5)')
1605! . "traitement ib=", ib, "itask,itaskb=", itask, itaskb," N_UNLINKED_L(IT)",N_UNLINKED_L(IT)
1606 IF(itaskb/=itask)cycle
1607 !Link Unlinked secnd cell to main IB,ICELL
1608 volcell = brick_list(nin,ib)%POLY(icell)%Vnew
1609 id = brick_list(nin,ib)%ID
1610 j = brick_list(nin,ib)%POLY(icell)%WhereIsMain(1)
1611 j1 = j/10
1612 j2 = mod(j,10)
1613 ibv = brick_list(nin,ib)%Adjacent_Brick(j1,4) !cell directly linked to a main
1614 IF(ibv==0)THEN
1615 print *, "inter22 :Error lagrangian surface is escaping eulerian domain."
1616 !with CALL ARRET() : truss might be with len=0 => time step issue.
1617 stop
1618 ENDIF
1619 fm = brick_list(nin,ibv)%Adjacent_Brick(j2,5)
1620 numsecnd = brick_list(nin,ibm)%SecndList%Num
1621 nsecndnod = brick_list(nin,ibm)%SecndList%NumSecndNodes
1622 numsecnd = numsecnd + 1
1623 nsecndnod = nsecndnod + brick_list(nin,ib)%POLY(icell)%NumNOD
1624 brick_list(nin,ibm)%SecndList%Num = numsecnd
1625 brick_list(nin,ibm)%SecndList%NumSecndNodes = nsecndnod
1626 brick_list(nin,ibm)%SecndList%FM(numsecnd) = fm
1627 brick_list(nin,ibm)%SecndList%FV(numsecnd) = j !j>10
1628 brick_list(nin,ibm)%SecndList%IV(numsecnd) = id
1629 brick_list(nin,ibm)%SecndList%IBV(numsecnd) = ib
1630 brick_list(nin,ibm)%SecndList%ICELLv(numsecnd) = icell
1631 brick_list(nin,ibm)%SecndList%VOL(numsecnd) = volcell
1632 brick_list(nin,ibm)%SecndList%NumNOD_cell(numsecnd) = brick_list(nin,ib)%POLY(icell)%NumNOD
1633 brick_list(nin,ibm)%SecndList%ListNodID(numsecnd,1:8) = brick_list(nin,ib)%POLY(icell)%ListNodID(1:8)
1634 brick_list(nin,ibm)%SecndList%SURF_v(numsecnd) = brick_list(nin,ib)%POLY(icell)%FACE(j1)%SURF
1635 enddo!next IUNLINK
1636 enddo!next IT
1637
1638 !augmenter le volume main
1639 !done later
1640
1641
1642 !------------------------------------------------------------!
1643 ! @23. TAG SECND CELLS WHICH ARE BASED ON FORMER main NODE !
1644 !------------------------------------------------------------!
1645! TARGET_DATA(:,:,:) = 0 !tag main cells linked to these secnd ones.
1646! Ntarget(NBF:NBL) = 0
1647 DO ib=nbf,nbl
1648 brick_list(nin,ib)%Ntarget = 0
1649 nbcut = brick_list(nin,ib)%NBCUT
1650 ncell = nbcut
1651 icell = 0
1652 pnodwasmain(1:8) => brick_list(nin,ib)%NODE(1:8)!%NodWasMain
1653 plistnodid(1)%p(1:8) => brick_list(nin,ib)%POLY(1)%ListNodID(1:8)
1654 plistnodid(2)%p(1:8) => brick_list(nin,ib)%POLY(2)%ListNodID(1:8)
1655 plistnodid(3)%p(1:8) => brick_list(nin,ib)%POLY(3)%ListNodID(1:8)
1656 plistnodid(4)%p(1:8) => brick_list(nin,ib)%POLY(4)%ListNodID(1:8)
1657 plistnodid(5)%p(1:8) => brick_list(nin,ib)%POLY(5)%ListNodID(1:8)
1658 plistnodid(6)%p(1:8) => brick_list(nin,ib)%POLY(6)%ListNodID(1:8)
1659 plistnodid(7)%p(1:8) => brick_list(nin,ib)%POLY(7)%ListNodID(1:8)
1660 plistnodid(8)%p(1:8) => brick_list(nin,ib)%POLY(8)%ListNodID(1:8)
1661 plistnodid(9)%p(1:8) => brick_list(nin,ib)%POLY(9)%ListNodID(1:8)
1662 mcell = brick_list(nin,ib)%MainID
1663 ntag = 0
1664 main_tag(:,:) = 0
1665
1666 brick_list(nin,ib)%MergeTarget(1:3,1:5) = 0
1667
1668 inod_w = brick_list(nin,ib)%OldMainStrongNode
1669 IF( inod_w*dt1 > zero)THEN
1670 IF (brick_list(nin,ib)%NODE(inod_w)%WhichCell==mcell)ncell = -1
1671 ENDIF
1672 !Finding secnd cell inside current brick element, which are based on former main cell nodes
1673 !------LOOP--------!
1674 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
1675 icell = icell +1
1676 IF (icell>ncell .AND. ncell/=0)icell=9
1677 IF(icell==mcell)cycle
1678 mnod = brick_list(nin,ib)%POLY(icell)%NumNOD
1679 CALL weighting_cell_nodes(nin,ib,icell,inod_w, idemerge)
1680 IF(sum(pnodwasmain(plistnodid(icell)%p(1:mnod))%NodWasMain)==0)cycle !if none of the polyhedra nodes is based on a former main cell node : nothing to do
1681 IF(pnodwasmain(inod_w)%NodWasMain<=0) cycle !(negatif en debug si IB nouveau dans le buffer)
1682 !ICELL has at least one node from the previous main cell
1683 ! lets get its main cell
1684
1685 !Continuity Criteria To avoid unexpected merging.
1686 !Consistent only with kinematic time step "several cycle to fully cross an element".
1687 vol_s = brick_list(nin,ib)%POLY(icell)%vnew
1688 uncutvol = brick_list(nin,ib)%UncutVol
1689 var = vol_s / (uncutvol)
1690 IF (var < critdvol22)cycle
1691 inod_w = brick_list(nin,ib)%OldMainStrongNode
1692 IF(inod_w==0)cycle !nothing to merge if no main or no cut cell in previous cycle.
1693
1694 ipos = brick_list(nin,ib)%POLY(icell)%WhereIsMain(1)
1695 icv = brick_list(nin,ib)%POLY(icell)%WhereIsMain(2)
1696 ibm = brick_list(nin,ib)%POLY(icell)%WhereIsMain(4)
1697 ntag = ntag + 1 !number of secnd cells wet by former main cell
1698c TARGET_DATA(IB,NTAG,1) = IPOS !facette localisant le plus grand main
1699c TARGET_DATA(IB,NTAG,2) = ICV !numero de cellule de ce main
1700 brick_list(nin,ib)%MergeTarget(1,ntag) = ipos
1701 brick_list(nin,ib)%MergeTarget(2,ntag) = icv
1702 brick_list(nin,ib)%MergeTarget(3,ntag) = ibm
1703 enddo! next ICELL
1704 !----END-LOOP------!
1705c Ntarget(IB) = NTAG ! number of secnd cell whose nodes were wet by the former main cell.= number of main to do material deportation (if one main per secnd cell)
1706 brick_list(nin,ib)%NTarget = ntag
1707 ! Ntarget always lower or equal to 9 (max number of secnd cell if one main per secnd)
1708 ! enumeration of the target
1709 !Warning : what happen if 2 secnd cells have the same main cell. Duplicates are not possible by tagging main target
1710 enddo!next IB
1711
1712 !------------------------------------------------------------!
1713 ! @24. VOLUME MERGING BETWEEN SECND CELL & ITS main ONE !
1714 ! UNCONDITIONAL !
1715 ! for law51 old volume need also to be merged !
1716 !------------------------------------------------------------!
1717 !loop on intersected bricks.
1718 DO ib=nbf,nbl
1719 ncell = brick_list(nin,ib)%NBCUT
1720 !pAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%Adjacent_Cell(1:5)
1721 !pNAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%NAdjCell
1722 padjbrick => brick_list(nin,ib)%Adjacent_Brick(1:6,1:5)
1723 psubvol(1:9) => brick_list(nin,ib)%POLY(1:9)
1724 mbuf => elbuf_tab(ngb(ib))%BUFLY(1)%MAT(1,1,1)
1725 mcell = brick_list(nin,ib)%MainID
1726 id = brick_list(nin,ib)%ID
1727 ng = brick_list(nin,ib)%NG
1728 idloc = brick_list(nin,ib)%IDLOC
1729 icell = 0
1730 brick_list(nin,ib)%Vnew_SCell = brick_list(nin,ib)%POLY(mcell)%Vnew
1731 !------LOOP--------!
1732 numsecnd = brick_list(nin,ib)%SecndList%Num
1733 DO ic=1,numsecnd
1734 icellv = brick_list(nin,ib)%SecndList%ICELLv(ic)
1735 ibv = brick_list(nin,ib)%SecndList%IBv(ic)
1736 ngv = brick_list(nin,ibv)%NG
1737 idlocv = brick_list(nin,ibv)%IDLOC
1738 iv = brick_list(nin,ibv)%ID
1739 !---MONO-MATERIAL---!
1740 psubvolv(1:9) => brick_list(nin,ibv)%POLY(1:9)
1741 volv = psubvolv(icellv)%Vnew
1742 !pSUBVOL(MCELL) = pSUBVOL(MCELL) + VOLv !no longer cumulated here
1743 brick_list(nin,ib)%Vnew_SCell = brick_list(nin,ib)%Vnew_SCell + volv
1744 enddo!next IC
1745 !initialisation de V1OLD
1746 IF(dt1==zero)THEN
1747 !---MULTI-MATERIAL----!
1748 mtn_ = iparg(1,ngb(ib))
1749 IF(mtn_==51)THEN
1750 llt_ = iparg(2,ngb(ib))
1751 !in sigesp51.F:
1752 !V1OLD = V1OLD - TIMESTEP*DDVOL1
1753 !Volume change correspond to DDVOL=DV/DT = f(fluxes). Gird deformation influences relative material velocity ans also fluxes
1754 DO itrimat=1,trimat
1755 ipos = 1
1756 k1 = m51_n0phas + (itrimat-1)*m51_nvphas +ipos-1 ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
1757 k = k1 * llt_
1758 vfrac(itrimat) = mbuf%VAR(k+idlocb(ib)) !backup for material merging
1759 ipos = 11
1760 k1 = m51_n0phas + (itrimat-1)*m51_nvphas +ipos-1 ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
1761 k = k1 * llt_
1762 var = brick_list(nin,ib)%Vnew_SCell*vfrac(itrimat)
1763 mbuf%VAR(k+idlocb(ib)) = var
1764 enddo!next ITRIMAT
1765 endif!(MTN_==51)
1766 endif!TIME==ZERO
1767 enddo!next IB
1768
1769
1770 !------------------------------------------------------------!
1771 ! @25. MATERIAL MERGING (IF main HAS TO BE DEPORTED) !
1772 !------------------------------------------------------------!
1773 ! if Ntarget > 0 then material deporation has to be done from former main cell buffer (GBUF)
1774 ! to the N=Ntarget target main cells.
1775 nbl1 = nbl
1776 IF(dt1==zero)nbl1 = 0
1777 if(ibug22_merge/=0)then
1778 if(itask==0)allocate(tmp22array(7,nb))
1779 call my_barrier
1780 tmp22array(1:7,nbf:nbl1)=zero
1781 endif
1782 DO ib=nbf,nbl1
1783 nbcut = brick_list(nin,ib)%NBCUT
1784 ncell = nbcut
1785 icell = 0
1786 gbuf => elbuf_tab(ngb(ib))%GBUF
1787 !pAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%Adjacent_Cell(1:5)
1788 padjbrick => brick_list(nin,ib)%Adjacent_Brick(1:6,1:5)
1789 mbuf => elbuf_tab(ngb(ib))%BUFLY(1)%MAT(1,1,1)
1790 idloc = brick_list(nin,ib)%IDLOC
1791 vol = zero
1792 volv = zero
1793 ntar = brick_list(nin,ib)%Ntarget
1794 llt_ = iparg(2,ngb(ib))
1795 !------LOOP--------!
1796 !loop on a main target to distribute former main data
1797 DO itar = 1,ntar
1798 j = brick_list(nin,ib)%MergeTarget(1,itar)
1799 icv = brick_list(nin,ib)%MergeTarget(2,itar)
1800 IF(j>=10)THEN
1801 j1 = j/10
1802 j2 = mod(j,10)
1803 ibv_i = brick_list(nin,ib )%Adjacent_Brick(j1,4)
1804 ngv = brick_list(nin,ibv_i)%Adjacent_Brick(j2,2)
1805 idlocv = brick_list(nin,ibv_i)%Adjacent_Brick(j2,3)
1806 ibv = brick_list(nin,ibv_i)%Adjacent_Brick(j2,4)
1807 ELSE
1808 ngv = brick_list(nin,ib )%Adjacent_Brick(j,2)
1809 idlocv = brick_list(nin,ib )%Adjacent_Brick(j,3)
1810 ibv = brick_list(nin,ib )%Adjacent_Brick(j,4) !TARGET CELL DATA IS STORED IN EINT_L(...,IBv), RHOL, UVARL...
1811 ENDIF
1812 gbufv => elbuf_tab(ngv)%GBUF
1813 ratio = one/ntar !in the loop but avoid conditional tests "IF(Ntar>0)" for NBF:NBL if not applicable
1814 itaskv = brick_list(nin,ibv)%ITASK
1815 imergel(itaskv) = 1
1816 brick_list(nin,ibv)%IsMergeTarget = 1
1817 !print *, "set isMergetarget ", ixs(11,brick_list(nin,ibv)%id)
1818 !---debug---!
1819 cod1=ixs(11,brick_list(nin,ib)%id)
1820 cod2=ixs(11,brick_list(nin,ibv)%id)
1821 if (ibug22_merge==-1)then
1822 !!!write(*,FMT='(A,I8,A,I8)') " MERGING : id=", Cod1, " to idv:", Cod2
1823 tmp22array(1,ib)=cod1
1824 tmp22array(2,ib)=cod2
1825 endif
1826 !-----------!
1827 !---------------------!
1828 ! MATERIAL MERGING ! !only if main cell loose all its nodes, main (which became secnd) contents is sent to the adjacent main.
1829 !---------------------!
1830 ! example:
1831 ! EINTv : target main on which secnd cell is linked
1832 ! EINT : old main cell buffer
1833 !-----------------------MONO-MATERIAL------------------------!
1834 vol = brick_list(nin,ib)%Vold_SCell
1835 supercellvol_l(itask,0,ibv) = supercellvol_l(itask,0,ibv) + ratio * vol
1836 !------ENERGIE--------!
1837 eint = ratio * vol * brick_list(nin,ib)%bakMAT%rhoE
1838 eint_l(itask,ibv) = eint_l(itask,ibv) + eint
1839 !------MASS----------!
1840 if (ibug22_merge==-1)then
1841 !!!write(*,FMT='(A,E30.16)') " RATIO=", RATIO
1842 !!!write(*,FMT='(A,E30.16)') " +rho=", GBUF%RHO(IDLOCB(IB))
1843 !!!write(*,FMT='(A,E30.16)') " +rhoUx=",GBUF%MOM(LLT_*(1-1) + IDLOCB(IB))
1844 !!!write(*,FMT='(A,E30.16)') " +Vol=", VOL
1845 tmp22array(3,ib)=ratio
1846 tmp22array(4,ib)=gbuf%RHO(idlocb(ib))
1847 tmp22array(5,ib)=gbuf%MOM(llt_*(1-1) + idlocb(ib))
1848 tmp22array(6,ib)=vol
1849 endif
1850 rho = ratio * vol * brick_list(nin,ib)%bakMAT%RHO
1851 rho_l(itask,ibv) = rho_l(itask,ibv) + rho
1852 !------TENSOR--------!
1853 DO j=1,nv46
1854 sig(j) = ratio * vol * brick_list(nin,ib)%bakMAT%SIG(j)
1855 sig_l(itask,j,ibv) = sig_l(itask,j,ibv) + sig(j)
1856 ENDDO
1857 !----MOMENTUM--------!
1858 !momentum is (rho.Ux, rho.Uy, rho.Uz)
1859 mom(1) = ratio * vol * brick_list(nin,ib)%bakMAT%rhoU(1)
1860 mom_l(1,itask,ibv) = mom_l(1,itask,ibv) + mom(1)
1861 mom(2) = ratio * vol * brick_list(nin,ib)%bakMAT%rhoU(2)
1862 mom_l(2,itask,ibv) = mom_l(2,itask,ibv) + mom(2)
1863 mom(3) = ratio * vol * brick_list(nin,ib)%bakMAT%rhoU(3)
1864 mom_l(3,itask,ibv) = mom_l(3,itask,ibv) + mom(3)
1865 !-----------------------MULTI-MATERIAL-----------------------!
1866 mtn_ = iparg(1,ngb(ib))
1867 IF(mtn_==37)THEN
1868 !UVAR(I,1) : massic percentage of liquid * global density (rho1*V1/V : it needs to give liquid mass multiplying by element volume in aleconve.F)
1869 !UVAR(I,2) : density of gas
1870 !UVAR(I,3) : density of liquid
1871 !UVAR(I,4) : volumetric fraction of liquid
1872 !UVAR(I,5) : volumetric fraction of gas
1873 llt_ = iparg(2,ngb(ib))
1874c print *, "merge-adding ratio uvar5=" , BRICK_LIST(NIN,IB)%bakMAT%UVAR(5)
1875c print *, "merge-adding ratio uvar4=" , BRICK_LIST(NIN,IB)%bakMAT%UVAR(4)
1876c print *, "merge-adding ratio uvar3=" , BRICK_LIST(NIN,IB)%bakMAT%UVAR(3)
1877c print *, "merge-adding ratio uvar2=" , BRICK_LIST(NIN,IB)%bakMAT%UVAR(2)
1878c print *, "merge-adding ratio uvar1=" , BRICK_LIST(NIN,IB)%bakMAT%UVAR(1)
1879
1880 uvarl(itask,ibv,5) = uvarl(itask,ibv,5)+ratio*vol* brick_list(nin,ib)%bakMAT%UVAR(5)
1881 uvarl(itask,ibv,4) = uvarl(itask,ibv,4)+ratio*vol* brick_list(nin,ib)%bakMAT%UVAR(4)
1882 uvarl(itask,ibv,3) = uvarl(itask,ibv,3)+ratio*vol* brick_list(nin,ib)%bakMAT%UVAR(3)*brick_list(nin,ib)%bakMAT%UVAR(4)
1883 uvarl(itask,ibv,2) = uvarl(itask,ibv,2)+ratio*vol* brick_list(nin,ib)%bakMAT%UVAR(2)*brick_list(nin,ib)%bakMAT%UVAR(5)
1884 uvarl(itask,ibv,1) = uvarl(itask,ibv,1)+ratio*vol* brick_list(nin,ib)%bakMAT%UVAR(1)
1885 supercellvol_l(itask,1,ibv) = supercellvol_l(itask,1,ibv) + uvarl(itask,ibv,4)
1886 supercellvol_l(itask,2,ibv) = supercellvol_l(itask,2,ibv) + uvarl(itask,ibv,5)
1887
1888 ELSEIF(mtn_==51)THEN
1889 llt_ = iparg(2,ngb(ib))
1890 !in sigesp51.F:
1891 !V1OLD = V1OLD - TIMESTEP*DDVOL1
1892 !Volume change correspond to DDVOL=DV/DT = f(fluxes). Gird deformation influences relative material velocity ans also fluxes
1893 DO itrimat=1,trimat
1894 ipos = 11
1895 k = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
1896 vol51(ib,itrimat) = brick_list(nin,ib)%bakMAT%UVAR(k) !backup for material merging
1897 var = ratio * vol51(ib,itrimat)
1898 IF(var/=zero) supercellvol_l(itask,itrimat,ibv) = supercellvol_l(itask,itrimat,ibv) + var
1899 enddo!next ITRIMAT
1900 uvarl(itask,ibv,1) = uvarl(itask,ibv,1) + brick_list(nin,ib)%bakMAT%UVAR(1)
1901 uvarl(itask,ibv,2) = uvarl(itask,ibv,2) + brick_list(nin,ib)%bakMAT%UVAR(2)
1902 uvarl(itask,ibv,3) = uvarl(itask,ibv,3) + brick_list(nin,ib)%bakMAT%UVAR(3)
1903c print *, "adding 1 " , BRICK_LIST(NIN,IB)%bakMAT%UVAR(1), UVARL(ITASK,IBv,1)
1904c print *, "adding 2 " , BRICK_LIST(NIN,IB)%bakMAT%UVAR(2), UVARL(ITASK,IBv,2)
1905c print *, "adding 3 " , BRICK_LIST(NIN,IB)%bakMAT%UVAR(3), UVARL(ITASK,IBv,3)
1906 DO itrimat=1,trimat
1907 var = ratio*vol51(ib,itrimat)
1908 IF(var==zero)cycle
1909 DO ipos = 1,m51_nvphas
1910 IF(ipos==11)cycle !volume already treated merged during geomeetrical merging.
1911 k = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
1912 pvar => brick_list(nin,ib)%bakMAT%UVAR(k)
1913 uvarl(itask,ibv,k) = uvarl(itask,ibv,k) + pvar * var
1914c print *, "adding j " , BRICK_LIST(NIN,IB)%bakMAT%UVAR(K) ,VAR, UVARL(ITASK,IBv,K)
1915 enddo!next IPOS
1916 enddo!next ITRIMAT
1917 endif!(MTN_==51)
1918 !------------------------------------------------------------!
1919 enddo! next ITAR
1920 !----END-LOOP------!
1921 enddo!next IB
1922
1923 !--------------!
1924 CALL my_barrier
1925 !--------------!
1926
1927
1928 if(ibug22_merge/=0 .and.dt1>zero)then
1929 if(itask==0)then
1930 do ib=1,nb
1931 cod1=nint(tmp22array(1,ib))
1932 cod2=nint(tmp22array(2,ib))
1933 if(cod1<=0)cycle
1934 write(*,fmt='(A,I8,A,I8)') " MERGING : id=", cod1, " to idv:", cod2
1935 write(*,fmt='(A,E30.16)') " RATIO=", tmp22array(3,ib)
1936 write(*,fmt='(A,E30.16)') " +rho=", tmp22array(4,ib)
1937 write(*,fmt='(A,E30.16)') " +rhoUx=", tmp22array(5,ib)
1938 write(*,fmt='(a,e30.16)') " +Vol=", TMP22ARRAY(6,IB)
1939 TMP22ARRAY(:,IB) = ZERO
1940 enddo
1941 endif
1942 call my_barrier !otherwise TMP22ARRAY IS MODIFIED
1943 endif
1944
1945
1946 !OPENMP STACKING
1947
1948! IF(IMERGEL(ITASK)==0)THEN
1949! NBL1=0
1950! ELSE
1951! IF(TRIMAT<0)THEN
1952! ALLOCATE(UVAR (I22LAW37))
1953! ELSE
1954! ALLOCATE(UVAR (M51_N0PHAS+TRIMAT*M51_NVPHAS))
1955! ENDIF
1956! ENDIF
1957
1958 DO IB=NBF,NBL1
1959 NBCUT = BRICK_LIST(NIN,IB)%NBCUT
1960 NCELL = NBCUT
1961 ICELL = 0
1962 GBUF => ELBUF_TAB(NGB(IB))%GBUF
1963 LBUF => ELBUF_TAB(NGB(IB))%BUFLY(1)%LBUF(1,1,1)
1964 !pAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%Adjacent_Cell(1:5)
1965 pAdjBRICK => BRICK_LIST(NIN,IB)%Adjacent_Brick(1:6,1:5)
1966 MBUF => ELBUF_TAB(NGB(IB))%BUFLY(1)%MAT(1,1,1)
1967 MLW = BRICK_LIST(NIN,IB)%MLW
1968 SOM = SUM(SuperCellVOL_L(0:NTHREAD-1,0,IB))
1969 LLT_ = IPARG(2,NGB(IB))
1970.OR. IF(MLW==37 MLW==51)THEN
1971 SOM_(1) = SUM(SuperCellVOL_L(0:NTHREAD-1,1,IB))
1972 SOM_(2) = SUM(SuperCellVOL_L(0:NTHREAD-1,2,IB))
1973 SOM_(3) = SUM(SuperCellVOL_L(0:NTHREAD-1,3,IB))
1974 SOM_(4) = SUM(SuperCellVOL_L(0:NTHREAD-1,4,IB))
1975 ENDIF
1976 IF(SOM==ZERO)CYCLE
1977 NG = BRICK_LIST(NIN,IB)%NG
1978 IDLOC = BRICK_LIST(NIN,IB)%IDLOC
1979 brickID = BRICK_LIST(NIN,IB)%ID
1980
1981 IF(BRICK_LIST(NIN,IB)%Ntarget>0)THEN
1982 DELTA = ZERO
1983 ELSE
1984 DELTA = ONE
1985 ENDIF
1986
1987 if (IBUG22_merge==-1)then
1988 !!!write(*,FMT='(a,e30.16)') " TARGET=", IXS(11,brick_list(nin,ib)%id)
1989 !!!write(*,FMT='(a,e30.16)') " DELTA=", DELTA
1990 !!!write(*,FMT='(a,e30.16)') " rho=", BRICK_LIST(NIN,IB)%bakMAT%RHO
1991 !!!write(*,FMT='(a,3e30.16)') " +rhoUx=", BRICK_LIST(NIN,IB)%bakMAT%rhoU(1:3)
1992 !!!write(*,FMT='(a,e30.16)') " +Vol=", BRICK_LIST(NIN,IB)%Vold_SCell
1993 TMP22ARRAY(1,IB) = IXS(11,brick_list(nin,ib)%id)
1994 TMP22ARRAY(2,IB) = DELTA
1995 TMP22ARRAY(3,IB) = BRICK_LIST(NIN,IB)%bakMAT%RHO
1996 TMP22ARRAY(4:6,IB) = BRICK_LIST(NIN,IB)%bakMAT%rhoU(1:3)
1997 TMP22ARRAY(7,IB) = BRICK_LIST(NIN,IB)%Vold_SCell
1998 endif
1999
2000 SOM = SOM + DELTA * BRICK_LIST(NIN,IB)%Vold_SCell
2001 SuperCellVOL_L(0:NTHREAD-1,0,IB) = ZERO
2002
2003 EINT = DELTA*BRICK_LIST(NIN,IB)%bakMAT%rhoE*BRICK_LIST(NIN,IB)%Vold_SCell +SUM(EINT_L(0:NTHREAD-1,IB))
2004 EINT = EINT / SOM
2005 EINT_L(0:NTHREAD-1,IB) = ZERO !to avoid full reinit
2006
2007 RHO = DELTA*BRICK_LIST(NIN,IB)%bakMAT%RHO*BRICK_LIST(NIN,IB)%Vold_SCell +SUM(RHO_L(0:NTHREAD-1,IB))
2008 RHO = RHO / SOM
2009 RHO_L(0:NTHREAD-1,IB) = ZERO !to avoid full reinit
2010
2011 !MomX
2012 MOM(1) = DELTA*BRICK_LIST(NIN,IB)%bakMAT%rhoU(1)*BRICK_LIST(NIN,IB)%Vold_SCell+SUM(MOM_L(1,0:NTHREAD-1,IB))
2013 MOM(1) = MOM(1) / SOM
2014 MOM_L(1,0:NTHREAD-1,IB)= ZERO !to avoid full reinit
2015 !MomY
2016 MOM(2) = DELTA*BRICK_LIST(NIN,IB)%bakMAT%rhoU(2)*BRICK_LIST(NIN,IB)%Vold_SCell+SUM(MOM_L(2,0:NTHREAD-1,IB))
2017 MOM(2) = MOM(2) / SOM
2018 MOM_L(2,0:NTHREAD-1,IB)= ZERO !to avoid full reinit
2019 !MomZ
2020 MOM(3) = DELTA*BRICK_LIST(NIN,IB)%bakMAT%rhoU(3)*BRICK_LIST(NIN,IB)%Vold_SCell+SUM(MOM_L(3,0:NTHREAD-1,IB))
2021 MOM(3) = MOM(3) / SOM
2022 MOM_L(3,0:NTHREAD-1,IB)= ZERO !to avoid full reinit
2023
2024 DO J=1,NV46
2025 SIG(J) = DELTA*BRICK_LIST(NIN,IB)%bakMAT%SIG(J)*BRICK_LIST(NIN,IB)%Vold_SCell+SUM(SIG_L(0:NTHREAD-1,J,IB))
2026 SIG(J) = SIG(J) / SOM
2027 SIG_L(0:NTHREAD-1,J,IB)= ZERO !to avoid full reinit
2028 ENDDO
2029
2030 GBUF%EINT(IDLOC) = EINT
2031 GBUF%RHO(IDLOC) = RHO ; ALEFVM_Buffer%FCELL(4, BRICK_LIST(NIN,IB)%ID ) = RHO
2032 GBUF%MOM(LLT_*(1-1) + IDLOC) = MOM(1) ; ALEFVM_Buffer%FCELL(1, BRICK_LIST(NIN,IB)%ID ) = MOM(1)
2033 GBUF%MOM(LLT_*(2-1) + IDLOC) = MOM(2) ; ALEFVM_Buffer%FCELL(2, BRICK_LIST(NIN,IB)%ID ) = MOM(2)
2034 GBUF%MOM(LLT_*(3-1) + IDLOC) = MOM(3) ; ALEFVM_Buffer%FCELL(3, BRICK_LIST(NIN,IB)%ID ) = MOM(3)
2035 !FCELL(5,.) voir bloc conditionnel ci-dessous
2036
2037 DO J=1,NV46
2038 GBUF%SIG(LLT_*(J-1)+IDLOC) = SIG(J)
2039 ENDDO
2040
2041 ! Cod1=ixs(11,brick_list(nin,ib)%id)
2042.or. ! if (IBUG22_merge==Cod1 IBUG22_merge==-1)then
2043 ! write(*,FMT='(a,i8,a,i8)') " AFTER MERGING : id=", Cod1, " is going to be updated"
2044 ! write(*,FMT='(a,e30.16)') " rho=", RHO
2045 ! write(*,FMT='(a,e30.16)') " rho.Ux=", MOM(1)
2046 ! print *,"-------------------------"
2047 ! endif
2048
2049 MTN_ = IPARG(1,NGB(IB))
2050 IF(MTN_==37)THEN
2051 !UVAR(I,1) : massic percentage of liquid * global density (rho1*V1/V : it needs to give liquid mass multiplying by element volume in aleconve.F)
2052 !UVAR(I,2) : density of gas
2053 !UVAR(I,3) : density of liquid
2054 !UVAR(I,4) : volumetric fraction of liquid
2055 !UVAR(I,5) : volumetric fraction of gas
2056 LLT_ = IPARG(2,NGB(IB))
2057
2058 SOM_(1) = SOM_(1) + DELTA*BRICK_LIST(NIN,IB)%bakMAT%UVAR(4)*BRICK_LIST(NIN,IB)%Vold_SCell
2059 SOM_(2) = SOM_(2) + DELTA*BRICK_LIST(NIN,IB)%bakMAT%UVAR(5)*BRICK_LIST(NIN,IB)%Vold_SCell
2060
2061
2062 UVAR(5) = DELTA*BRICK_LIST(NIN,IB)%bakMAT%UVAR(5)*BRICK_LIST(NIN,IB)%Vold_SCell + SUM(UVARL(0:NTHREAD-1,IB,5))
2063 UVAR(4) = DELTA*BRICK_LIST(NIN,IB)%bakMAT%UVAR(4)*BRICK_LIST(NIN,IB)%Vold_SCell + SUM(UVARL(0:NTHREAD-1,IB,4))
2064 UVAR(3) = DELTA*BRICK_LIST(NIN,IB)%bakMAT%UVAR(3)*BRICK_LIST(NIN,IB)%bakMAT%UVAR(4)*BRICK_LIST(NIN,IB)%Vold_SCell
2065 UVAR(3) = UVAR(3) + SUM(UVARL(0:NTHREAD-1,IB,3))
2066 UVAR(2) = DELTA*BRICK_LIST(NIN,IB)%bakMAT%UVAR(2)*BRICK_LIST(NIN,IB)%bakMAT%UVAR(5)*BRICK_LIST(NIN,IB)%Vold_SCell
2067 UVAR(2) = UVAR(2) + SUM(UVARL(0:NTHREAD-1,IB,2))
2068 UVAR(1) = DELTA*BRICK_LIST(NIN,IB)%bakMAT%UVAR(1)*BRICK_LIST(NIN,IB)%Vold_SCell + SUM(UVARL(0:NTHREAD-1,IB,1))
2069 UVARL(0:NTHREAD-1,IB,1:5) = ZERO
2070
2071c print *, "merge uvar-5", BRICK_LIST(NIN,IB)%bakMAT%UVAR(5),UVAR(1) / SOM
2072c print *, "merge uvar-4", BRICK_LIST(NIN,IB)%bakMAT%UVAR(4),UVAR(2) / SOM
2073c print *, "merge uvar-3", BRICK_LIST(NIN,IB)%bakMAT%UVAR(3),UVAR(3) / SOM
2074c print *, "merge uvar-2", BRICK_LIST(NIN,IB)%bakMAT%UVAR(2),UVAR(4) / SOM
2075c print *, "merge uvar-1", BRICK_LIST(NIN,IB)%bakMAT%UVAR(1),UVAR(5) / SOM
2076c print *, ""
2077
2078 MBUF%VAR((0*LLT_ + IDLOC)) = UVAR(1) / SOM !SOM=VOL
2079 IF(SOM_(2)>EM20)THEN
2080 MBUF%VAR((1*LLT_ + IDLOC)) = UVAR(2) / SOM_(2) !SOM2 = VOL_G =sum (vol_g , cell)
2081 MBUF%VAR((4*LLT_ + IDLOC)) = UVAR(5) / SOM
2082 ENDIF
2083 IF(SOM_(1)>EM20)THEN
2084 MBUF%VAR((2*LLT_ + IDLOC)) = UVAR(3) / SOM_(1) !SOM3 = VOL_L =sum (vol_l , cell)
2085 MBUF%VAR((3*LLT_ + IDLOC)) = UVAR(4) / SOM
2086 ENDIF
2087 SuperCellVOL_L(0:NTHREAD-1,1:4,IB) = ZERO
2088
2089
2090 !APPEL SIGEPS37 POUR EQUILIBRE PRESSION
2091 brickID = BRICK_LIST(NIN,IB)%ID
2092 IDLOC = IDLOCB(IB)
2093 NG = NGB(IB)
2094 VOL = SOM
2095 CALL SIGEPS37_SINGLE_CELL (
2096 1 ELBUF_TAB, IXS, BUFMAT, IPARG, IPM,
2097 2 IDLOC , NG , brickID, VOL
2098 . )
2099
2100
2101
2102 ELSEIF(MTN_==51)THEN
2103 LLT_ = IPARG(2,NGB(IB))
2104 UVAR(1) = DELTA*BRICK_LIST(NIN,IB)%bakMAT%UVAR(1) + SUM(UVARL(0:NTHREAD-1,IB,1))
2105 UVAR(2) = DELTA*BRICK_LIST(NIN,IB)%bakMAT%UVAR(2) + SUM(UVARL(0:NTHREAD-1,IB,2))
2106 UVAR(3) = DELTA*BRICK_LIST(NIN,IB)%bakMAT%UVAR(3) + SUM(UVARL(0:NTHREAD-1,IB,3))
2107 UVARL(0:NTHREAD-1,IB,1:3) = ZERO
2108 MBUF%VAR((0*LLT_ + IDLOC))= UVAR(1) / SOM
2109 MBUF%VAR((1*LLT_ + IDLOC))= UVAR(2) / SOM
2110 MBUF%VAR((2*LLT_ + IDLOC))= UVAR(3) / SOM
2111 DO ITRIMAT=1,TRIMAT
2112 SOMi = SUM(SuperCellVOL_L(0:NTHREAD-1,ITRIMAT,IB))
2113 SuperCellVOL_L(0:NTHREAD-1,ITRIMAT,IB) = ZERO
2114 IPOS = 11
2115 K = ((M51_N0PHAS + (ITRIMAT-1)*M51_NVPHAS )+IPOS)
2116 SOMi = SOMi + DELTA*BRICK_LIST(NIN,IB)%bakMAT%UVAR(K)
2117 VOL51(IB,ITRIMAT) = BRICK_LIST(NIN,IB)%bakMAT%UVAR(K)
2118 UVAR(:) = ZERO
2119 IF(SOMi>ZERO)THEN
2120 DO IPOS=1,M51_NVPHAS
2121 IF(IPOS==11)THEN
2122 K1 = ((M51_N0PHAS + (ITRIMAT-1)*M51_NVPHAS )+IPOS-1) ! ici UVAR(1:M51_N0PHAS+TRIMAT*M51_NVPHAS)
2123 K = K1 * LLT_
2124 MBUF%VAR(K+IDLOCB(IB)) = SOMi
2125 CYCLE
2126 ELSEIF(IPOS==1)THEN
2127 K1 = ((M51_N0PHAS + (ITRIMAT-1)*M51_NVPHAS )+IPOS-1) ! ici UVAR(1:M51_N0PHAS+TRIMAT*M51_NVPHAS)
2128 K = K1 * LLT_
2129 MBUF%VAR(K+IDLOCB(IB)) = SOMi / SOM
2130 CYCLE
2131 ENDIF
2132 K = ((M51_N0PHAS + (ITRIMAT-1)*M51_NVPHAS )+IPOS)
2133 K1 = (K-1)*LLT_
2134 pVAR => MBUF%VAR(K1+IDLOCB(IB))
2135 UVAR(K) = SUM(UVARL(0:NTHREAD-1,IB,K)) + pVAR * VOL51(IB,ITRIMAT) * DELTA
2136 UVARL(0:NTHREAD-1,IB,K) = ZERO
2137 UVAR(K) = UVAR(K) / SOMi
2138 IF(IPOS/=1)pVar = UVAR(K)
2139 ENDDO!next IPOS
2140 ELSE
2141 !vfrac_i=0 Vol_i=0
2142 IPOS = 1
2143 K1 = ((M51_N0PHAS + (ITRIMAT-1)*M51_NVPHAS )+IPOS-1)
2144 K = K1 * LLT_
2145 MBUF%VAR(K+IDLOCB(IB)) = ZERO
2146 IPOS = 11
2147 K1 = ((M51_N0PHAS + (ITRIMAT-1)*M51_NVPHAS )+IPOS-1)
2148 K = K1 * LLT_
2149 MBUF%VAR(K+IDLOCB(IB)) = ZERO
2150 ENDIF!(SOMi>ZERO)
2151 ENDDO!next ITRIMAT
2152
2153 ENDIF! (MTN_==51)
2154
2155 BRICK_LIST(NIN,IB)%Vold_Scell= SOM
2156
2157
2158 ENDDO!next IB
2159
2160
2161
2162
2163 if(IBUG22_merge/=0 )then
2164 call my_barrier
2165 if(itask==0)then
2166 if(DT1>ZERO)then
2167 do ib=1,nb
2168 if (TMP22ARRAY(1,IB)==zero )cycle
2169 write(*,FMT='(a,e30.16)') " TARGET=", NINT(TMP22ARRAY(1,IB))
2170 write(*,FMT='(a,e30.16)') " DELTA=", TMP22ARRAY(2,IB)
2171 write(*,FMT='(a,e30.16)') " rho=", TMP22ARRAY(3,IB)
2172 write(*,FMT='(a,3e30.16)') " +rhoUx=", TMP22ARRAY(4,IB)
2173 write(*,FMT='(a,e30.16)') " +Vol=", TMP22ARRAY(5,IB)
2174 enddo
2175 endif!if(DT1>ZERO)
2176 !print *, "deallocate"
2177 deallocate(TMP22ARRAY)
2178 endif
2179 endif
2180
2181 !------------------------------------------------------------!
2182 ! @26. DEMERGING CRITERIA !
2183 ! (IF main HAS TO BE DEPORTED THEN !
2184 ! DEMERGING WILL OCCURS IN NEW main CELL) !
2185 ! STEP 1/2 - listing origin cell !
2186 !------------------------------------------------------------!
2187 ! if Ntar > 0 then material deporation has to be done from former main cell buffer (GBUF)
2188 ! to the N=Ntar target main cells.
2189 NBL1 = NBL
2190 IF(DT1==ZERO)NBL1 = 0
2191 ORIGIN_DATA(:,:,:) = 0 !tag main cells linked to these secnd ones.
2192 Norigin(NBF:NBL) = 0
2193 DO IB=NBF,NBL1
2194 NBCUT = BRICK_LIST(NIN,IB)%NBCUT
2195 NCELL = NBCUT
2196 ICELL = 0
2197 MCELL = BRICK_LIST(NIN,IB)%mainID
2198 GBUF => ELBUF_TAB(NGB(IB))%GBUF
2199 pNodWasMain(1:8) => BRICK_LIST(NIN,IB)%NODE(1:8)!%NodWasMain
2200 pWhereWasMain(1:8) => BRICK_LIST(NIN,IB)%NODE(1:8)!%WhereWasMain
2201 brickID = BRICK_LIST(NIN,IB)%ID
2202 MBUF => ELBUF_TAB(NGB(IB))%BUFLY(1)%MAT(1,1,1)
2203 NTAG = 0
2204 MTN_ = IPARG(1,NGB(IB))
2205 NewMnod(1:8) = 0
2206 ITAG(:) = 0
2207 IF(MCELL == 0)CYCLE
2208 MNOD = BRICK_LIST(NIN,IB)%POLY(MCELL)%NumNOD
2209ccc print *, "demerging criteria for id=", ixs(11,brick_list(nin,ib)%id)
2210 !------CHECK NODES--------!
2211 DO K=1,MNOD
2212 INOD = BRICK_LIST(NIN,IB)%POLY(MCELL)%ListNodID(K)
2213 IF(pNodWasMain(INOD)%NodWasMain==1)CYCLE !node not concerned
2214 !Here node was secnd and just became main
2215 NewMnod(INOD) = 1
2216 ENDDO
2217
2218 !Continuity Criteria To avoid unexpected demerging.
2219 !Consistent only with kinematic time step "several cycle to fully cross an element".
2220 VOL_M = BRICK_LIST(NIN,IB)%secndlist%VOL_unmerged !BRICK_LIST(NIN,IB)%POLY(MCELL)vnew
2221 UncutVOL = BRICK_LIST(NIN,IB)%UncutVol
2222 VAR = VOL_M / UncutVOL
2223ccc print *, " VAR > 1-CritDVol22 : CYCLE", VAR > 1-CritDVol22, VAR, VOL_M
2224ccc print *, " MNOD==0) : CYCLE", MNOD==0
2225!not correct with sharp edges
2226! IF (VAR > 1-CritDVol22)CYCLE
2227 IF(MNOD==0) CYCLE
2228 CALL weighting_Cell_Nodes(NIN,IB,MCELL,INod_W,IDEMERGE)
2229
2230 ! IF(pNodWasMain(BRICK_LIST(NIN,IB)%OldMainStrongNode) == 1) IDEMERGE = 0
2231
2232 INod_W_old = BRICK_LIST(NIN,IB)%OldMainStrongNode
2233ccc print *, " INOD_W_old==0 : CYCLE", INOD_W_old==0
2234 IF(INOD_W_old<=0)CYCLE
2235 IF(BRICK_LIST(NIN,IB)%NODE(INod_W_old)%WhichCell == MCELL) IDEMERGE = 0
2236ccc print *, " IDEMERGE= : ", IDEMERGE
2237
2238 IF( IDEMERGE==1 ) THEN !(negatif en debug si IB nouveau dans le buffer)
2239! IF(pNodWasMain(INod_W)<=0 ) THEN !(negatif en debug si IB nouveau dans le buffer)
2240.AND.! IF(SUM(NewMnod)>=MNOD/2 MNOD/=0)THEN !DEMERGE criteria : all main nodes were secndo ones
2241 !get related main cell location (through its face)
2242 DO K=1,MNOD
2243 INOD = BRICK_LIST(NIN,IB)%POLY(MCELL)%ListNodID(K)
2244 J = pWhereWasMain(INOD)%WhereWasMain
2245 IF(J==0)CYCLE !means node was already main
2246 IF(ITAG(J)==1)CYCLE !list without doublon
2247 ITAG(J) = 1
2248 !Check if main cell was moved (merged)
2249 !Get main cut cell address
2250 IF(J>=10)THEN
2251 J1 = J/10
2252 J2 = MOD(J,10)
2253 IBv_i = BRICK_LIST(NIN,IB )%Adjacent_Brick(J1,4)
2254 IBv = BRICK_LISt(NIN,IBv_i)%Adjacent_Brick(J2,4)
2255 ELSE
2256 IBv = BRICK_LIST(NIN,IB )%Adjacent_Brick(J,4)
2257 ENDIF
2258 Ntar = BRICK_LIST(NIN,IBv)%Ntarget
2259 IF(Ntar == 0)THEN
2260 NTAG = NTAG + 1
2261 ORIGIN_DATA(IB,NTAG,1) = IBv
2262 ORIGIN_DATA(IB,NTAG,2) = J
2263 ORIGIN_DATA(IB,NTAG,3) = IBv
2264 ELSE
2265 DO Itar=1,NTar
2266 NTAG = NTAG + 1
2267 ORIGIN_DATA(IB,NTAG,1) = BRICK_LIST(NIN,IBv)%MergeTarget(3,ITAR)
2268 ORIGIN_DATA(IB,NTAG,2) = J
2269 ORIGIN_DATA(IB,NTAG,3) = IBv
2270 ENDDO
2271 ENDIF
2272 ENDDO
2273 Norigin(IB) = NTAG
2274 IF(NTAG==0)THEN
2275 print *, "**error : inter22, topology issue."
2276 ENDIF
2277 ENDIF
2278 ENDDO! next IB
2279
2280 !--------------!
2281 CALL MY_BARRIER
2282 !--------------!
2283
2284 !------------------------------------------------------------!
2285 ! @27. MATERIAL DEMERGING !
2286 ! STEP 2 : filling main cell with data from origin cells !
2287 !------------------------------------------------------------!
2288 ! if Ntar > 0 then material deporation has to be done from former main cell buffer (GBUF)
2289 ! to the N=Ntar target main cells.
2290 NBL1 = NBL
2291 IF(DT1==ZERO)NBL1 = 0
2292
2293 if(IBUG22_merge/=0)then
2294 if(itask==0)ALLOCATE (TMP22ARRAY(8,NB))
2295 IN22 = 0
2296 call my_barrier
2297 TMP22ARRAY(1:8,NBF:NBL1)=ZERO
2298 endif
2299
2300
2301 DO IB=NBF,NBL1
2302 IF(Norigin(IB)==0)CYCLE
2303 NBCUT = BRICK_LIST(NIN,IB)%NBCUT
2304 IF(NBCUT==0)CYCLE
2305 NCELL = NBCUT
2306 ICELL = 0
2307 MCELL = BRICK_LIST(NIN,IB)%mainID
2308 GBUF => ELBUF_TAB(NGB(IB))%GBUF
2309 LBUF => ELBUF_TAB(NGB(IB))%BUFLY(1)%LBUF(1,1,1)
2310 !pAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%Adjacent_Cell(1:5)
2311 pAdjBRICK => BRICK_LIST(NIN,IB)%Adjacent_Brick(1:6,1:3)
2312 pNodWasMain(1:8) => BRICK_LIST(NIN,IB)%NODE(1:8)!%NodWasMain
2313 pWhereWasMain(1:8) => BRICK_LIST(NIN,IB)%NODE(1:8)!%WhereWasMain
2314 pSUBVOL(1:9) => BRICK_LIST(NIN,IB)%POLY(1:9)
2315 brickID = BRICK_LIST(NIN,IB)%ID
2316 IDLOC = BRICK_LIST(NIN,IB)%IDLOC
2317 MTN_ = IPARG(1,NGB(IB))
2318 LLT_ = IPARG(2,NGB(IB))
2319 MBUF => ELBUF_TAB(NGB(IB))%BUFLY(1)%MAT(1,1,1)
2320 NTAG = 0
2321 VOLorig(1:24) = ZERO
2322 VOLcell = ZERO
2323 VOLcell51(:) = ZERO
2324 VolSECND51(:) = ZERO
2325 VOL = ZERO
2326 RHO = ZERO
2327 EINT = ZERO
2328 SIG(1:6) = ZERO
2329 MOM(1:3) = ZERO
2330 SSP = ZERO
2331 UVAR(:) = ZERO
2332 RHO_ADV = ZERO
2333 EINT_ADV = ZERO
2334 SIG_ADV(1:6) = ZERO
2335 MOM_ADV(1:3) = ZERO
2336 UVAR_ADV(:) = ZERO
2337
2338 !-----------------------------------------------!
2339 ! INIT CELL CONTENT IF IT CAME FROM A MERGE !
2340 !-----------------------------------------------!
2341 IF(Brick_list(nin,ib)%IsMergeTarget==1)THEN
2342 ! pause
2343 !!!print *, "ismergeTarget , Norigin", ixs(11,brickID), Norigin(IB), BRICK_LIST(NIN,IB)%Ntarget
2344 ! This cell received material from a merge.it also has to be taken into account
2345 VolCell = BRICK_LIST(NIN,IB)%Vold_Scell
2346 RHO = VolCell* GBUF%RHO(IDLOC)
2347 EINT = VolCell* GBUF%EINT(IDLOC)
2348 MOM(1) = VolCell* GBUF%MOM(LLT_*(1-1) + IDLOC)
2349 MOM(2) = VolCell* GBUF%MOM(LLT_*(2-1) + IDLOC)
2350 MOM(3) = VolCell* GBUF%MOM(LLT_*(3-1) + IDLOC)
2351 SIG(1) = VolCell* GBUF%SIG(LLT_*(1-1)+IDLOC)
2352 SIG(2) = VolCell* GBUF%SIG(LLT_*(2-1)+IDLOC)
2353 SIG(3) = VolCell* GBUF%SIG(LLT_*(3-1)+IDLOC)
2354 SIG(4) = VolCell* GBUF%SIG(LLT_*(4-1)+IDLOC)
2355 SIG(5) = VolCell* GBUF%SIG(LLT_*(5-1)+IDLOC)
2356 SIG(6) = VolCell* GBUF%SIG(LLT_*(6-1)+IDLOC)
2357 SSP = VolCell* LBUF%SSP(IDLOC)
2358 VolCell = ZERO
2359 IF(MTN_==37)THEN
2360 UVAR(5) = VolCell*MBUF%VAR(((5-1)*LLT_ + IDLOC))
2361 UVAR(4) = VolCell*MBUF%VAR(((4-1)*LLT_ + IDLOC))
2362 UVAR(3) = VolCell*MBUF%VAR(((3-1)*LLT_ + IDLOC))*MBUF%VAR(((4-1)*LLT_ + IDLOC))
2363 UVAR(2) = VolCell*MBUF%VAR(((2-1)*LLT_ + IDLOC))*MBUF%VAR(((5-1)*LLT_ + IDLOC))
2364 UVAR(1) = VolCell*MBUF%VAR(((1-1)*LLT_ + IDLOC))
2365 ELSEIF(MTN_==51)THEN
2366 IPOS = 11
2367 volSECND51(1) = MBUF%VAR(IDLOC + ((M51_N0PHAS + (1-1)*M51_NVPHAS )+IPOS-1)*LLT_)
2368 volSECND51(2) = MBUF%VAR(IDLOC + ((M51_N0PHAS + (2-1)*M51_NVPHAS )+IPOS-1)*LLT_)
2369 volSECND51(3) = MBUF%VAR(IDLOC + ((M51_N0PHAS + (3-1)*M51_NVPHAS )+IPOS-1)*LLT_)
2370 volSECND51(4) = MBUF%VAR(IDLOC + ((M51_N0PHAS + (4-1)*M51_NVPHAS )+IPOS-1)*LLT_)
2371 DO ITRIMAT=1,TRIMAT
2372 DO IPOS = 1 , M51_NVPHAS
2373 IF(IPOS==11) CYCLE
2374 K = ((M51_N0PHAS + (ITRIMAT-1)*M51_NVPHAS )+IPOS-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
2375 K1 = K * LLT_
2376 UVAR(K+1) = volSECND51(ITRIMAT) * MBUF%VAR(K1+IDLOC)
2377 ENDDO!next IPOS
2378 ENDDO!next ITRIMAT
2379 volSECND51(1:4) = ZERO
2380 ENDIF
2381 ENDIF
2382
2383 ! Velocity of main is velocity wherever inside the supercell. It will be used to make an advection step get pressure gradient which did not occur on previous step due to homogeneity.
2384 VEL_M(1) = GBUF%MOM(3*(IDLOC-1)+1) / GBUF%RHO(IDLOC)
2385 VEL_M(2) = GBUF%MOM(3*(IDLOC-1)+2) / GBUF%RHO(IDLOC)
2386 VEL_M(3) = GBUF%MOM(3*(IDLOC-1)+3) / GBUF%RHO(IDLOC)
2387
2388 NTAR = BRICK_LIST(NIN,IB)%Ntarget
2389 !-----------------------------------------------!
2390 ! LOOP ON ORIGIN CELLS !
2391 !-----------------------------------------------!
2392 DO IORIG = 1, Norigin(IB)
2393 !IORIG data regarding IB
2394 IBm = ORIGIN_DATA(IB,IORIG,1) !real target
2395 IBo = ORIGIN_DATA(IB,IORIG,3) !original brick before knowing it was burged below
2396 IV = BRICK_LIST(NIN,IBm)%ID
2397 NGv = BRICK_LIST(NIN,IBm)%NG
2398 IDLOCv = BRICK_LIST(NIN,IBm)%IDLOC
2399 J = ORIGIN_DATA(IB,IORIG,2)
2400 IF(IBm==IB)CYCLE
2401 IF(J>=10)THEN
2402 J1 = J/10
2403 J2 = MOD(J,10)
2404 IBv = BRICK_LIST(NIN,IB )%Adjacent_Brick(J1,4)
2405 JV = BRICK_LIST(NIN,IBv)%Adjacent_Brick(J2,5)
2406 ELSE
2407 JV = BRICK_LIST(NIN,IB)%Adjacent_Brick(J,5)
2408 ENDIF
2409
2410 NumSECND = OLD_SecndList(NIN,IBo)%Num !this is list from previous cycle
2411 !-----------------------------------------------!
2412 ! GET SECND CONNECTED TO IT !
2413 !-----------------------------------------------!
2414 DO K=1,NumSECND
2415 IF (OLD_SecndList(NIN,IBo)%IBv(K) /= IB)CYCLE
2416 IF (OLD_SecndList(NIN,IBo)%FM(K) /= JV)CYCLE !next secnd : this one is not connecter through main face id %FM
2417
2418 !---debug---!
2419 cod1 = ixs(11,brick_list(nin,ib)%id)
2420 cod2 = ixs(11,brick_list(nin,ibo)%id)
2421 if (IBUG22_merge==-1)then
2422#include "lockon.inc"
2423 if(IBm==ibo)then
2424 ! write(*,FMT='(a,i8,a,i8)') "DEMERGING : id=", Cod1,
2425 ! . " from idv:", ixs(11,brick_list(nin,ibo)%id)
2426
2427 IN22 = IN22 + 1
2428 tmp22array(1,IN22)=Cod1
2429 tmp22array(2,IN22)=ixs(11,brick_list(nin,ibo)%id)
2430 tmp22array(3,IN22)=0
2431
2432 else
2433 ! write(*,FMT='(a,i8,a,i8,a,i8)') "DEMERGING : id=", Cod1,
2434 ! . " from idv:", ixs(11,brick_list(nin,ibo)%id)," moved to", ixs(11,brick_list(nin,IBm)%id)
2435
2436 IN22 = IN22 + 1
2437 tmp22array(1,IN22)=Cod1
2438 tmp22array(2,IN22)=ixs(11,brick_list(nin,ibo)%id)
2439 tmp22array(3,IN22)=ixs(11,brick_list(nin,IBm)%id)
2440
2441 endif
2442 ! write (*,FMT='(i10,a,f30.16,a,f30.16)')
2443 ! . Cod1," Vold=",brick_list(nin,ib)%Vold_scell, " Vnew=",brick_list(nin,ib)%Vnew_scell
2444 ! write (*,FMT='(i10,a,f30.16,a,f30.16)')
2445 ! . Cod2," Vold=",brick_list(nin,ibo)%Vold_scell," Vnew=",brick_list(nin,ibo)%Vnew_scell
2446
2447 tmp22array(4,IN22)= zero !brick_list(nin,ib)%Vold_scell !value not consistent, old main was merge elsewhere so scell has change of side since a demerge process is emerging a new scell
2448 tmp22array(5,IN22)= brick_list(nin,ib)%Vnew_scell
2449 tmp22array(6,IN22)= zero !brick_list(nin,ibo)%Vold_scell
2450 tmp22array(7,IN22)= brick_list(nin,ibo)%Vnew_scell
2451 tmp22array(8,IN22)=tmp22array(8,IN22)+1
2452#include "lockoff.inc"
2453 endif
2454 !-----------!
2455
2456 !here secnd is merged on the correct face, and also to the correct target main (IB)
2457 !!!!!!!!!!!!! VOLsecnd = OLD_SecndList(NIN,IBo)%VOL(K)
2458 !-----------------------------------------------!
2459 ! START TO CALCULATE SECND NEW QUANTITIES !
2460 !-----------------------------------------------!
2461 VOLsecnd = OLD_SecndList(NIN,IBo)%VOL(K)
2462 !!! VolsecndOLD = OLD_SecndList(NIN,IBo)%VOL(K)
2463 VOLcell = VOLcell + VOLsecnd !cumulative salve volume from previous cycle = old volume for new main cell
2464 !!! VOLcellOLD = VOLcell + VOLsecndOLD
2465 GBUF => ELBUF_TAB(NGv)%GBUF
2466 LBUF => ELBUF_TAB(NGv)%BUFLY(1)%LBUF(1,1,1)
2467 LLT_v = IPARG(2,NGv)
2468 !secnd cell surf
2469 SURF_S = OLD_SecndList(NIN,IBo)%SURF_V(K)
2470 FV = OLD_SecndList(NIN,IBo)%FV(K)
2471 IF(FV<10)THEN
2472 NORM_S(1:3) = BRICK_LIST(NIN,IB)%N(FV,1:3)
2473 ELSE
2474 !FV = FV/10
2475 !NORM_S(1:3) = BRICK_LIST(NIN,IB)%N(FV,1:3)
2476 !IL Y A PLUSIEURS FACETTES POSSIBLE POUR LES VOISINS INDIRECT. IL FAUDRAIT SAUVEGARGER J>10 pour chaque chemin, puis advection
2477 !ou alors changer la normale en prenant le vecteur centroid_M-centroid_S
2478 ZM(1:3) = BRICK_LIST(NIN,IBo)%SCellCenter(1:3)
2479 ZS(1:3) = BRICK_LIST(NIN,IB)%POLY(MCELL)%CellCenter(1:3)
2480 NORM_S(1) = ZM(1)-ZS(1)
2481 NORM_S(2) = ZM(2)-ZS(2)
2482 NORM_S(3) = ZM(3)-ZS(3)
2483 NORM = SQRT(NORM_S(1)*NORM_S(1)+NORM_S(2)*NORM_S(2)+NORM_S(3)*NORM_S(3))
2484 NORM_S(1) = NORM_S(1) / NORM
2485 NORM_S(2) = NORM_S(2) / NORM
2486 NORM_S(3) = NORM_S(3) / NORM
2487 ENDIF
2488
2489
2490 !----------------------------------!
2491 ! ADVECTION STEP IN DEMERGED CELL !
2492 !----------------------------------!
2493 ADV = ZERO
2494 IF(IADV==1)THEN
2495 !MCELL = BRICK_LIST(NIN,IB)%MainID
2496 VM = ONE ! BRICK_LIST(NIN,IBm)%POLY(MCELL)%Vnew
2497 ADV = -HALF* HALF* DT1*VM*SURF_S* (VEL_M(1)*NORM_S(1)+VEL_M(2)*NORM_S(2)+VEL_M(3)*NORM_S(3))
2498 RHO_ADV = RHO_ADV + ADV * GBUF%RHO(IDLOCv)
2499 EINT_ADV = EINT + ADV * GBUF%EINT(IDLOCv)
2500 MOM_ADV(1) = MOM_ADV(1) + ADV * GBUF%MOM(3*(IDLOCv-1)+1)
2501 MOM_ADV(2) = MOM_ADV(2) + ADV * GBUF%MOM(3*(IDLOCv-1)+2)
2502 MOM_ADV(3) = MOM_ADV(3) + ADV * GBUF%MOM(3*(IDLOCv-1)+3)
2503 SIG_ADV(1) = SIG_ADV(1) + ADV * GBUF%SIG(LLT_v*(1-1)+IDLOCv)
2504 SIG_ADV(2) = SIG_ADV(2) + ADV * GBUF%SIG(LLT_v*(2-1)+IDLOCv)
2505 SIG_ADV(3) = SIG_ADV(3) + ADV * GBUF%SIG(LLT_v*(3-1)+IDLOCv)
2506 SIG_ADV(4) = SIG_ADV(4) + ADV * GBUF%SIG(LLT_v*(4-1)+IDLOCv)
2507 SIG_ADV(5) = SIG_ADV(5) + ADV * GBUF%SIG(LLT_v*(5-1)+IDLOCv)
2508 SIG_ADV(6) = SIG_ADV(6) + ADV * GBUF%SIG(LLT_v*(6-1)+IDLOCv)
2509 ENDIF
2510 !-------------------------------!
2511 ! MONO-MATERIAL CUMULATION !
2512 !-------------------------------!
2513 RHO = RHO + VOLsecnd * GBUF%RHO(IDLOCv)
2514 EINT = EINT + VOLsecnd * GBUF%EINT(IDLOCv)
2515 MOM(1) = MOM(1) + VOLsecnd * GBUF%MOM(LLT_v*(1-1)+IDLOCv)
2516 MOM(2) = MOM(2) + VOLsecnd * GBUF%MOM(LLT_v*(2-1)+IDLOCv)
2517 MOM(3) = MOM(3) + VOLsecnd * GBUF%MOM(LLT_v*(3-1)+IDLOCv)
2518 SIG(1) = SIG(1) + VOLsecnd * GBUF%SIG(LLT_v*(1-1)+IDLOCv)
2519 SIG(2) = SIG(2) + VOLsecnd * GBUF%SIG(LLT_v*(2-1)+IDLOCv)
2520 SIG(3) = SIG(3) + VOLsecnd * GBUF%SIG(LLT_v*(3-1)+IDLOCv)
2521 SIG(4) = SIG(4) + VOLsecnd * GBUF%SIG(LLT_v*(4-1)+IDLOCv)
2522 SIG(5) = SIG(5) + VOLsecnd * GBUF%SIG(LLT_v*(5-1)+IDLOCv)
2523 SIG(6) = SIG(6) + VOLsecnd * GBUF%SIG(LLT_v*(6-1)+IDLOCv)
2524 SSP = SSP + VOLsecnd * LBUF%SSP(IDLOCv)
2525 !-------------------------------!
2526 ! MULTI-MATERIAL CUMULULATION !
2527 !-------------------------------!
2528 IF(MTN_==37)THEN
2529 !UVAR(I,1) : massic percentage of liquid * global density (rho1*V1/V : it needs to give liquid mass multiplying by element volume in aleconve.F)
2530 !UVAR(I,2) : density of gas
2531 !UVAR(I,3) : density of liquid
2532 !UVAR(I,4) : volumetric fraction of liquid
2533 !UVAR(I,5) : volumetric fraction of gas
2534 MBUFv => ELBUF_TAB(NGv)%BUFLY(1)%MAT(1,1,1)
2535 LLT_v = IPARG(2,NGv)
2536 UVAR(5) = UVAR(5) +VOLsecnd*MBUFv%VAR(((5-1)*LLT_v + IDLOCv))
2537 UVAR(4) = UVAR(4) +VOLsecnd*MBUFv%VAR(((4-1)*LLT_v + IDLOCv))
2538 UVAR(3) = UVAR(3) +VOLsecnd*MBUFv%VAR(((3-1)*LLT_v + IDLOCv))*MBUFv%VAR(((4-1)*LLT_v + IDLOCv))
2539 UVAR(2) = UVAR(2) +VOLsecnd*MBUFv%VAR(((2-1)*LLT_v + IDLOCv))*MBUFv%VAR(((5-1)*LLT_v + IDLOCv))
2540 UVAR(1) = UVAR(1) +VOLsecnd*MBUFv%VAR(((1-1)*LLT_v + IDLOCv))
2541 UVAR_ADV(1) = UVAR_ADV(1) +ADV * MBUFv%VAR(((1-1)*LLT_v + IDLOCv))
2542
2543 ELSEIF(MTN_==51)THEN
2544 MBUFv => ELBUF_TAB(NGv)%BUFLY(1)%MAT(1,1,1)
2545 LLT_ = IPARG(2,NGB(IB))
2546 LLT_v = IPARG(2,NGv)
2547 !in sigesp51.F:
2548 !V1OLD = V1OLD - TIMESTEP*DDVOL1
2549 !Volume change correspond to DDVOL=DV/DT = f(fluxes). Gird deformation influences relative material velocity ans also fluxes
2550 DO ITRIMAT=1,TRIMAT
2551 !get,store Volumetric Fraction
2552 IPOS = 1
2553 K1 = ((M51_N0PHAS + (ITRIMAT-1)*M51_NVPHAS )+IPOS-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
2554 Kv = K1 * LLT_v
2555 Vfrac(ITRIMAT) = MBUFv%VAR(Kv+IDLOCv)
2556 VolSECND51(ITRIMAT) = Vfrac(ITRIMAT) * VOLsecnd !submaterial volume in secnd cell
2557 VolCELL51(ITRIMAT) = VolCELL51(ITRIMAT) + VolSECND51(ITRIMAT) !submaterial volume in ALL secnd cell
2558 ENDDO!next ITRIMAT
2559 UVAR(1) = UVAR(1) + MBUFv%VAR((0*LLT_v + IDLOCv)) * VOLsecnd
2560 UVAR(2) = UVAR(2) + MBUFv%VAR((1*LLT_v + IDLOCv)) * VOLsecnd
2561 UVAR(3) = UVAR(3) + MBUFv%VAR((2*LLT_v + IDLOCv)) * VOLsecnd
2562 UVAR_ADV(1) = UVAR_ADV(1) + MBUFv%VAR((0*LLT_v + IDLOCv)) * ADV
2563 UVAR_ADV(2) = UVAR_ADV(2) + MBUFv%VAR((1*LLT_v + IDLOCv)) * ADV
2564 UVAR_ADV(3) = UVAR_ADV(3) + MBUFv%VAR((2*LLT_v + IDLOCv)) * ADV
2565 IF(IADV==1)print *, "to do compute/get vfrac on face" !; pause
2566 DO ITRIMAT=1,TRIMAT
2567 DO IPOS = 1 , M51_NVPHAS
2568 Ko = ((M51_N0PHAS + (ITRIMAT-1)*M51_NVPHAS ))
2569 K1 = ((M51_N0PHAS + (ITRIMAT-1)*M51_NVPHAS )+IPOS-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
2570 Kv = K1 * LLT_v
2571 IF(IPOS==11)THEN
2572 UVAR(K1+1) = UVAR(K1+1) + volSECND51(ITRIMAT)
2573 ELSE
2574 UVAR(K1+1) = UVAR(K1+1) + volSECND51(ITRIMAT) * MBUFv%VAR(Kv+IDLOCv)
2575 UVAR_ADV(K1+1) = UVAR_ADV(K1+1) + ADV * MBUFv%VAR(Kv+IDLOCv)
2576 IF(IADV==1)print *, "to do compute/get vfrac on face" !; pause
2577 ENDIF
2578 ENDDO!next IPOS
2579 ENDDO!next ITRIMAT
2580 ENDIF!(MTN_==51)
2581 !---------------------------
2582#include "lockon.inc"
2583 !lock method : to be update later can be removed
2584 !SUBSTRACTING SECND VOLUME FROM SUPERCELL TO DEMERGE
2585 !--MONOMAT
2586 BRICK_LIST(NIN,IBm)%Vold_Scell = BRICK_LIST(NIN,IBm)%Vold_Scell - VOLsecnd
2587 !--MULTIMAT
2588 IF(MTN_==51)THEN
2589 DO ITRIMAT=1,TRIMAT
2590 !remove submaterial volume in detached secnd cell
2591 IPOS = 11
2592 K1 = ((M51_N0PHAS + (ITRIMAT-1)*M51_NVPHAS )+IPOS-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
2593 Kv = K1 * LLT_v
2594 MBUFv%VAR(Kv+IDLOCv) = MBUFv%VAR(Kv+IDLOCv) - VolSECND51(ITRIMAT)
2595 ENDDO!next ITRIMAT
2596 ENDIF
2597#include "lockoff.inc"
2598 ENDDO!next K
2599 ENDDO!next IORIG
2600
2601 IF(VOLcell==ZERO)CYCLE
2602
2603 !-----------------------------------------------!
2604 ! GET VALUE FROM CUMULATIVE ONE !
2605 !---MONOMAT-------------------------------------!
2606 RHO = (RHO +IADV*RHO_ADV )/ VOLcell
2607 EINT = (EINT +IADV*EINT_ADV )/ VOLcell
2608 MOM(1) = (MOM(1) +IADV*MOM_ADV(1))/ VOLcell
2609 MOM(2) = (MOM(2) +IADV*MOM_ADV(2))/ VOLcell
2610 MOM(3) = (MOM(3) +IADV*MOM_ADV(3))/ VOLcell
2611 SIG(:) = (SIG(:) +IADV*SIG_ADV(:))/ VOLcell
2612 SSP = SSP / VOLcell ! -> run EOS
2613 !---MULTIMAT -----------------------------------!
2614 IF(MTN_==37)THEN
2615 UVAR(1) = (UVAR(1)+ IADV*UVAR_ADV(1)) / VOLcell
2616 UVAR(2) = UVAR(2) / VOLcell
2617 UVAR(3) = UVAR(3) / VOLcell
2618 UVAR(4) = UVAR(4) / VOLcell
2619 UVAR(5) = UVAR(5) / VOLcell
2620 !call sigeps37_single_cell before storing in MBUF%VAR
2621 ELSEIF(MTN_==51)THEN
2622 UVAR(1) = (UVAR(1) + IADV*UVAR_ADV(1)) / VOLcell !sum of full volume because UVAR(1:3) are homogenized data, not submaterial one.
2623 UVAR(2) = (UVAR(2) + IADV*UVAR_ADV(2)) / VOLcell
2624 UVAR(3) = (UVAR(3) + IADV*UVAR_ADV(3)) / VOLcell
2625 DO ITRIMAT=1,TRIMAT
2626 Ko = M51_N0PHAS + (ITRIMAT-1)*M51_NVPHAS
2627 IF(VOLcell51(ITRIMAT)/=ZERO)THEN
2628 UVAR(Ko+01:Ko+10) = (UVAR(Ko+01:Ko+10) + IADV*UVAR_ADV(Ko+01:Ko+10)) / VOLcell51(ITRIMAT) !submaterial data are averaged with sum of submaterial secnd volume
2629 UVAR(Ko+12:Ko+M51_NVPHAS) = (UVAR(Ko+12:Ko+M51_NVPHAS) + IADV*UVAR_ADV(Ko+12:Ko+M51_NVPHAS)) / VOLcell51(ITRIMAT) !except ipos=11 (volume)
2630 ELSE
2631 UVAR(Ko+01) = ZERO
2632 UVAR(Ko+11) = ZERO
2633 ENDIF
2634 ENDDO
2635 ENDIF!(MTN_==51)
2636 !-----------------------------------------------!
2637
2638 !-----------------------------------------------!
2639 ! STORE VALUES IN DEMERGED CELLS !
2640 !---MONOMAT-------------------------------------!
2641 NG = BRICK_LIST(NIN,IB)%NG
2642 IDLOC = BRICK_LIST(NIN,IB)%IDLOC
2643 LLT_ = IPARG(2,NG)
2644 GBUF => ELBUF_TAB(NG)%GBUF
2645 LBUF => ELBUF_TAB(NG)%BUFLY(1)%LBUF(1,1,1)
2646 GBUF%RHO(IDLOC) = RHO
2647 GBUF%EINT(IDLOC) = EINT
2648 GBUF%MOM(LLT_*(1-1)+IDLOC) = MOM(1)
2649 GBUF%MOM(LLT_*(2-1)+IDLOC) = MOM(2)
2650 GBUF%MOM(LLT_*(3-1)+IDLOC) = MOM(3)
2651 GBUF%SIG(LLT_*(1-1)+IDLOC) = SIG(1)
2652 GBUF%SIG(LLT_*(2-1)+IDLOC) = SIG(2)
2653 GBUF%SIG(LLT_*(3-1)+IDLOC) = SIG(3)
2654 GBUF%SIG(LLT_*(4-1)+IDLOC) = SIG(4)
2655 GBUF%SIG(LLT_*(5-1)+IDLOC) = SIG(5)
2656 GBUF%SIG(LLT_*(6-1)+IDLOC) = SIG(6)
2657 LBUF%SSP(IDLOC) = SSP
2658 ALEFVM_Buffer%FCELL(4, BRICK_LIST(NIN,IB)%ID ) = RHO
2659 ALEFVM_Buffer%FCELL(1, BRICK_LIST(NIN,IB)%ID ) = MOM(1)
2660 ALEFVM_Buffer%FCELL(2, BRICK_LIST(NIN,IB)%ID ) = MOM(2)
2661 ALEFVM_Buffer%FCELL(3, BRICK_LIST(NIN,IB)%ID ) = MOM(3)
2662 ALEFVM_Buffer%FCELL(6, BRICK_LIST(NIN,IB)%ID ) = -THIRD*(SIG(1)+SIG(2)+SIG(3))
2663 ALEFVM_Buffer%FCELL(5, BRICK_LIST(NIN,IB)%ID ) = SSP
2664
2665 !---MUKLTIMAT-------------------------------------!
2666 IF(MTN_==37)THEN
2667 LLT_ = IPARG(2,NGB(IB))
2668 IDLOC = IDLOCB(IB)
2669 MBUF => ELBUF_TAB(NGB(IB))%BUFLY(1)%MAT(1,1,1)
2670 MT = IXS(1,brickID)
2671 IADBUF = IPM(7,MT)
2672 RHO10 = BUFMAT(IADBUF-1+11)
2673 RHO20 = BUFMAT(IADBUF-1+12)
2674c print *, ""
2675c print *, "demerge-5", MBUF%VAR((5-1)*LLT_+IDLOC) , UVAR(5)
2676c print *, "demerge-4", MBUF%VAR((4-1)*LLT_+IDLOC) , UVAR(4)
2677c print *, "demerge-3", MBUF%VAR((3-1)*LLT_+IDLOC) , UVAR(3)
2678c print *, "demerge-2", MBUF%VAR((2-1)*LLT_+IDLOC) , UVAR(2)
2679c print *, "demerge-1", MBUF%VAR((1-1)*LLT_+IDLOC) , UVAR(1)
2680 IF(UVAR(3) == ZERO) UVAR(3) = RHO10
2681 IF(UVAR(2) == ZERO) UVAR(2) = RHO20
2682 MBUF%VAR((5-1)*LLT_+IDLOC) = UVAR(5)
2683 MBUF%VAR((4-1)*LLT_+IDLOC) = UVAR(4)
2684 MBUF%VAR((3-1)*LLT_+IDLOC) = UVAR(3)
2685 MBUF%VAR((2-1)*LLT_+IDLOC) = UVAR(2)
2686 MBUF%VAR((1-1)*LLT_+IDLOC) = UVAR(1)
2687 ELSEIF(MTN_==51)THEN
2688 LLT_ = IPARG(2,NGB(IB))
2689 IDLOC = IDLOCB(IB)
2690 MBUF => ELBUF_TAB(NGB(IB))%BUFLY(1)%MAT(1,1,1)
2691 DO ITRIMAT=1,TRIMAT
2692 IF(VOLcell51(ITRIMAT)/=ZERO)THEN
2693 !get,store Volumetric Fraction
2694 DO IPOS=1,M51_NVPHAS
2695 K0 = ((M51_N0PHAS + (ITRIMAT-1)*M51_NVPHAS )+IPOS-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
2696 K1 = K0 * LLT_
2697 MBUF%VAR(K1+IDLOC) = UVAR(M51_N0PHAS+(ITRIMAT-1)*M51_NVPHAS+IPOS)
2698 ENDDO
2699
2700 IPOS=1
2701 K0 = ((M51_N0PHAS + (ITRIMAT-1)*M51_NVPHAS )+IPOS-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
2702 K1 = K0 * LLT_
2703 MBUF%VAR(K1+IDLOC) = VolCELL51(ITRIMAT)/VolCELL
2704 IPOS=11
2705 K0 = ((M51_N0PHAS + (ITRIMAT-1)*M51_NVPHAS )+IPOS-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
2706 K1 = K0 * LLT_
2707 MBUF%VAR(K1+IDLOC) = VolCELL51(ITRIMAT)
2708 ELSE
2709 IPOS=1
2710 K0 = ((M51_N0PHAS + (ITRIMAT-1)*M51_NVPHAS )+IPOS-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
2711 K1 = K0 * LLT_
2712 MBUF%VAR(K1+IDLOC) = ZERO
2713 IPOS=11
2714 K0 = ((M51_N0PHAS + (ITRIMAT-1)*M51_NVPHAS )+IPOS-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
2715 K1 = K0 * LLT_
2716 MBUF%VAR(K1+IDLOC) = ZERO
2717 ENDIF
2718 ENDDO
2719 ENDIF!(MTN_==51)
2720 !-----------------------------------------------!
2721
2722 !-----------------------------------------------!
2723 ! NEW main CELL UPDATE !
2724 !-----------------------------------------------!
2725 BRICK_LIST(NIN,IB)%Vold_SCell = VOLcell
2726
2727 !APPEL SIGEPS37 POUR EQUILIBRE PRESSION
2728 brickID = BRICK_LIST(NIN,IB)%ID
2729 IDLOC = IDLOCB(IB)
2730 NG = NGB(IB)
2731 VOL = VOLcell
2732 MTN_ = IPARG(1,NG)
2733 IF(MTN_==37)THEN
2734 CALL SIGEPS37_SINGLE_CELL (
2735 1 ELBUF_TAB, IXS, BUFMAT, IPARG, IPM,
2736 2 IDLOC , NG , brickID, VOL
2737 . )
2738 ENDIF
2739
2740 !
2741 !VOLcell = sum (VOLslacve) : each volsecnd is removed from
2742 !its former main cell and cumulative value is set in new
2743 !main cell to demerge.
2744 !-----------------------------------------------------!
2745 ! ADVECTION STEP UPDATE ORIGIN CELL (PREVIOUS MAIN) !
2746 !-----------------------------------------------------!
2747 IF(IADV==1)THEN
2748 NG = BRICK_LIST(NIN,IBm)%NG
2749 IDLOC = BRICK_LIST(NIN,IBm)%IDLOC
2750 GBUF => ELBUF_TAB(NG)%GBUF
2751 LBUF => ELBUF_TAB(NG)%BUFLY(1)%LBUF(1,1,1)
2752 GBUF%RHO(IDLOC) = GBUF%RHO(IDLOC) - RHO_ADV/VolCEll
2753 GBUF%EINT(IDLOC) = GBUF%EINT(IDLOC) - EINT_ADV/VolCEll
2754 GBUF%MOM(3*(IDLOC-1) +1) = GBUF%MOM(3*(IDLOC-1) +1) - MOM_ADV(1)/VolCELL
2755 GBUF%MOM(3*(IDLOC-1) +2) = GBUF%MOM(3*(IDLOC-1) +2) - MOM_ADV(2)/VolCELL
2756 GBUF%MOM(3*(IDLOC-1) +3) = GBUF%MOM(3*(IDLOC-1) +3) - MOM_ADV(3)/VolCELL
2757 GBUF%SIG(LLT_*(1-1) +IDLOC) = GBUF%SIG(LLT_*(1-1) +IDLOC) - SIG_ADV(1)/VolCELL
2758 GBUF%SIG(LLT_*(2-1) +IDLOC) = GBUF%SIG(LLT_*(2-1) +IDLOC) - SIG_ADV(2)/VolCELL
2759 GBUF%SIG(LLT_*(3-1) +IDLOC) = GBUF%SIG(LLT_*(3-1) +IDLOC) - SIG_ADV(3)/VolCELL
2760 GBUF%SIG(LLT_*(4-1) +IDLOC) = GBUF%SIG(LLT_*(4-1) +IDLOC) - SIG_ADV(4)/VolCELL
2761 GBUF%SIG(LLT_*(5-1) +IDLOC) = GBUF%SIG(LLT_*(5-1) +IDLOC) - SIG_ADV(5)/VolCELL
2762 GBUF%SIG(LLT_*(6-1) +IDLOC) = GBUF%SIG(LLT_*(6-1) +IDLOC) - SIG_ADV(6)/VolCELL
2763 MTN_ = IPARG(1,NG)
2764 IF(MTN_==37)THEN
2765 LLT_ = IPARG(2,NG)
2766 IDLOC = IDLOCB(IBm)
2767 MBUF => ELBUF_TAB(NG)%BUFLY(1)%MAT(1,1,1)
2768 brickID = BRICK_LIST(NIN,IBm)%ID
2769 MT = IXS(1,brickID)
2770 IADBUF = IPM(7,MT)
2771 RHO10 = BUFMAT(IADBUF-1+11)
2772 RHO20 = BUFMAT(IADBUF-1+12)
2773 IF(UVAR(3) == ZERO) UVAR(3) = RHO10
2774 IF(UVAR(2) == ZERO) UVAR(2) = RHO20
2775 !UVAR(I,1) : massic percentage of liquid * global density (rho1*V1/V : it needs to give liquid mass multiplying by element volume in aleconve.F)
2776 !UVAR(I,2) : density of gas
2777 !UVAR(I,3) : density of liquid
2778 !UVAR(I,4) : volumetric fraction of liquid
2779 !UVAR(I,5) : volumetric fraction of gas
2780 MBUF%VAR((1-1)*LLT_+IDLOC) = UVAR(1) - UVAR_ADV(1)/VolCELL
2781
2782 !APPEL SIGEPS37 POUR EQUILIBRE PRESSION
2783 brickID = BRICK_LIST(NIN,IBm)%ID
2784 IDLOC = IDLOCB(IBm)
2785 NG = NGB(IBm)
2786 VOL = BRICK_LIST(NIN,IBm)%Vold_Scell
2787 MTN_ = IPARG(1,NG)
2788 IF(MTN_==37)THEN
2789 CALL SIGEPS37_SINGLE_CELL (
2790 1 ELBUF_TAB, IXS, BUFMAT, IPARG, IPM,
2791 2 IDLOC , NG , brickID, VOL
2792 . )
2793 ENDIF
2794
2795 ELSEIF(MTN_==51)THEN
2796 ! print *, "todo" ; pause
2797 ENDIF
2798 ENDIF
2799
2800
2801 !-----------------------------------------------!
2802 ENDDO!next IB
2803
2804 if (IBUG22_merge==-1)then
2805 call my_barrier
2806.and. if(itask==0 DT1>ZERO)then
2807 allocate(order(IN22),value(IN22))
2808 DO Ib=1,in22
2809 VALUE(ib)=tmp22array(1,ib)
2810 ENDDO
2811 order=0
2812 !CALL QUICKSORT_I2 !(ORDER,IN22,value)
2813 do ib2=1,IN22
2814 ib=order(ib2)
2815 cod1=tmp22array(1,ib)
2816 cod2=tmp22array(2,ib)
2817 cod3=tmp22array(3,ib)
2818 !if(cod1==0)cycle
2819 if(cod3==0)then
2820 write(*,FMT='(a,i8,a,i8)') "DEMERGING : id=", Cod1, " from idv:", cod2
2821 else
2822 write(*,FMT='(a,i8,a,i8,a,i8)') "DEMERGING : id=", Cod1, " from idv:", cod2," moved to", cod3
2823 endif
2824 write (*,FMT='(i10,a,f30.16,a,f30.16)') Cod1," Vold=", tmp22array(4,ib), " Vnew=", tmp22array(5,ib)
2825 write (*,FMT='(i10,a,f30.16,a,f30.16)') Cod2," Vold=", tmp22array(6,ib), " Vnew=", tmp22array(7,ib)
2826 write (*,FMT='(a,f30.16)') " Number of access =", tmp22array(8,ib)
2827 enddo
2828 endif
2829 call my_barrier
2830 if(itask==0)then
2831 if(allocated(tmp22array))deallocate(tmp22array)
2832.and. if(DT1>ZERO allocated(order) )deallocate(order)
2833 endif
2834 endif
2835
2836
2837 !------------------------------------------------------------!
2838 ! @28. COMPUTE VOLD_Vnew_Cell for newly cell in buffer !
2839 !------------------------------------------------------------!
2840 !OLD_Vnew_Cell(1:9) are old cell volumes for a cut cell.
2841 !It is taken from previous cycle.
2842 !If brick was not cut on previous cycle then it must be initialized here
2843 IF(DT1==ZERO)NBL1 = 0
2844 DO IB=NBF,NBL1
2845 NBCUT = BRICK_LIST(NIN,IB)%NBCUT
2846 WASCUT = BRICK_LIST(NIN,IB)%WasCut
2847 NewInBuffer = BRICK_LIST(NIN,IB)%NewInBuffer
2848CCC IF(WASCUT==0.and. NBCUT==0)THEN
2849 IF( NewInBuffer == 1 )THEN
2850
2851 MCELL = 1
2852 BRICK_LIST(NIN,IB)%POLY(MCELL)%OLD_Vnew = BRICK_LIST(NIN,IB)%POLY(MCELL)%Vnew
2853 BRICK_LIST(NIN,IB)%POLY(2:9)%OLD_Vnew = ZERO
2854! print *, "new brick in cut cell buffer:", IXS(11,brick_list(nin,ib)%id)
2855! print *, "new brick in cut cell buffer: BRICK_LIST(NIN,IB)%POLY(MCELL)%OLD_Vnew=",BRICK_LIST(NIN,IB)%POLY(MCELL)%Vnew
2856
2857 NG = BRICK_LIST(NIN,IB)%NG
2858 IDLOC = BRICK_LIST(NIN,IB)%IDLOC
2859 ID = BRICK_LIST(NIN,IB)%ID
2860 GBUF => ELBUF_TAB(NG)%GBUF
2861 BRICK_LIST(NIN,IB)%Vold_Scell = GBUF%VOL(IDLOC)
2862
2863c print *, "NEW IN BUFFER : ", ixs(11,id)
2864c print *, " vold_scell : ", brick_list(nin,ib)%vold_scell
2865c pause
2866
2867 ENDIF
2868 ENDDO!next IB
2869
2870 !------------------------------------------------------------!
2871 ! @29. COMPUTE VOLUME CHANGE FOR EACH CELL !
2872 !------------------------------------------------------------!
2873.OR..AND. IF(IBUG22_TRACKING>0 IBUG22_TRACKING==-1ITASK==0)print *, "=== TOPOLOGICAL TRACKING ==="
2874 !topology is changing so this is not trivial and must be taken into account
2875 !example : a single secnd hexa can be split into {1 penta + 2 tetra} so
2876 !it is not possible to just make difference between two polyhedra volumes.
2877 NBL1 = NBL
2878 IF(DT1==ZERO)THEN
2879 NBL1 = 0
2880 DO IB=NBF,NBL
2881 BRICK_LIST(NIN,IB)%POLY(1:9)%DVOL(1) = ZERO
2882 BRICK_LIST(NIN,IB)%DVOL = ZERO
2883 ENDDO
2884 ENDIF
2885 DO IB=NBF,NBL1
2886 ID = BRICK_LIST(NIN,IB)%ID
2887 MCELL = BRICK_LIST(NIN,IB)%mainID
2888 NCELL = BRICK_LIST(NIN,IB)%NBCUT
2889 ICELL = 0
2890 M(:,:) = ZERO
2891 WasCut = BRICK_LIST(NIN,IB)%WasCut
2892 MTN_ = IPARG(1,NGB(IB))
2893
2894 IF(WasCut>0)THEN !forcemment M(9,JJ)>ZERO, sorti du test conditionel
2895 !filling relation matrix
2896 DO INOD = 1,8
2897 JJ = BRICK_LIST(NIN,IB)%NODE(INOD)%WhichCell
2898 II = BRICK_LIST(NIN,IB)%NODE(INOD)%OLD_WhichCell
2899 IF(II <0) II = 1 !just being cut on this cycle so whichcell node was not set on previous cycle
2900 M(II,JJ) = BRICK_LIST(NIN,IB)%POLY(JJ)%Vnew !instead of setting kronecker symbol, volume is directly stored. V=0 <=> delta_ij=0
2901 ENDDO!next INOD
2902.AND. IF(NBCUT>2 M(9,9)==ZERO)M(9,9) = BRICK_LIST(NIN,IB)%POLY(JJ)%Vnew
2903 !cas particulier cellule 9 et cellule i (sigma_col>1)
2904 !les 2 ne peuvent etre liees car pas le meme cote de la surface.
2905 ICELL = 9
2906 JJ = ICELL
2907 VAR = ZERO
2908 DO II=1,9
2909 IF(M(II,JJ)>ZERO)VAR=VAR+ONE
2910 ENDDO!next II
2911 IF(NCELL==0)THEN
2912 II=MAXLOC(BRICK_LIST(NIN,IB)%POLY(1:9)%OLD_Vnew,1)
2913 DO I=1,9
2914 IF (I==II) CYCLE
2915 M(I,JJ) = ZERO
2916 ENDDO
2917 ELSE
2918 IF(VAR>TWO )THEN
2919 M(1:8,JJ) = ZERO !plusieurs polyedres disparraissent, on ne les compte plus.
2920 ELSEIF(VAR==TWO )THEN
2921 !garder le volume le plus gros
2922cc DO II=1,8 ! 1st value is between 1,8 ; last one is at index 9
2923cc VAR = BRICK_LIST(NIN,IB)%POLY(II)%OLD_Vnew
2924cc IF(VAR>ZERO)EXIT
2925cc ENDDO
2926cc IF(VAR>BRICK_LIST(NIN,IB)%POLY(9)%OLD_Vnew)THEN
2927cc M(9,JJ) = ZERO
2928cc ELSE
2929cc M(II,JJ) = ZERO
2930cc ENDIF
2931
2932 !essai correction 20 avril 2010
2933 VAR = BRICK_LIST(NIN,IB)%POLY(1)%OLD_Vnew
2934 VAR_ = BRICK_LIST(NIN,IB)%POLY(9)%OLD_Vnew
2935 lCOND1 = VAR<VAR_
2936 VAR = BRICK_LIST(NIN,IB)%POLY(1)%Vnew
2937 VAR_ = BRICK_LIST(NIN,IB)%POLY(9)%Vnew
2938 lCOND2 = VAR<VAR_
2939
2940 M(9,9) = ZERO
2941
2942! IF(lCOND1/=lCOND2)THEN
2943! VAR=M(9,1)
2944! M(9,1) = M(1,9)
2945! M(1,9) = VAR
2946! ENDIF
2947
2948
2949
2950 ENDIF
2951 ENDIF!NCELL==0
2952
2953 ELSEIF(WASCUT == 0)THEN
2954 MCELL = BRICK_LIST(NIN,IB)%MainID
2955 M(1:9,1:9) = ZERO
2956 M(1,MCELL) = BRICK_LIST(NIN,IB)%POLY(MCELL)%Vnew
2957
2958 ENDIF
2959
2960c IF(IBUG22_TRACKING>0 .OR. IBUG22_TRACKING==-1)then
2961c !display matrix
2962c if (ncell > 0)then
2963c write(*,FMT=('(A,9I15 )')) " ", 1,2,3,4,5,6,7,8,9
2964c do II=1,9
2965c write(*,FMT=('(A,9F15.5 )')) " |", M(II,1:9 )
2966c enddo
2967c endif
2968c ENDIF
2969
2970
2971
2972 !calculating old_volume for each cell
2973 ICELL = 0
2974 DO WHILE (ICELL<=NCELL) ! loop on polyhedron {1:NCELL} U {9}
2975 ICELL = ICELL + 1
2976.AND. IF (ICELL>NCELL NCELL/=0)ICELL=9
2977 VAR = ZERO
2978 VAR_ = ZERO
2979 JJ = ICELL
2980 VFRAC(1:4) = ZERO
2981 VAR_VF(1:4) = ZERO
2982 VOLD_phase(1:4) = ZERO
2983 Vj = BRICK_LIST(NIN,IB)%POLY(JJ)%Vnew
2984 DO II=1,9
2985 IF(M(II,JJ)==ZERO)CYCLE
2986 Vi = BRICK_LIST(NIN,IB)%POLY(II)%OLD_Vnew
2987 !Vi : Vprev, Vj: Vnew for ratio
2988 VAR = VAR + Vi *Vj / SUM(M(II,:)) !this is not a sum of kronecker symbol since volume were stored instead
2989 ENDDO
2990 BRICK_LIST(NIN,IB)%POLY(ICELL)%Vold = VAR
2991
2992c IF(IBUG22_TRACKING==-1 .OR. IBUG22_TRACKING==IXS(11,brick_list(nin,ib)%id))then
2993c write(*,FMT=('(A,I8,A1,I1 )')) " brickID =", ixs(11,id),".",ICELL
2994c write(*,FMT=('(A,I8 )')) " NBCUT =", brick_list(nin,ib)%nbcut
2995c write(*,FMT=('(A,I8 )')) " WASCUT =", brick_list(nin,ib)%wascut
2996c write(*,FMT=('(A,F30.16 )')) " VOLD =", BRICK_LIST(NIN,IB)%POLY(ICELL)%Vold
2997c write(*,FMT=('(A,F30.16 )')) " VNEW =", BRICK_LIST(NIN,IB)%POLY(ICELL)%Vnew
2998c var_ = zero
2999c var__ = zero
3000c icellTag(1:9) = 0
3001c icellTag_old(1:9) = 0
3002c do ii=1,9
3003c if(M(ii,ICELL)>ZERO)then
3004c var_=var_+ONE
3005c icellTag(ii) = 1
3006c endif
3007c enddo
3008c if(var_==ONE)then
3009c var__=zero
3010c do ii=1,9
3011c K = MAXLOC(icellTag,1)
3012c if( M( K ,ii) > ZERO )then
3013c var__=var__+ONE
3014c icelltag_old(ii) = 1
3015c endif
3016c enddo
3017c endif
3018c if ( var_ >= TWO )then
3019c SOM = ZERO
3020c VFRAC(1:4) = ZERO
3021c DO k=1,9
3022c if(icelltag(k)==0)CYCLE
3023c Vi = BRICK_LIST(NIN,IB)%POLY(K)%OLD_Vnew
3024c SOM = SOM + Vi
3025c enddo
3026c write(*,FMT=('(A,I8)')) " merge from poly:", IcellTag(1:9)
3027c elseif (var_ == ONE .AND. var__==ONE)then
3028c write(*,FMT=('(A,I8)')) " cell is stable (same number as previous cycle)", Icelltag(1)
3029c elseif (var__ >= TWO)then
3030c K = MAXLOC(icellTag,1)
3031c write(*,FMT=('(A,I8)')) " cell is a fragment of a previous polyhedron ",K !1seul dans la colonne (var_) et sur la ligne (var__)
3032c endif
3033c ENDIF
3034
3035 ENDDO!next ICELL
3036 ENDDO !next IB
3037
3038
3039
3040 !------------------------------------------------------------!
3041 ! @30. LINK SWITCHING !
3042 ! link has change then transportation must occur !
3043 !------------------------------------------------------------!
3044 NBL1 = NBL
3045 IN = 0
3046 debug_outp2 = .false.
3047 if(IBUG22_LINK_SWITCH/=0)then
3048 if(itask==0)ALLOCATE (TMP22ARRAY(6,NB))
3049 call my_barrier
3050 TMP22ARRAY(1:6,NBF:NBL)=ZERO
3051 do ib=nbf,nbl
3052 ie = brick_list(nin,ib)%id
3053.OR. if(ixs(11,ie)==IBUG22_LINK_SWITCH IBUG22_LINK_SWITCH==-1)then
3054 debug_outp2 = .true.
3055 exit
3056 endif
3057 enddo
3058 endif!(IBUG22_LINK_SWITCH/=0)
3059 IF(DT1==ZERO)NBL1=0
3060 DO IB=NBF,NBL1
3061 NCELL = BRICK_LIST(NIN,IB)%NBCUT
3062 !pAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%Adjacent_Cell(1:5)
3063 !pNAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%NAdjCell
3064 pAdjBRICK => BRICK_LIST(NIN,IB)%Adjacent_Brick(1:6,1:5)
3065 pSUBVOL(1:9) => BRICK_LIST(NIN,IB)%POLY(1:9)
3066 MBUF => ELBUF_TAB(NGB(IB))%BUFLY(1)%MAT(1,1,1)
3067 MCELL = BRICK_LIST(NIN,IB)%MainID
3068 ID = BRICK_LIST(NIN,IB)%ID
3069 NG = BRICK_LIST(NIN,IB)%NG
3070 IDLOC = BRICK_LIST(NIN,IB)%IDLOC
3071 GBUF => ELBUF_TAB(NGB(IB))%GBUF
3072 LBUF => ELBUF_TAB(NGB(IB))%BUFLY(1)%LBUF(1,1,1)
3073 ICELL = 0
3074c print *, "check switch in brickID=", ixs(11,brick_list(nin,ib)%id)
3075 !------LOOP--------!
3076 NumSECND = BRICK_LIST(NIN,IB)%SecndList%Num
3077 !BRICK_LIST(NIN,IB)%SecndList%ISWITCH(1:NumSECND) = 0
3078 DO IC=1,NumSECND
3079 ICELLv = BRICK_LIST(NIN,IB)%SecndList%ICELLv(IC)
3080 IBv = BRICK_LIST(NIN,IB)%SecndList%IBv(IC)
3081 FM = BRICK_LIST(NIN,IB)%SecndList%FM(IC)
3082 FV = BRICK_LIST(NIN,IB)%SecndList%FV(IC)
3083 brickID = BRICK_LIST(NIN,IB)%ID
3084 IBv = IBv
3085 NGv = BRICK_LIST(NIN,IBv)%NG
3086 IDLOCv = BRICK_LIST(NIN,IBv)%IDLOC
3087 IV = BRICK_LIST(NIN,IBv)%ID
3088 ITASKv = BRICK_LIST(NIN,IBv)%ITASK
3089 ! loop over the nodes of the adjacent secondary cell
3090 NNODES = BRICK_LIST(NIN,IBv)%POLY(ICELLV)%NumNOD
3091 J1 = 0
3092 J2 = 0
3093 DO K1=1,NNODES !Nodes attached to IB
3094 INOD = BRICK_LIST(NIN,IBv)%POLY(ICELLv)%ListNodID(K1)
3095 FV_old = BRICK_LIST(NIN,IBv)%NODE(INOD)%WhereWasMain !FACE
3096 IF (FV_old == 0)CYCLE !the node was wet by main cell : MERGE PROCESS APPLIES (see related chapter)
3097 !OLD_M
3098 IF (FV_old == FV)CYCLE ! NO LINK CHANGE : this is still linked to the same direction
3099 !else FV_old /= FV : THERE WAS A SWITCH
3100 ! GET PREVIOUS FACE TO GET PREVIOUS main CELL
3101 IF(FV_old<=NV46)THEN
3102 IBmo = BRICK_LIST(NIN,IBv)%Adjacent_Brick(FV_old,4)
3103 ELSE
3104 J = FV_old
3105 J1 = J/10
3106 J2 = MOD(J,10)
3107 IBv = BRICK_LIST(NIN,IBv)%Adjacent_Brick(J1,4)
3108 IBmo = BRICK_LIST(NIN,IBv)%Adjacent_Brick(J2,4)
3109 ENDIF
3110
3111 IBm = IB
3112 LLTm = IPARG(2,NGB(IB))
3113 ! IF(IBm==IBmo)CYCLE : impossible it meansFV_old=FV
3114 lCYCLE = .FALSE.
3115 NumTarget = BRICK_LISt(NIN,IBmo)%Ntarget
3116 IF(NumTarget>=2)print *, "**inter22 - Warning Multiple targets",ixs(11,BRICK_LIST(NIN,IBv)%ID), ICELLv
3117 DO ITARGET=1,NumTarget
3118 IF (BRICK_LISt(NIN,IBmo)%MergeTarget(3,ITarget)==IBm)THEN
3119 lCYCLE = .TRUE.
3120 EXIT
3121 ENDIF
3122 IBmo = BRICK_LIST(NIN,IBmo)%MergeTarget(3,ITarget) !en tout rigueur il faudrait prelever dans toutes les target, ici on prend le premier
3123 print *, "**inter22 : multiple targets", ixs(11,BRICK_LIST(NIN,IBv)%ID), ICELLv
3124 ENDDO
3125 if(lCYCLE)CYCLE
3126 lCYCLE = .FALSE.
3127 !remove salve cell volume from IBMo(Vold) and add it in IBo(Vold)
3128 !addition is treated now since only current thread is handling this cell IBm \in [NBF,NBL]
3129 !search corresponding secnd cell in IBmo%old_secnds_lit
3130 lFOUND = .FALSE. !if found old salve cell attached to INOD
3131 NumSecnd2 = OLD_SecndList(NIN,IBmo)%Num
3132 NGo = BRICK_LIST(NIN,IBmo)%NG
3133 LLTo = IPARG(2,NGo)
3134 DO IC2 = 1, NumSecnd2
3135 !get previous secnd cell
3136 FV2 = OLD_SecndList(NIN,IBmo)%FV(IC2)
3137 ICELLv2 = OLD_SecndList(NIN,IBmo)%ICELLv(IC2)
3138 IBv2 = OLD_SecndList(NIN,IBmo)%IBv(IC2)
3139 !Si la cellule escalve etait liee par une autre face que FV_old alors ce n'etait pas la bonne cellule main et on passe outres
3140 IF(fv2 /= fv_old )cycle
3141 !la cellule second. est attachee a la cellule main qui nous interesse IBmo
3142 !si on arrive ici c'est que la cellule second. (0, en partie, ou en totalite) change de liaison
3143 !on verifie si un de ses noeuds correspond au noeud INOD de la nouvelle cellule escalve dans IBv
3144 !on boucle sur les noeuds sommets de l'ancienne cellule pour voir quel noeud sommet a commute
3145 volcell = zero
3146 volsecnd = zero
3147 nnodes2 = old_secndlist(nin,ibmo)%NumNOD_Cell(ic2)
3148 DO k2=1,nnodes2
3149 inod2 = old_secndlist(nin,ibmo)%ListNodID(ic2,k2)
3150 icell2 = brick_list(nin,ibv2)%NODE(inod2)%WhichCell
3151 ibmcur = brick_list(nin,ibv2)%POLY(icell2)%WhereIsMain(4)
3152 IF(ibmcur == ib) volsecnd = volsecnd + one/nnodes*brick_list(nin,ibv2)%POLY(icell2)%Vnew !example brick 20 devient main dans cas 15.14 : le penta dans 41 a deux noeuds, si on boucle 2 fois on prend 2 fois son volume pour l'injecter dans bric 30.
3153 volcell = volcell + one/nnodes*brick_list(nin,ibv2)%POLY(icell2)%Vnew
3154 ENDDO
3155 if (volsecnd == zero)cycle
3156 !BRICK_LIST(NIN,IB)%SecndList%ISWITCH(IC) = 1
3157 !print *, ""
3158 !print *, "VolCell =", VolCell
3159 !print *, "VolSECND =", VolSECND
3160 ratio = volsecnd/volcell !est le ratio de prelevement dans le volume de l'ancienne cellule second.
3161 !here IC2(cycle-1) corresponds to IC(cycle+0)
3162 volcell = ratio * old_secndlist(nin,ibmo)%VOL(ic2)
3163 vold = brick_list(nin,ibm)%Vold_SCell
3164 gbufo => elbuf_tab(ngb(ibmo))%GBUF
3165 !---A.MATERIAL MERGING---!
3166 !Ratio is <1 in case 1 secnd split into several secnd cells which different main links, then collect only relevant quantity not all from old secnd cell.
3167 eint = volcell * gbufo%EINT(idlocb(ibmo))
3168 rho = volcell * gbufo%RHO(idlocb(ibmo))
3169 mom(1) = volcell * gbufo%MOM(llto*(1-1) + idlocb(ibmo))
3170 mom(2) = volcell * gbufo%MOM(llto*(2-1) + idlocb(ibmo))
3171 mom(3) = volcell * gbufo%MOM(llto*(3-1) + idlocb(ibmo))
3172 DO j=1,nv46
3173 sig(j) = volcell * gbufo%SIG(llto*(j-1)+idlocb(ibmo))
3174 ENDDO
3175 gbuf%EINT(idlocb(ibm)) = (gbuf%EINT(idlocb(ibm)) * vold + eint) / (vold+volcell)
3176 gbuf%RHO(idlocb(ibm)) = (gbuf%RHO(idlocb(ibm)) * vold + rho) / (vold+volcell)
3177 gbuf%MOM(lltm*(1-1) + idlocb(ibm)) = (gbuf%MOM(lltm*(1-1) + idlocb(ibm)) * vold + mom(1)) / (vold+volcell)
3178 gbuf%MOM(lltm*(2-1) + idlocb(ibm)) = (gbuf%MOM(lltm*(2-1) + idlocb(ibm)) * vold + mom(2)) / (vold+volcell)
3179 gbuf%MOM(lltm*(3-1) + idlocb(ibm)) = (gbuf%MOM(lltm*(3-1) + idlocb(ibm)) * vold + mom(3)) / (vold+volcell)
3180 alefvm_buffer%FCELL(4, brick_list(nin,ibm)%id )= gbuf%RHO(idlocb(ibm))
3181 alefvm_buffer%FCELL(1, brick_list(nin,ibm)%id )= gbuf%MOM(lltm*(1-1) + idlocb(ibm))
3182 alefvm_buffer%FCELL(2, brick_list(nin,ibm)%id )= gbuf%MOM(lltm*(2-1) + idlocb(ibm))
3183 alefvm_buffer%FCELL(3, brick_list(nin,ibm)%id )= gbuf%MOM(lltm*(3-1) + idlocb(ibm))
3184
3185 !FCELL(5, brick_list(nin,ibm)%id) = ! see mtn 37 conditional block below
3186 !FCELL(6, brick_list(nin,ibm)%id) = ! idem
3187
3188 DO j=1,nv46
3189 gbuf%SIG(lltm*(j-1)+idlocb(ibm)) = (gbuf%SIG(lltm*(j-1)+idlocb(ibm)) * vold + sig(j)) / (vold+volcell)
3190 ENDDO
3191
3192 if(debug_outp2)then
3193 !!!write (*,FMT='(A)') " === LINK SWITCH ==="
3194 in = in + 1
3195 endif
3196
3197 if(ibug22_link_switch==ixs(11,brick_list(nin,ib)%id) .or. ibug22_link_switch==-1 )then
3198 ! print *, "brick target =", ixs(11,brick_list(nin,ib)%id)
3199 ! print *, "brick origin =", ixs(11,brick_list(nin,ibmo)%id)
3200 ! print *, "adding",VOLcell ,'to',BRICK_LIST(NIN,IBm)%Vold_SCell
3201 ! print *, "updated target -old volume- =",BRICK_LIST(NIN,IBm)%Vold_SCell +VOLcell
3202 tmp22array(1,ibm)= ixs(11,brick_list(nin,ib)%id)
3203 tmp22array(2,ibm)= ixs(11,brick_list(nin,ibmo)%id)
3204 tmp22array(3,ibm)= ixs(11,brick_list(nin,ibm)%id)
3205 tmp22array(4,ibm)= brick_list(nin,ibm)%Vold_SCell
3206 tmp22array(5,ibm)= brick_list(nin,ibm)%Vold_SCell +volcell
3207 tmp22array(6,ibm)= volcell
3208 endif
3209
3210 !GBUF%VOL(IDLOCB(IBm)) = GBUF%VOL(IDLOCB(IBm)) + VOLcell
3211 brick_list(nin,ibm)%Vold_SCell = brick_list(nin,ibm)%Vold_SCell + volcell
3212
3213 mtn_ = iparg(1,ngb(ibmo))
3214 volcell51(1:trimat) = zero
3215 IF(mtn_==37)THEN
3216 mbuf => elbuf_tab(ngb(ib))%BUFLY(1)%MAT(1,1,1)
3217 mt = ixs(1,brickid)
3218 iadbuf = ipm(7,mt)
3219 rho10 = bufmat(iadbuf-1+11)
3220 rho20 = bufmat(iadbuf-1+12)
3221 !UVAR(I,1) : massic percentage of liquid * global density (rho1*V1/V : it needs to give liquid mass multiplying by element volume in aleconve.F)
3222 !UVAR(I,2) : density of gas
3223 !UVAR(I,3) : density of liquid
3224 !UVAR(I,4) : volumetric fraction of liquid
3225 !UVAR(I,5) : volumetric fraction of gas
3226 mbuf => elbuf_tab(ngb(ibm))%BUFLY(1)%MAT(1,1,1)
3227 mbufo => elbuf_tab(ngb(ibmo))%BUFLY(1)%MAT(1,1,1)
3228 llt_ = iparg(2,ngb(ibm))
3229 llt_o = iparg(2,ngb(ibmo))
3230
3231 ! UVAR(5) - alpha_gas
3232 pvar => mbuf%VAR((5-1)*llt_ +idlocb(ibm))
3233 pvaro => mbufo%VAR((5-1)*llt_o+idlocb(ibmo))
3234C print *, "link switch 5", pVAR, (pVAR * VOLD + VOLcell * pVARo)/ (VOLD+VOLcell)
3235 pvar = (pvar * vold + volcell * pvaro)/ (vold+volcell)
3236 pvar = max(pvar,zero)
3237
3238 ! UVAR(4) - alpha_liquid
3239 pvar => mbuf%VAR((4-1)*llt_ +idlocb(ibm))
3240 pvaro => mbufo%VAR((4-1)*llt_o+idlocb(ibmo))
3241C print *, "link switch 4", pVAR, (pVAR * VOLD + VOLcell * pVARo)/ (VOLD+VOLcell)
3242 pvar = (pvar * vold + volcell * pvaro)/ (vold+volcell)
3243 pvar = max(pvar,zero)
3244
3245 ! UVAR(3) - rho_liquid
3246 pvar => mbuf%VAR((3-1)*llt_ +idlocb(ibm))
3247 pvaro => mbufo%VAR((3-1)*llt_o+idlocb(ibmo))
3248 alp = mbuf%VAR((4-1)*llt_ +idlocb(ibm))
3249 alpo = mbufo%VAR((4-1)*llt_o+idlocb(ibmo))
3250C print *, "link switch 3", pVAR, (pVAR*ALP*VOLD+VOLcell*pVARo*ALPo /(VOLD+VOLcell)
3251 pvar = (pvar*alp*vold+volcell*pvaro*alpo)
3252 IF(pvar>zero)pvar=pvar/(alp*vold+alpo*volcell)
3253 pvar = max(pvar,zero)
3254 !cts14june2017!
3255 IF( pvar == zero) pvar = rho10
3256
3257 ! UVAR(2) - rho_gas
3258 pvar => mbuf%VAR((2-1)*llt_ +idlocb(ibm))
3259 pvaro => mbufo%VAR((2-1)*llt_o+idlocb(ibmo))
3260 alp = mbuf%VAR((5-1)*llt_ +idlocb(ibm))
3261 alpo = mbufo%VAR((5-1)*llt_o+idlocb(ibmo))
3262C print *, "link switch 2", pVAR, (pVAR*ALP*VOLD+VOLcell*pVARo*ALPo /(VOLD+VOLcell)
3263 pvar = (pvar*alp*vold+volcell*pvaro*alpo)
3264 IF(pvar>zero)pvar=pvar/(alp*vold+alpo*volcell)
3265 pvar = max(pvar,zero)
3266 !cts14june2017!
3267 IF( pvar == zero) pvar = rho20
3268
3269 ! UVAR(1) - rho1.V1/V = rho1.alpha_liquid
3270 pvar => mbuf%VAR((1-1)*llt_ +idlocb(ibm))
3271 pvaro => mbufo%VAR((1-1)*llt_o+idlocb(ibmo))
3272C print *, "link switch 1", pVAR, (pVAR * VOLD + VOLcell * pVARo)/ (VOLD+VOLcell)
3273 pvar = (pvar * vold + volcell * pvaro)/ (vold+volcell)
3274 pvar = max(pvar,zero)
3275
3276 !GBUF%RHO(IDLOCB(IBm)) = (GBUF%RHO(IDLOCB(IBm)) * VOLD + RHO) / (VOLD+VOLcell)
3277
3278
3279
3280 !APPEL SIGEPS37 POUR EQUILIBRE PRESSION
3281 brickid = brick_list(nin,ibm)%ID
3282 idloc = idlocb(ibm)
3283 ng = ngb(ibm)
3284 vol = vold+volcell
3285 CALL sigeps37_single_cell (
3286 1 elbuf_tab, ixs, bufmat, iparg, ipm,
3287 2 idloc , ng , brickid, vol
3288 . )
3289
3290
3291 ELSEIF(mtn_==51)THEN
3292 mbuf => elbuf_tab(ngb(ibm))%BUFLY(1)%MAT(1,1,1)
3293 llt_ = iparg(2,ngb(ibm))
3294 mbufo => elbuf_tab(ngb(ibmo))%BUFLY(1)%MAT(1,1,1)
3295 llt_o = iparg(2,ngb(ibmo))
3296
3297 DO itrimat=1,trimat
3298 !IBm_old - volume fraction!
3299 ipos = 1
3300 k0 = m51_n0phas + (itrimat-1)*m51_nvphas +ipos-1 ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
3301 k = k0 * llt_o
3302 vfrac(itrimat) = mbufo%VAR(k+idlocb(ibmo))
3303 enddo!next ITRIMAT
3304
3305 DO itrimat=1,trimat
3306 ipos = 11
3307 k0 = m51_n0phas + (itrimat-1)*m51_nvphas +ipos-1 ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
3308 k = k0 * llt_
3309 volcell51(itrimat) = volcell * vfrac(itrimat)
3310 mbuf%VAR(k+idlocb(ibm)) = mbuf%VAR(k+idlocb(ibm)) + volcell51(itrimat)
3311 enddo!next ITRIMAT
3312
3313 DO itrimat=1,trimat
3314 !IBm - volume fraction!
3315 ipos = 1
3316 k0 = m51_n0phas + (itrimat-1)*m51_nvphas +ipos-1 ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
3317 k = k0 * llt_
3318 vfrac(itrimat) = mbuf%VAR(k+idlocb(ibm))
3319 IF(volcell51(itrimat)<=zero)cycle !unplug if nothiong to add (otherwise div by zero)
3320 DO ipos = 1,m51_nvphas
3321 IF(ipos==11)cycle !volume already treated
3322 k2 = ((m51_n0phas + (itrimat-1)*m51_nvphas )+ipos-1) ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
3323 k = llt_ * k2
3324 ko = llt_o * k2
3325 pvar => mbuf%VAR(k+idlocb(ibm))
3326 pvar = ( pvar * vfrac(itrimat)*vold + volcell51(itrimat)*mbufo%VAR(ko+idlocb(ibmo)) )
3327 pvar = pvar / (vfrac(itrimat)*vold + volcell51(itrimat))
3328 enddo!next IPOS
3329 enddo!next ITRIMAT
3330
3331 endif! law37 | law51
3332 !---THREAD TAG---!
3333C print *, " --- set VOLD_L += VOLcell, IB, ID ", VOLcell , IBMo, ixs(11,brick_list(nin,ibmo)%ID)
3334 vold_l(itask,0,ibmo) = vold_l(itask,0,ibmo) + volcell
3335 vold_l(itask,1:trimat,ibmo) = vold_l(itask,1:trimat,ibmo) + volcell51(1:trimat)
3336 enddo!next IC2
3337 EXIT !next IC (since this one is treated)
3338 enddo!next K1
3339 enddo!next IC
3340 enddo!next IB
3341
3342 CALL my_barrier
3343
3344 if(ibug22_link_switch/=0 .AND. itask==0 .and. dt1>zero)THEN
3345 write (*,fmt='(A)') " === LINK SWITCH ==="
3346 DO ib=1,nb
3347 IF(tmp22array(1,ib)==0)cycle
3348 print *, "brick target =", tmp22array(1,ib)
3349 print *, "brick origin =", tmp22array(2,ib)
3350 print *, "brick main =", tmp22array(3,ib)
3351 print *, "adding",tmp22array(6,ib) ,'to', tmp22array(4,ib)
3352 print *, "updated target -old volume- =", tmp22array(5,ib)
3353 ENDDO
3354 endif
3355 call my_barrier
3356
3357 !print
3358 if(itask==0)then
3359 if(ibug22_link_switch/=0)then
3360 DO ib=1,nb
3361 DO it = 0,nthread-1
3362 IF (vold_l(it,0,ib) == zero)cycle
3363 print *, ""
3364 print *, " brick ID =", ixs(11,brick_list(nin,ib)%id)
3365 print *, " removing ",vold_l(it,0,ib),'from',brick_list(nin,ib)%Vold_SCell
3366 print *, " new origin volume =", ixs(11,brick_list(nin,ib)%id)
3367 print *, " %vold, %vnew " , brick_list(nin,ib)%Vold_SCell- vold_l(it,0,ib), brick_list(nin,ib)%Vnew_SCell
3368 print *, ""
3369 ENDDO
3370 ENDDO
3371 endif
3372 endif
3373 call my_barrier
3374
3375 !---B.OLD VOLUME UPDATE---!
3376 !substract only on current thread
3377 DO ib=nbf,nbl1
3378 !GBUF => ELBUF_TAB(NGB(IB))%GBUF
3379 DO it = 0,nthread-1
3380 IF (vold_l(it,0,ib) == zero)cycle
3381 if(ibug22_link_switch==ixs(11,brick_list(nin,ib)%id) .or. ibug22_link_switch==-1 )then
3382 !print *, ""
3383 !print *, " brick ID =", ixs(11,brick_list(nin,ib)%id)
3384 !print *, " removing ",VOLD_L(IT,0,IB),'from',BRICK_LIST(NIN,IB)%Vold_SCell
3385 !print *, " new origin volume =", ixs(11,brick_list(nin,ib)%id)
3386 !print *, " %vold, %vnew " , BRICK_LIST(NIN,IB)%Vold_SCell- VOLD_L(IT,0,IB), BRICK_LIST(NIN,IB)%Vnew_SCell
3387 !print *, ""
3388 endif
3389
3390 brick_list(nin,ib)%Vold_SCell = brick_list(nin,ib)%Vold_SCell - vold_l(it,0,ib)
3391
3392 !!!!!!!!!!!!BRICK_LIST(NIN,IB)%Vold_SCell = BRICK_LIST(NIN,IB)%Vold_SCell - VOLD_L(IT,0,IB)
3393 mtn_ = iparg(1,ngb(ib))
3394 IF(mtn_==51)THEN
3395 mbuf => elbuf_tab(ngb(ib))%BUFLY(1)%MAT(1,1,1)
3396 llt_ = iparg(2,ngb(ib))
3397 DO itrimat=1,trimat
3398 ipos = 11
3399 k0 = m51_n0phas + (itrimat-1)*m51_nvphas +ipos-1 ! example : IPOS=1 => VFRAC {UVAR(I,ADD)=UVAR(K+I)}
3400 k = k0 * llt_
3401 mbuf%VAR(k+idlocb(ib)) = mbuf%VAR(k+idlocb(ib)) - vold_l(it,itrimat,ib)
3402 enddo!next ITRIMAT
3403 endif!(MTN_==51)
3404
3405
3406ccc if(IBUG22_LINK_SWITCH==ixs(11,brick_list(nin,ib)%id) .or. IBUG22_LINK_SWITCH==-1 )then
3407ccc print *, " --reset VOLD_L = ZERO, IB, ID ", ZERO , IB, ixs(11,brick_list(nin,ib)%ID)
3408ccc endif
3409 vold_l(it,0:max(0,trimat),ib) = zero !reset VOLD_L for next cycle
3410 enddo!next IT
3411 enddo!next IB
3412
3413
3414
3415
3416 if(ibug22_link_switch/=0)THEN
3417 call my_barrier
3418 if(itask==0)DEALLOCATE (tmp22array)
3419 endif
3420
3421 !------------------------------------------------------------!
3422 ! @32. LEVEL2 NEIGHBORHOOD (FOR FLUXES & CONVECTION) !
3423 !------------------------------------------------------------!
3424 DO ib=nbf,nbl
3425 nbcut = brick_list(nin,ib)%NBCUT
3426 pismain(1:9) => brick_list(nin,ib)%POLY(1:9)
3427 !pAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%Adjacent_Cell(1:5)
3428 padjbrick => brick_list(nin,ib)%Adjacent_Brick(1:6,1:5)
3429 pmainid => brick_list(nin,ib)%MainID
3430 ncell = nbcut
3431 icell = 0
3432 ie = brick_list(nin,ib)%ID
3433 !on traite les cellules mains pour eviter les lockon/lockoff (en evitant boucle sur les escalves)
3434 mcell = pmainid
3435 IF(mcell == 0)cycle ! next hexahedron cell
3436 brick_list(nin,ib)%ADJ_ELEMS%Num = 0
3437 DO j=1,6
3438 iv = padjbrick(j,1)
3439 ngv = padjbrick(j,2)
3440 idlocv = padjbrick(j,3)
3441 ibv = padjbrick(j,4)
3442 IF(ibv>0)THEN !if adjacent brick to face J is in cut cell buffer
3443 nbcutv = brick_list(nin,ibv)%NBCUT
3444 IF(nbcutv==0)cycle
3445 nadjcell = brick_list(nin,ib)%POLY(mcell)%FACE(j)%NAdjCell
3446 DO ic=1,nadjcell
3447 icellv = brick_list(nin,ib)%POLY(mcell)%FACE(j)%Adjacent_Cell(ic)
3448 mcellv = brick_list(nin,ibv)%MainID
3449 IF(icellv == mcellv)cycle
3450 !this adjacent cell is a secnd one, use it as bridge to a new adjacent main cell.
3451 !pWhereIsMainV => BRICK_LIST(NIN,IBv)%POLY(1:9)%WhereIsMain(1:4)
3452 jv = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(1) ! pWhereIsMainV(ICELLv,1)
3453 idm = brick_list(nin,ibv)%POLY(icellv)%WhereIsMain(3) ! pwhereismainv(icellv,3)
3454 IF(ie/=idm)THEN !secnd cell merged to the main IB
3455 !prendre tous les voisins main de la cellule salve sauf face JV
3456 !ajouter un seul voisin celui au quel la secnd cell est mergee
3457 brick_list(nin,ib)%ADJ_ELEMS%Num = brick_list(nin,ib)%ADJ_ELEMS%Num + 1
3458 k = brick_list(nin,ib)%ADJ_ELEMS%Num
3459 brick_list(nin,ib)%ADJ_ELEMS%IV(k) = idm
3460 brick_list(nin,ib)%ADJ_ELEMS%IB(k) = ibm
3461 brick_list(nin,ib)%ADJ_ELEMS%ICELL(k) = icellv
3462 brick_list(nin,ib)%ADJ_ELEMS%SecndFACE(k) = jv !to retrieve cut face area
3463 else! !secnd cell merged to main /= IB
3464 padjbrickv => brick_list(nin,ibv)%Adjacent_Brick(1:6,1:5)
3465 DO j2=1,6 !loop on IV.icellv (secnd) adjacent bricks
3466 IF (j2==jv)cycle !behind JV there is IB itself, so skip it
3467 ivv = padjbrick(j2,1)
3468 ngvv = padjbrick(j2,2)
3469 idlocvv = padjbrick(j2,3)
3470 ibvv = padjbrick(j2,4)
3471 ifvv = padjbrick(j2,5)
3472 IF(ivv == 0)cycle
3473
3474 IF(ibvv > 0)THEN
3475 mcellvv = brick_list(nin,ibvv)%MainID
3476 !if main IVV.MCELvv has connexion with IV.ICELLv then add it
3477 nadjcellv = brick_list(nin,ibvv)%POLY(mcellvv)%FACE(ifvv)%NAdjCell
3478 DO k=1,nadjcellv
3479 icv = brick_list(nin,ibvv)%POLY(mcellvv)%FACE(ifvv)%Adjacent_Cell(k)
3480 IF(icv/=icellv)cycle
3481 brick_list(nin,ib)%ADJ_ELEMS%Num = brick_list(nin,ib)%ADJ_ELEMS%Num + 1
3482 k2 = brick_list(nin,ib)%ADJ_ELEMS%Num
3483 brick_list(nin,ib)%ADJ_ELEMS%IV(k2) = ivv
3484 brick_list(nin,ib)%ADJ_ELEMS%IB(k2) = ibvv
3485 brick_list(nin,ib)%ADJ_ELEMS%ICELL(k2) = icellv !=ICv
3486 brick_list(nin,ib)%ADJ_ELEMS%SecndFACE(k2) = j2 !to retrieve cut face area
3487 EXIT
3488 enddo!next K
3489 ELSEIF(ibvv == 0)THEN
3490 brick_list(nin,ib)%ADJ_ELEMS%Num = brick_list(nin,ib)%ADJ_ELEMS%Num + 1
3491 k2 = brick_list(nin,ib)%ADJ_ELEMS%Num
3492 brick_list(nin,ib)%ADJ_ELEMS%IV(k2) = ivv
3493 brick_list(nin,ib)%ADJ_ELEMS%IB(k2) = ibvv
3494 brick_list(nin,ib)%ADJ_ELEMS%ICELL(k2) = 1
3495 brick_list(nin,ib)%ADJ_ELEMS%SecndFACE(k2) = j2 !to retrieve cut face area
3496 ENDIF
3497 enddo!next J2
3498 endif!(IE==IDM)
3499 enddo!next IC
3500 ENDIF
3501 enddo!next J
3502 ENDDO
3503
3504
3505
3506 !------------------------------------------------------------!
3507 ! @32. COUNT POINTS ON A GIVEN POLYHEDRON FACE !
3508 !------------------------------------------------------------!
3509 DO ib=nbf,nbl
3510 nbcut = brick_list(nin,ib)%NBCUT
3511 mcell = brick_list(nin,ib)%MainID
3512 ncell = brick_list(nin,ib)%NBCUT
3513 icell = 0
3514 nintp(1:6,1:9) = 0
3515 nnod(1:6,1:9) = 0
3516 IF(nbcut == 0) THEN
3517 brick_list(nin,ib)%POLY(1)%FACE(1:6)%NumPOINT=(/4,4,4,4,4,4/)
3518 cycle
3519 ENDIF
3520 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
3521 icell = icell + 1
3522 IF (icell>ncell .AND. ncell/=0)THEN
3523 icell=9
3524 cycle
3525 ENDIF
3526 secid = brick_list(nin,ib)%SECID_Cell(icell)
3527 brick_list(nin,ib)%POLY(icell)%FACE(1:6)%NumPOINT = 0
3528 DO j=1,6
3529 jj = gface(j,secid)
3530 IF(jj==0)EXIT
3531 np = gnpt(j,secid)
3532 nn = gnnod(j,secid)
3533 brick_list(nin,ib)%POLY(icell)%FACE(jj)%NumPOINT = np
3534 nintp(jj,icell) = np-nn !only intersection point
3535 nnod(jj,icell) = nn
3536 ENDDO
3537 enddo!next ICELL
3538 IF(icell==9)THEN
3539 DO j=1,6
3540 !sommets
3541 nn = sum( nnod(j,1:ncell) )
3542 nn = 4 - nn !complementary
3543 !points dintersection
3544 np = sum( nintp(j,1:ncell) )
3545 !nombre total de point = sommets + intersection pt
3546 brick_list(nin,ib)%POLY(9)%FACE(j)%NumPOINT = np + nn
3547 ENDDO
3548 ENDIF
3549 ENDDO
3550
3551 !------------------------------------------------------------!
3552 ! @33. LOCATE WHERE WAS main FOR NEXT CYCLE !
3553 ! & TAG ITS BRICK NODES (update in i22intersect.F) !
3554 !------------------------------------------------------------!
3555 DO ib=nbf,nbl
3556 pwherewasmain(1:8) => brick_list(nin,ib)%NODE(1:8)!%WhereWasMain
3557 !pWhereIsMain => BRICK_LIST(NIN,IB)%POLY(1:9)%WhereIsMain(1:4)
3558 plistnodid(1)%p(1:8) => brick_list(nin,ib)%POLY(1)%ListNodID(1:8)
3559 plistnodid(2)%p(1:8) => brick_list(nin,ib)%POLY(2)%ListNodID(1:8)
3560 plistnodid(3)%p(1:8) => brick_list(nin,ib)%POLY(3)%ListNodID(1:8)
3561 plistnodid(4)%p(1:8) => brick_list(nin,ib)%POLY(4)%ListNodID(1:8)
3562 plistnodid(5)%p(1:8) => brick_list(nin,ib)%POLY(5)%ListNodID(1:8)
3563 plistnodid(6)%p(1:8) => brick_list(nin,ib)%POLY(6)%ListNodID(1:8)
3564 plistnodid(7)%p(1:8) => brick_list(nin,ib)%POLY(7)%ListNodID(1:8)
3565 plistnodid(8)%p(1:8) => brick_list(nin,ib)%POLY(8)%ListNodID(1:8)
3566 plistnodid(9)%p(1:8) => brick_list(nin,ib)%POLY(9)%ListNodID(1:8)
3567 pnodwasmain(1:8) => brick_list(nin,ib)%NODE(1:8)!%NodWasMain
3568 ncell = brick_list(nin,ib)%NBCUT
3569 icell = 0
3570 mcell = brick_list(nin,ib)%MainID
3571 ! nodes wet by main cell
3572 mnod = brick_list(nin,ib)%POLY(mcell)%NumNOD
3573 pnodwasmain(1:8)%NodWasMain = 0
3574 pwherewasmain(1:8)%WhereWasMain = 0
3575 DO j=1,mnod
3576 pnodwasmain(plistnodid(mcell)%p(j))%NodWasMain = 1
3577 ENDDO
3578 ! locate main cell for nodes of a given secnd cell
3579 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
3580 icell = icell +1
3581 IF (icell>ncell .AND. ncell/=0)icell=9
3582 IF(icell == mcell)cycle
3583 ipos = brick_list(nin,ib)%POLY(icell)%WhereIsMain(1) !pwhereismain(icell,1)
3584 mnod = brick_list(nin,ib)%POLY(icell)%NumNOD
3585 !All polyhedron nodes have the same main cell : face=IPOS
3586 DO j=1,mnod
3587 pwherewasmain(plistnodid(icell)%p(j))%WhereWasMain = ipos
3588 enddo!next J
3589 enddo!next ICELL
3590 ENDDO
3591
3592
3593
3594 !------------------------------------------------------------!
3595 ! @34. SET IIAD22(1:NUMELS) = %ELBUG_TAB(NG)%TAG22(I) !
3596 !------------------------------------------------------------!
3597 DO ib=nbf,nbl
3598 ie = brick_list(nin,ib)%ID
3599 iiad22(nin, ie) = ib
3600 ENDDO
3601
3602 !------------------------------------------------------------!
3603 ! @35. SET MLW !
3604 !------------------------------------------------------------!
3605 DO ib=nbf,nbl
3606 mlw = iparg(1,ngb(ib))
3607 brick_list(nin,ib)%MLW = mlw
3608 ENDDO
3609
3610 !------------------------------------------------------------!
3611 ! @36. RESET OLD SECND LIST CONNECTIVITY !
3612 ! WILL BE DEFINED IN I22_GET_PRV_DATA() !
3613 !------------------------------------------------------------!
3614 !RESET ONLY NON ZERO VALUES
3615 DO ib=nbf,nbl
3616 numsecnd = old_secndlist(nin,ib)%Num
3617 IF (numsecnd==0)cycle
3618 old_secndlist(nin,ib)%Num = 0
3619 old_secndlist(nin,ib)%NumSecndNodes = 0
3620 DO j=1,24
3621 old_secndlist(nin,ib)%FM(j) = 0
3622 old_secndlist(nin,ib)%FV(j) = 0
3623 old_secndlist(nin,ib)%IV(j) = 0
3624 old_secndlist(nin,ib)%IBV(j) = 0
3625 old_secndlist(nin,ib)%ICELLv(j) = 0
3626 old_secndlist(nin,ib)%VOL(j) = zero
3627 old_secndlist(nin,ib)%NumNOD_Cell(j) = 0
3628 old_secndlist(nin,ib)%ListNodID(j,1:8) = 0
3629 ENDDO
3630 enddo!next IB
3631
3632 !------------------------------------------------------------!
3633 ! @37. SET main STRONG NODE !
3634 !------------------------------------------------------------!
3635 DO ib=nbf,nbl
3636 mcell = brick_list(nin,ib)%mainID
3637 CALL weighting_cell_nodes(nin,ib,mcell,inod_w, idemerge)
3638 IF(mcell==0)THEN
3639 brick_list(nin,ib)%OldMainStrongNode = 0
3640 ELSE
3641 brick_list(nin,ib)%OldMainStrongNode = inod_w
3642 ENDIF
3643 ENDDO
3644
3645 !------------------------------------------------------------!
3646 ! @38. DRAW SUPERCELL CONNEXION WITH 1D TRUSSES !
3647 !------------------------------------------------------------!
3648 lstilltruss = .true.
3649 igr = ipari(52) !reserved GRTRUS ---> truss group identifier
3650 ntrus = ipari(53 )!reserved GRTRUS ---> truss group NEL
3651 IF( itask==0 .AND. ntrus/=0 )THEN
3652 ii = 1 ! first node of group node
3653 DO ib=nbf,nbl
3654 mcell = brick_list(nin,ib)%mainID
3655 IF (mcell==0 ) cycle
3656 IF (.NOT.lstilltruss) cycle
3657 point0(1:3) = brick_list(nin,ib)%POLY(mcell)%CellCenter(1:3)
3658
3659 !LOOP ON SECND CELLS TO DRAW CONNEXION
3660 numsecnd = brick_list(nin,ib)%SecndList%Num
3661 DO isecnd=1,numsecnd
3662 ibv = brick_list(nin,ib)%SecndList%IBV(isecnd)
3663 icellv = brick_list(nin,ib)%SecndList%ICELLv(isecnd)
3664 pointtmp(1:3) = brick_list(nin,ibv)%POLY(icellv)%CellCenter(1:3)
3665 !check if group has still one truss
3666 IF (ii>ntrus) THEN
3667 lstilltruss=.false.
3668 EXIT
3669 print *, "** Warning inter22 : no more truss in group to mark cell links"
3670 ENDIF
3671 x(1:3,ixt(2,igrtruss(igr)%ENTITY(ii))) = point0(1:3)
3672 x(1:3,ixt(3,igrtruss(igr)%ENTITY(ii))) = pointtmp(1:3)
3673 !write (*,*) "set TRUS_id=", IXT(NIXT,IGRTRUSS(IGR)%ENTITY(II)) ," corr=(", Point0(1:3), ")-(" ,PointTMP(1:3) ,")"
3674 if(ibug22_truss==1)then
3675 write (*,fmt='(A,I10,A,I10,A,I10)') "set TRUS_id=", ixt(nixt,igrtruss(igr)%ENTITY(ii)) ,
3676 . " main=", ixs(11,brick_list(nin,ib)%id) ," secnd=", ixs(11,brick_list(nin,ibv)%id)
3677 endif
3678 ii = ii + 1 !next truss
3679 ENDDO
3680 ENDDO !next IB
3681 DO ii = 1,ntrus
3682 !print *, "reset TRUS_id=", IXT(NIXT,IGRTRUSS(IGR)%ENTITY(II))
3683 x(1:3,ixt(2,igrtruss(igr)%ENTITY(ii))) = (/zero, zero, zero/)
3684 x(1:3,ixt(3,igrtruss(igr)%ENTITY(ii))) = (/ one, zero, zero/)
3685 ENDDO
3686 ENDIF !ITASK==0
3687
3688
3689 !------------------------------------------------------------!
3690 ! @39. COMPUTE DVOL PREDICTION FOR EACH CELL !
3691 !------------------------------------------------------------!
3692 DO ib=nbf,nbl
3693 nbcut = brick_list(nin,ib)%NBCUT
3694 icell = 0
3695 ncell = nbcut
3696 IF(nbcut == 0)THEN
3697 icell = 1
3698 brick_list(nin,ib)%POLY(icell)%DVOL(1) = zero
3699 cycle
3700 ENDIF
3701 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
3702 icell = icell +1
3703 IF (icell>ncell .AND. ncell/=0)icell=9
3704 IF(icell == 9)THEN
3705 brick_list(nin,ib)%POLY(9)%DVOL(1) = -sum(brick_list(nin,ib)%POLY(1:nbcut)%DVOL(1))
3706 ELSE
3707 ! V_face
3708 vsum(1) = brick_list(nin,ib)%PCUT(icell)%Vel(1)
3709 vsum(2) = brick_list(nin,ib)%PCUT(icell)%Vel(2)
3710 vsum(3) = brick_list(nin,ib)%PCUT(icell)%Vel(3)
3711 !VSUM(1) = BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE0%VEL(1)
3712 !VSUM(2) = BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE0%VEL(2)
3713 !VSUM(3) = BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE0%VEL(3)
3714 !not unitary : NN= 2Sn
3715 n_(1) = brick_list(nin,ib)%PCUT(icell)%N(1)
3716 n_(2) = brick_list(nin,ib)%PCUT(icell)%N(2)
3717 n_(3) = brick_list(nin,ib)%PCUT(icell)%N(3)
3718 dvol_predic = dt1*(vsum(1)*n_(1) + vsum(2)*n_(2) + vsum(3)*n_(3))
3719 brick_list(nin,ib)%POLY(icell)%DVOL(1) = dvol_predic !for next CYCLE
3720 ENDIF
3721 enddo!next ICELL
3722 enddo!next ib
3723 !!!!!!!!
3724 !output!
3725 !!!!!!!!
3726 ! IF(IBUG22_PREDICTION/=0)THEN
3727 ! print *, "==PREDICTION DVOL"
3728 ! DO IB=NBF,NBL
3729 ! IF(IBUG22_PREDICTION>0 .AND. IBUG22_PREDICTION/=IXS(11,brick_list(nin,ib)%id))cycle
3730 ! NCELL = BRICK_LIST(NIN,IB)%NBCUT
3731 ! ICELL = 0
3732 ! DO WHILE (ICELL<=NCELL) ! loop on polyhedron {1:NCELL} U {9}
3733 ! ICELL = ICELL +1
3734 ! IF (ICELL>NCELL .AND. NCELL/=0)ICELL = 9
3735 ! write(*,FMT=('(A,I8,A,I1,A,F30.16)'))"brickID=", ixs(11,brick_list(nin,ib)%id),".",ICELL,"=",
3736 ! . BRICK_LIST(NIN,IB)%POLY(ICELL)%DVOL(1)
3737 ! enddo!next ICELL
3738 ! ENDDO!next IB
3739 ! ENDIF
3740 !------------------------------------------------------------!
3741 ! @40. BUILD SUPERCELL DVOL (PREDICTION) !
3742 !------------------------------------------------------------!
3743 !--------------!
3744 CALL my_barrier
3745 !--------------!
3746 !THEN COMPUTE SUPERCELL GEOMETRICAL DVOL
3747 DO ib=nbf,nbl
3748 ncell = brick_list(nin,ib)%NBCUT
3749 mcell = brick_list(nin,ib)%MainID
3750 brick_list(nin,ib)%DVOL = brick_list(nin,ib)%POLY(mcell)%DVOL(1)
3751 !------LOOP--------!
3752 numsecnd = brick_list(nin,ib)%SecndList%Num
3753 DO ic=1,numsecnd
3754 icellv = brick_list(nin,ib)%SecndList%ICELLv(ic)
3755 ibv = brick_list(nin,ib)%SecndList%IBv(ic)
3756 ngv = brick_list(nin,ibv)%NG
3757 idlocv = brick_list(nin,ibv)%IDLOC
3758 iv = brick_list(nin,ibv)%ID !'0' PREDICTION with old conformation
3759 brick_list(nin,ib)%DVOL = brick_list(nin,ib)%DVOL + brick_list(nin,ibv)%POLY(icellv)%DVOL(1) !'1' PREDICTION with new conformation
3760 enddo!next IC
3761 enddo!next IB
3762
3763
3764 !attention faire traitemtn identique en loi 51 pour V1OLD,V2OLD..
3765
3766!!!!!!!!!!!!!!!!!!!!!!!! attention en modifiant ici les vold_scell, on perturbe le critre predictif ci-dessous.
3767
3768
3769
3770
3771 !------------------------------------------------------------!
3772 ! @41. FIX SUPERCELL DVOL (CORRECTION) !
3773 !------------------------------------------------------------!
3774
3775 lskip = .false.
3776
3777c if(IBUG22_dvol /= 0)then
3778c print *, "+---Sinit22_fvm : DVOL majoration"
3779c endif
3780 if(dt1==zero)then
3781 do ib=nbf,nbl
3782 brick_list(nin,ib)%Vold_SCell = brick_list(nin,ib)%Vnew_SCell
3783 enddo
3784 endif
3785CCC RATIO22 = 1.02D00
3786
3787 if(.NOT.lskip)then
3788 DO ib=nbf,1*nbl
3789 icode = brick_list(nin,ib)%ICODE
3790 old_icode = brick_list(nin,ib)%OLD_ICODE
3791 if (brick_list(nin,ib)%Vold_SCell==-111.d00) brick_list(nin,ib)%Vold_SCell = brick_list(nin,ib)%Vnew_SCell !never init.
3792 vnew = brick_list(nin,ib)%Vnew_SCell
3793 vold = brick_list(nin,ib)%Vold_SCell
3794 dvol_numeric = vnew-vold
3795 dvol_predic = brick_list(nin,ib)%dvol
3796 uncutvol = brick_list(nin,ib)%UncutVOL
3797 !IF GEOMETRICAL VALUE IS RISEN ABOVE ITS PREDICTION : CORRECT
3798c if((IBUG22_prediction ==-1 .OR. IBUG22_prediction==ixs(11,brick_list(nin,ib)%id)) )then
3799c print *, "+"
3800c print *, "+------elem_id =",ixs(11,brick_list(nin,ib)%id)
3801c print *, "+--------old_icode =",OLD_ICODE
3802c print *, "+--------icode =",ICODE
3803c print *, "+--------dvol_prediction =",brick_list(nin,ib)%dvol
3804c print *, "+--------dvol_numeric =",VNEW-VOLD
3805c print *, "+--------vnew =",VNEW
3806c print *, "+--------ratio22 =",RATIO22
3807c print *, "+--------DVOL_NUM/V =",ABS(DVOL_NUMERIC)/UncutVOL
3808c endif
3809
3810 IF(abs(dvol_numeric) > ratio22*abs(dvol_predic) .AND. dvol_predic /= zero .AND. ratio22 < 1000 )THEN
3811 IF((icode /= old_icode ) )THEN
3812! IF( (ICODE /= OLD_ICODE ) .OR. (ICODE==OLD_ICODE>ZERO .AND. ABS(DVOL_NUMERIC/DVOL_PREDIC)>1.5d00 ))THEN
3813c if((IBUG22_prediction ==-1 .OR. IBUG22_prediction==ixs(11,brick_list(nin,ib)%id)) )then
3814c print *, "+"
3815c print *, "+------elem_id =",ixs(11,brick_list(nin,ib)%id)
3816c print *, "+--------old_icode =",OLD_ICODE
3817c print *, "+--------icode =",ICODE
3818c print *, "+--------dvol_prediction =",brick_list(nin,ib)%dvol
3819c print *, "+--------dvol_numeric =",VNEW-VOLD
3820c print *, "+--------vnew =",VNEW
3821c endif
3822 brick_list(nin,ib)%Vold_SCell = brick_list(nin,ib)%Vnew_SCell - dvol_predic
3823 ng = brick_list(nin,ib)%NG
3824 idloc = brick_list(nin,ib)%IDLOC
3825 gbuf => elbuf_tab(ng)%GBUF
3826C write (*,FMT='(A,F30.16)' ) " UncutVOL =", brick_list(nin,ib)%uncutvol
3827C write (*,FMT='(A,F30.16)' ) " ABS(DVOL_NUMERIC)/UncutVOL =",ABS(DVOL_NUMERIC)/UncutVOL
3828 !################################################################!
3829 !# TRAITEMENT MULTIMATERIAU LAW51 #!
3830 !################################################################!
3831 mlw = brick_list(nin,ib)%MLW
3832 ng = brick_list(nin,ib)%NG
3833 idloc = brick_list(nin,ib)%IDLOC
3834 IF(mlw==51)THEN
3835 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
3836 llt = iparg(2,ng)
3837 DO itrimat=1,trimat
3838 !IPOS=01 : fractions volumiques
3839 !IPOS=11 : anciens volumes de chaque phase
3840 ipos = 1
3841 k = llt * (m51_n0phas + (itrimat-1)*m51_nvphas +ipos-1)
3842 vfrac(itrimat) = mbuf%VAR(k+idloc)
3843 ipos = 11
3844 k = llt * (m51_n0phas + (itrimat-1)*m51_nvphas +ipos-1)
3845 volcell51_old(itrimat) = mbuf%VAR(k+idloc)
3846 mbuf%VAR(k+idloc) = max(zero, (vnew-dvol_predic)*vfrac(itrimat) )
3847 volcell51(itrimat) = mbuf%VAR(k+idloc)
3848 enddo!next ITRIMAT
3849 ENDIF
3850 if(ibug22_prediction ==-1 .OR. ibug22_prediction==ixs(11,brick_list(nin,ib)%id))then
3851c print *, "+--------vold =",VOLD ,"->", BRICK_LIST(NIN,IB)%Vold_SCell
3852! print *, "+--------vold_1 =",VolCELL51_OLD(1) ,"->", VolCELL51(1)
3853! print *, "+--------vold_2 =",VolCELL51_OLD(2) ,"->", VolCELL51(2)
3854! print *, "+--------vold_3 =",VolCELL51_OLD(3) ,"->", VolCELL51(3)
3855! print *, "+--------vold_4 =",volcell51_old(4) ,"->", volcell51(4)
3856 endif
3857 endif!IF(ICODE /= OLD_ICODE .AND. OLD_ICODE /= 0)
3858 endif!IF(ABS(DVOL_NUMERIC) > RATIO22*ABS(DVOL_PREDIC) .AND. RATIO22 < 1000 .AND. ABS(DVOL_NUMERIC)/UncutVOL>0.05)
3859 enddo!next IB
3860 endif
3861
3862 IF( (ibug22_prediction ==-1 ))THEN
3863 call my_barrier
3864 IF(itask==0)THEN
3865 DO ib=1,nb
3866 icode = brick_list(nin,ib)%ICODE
3867 old_icode = brick_list(nin,ib)%OLD_ICODE
3868 vnew = brick_list(nin,ib)%Vnew_SCell
3869 vold = brick_list(nin,ib)%Vold_SCell
3870 dvol_numeric = vnew-vold
3871 dvol_predic = brick_list(nin,ib)%dvol
3872 uncutvol = brick_list(nin,ib)%UncutVOL
3873 print *, "+------elem_id =",ixs(11,brick_list(nin,ib)%id)
3874 print *, "+--------old_icode =",old_icode
3875 print *, "+--------icode =",icode
3876 print *, "+--------dvol_prediction =",brick_list(nin,ib)%dvol
3877 print *, "+--------dvol_numeric =",vnew-vold
3878 print *, "+--------vnew =",vnew
3879 print *, "+--------vold =",vold ,"->", brick_list(nin,ib)%Vold_SCell
3880 ENDDO
3881 ENDIF
3882 ENDIF
3883
3884
3885 !------------------------------------------------------------!
3886 ! @42. TIME ZERO VOLUME INIT (NOT DONE IN STARTER) !
3887 ! needed for convection subroutines !
3888 ! sforc3>voln22 not called before convection !
3889 !------------------------------------------------------------!
3890 IF(dt1==zero)THEN
3891 DO ib=nbf,nbl
3892 psubvol(1:9) => brick_list(nin,ib)%POLY(1:9)
3893 psubvold => brick_list(nin,ib)%Vold_SCell
3894 mcell = brick_list(nin,ib)%MainID
3895 IF(mcell==0)cycle
3896 psubvold = brick_list(nin,ib)%Vnew_Scell
3897 ENDDO
3898 ENDIF
3899 DO ib=nbf,nbl
3900 psubvol(1:9) => brick_list(nin,ib)%POLY(1:9)
3901 psubvold => brick_list(nin,ib)%Vold_SCell
3902 mcell = brick_list(nin,ib)%MainID
3903 !IF(MCELL==0)CYCLE
3904 gbuf => elbuf_tab(ngb(ib))%GBUF
3905 IF(psubvold>zero)gbuf%VOL(idlocb(ib)) = psubvold
3906C write (*,FMT='(A,I10,A,F30.16,A,F30.16)')"brickID=", IXS(11,brick_list(nin,ib)%id),
3907C . " vold=",brick_list(nin,ib)%vold_scell," vnew=",brick_list(nin,ib)%vnew_scell
3908 ENDDO
3909
3910
3911 !------------------------------------------------------------!
3912 ! @43. DEBUG OUTPUT !
3913 !------------------------------------------------------------!
3914 if (debug_outp .AND. ibug22_sinit/=0)then
3915 do ib=nbf,nbl
3916 psubvol(1:9) => brick_list(nin,ib)%POLY(1:9)
3917 plistnodid(1)%p(1:8) => brick_list(nin,ib)%POLY(1)%ListNodID(1:8)
3918 plistnodid(2)%p(1:8) => brick_list(nin,ib)%POLY(2)%ListNodID(1:8)
3919 plistnodid(3)%p(1:8) => brick_list(nin,ib)%POLY(3)%ListNodID(1:8)
3920 plistnodid(4)%p(1:8) => brick_list(nin,ib)%POLY(4)%ListNodID(1:8)
3921 plistnodid(5)%p(1:8) => brick_list(nin,ib)%POLY(5)%ListNodID(1:8)
3922 plistnodid(6)%p(1:8) => brick_list(nin,ib)%POLY(6)%ListNodID(1:8)
3923 plistnodid(7)%p(1:8) => brick_list(nin,ib)%POLY(7)%ListNodID(1:8)
3924 plistnodid(8)%p(1:8) => brick_list(nin,ib)%POLY(8)%ListNodID(1:8)
3925 plistnodid(9)%p(1:8) => brick_list(nin,ib)%POLY(9)%ListNodID(1:8)
3926 !pAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%Adjacent_Cell(1:5)
3927 !pNAdjCELL => BRICK_LIST(NIN,IB)%POLY(1:9)%FACE(1:6)%NAdjCell
3928 padjbrick => brick_list(nin,ib)%Adjacent_Brick(1:6,1:5)
3929 pismain(1:9) => brick_list(nin,ib)%POLY(1:9)
3930 pmainid => brick_list(nin,ib)%MainID
3931 id = brick_list(nin,ib)%ID
3932 icell = 0
3933 mcell = mcell
3934 ncell = brick_list(nin,ib)%NBCUT
3935 gbuf => elbuf_tab(ngb(ib))%GBUF
3936 vol = gbuf%VOL(idlocb(ib))
3937 brickid = idb(ib)
3938 if(ibug22_sinit/=ixs(11,brickid) .AND. ibug22_sinit/=-1) cycle
3939 write (*,fmt='(A,I12)') "+=== BRICK ID===", ixs(11,id)
3940 if(ncell==0)then
3941 write (*,fmt='(A )') "| uncut: "
3942 write (*,fmt='(A,1F30.20)') "| volume: ", vol
3943 write (*,fmt='(A,1F30.20)') "| ext. volume: ", brick_list(nin,ib)%Vnew_Scell
3944 write (*,fmt='(A,1F30.20)') "| masse: ", gbuf%VOL(idlocb(ib)) * gbuf%RHO(idlocb(ib))
3945 if(brick_list(nin,ib)%SecndList%Num>0)then
3946 do j=1,brick_list(nin,ib)%SecndList%Num
3947 write (*,fmt='(A )') "| secnd list : "
3948 write (*,fmt='(A,I10 )') "| + J : ", brick_list(nin,ib)%SecndList%FM(j)
3949 write (*,fmt='(A,I10 )') "| + IB : ", brick_list(nin,ib)%SecndList%IBv(j)
3950 write (*,fmt='(A,I10 )') "| + brickID : ", ixs(11,brick_list(nin,ib)%SecndList%IV(j))
3951 enddo
3952 endif
3953 cycle
3954 endif
3955 write (*,fmt='(a,1f30.20)') "| volume: ", vol
3956 write (*,FMT='(a,6f30.20)') "| faces: ", F(1:6,IB)
3957 write (*,FMT='(a,1f30.20)') "| masse: ", GBUF%VOL(IDLOCB(IB)) * GBUF%RHO(IDLOCB(IB))
3958 DO WHILE (ICELL<=NCELL) ! loop on polyhedron {1:NCELL} U {9}
3959 ICELL = ICELL +1
3960.AND. IF (ICELL>NCELL NCELL/=0)ICELL=9
3961 debugMAINSECND = '.........|'
3962 mnod = BRICK_LIST(NIN,IB)%POLY(icell)%NumNOD
3963 write (*,FMT='(a )') "|"
3964 if(icell/=9)then
3965 write (*,FMT='(a,i1,a,a,a1,i2,a,a6)')
3966 . "+== ICELL= ", ICELL , ", SecType=", BRICK_LIST(NIN,IB)%SECTYPE(ICELL) ,
3967 . "(", BRICK_LIST(NIN,IB)%SecID_Cell(ICELL) , ") - ", Char1
3968 else
3969 write (*,FMT='(a,a6)') "+== REMAINING POLYHEDRON - ", Char1
3970 endif
3971
3972 write (*,FMT='(a )') "| |"
3973 write (*,FMT='(a,i1)') "| +===Main/Secnd=", pIsMain(ICELL)%IsMain
3974 write (*,FMT='(a,f30.20)') "| +======SUVOLUMES=", pSUBVOL(ICELL)%Vnew
3975 write (*,FMT='(a,6f30.20)') "| +=======SUBFACES=", BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(1:6)%Surf
3976 write (*,FMT='(a,f30.20)') "| +=======CUT AERA=", brick_list(nin,ib)%PCUT(icell)%SCUT(1)
3977 write (*,fmt='(A,A,I2)') "| +======NUM POINT=", " ",brick_list(nin,ib)%POLY(icell)%NumPOINT
3978 write (*,fmt='(A,A,I1,A,8I12)') "| +======NODE LIST=",
3979 . " (",mnod,") ", plistnodid(icell)%p(1:mnod)
3980 write (*,fmt='(A,A,8I12)') "| | radIDs=",
3981 . " ", ixs(1+plistnodid(icell)%p(1:mnod),id)
3982 write (*,fmt='(A,A,8I12)') "| | userIDs=",
3983 . " ", itab(ixs(1+plistnodid(icell)%p(1:mnod),id))
3984 iad2 = ale_connectivity%ee_connect%iad_connect(brickid)
3985 lgth2 = ale_connectivity%ee_connect%iad_connect(brickid+1) -
3986 . ale_connectivity%ee_connect%iad_connect(brickid)
3987 If(sum(iabs(ale_connectivity%ee_connect%connected(iad2:iad2 + lgth2 - 1)))/=0)then
3988 write (*,fmt='(A,6I10)') "| +===Adj Brick(i)=", padjbrick(1:6,1)
3989 DO j=1,6
3990 IF( padjbrick(j,1)/=0 )THEN
3991 write (*,fmt='(A,6I10)') "| +===Adj Brick(u)=", ixs(11,padjbrick(j,1))
3992 ENDIF
3993 ENDDO
3994 DO j=1,6
3995 nadjcell = brick_list(nin,ib)%POLY(icell)%FACE(j)%NAdjCell
3996 IF(nadjcell/=0)THEN
3997 write (*,fmt='(A,I1,A,5I3)') "| +======Adj Cells, face=",j, " :",
3998 . brick_list(nin,ib)%POLY(icell)%FACE(j)%Adjacent_Cell(1:nadjcell)
3999 ENDIF
4000 ENDDO
4001 ENDIF
4002 enddo!next icell
4003 write (*,fmt='(A )') " "
4004 enddo!next IB
4005 if(itask==0.AND.debug_outp)then
4006 write (*,fmt='(A )') " ----sini22_end----"
4007 write (*,fmt='(A )') " "
4008 endif
4009 endif! (IBUG22_sinit/=0)
4010
4011 !------------------------------------------------------------!
4012 ! @44. DEBUG - WRITE CUT CELL BUFFER !
4013 !------------------------------------------------------------!
4014 ! CALL WRITE_CUT_CELL_BUFFER() !post/debug
4015
4016 !------------------------------------------------------------!
4017 ! @45. SUPERCELL CENTERS !
4018 !------------------------------------------------------------!
4019 DO ib=nbf,nbl
4020 ncell = brick_list(nin,ib)%NBCUT
4021 numsecnd = brick_list(nin,ib)%SecndList%Num
4022 mainid = brick_list(nin,ib)%mainID
4023 IF (mainid ==0) cycle
4024 volcell = brick_list(nin,ib)%POLY(mainid)%Vnew
4025 vol = volcell !cumul
4026 ppoint(1:3) = brick_list(nin,ib)%POLY(mainid)%CellCenter(1:3) * volcell
4027 DO ic=1,numsecnd
4028 icellv = brick_list(nin,ib)%SecndList%ICELLv(ic)
4029 ibv = brick_list(nin,ib)%SecndList%IBv(ic)
4030 volcell = brick_list(nin,ibv)%POLY(icellv)%Vnew
4031 point0(1:3) = brick_list(nin,ibv)%POLY(icellv)%CellCenter(1:3)
4032 ppoint(1) = ppoint(1) + volcell * point0(1)
4033 ppoint(2) = ppoint(2) + volcell * point0(2)
4034 ppoint(3) = ppoint(3) + volcell * point0(3)
4035 vol = vol + volcell
4036 enddo!next IC
4037 ppoint(1) = ppoint(1) / vol
4038 ppoint(2) = ppoint(2) / vol
4039 ppoint(3) = ppoint(3) / vol
4040 brick_list(nin,ib)%SCellCenter(1:3) = ppoint(1:3)
4041 enddo!next IB
4042
4043 !------------------------------------------------------------!
4044 ! @46. MARK SUPER-CELL CENTERS WITH ORPHAN NODES !
4045 !------------------------------------------------------------!
4046 lstillnode = .true.
4047 IF(ipari(70)/=0)THEN
4048 nnodes = igrnod(ipari(70))%NENTITY
4049 !!!WRITE(*,*)(ITAB(IBUFSSG(J)),J=IAD0,IAD0+NNODES-1)
4050 IF( itask==0 .AND. nnodes/=0 )THEN
4051 ii = 1 ! first node of group node
4052 DO ib=nbf,nbl
4053 mcell = brick_list(nin,ib)%mainID
4054 IF (mcell==0 ) cycle
4055 IF (.NOT.lstillnode) cycle
4056 point0(1:3) = brick_list(nin,ib)%SCellCenter(1:3)
4057 IF (ii>nnodes) THEN
4058 lstillnode = .false.
4059 print *, "** Warning inter22 : no more node in group to mark cell center. Last one was" ,
4060 . itab(igrnod(ipari(70))%ENTITY(nnodes))
4061 EXIT
4062 ENDIF
4063 x(1:3,igrnod(ipari(70))%ENTITY(ii)) = point0(1:3)
4064 if(ibug22_orphannodes == 1)then
4065 write (*,fmt='(A,I10,A,I10,A,I10)')"set orphan_node_id=",itab(igrnod(ipari(70))%ENTITY(ii))
4066 endif
4067 ii = ii + 1 !next node
4068 ENDDO !next IB
4069 DO ii = 1, nnodes
4070 !print *, "reset SPHCEL_id=", IXT(NIXT,IGRNOD(IPARI(70))%ENTITY(II))
4071 x(1:3,igrnod(ipari(70))%ENTITY(ii)) = (/zero, zero, zero/)
4072 ENDDO
4073 ENDIF !ITASK==0
4074 ENDIF
4075
4076 !------------------------------------------------------------!
4077 ! @47. UNCUT CELLS + POLY 9 : CENTERS !
4078 !------------------------------------------------------------!
4079 !FOR UNCUT BRICKS, Centers for polyehedra faces computed in i22subvol
4080 DO ib=nbf,nbl
4081 ncell = brick_list(nin,ib)%NBCUT
4082 IF(ncell==0)THEN
4083 ie = brick_list(nin,ib)%ID
4084 nc(1:8) = ixs(2:9,ie)
4085 DO j=1,6
4086 brick_list(nin,ib)%POLY(1)%FACE(j)%Center(1) = fourth * sum( x(1, nc(icf(1:4,j)) ) )
4087 brick_list(nin,ib)%POLY(1)%FACE(j)%Center(2) = fourth * sum( x(2, nc(icf(1:4,j)) ) )
4088 brick_list(nin,ib)%POLY(1)%FACE(j)%Center(3) = fourth * sum( x(3, nc(icf(1:4,j)) ) )
4089 ENDDO
4090 ELSE
4091 !simplification, car approximation precendente invalide lorsque limit(FACE9) = 0
4092 ie = brick_list(nin,ib)%ID
4093 nc(1:8) = ixs(2:9,ie)
4094 !
4095 !ICELL 9 only since polyhedra face centers were computed in i22subvol.
4096 ! Poly9 is complementary polyhedron, it is built without graph but deduced by boolean operation from full brick.
4097 ! its face centers must be computed to display face center (velocity post-treatment)
4098 ! face center are only known for cut polyhedra (A & C below), center for poly9 must be deduced (poly B below)
4099 !
4100 ! N4 N3
4101 ! !------------+-------!
4102 ! ! | I3 !
4103 ! ! | !
4104 ! ! | ! Ni : nodes
4105 ! ! | ! Ii : intersection points
4106 ! ! I4 B | C ! A,C : cut polyhedra
4107 ! !\ | ! B : complemntary polyhedron
4108 ! ! \ | ! Ca,Cb,Cc : face centers for each A,B,C polyhedra
4109 ! ! \ | ! Npa,Npb,Npc : number of points for polygon on face
4110 ! ! A \ I1 | I2 !
4111 ! !----+-------|-------!
4112 ! N1 N2
4113 !
4114 ! Ca = (N1+I1+I4)/3
4115 ! Cb = (N4+I1+I2+I3+I4)/5
4116 ! Cc = (N2+N3+I2+I3)/4
4117 ! 3Ca + 5Cb + 4Cc = Sum(Ni) + 2 Sum(Ii)
4118 ! => poly9 C9 = [-Npa.Ca -Npc.Cc + Sum(Ni) +2.Sum(Ii)] / Np9
4119 !
4120 !
4121 Do j=1,6
4122 face = brick_list(nin,ib)%POLY(9)%FACE(j)%Surf
4123 np_(9) = brick_list(nin,ib)%POLY(9)%FACE(j)%NumPOINT
4124 IF(abs(face)<=em10 .OR. np_(9)==0)THEN
4125 brick_list(nin,ib)%POLY(9)%FACE(j)%Center(1:3) = zero
4126 ELSE
4127 ncell = brick_list(nin,ib)%NBCUT
4128 DO icell = 1 ,ncell
4129 np_(icell) = brick_list(nin,ib)%POLY(icell)%FACE(j)%NumPOINT
4130 center(1:3,icell) = brick_list(nin,ib)%POLY(icell)%FACE(j)%Center(1:3)
4131 ENDDO
4132 !NP_ is number of points (intersection + nodes)
4133 ! pPoint is sum of nodes
4134 !CUT_point is sum of intersection points
4135 ! Center : are face center
4136 ! Point0 : is center of poly 9
4137 np_(ncell+1:9) = 0
4138 np_(9) = brick_list(nin,ib)%POLY(9)%FACE(j)%NumPOINT
4139 ppoint(1) = sum( x(1, nc(icf(1:4,j)) ) )
4140 ppoint(2) = sum( x(2, nc(icf(1:4,j)) ) )
4141 ppoint(3) = sum( x(3, nc(icf(1:4,j)) ) )
4142 cut_point(1:3) = brick_list(nin,ib)%POLY(9)%FACE(j)%Center(1:3) !stored there in i22subvol.F it is SUM(cutpoints on the face J)
4143 point0(1) = ppoint(1) + two*cut_point(1)
4144 . - np_(1)*center(1,1)- np_(2)*center(1,2)- np_(3)*center(1,3)- np_(4)*center(1,4)
4145 . - np_(5)*center(1,5)- np_(6)*center(1,6)- np_(7)*center(1,7)- np_(8)*center(1,8)
4146 point0(2) = ppoint(2) + two*cut_point(2)
4147 . - np_(1)*center(2,1)- np_(2)*center(2,2)- np_(3)*center(2,3)- np_(4)*center(2,4)
4148 . - np_(5)*center(2,5)- np_(6)*center(2,6)- np_(7)*center(2,7)- np_(8)*center(2,8)
4149 point0(3) = ppoint(3) + two*cut_point(3)
4150 . - np_(1)*center(3,1)- np_(2)*center(3,2)- np_(3)*center(3,3)- np_(4)*center(3,4)
4151 . - np_(5)*center(3,5)- np_(6)*center(3,6)- np_(7)*center(3,7)- np_(8)*center(3,8)
4152 point0(1:3) = point0(1:3) / np_(9)
4153 brick_list(nin,ib)%POLY(9)%FACE(j)%Center(1:3) = point0(1:3)
4154 !print *, "poly9 center, ixs(11,IE), face J ", ixs(11,IE), J
4155 !print *, BRICK_LIST(NIN,IB)%POLY(9)%FACE(J)%Center(1:3)
4156 ENDIF
4157 enddo!next J
4158 ENDIF
4159 enddo!next IB
4160
4161
4162
4163 !------------------------------------------------------------!
4164 ! @48. MARK CELL CENTERS WITH ORPHAN NODES !
4165 !------------------------------------------------------------!
4166 lstillnode = .true.
4167 IF(ipari(81)/=0)THEN
4168 nnodes = igrnod(ipari(81))%NENTITY
4169 !!!WRITE(*,*)(ITAB(IBUFSSG(J)),J=IAD0,IAD0+NNODES-1)
4170 IF( itask==0 .AND. nnodes/=0 )THEN
4171 ii = 1 ! first node of group node
4172 DO ib=nbf,nbl
4173 icell = 0
4174 ncell = brick_list(nin,ib)%NBCUT
4175 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
4176 icell = icell +1
4177 IF (icell>ncell .AND. ncell/=0)icell=9
4178 IF (.NOT.lstillnode) cycle
4179 point0(1:3) = brick_list(nin,ib)%POLY(icell)%CellCenter(1:3)
4180 IF(ii>nnodes)THEN
4181 lstillnode = .false.
4182 print *, "** Warning inter22 : no more node in group to mark cell center",
4183 . itab(igrnod(ipari(81))%ENTITY(nnodes))
4184 EXIT
4185 ENDIF
4186 brick_list(nin,ib)%POLY(icell)%ID_FREE_NODE = igrnod(ipari(81))%ENTITY(ii)
4187 x(1:3,igrnod(ipari(81))%ENTITY(ii)) = point0(1:3)
4188 if(ibug22_orphannodes == 1)then
4189 write (*,fmt='(A,I10,A,I10,A,I10)')"set orphan_node_id=",
4190 . itab(igrnod(ipari(81))%ENTITY(ii)),"in brick_id=",ixs(11,brick_list(nin,ib)%id)
4191 endif
4192 ii = ii + 1 !next node
4193 ENDDO
4194 ENDDO !next IB
4195 DO ii = 1,nnodes
4196 !print *, "reset SPHCEL_id=", IXT(NIXT,IGRNOD(IPARI(81))%ENTITY(II))
4197 x(1:3,igrnod(ipari(81))%ENTITY(ii)) = zero
4198 ENDDO
4199 ENDIF !ITASK==0
4200 ENDIF
4201
4202
4203 !------------------------------------------------------------!
4204 ! @49. MARK FACE CENTERS WITH ORPHAN NODES !
4205 !------------------------------------------------------------!
4206 IF(ipari(82)/=0)THEN
4207 nnodes = igrnod(ipari(82))%NENTITY
4208 IF(nnodes==0)RETURN
4209 ii = 1 ! first node of group node
4210 DO ib=nbf,nbl
4211 ie = brick_list(nin,ib)%ID
4212 icell = 0
4213 ncell = brick_list(nin,ib)%NBCUT
4214 DO WHILE (icell<=ncell) ! loop on polyhedron {1:NCELL} U {9}
4215 icell = icell +1
4216 IF (icell>ncell .AND. ncell/=0)icell=9
4217 IF(.NOT.lstillnode) cycle
4218 DO j=1, 6
4219 IF (ii>nnodes) THEN
4220 lstillnode = .false.
4221 print *, "** Warning inter22 : no more node in group to mark face centers." ,
4222 . itab(igrnod(ipari(82))%ENTITY(nnodes))
4223 EXIT
4224 ENDIF
4225 node_id = igrnod(ipari(82))%ENTITY(ii)
4226 x(1:3,node_id) = brick_list(nin,ib)%POLY(icell)%FACE(j)%Center(1:3)
4227 ii = ii + 1 !next node
4228 ENDDO
4229 ENDDO !next ICELL
4230 enddo!next IB
4231 ENDIF
4232
4233
4234
4235 !print *, "SINIT22_FVM : NB=", NBL
4236 !do ib=nbf,nbl
4237 ! id = brick_list(nin,ib)%id
4238 ! print *, "ixs =", ixs(11,id)
4239 ! print *, "nbcut = ", brick_list(1,ib)%nbcut
4240 ! print *, "sectype1 =", brick_list(1,ib)%sectype(1)
4241 ! print *, "sectype2 =", brick_list(1,ib)%sectype(2)
4242 ! print *, "---------------------"
4243 !enddo
4244
4245
4246 !------------------------------------------------------------!
4247 ! @50. DEALLOCATE !
4248 !------------------------------------------------------------!
4249 IF(ALLOCATED(debugmainsecndv))DEALLOCATE (debugmainsecndv)
4250 IF(ALLOCATED(ismainv))DEALLOCATE (ismainv)
4251 IF(ALLOCATED(f))DEALLOCATE (f)
4252 IF(ALLOCATED(vol51))DEALLOCATE (vol51,vol51v)
4253 IF(ALLOCATED(origin_data))DEALLOCATE (origin_data)
4254 IF(ALLOCATED(destroy_data))DEALLOCATE (destroy_data)
4255 IF(ALLOCATED(norigin))DEALLOCATE (norigin)
4256
4257
4258 RETURN
4259 END
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
i22aera
function i22aera(npts, p, c)
Definition i22subvol.F:2380
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
for
for(i8=*sizetab-1;i8 >=0;i8--)
Definition mumps_common.c:91
ale_connectivity_mod
Definition ale_connectivity_mod.F:223
alefvm_mod
Definition alefvm_mod.F:88
alefvm_mod::alefvm_buffer
type(alefvm_buffer_), target alefvm_buffer
Definition alefvm_mod.F:120
groupdef_mod
Definition groupdef_mod.F:662
i22bufbric_mod
Definition cut-cell-search_mod.F:95
i22bufbric_mod::brick_list
type(brick_entity), dimension(:,:), allocatable, target brick_list
Definition cut-cell-search_mod.F:105
i22tri_mod
Definition cut-cell-search_mod.F:212
i22tri_mod::ibug22_truss
integer ibug22_truss
Definition cut-cell-search_mod.F:372
i22tri_mod::nb
integer nb
Definition cut-cell-search_mod.F:487
i22tri_mod::ibug22_ccbuflist
integer ibug22_ccbuflist
Definition cut-cell-search_mod.F:428
i22tri_mod::ibug22_orphannodes
integer ibug22_orphannodes
Definition cut-cell-search_mod.F:373
i22tri_mod::ibug22_dtmin
integer ibug22_dtmin
Definition cut-cell-search_mod.F:383
i22tri_mod::ibug22_undirectlink
integer ibug22_undirectlink
Definition cut-cell-search_mod.F:357
i22tri_mod::ibug22_link_switch
integer ibug22_link_switch
Definition cut-cell-search_mod.F:408
i22tri_mod::ibug22_prediction
integer ibug22_prediction
Definition cut-cell-search_mod.F:418
i22tri_mod::ibug22_merge
integer ibug22_merge
Definition cut-cell-search_mod.F:352
i22tri_mod::ibug22_sinit
integer ibug22_sinit
Definition cut-cell-search_mod.F:316
initbuf_mod
Definition initbuf.F:160
sigeps37_single_cell
subroutine sigeps37_single_cell(elbuf_tab, ixs, bufmat, iparg, ipm, idloc, ng, brickid, vol)
Definition sigeps37_single_cell.F:36
sinit22_fvm
subroutine sinit22_fvm(ixs, elbuf_tab, iparg, itab, itask, bufbric, nbric_l, x, ale_connectivity, v, nv46, veul, igrnod, ipari, igrtruss, ixt, bufmat, ipm)
Definition sinit22_fvm.F:52
main
int main(int argc, char *argv[])
Definition CMakeCCompilerId.c:841
my_barrier
subroutine my_barrier
Definition machine.F:31
ale_connectivity_mod::t_ale_connectivity
Definition ale_connectivity_mod.F:254
groupdef_mod::group_
Definition groupdef_mod.F:700
weighting_cell_nodes
subroutine weighting_cell_nodes(nin, ib, icell, ires, idemerge)
Definition weighting_cell_nodes.F:32