OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i22wetsurf.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!|| i22wetsurf ../engine/source/interfaces/int22/i22wetsurf.F
25!||--- called by ------------------------------------------------------
26!|| i22mainf ../engine/source/interfaces/int22/i22mainf.F
27!||--- calls -----------------------------------------------------
28!|| getprojectedfacetype ../engine/source/interfaces/int22/i22wetsurf.F
29!|| polygonalclipping ../engine/source/interfaces/int22/i22clip_tools.F
30!|| setcounterclockwisepolyg ../engine/source/interfaces/int22/i22clip_tools.F
31!||--- uses -----------------------------------------------------
32!|| element_mod ../common_source/modules/elements/element_mod.F90
33!|| i22bufbric_mod ../common_source/modules/interfaces/cut-cell-search_mod.F
34!|| i22tri_mod ../common_source/modules/interfaces/cut-cell-search_mod.F
35!||====================================================================
36 SUBROUTINE i22wetsurf(
37 1 JLT ,CAND_B ,CAND_E ,CB_LOC ,CE_LOC ,
38 2 X1 ,X2 ,X3 ,X4 ,Y1 ,
39 3 Y2 ,Y3 ,Y4 ,Z1 ,Z2 ,
40 4 Z3 ,Z4 ,XI ,YI ,ZI ,
41 5 NX1 ,NX2 ,NX3 ,NX4 ,NY1 ,
42 6 NY2 ,NY3 ,NY4 ,NZ1 ,NZ2 ,
43 7 NZ3 ,NZ4 ,LB1 ,LB2 ,LB3 ,
44 8 LB4 ,LC1 ,LC2 ,LC3 ,LC4 ,
45 9 P1 ,P2 ,P3 ,P4 ,IX1 ,
46 A IX2 ,IX3 ,IX4 ,NSVG ,STIF ,
47 B JLT_NEW,GAPV ,INACTI ,CAND_P ,N_SCUT ,
48 C INDEX ,VXI ,VYI ,
49 D VZI ,MSI ,KINI ,SURF ,IBAG ,
50 E ITAB ,IRECT ,I_STOK ,IXS ,NFT ,
51 F CoG ,Seff ,Delta ,X)
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
58C-----------------------------------------------
59C M o d u l e s
60C-----------------------------------------------
61 USE i22bufbric_mod
62 USE i22tri_mod
63 use element_mod , only : nixs
64C-----------------------------------------------
65C I m p l i c i t T y p e s
66C-----------------------------------------------
67#include "implicit_f.inc"
68#include "comlock.inc"
69#include "inter22.inc"
70#include "com04_c.inc"
71C-----------------------------------------------
72C G l o b a l P a r a m e t e r s
73C-----------------------------------------------
74#include "mvsiz_p.inc"
75C-----------------------------------------------
76C D u m m y A r g u m e n t s
77C-----------------------------------------------
78 INTEGER JLT, JLT_NEW,INACTI, CAND_B(*),CB_LOC(MVSIZ),NFT,
79 . CAND_E(*),CE_LOC(MVSIZ), NSVG(MVSIZ), KINI(*), IXS(NIXS,*)
80 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
81 . INDEX(MVSIZ),IBAG, ITAB(*),IRECT(4,*),I_STOK
82 my_real
83 . NX1(MVSIZ), NX2(MVSIZ), NX3(MVSIZ), NX4(MVSIZ),
84 . NY1(MVSIZ), NY2(MVSIZ), NY3(MVSIZ), NY4(MVSIZ),
85 . NZ1(MVSIZ), NZ2(MVSIZ), NZ3(MVSIZ), NZ4(MVSIZ),
86 . LB1(MVSIZ), LB2(MVSIZ), LB3(MVSIZ), LB4(MVSIZ),
87 . LC1(MVSIZ), LC2(MVSIZ), LC3(MVSIZ), LC4(MVSIZ),
88 . X1(MVSIZ), X2(MVSIZ), X3(MVSIZ), X4(MVSIZ),
89 . Y1(MVSIZ), Y2(MVSIZ), Y3(MVSIZ), Y4(MVSIZ),
90 . Z1(MVSIZ), Z2(MVSIZ), Z3(MVSIZ), Z4(MVSIZ),
91 . XI(MVSIZ), YI(MVSIZ), ZI(MVSIZ), STIF(MVSIZ),
92 . P1(MVSIZ), P2(MVSIZ), P3(MVSIZ), P4(MVSIZ),
93 . gapv(mvsiz), cand_p(*),
94 . vxi(mvsiz), vyi(mvsiz), vzi(mvsiz), msi(mvsiz),
95 . surf(3,*), n_scut(1:3,nbcut_max,mvsiz),cog(3,nbcut_max,mvsiz),
96 . seff(nbcut_max,mvsiz),delta(4,nbcut_max,mvsiz),x(3,numnod)
97
98C-----------------------------------------------
99C L o c a l V a r i a b l e s
100C-----------------------------------------------
101
102 TYPE polygon
103 my_real :: node(16,2) !index1 : numpt, index2:coor
104 INTEGER :: NumNodes
105 END TYPE
106
107 INTEGER :: I, J, NIN, IB, IFT,ILT, K, Q, L
108 INTEGER :: NBCUT , IE, IEDG_LAG, IEDG_SCUT
109 INTEGER :: ICUT, IDX_LAG(7), IDX_SCUT(7)
110 INTEGER :: TAG_SCUT(4,NBCUT_MAX), TAG_NOD(4,NBCUT_MAX) !taging lagrangian node regarding Pcut (proj criteria), 4 nodes, 2= max(nbcut)
111 INTEGER :: IS_ON_SCUT(NBCUT_MAX), NP, Snod_On_Lag(6)
112C INTEGER :: TAG_EL(NBCUT_MAX)
113 INTEGER :: IADD(4,6), NCUT(NBCUT_MAX), NP_RECT(MVSIZ),ITMP, NN
114 INTEGER :: HULL_PT, HULL_IDX(24)
115 INTEGER :: NP_INTER, NP_LAG, NP_SURF, NODE_NUM,brickID, ListP(4), NBCUT_ADD
116 my_real :: NODE_XY(2,24), NODE_XYZ(3,24), PS, MARGE
117 my_real :: dist_pcut(nbcut_max, 4,i_stok) !2 pcuts, 4 points, ncande couples
118 my_real :: dmin
119 my_real :: normal(3), basis(3), ph(3,nbcut_max,4),p(3,4), a,b,c,d , a2,b2,c2 ,a2b2, l2, dist(4), pph(3,4)
120 my_real :: distance(4), edge(4), area_tria(4)
121 my_real :: sh(3,nbcut_max,8), s(3,8)
122 my_real :: scut_pts(3,6) !Scut Points
123 my_real :: k1(2),k2(2),k3(2),k4(2)
124 my_real :: point(3),ratio,cog3,p_grav(3,nbcut_max,mvsiz),ss2(4), s2, sel(mvsiz)
125 my_real
126 . x0(mvsiz), y0(mvsiz), z0(mvsiz),
127 . x01(mvsiz), x02(mvsiz), x03(mvsiz), x04(mvsiz),
128 . y01(mvsiz), y02(mvsiz), y03(mvsiz), y04(mvsiz),
129 . z01(mvsiz), z02(mvsiz), z03(mvsiz), z04(mvsiz)
130 my_real
131 . xi0,sx0,
132 . yi0,sy0,
133 . zi0,sz0
134 my_real
135 . nnx1(mvsiz), nnx2(mvsiz), nnx3(mvsiz), nnx4(mvsiz),
136 . nny1(mvsiz), nny2(mvsiz), nny3(mvsiz), nny4(mvsiz),
137 . nnz1(mvsiz), nnz2(mvsiz), nnz3(mvsiz), nnz4(mvsiz)
138 my_real
139 . val1,val2,val3, v1(3), v2(3), v3(3), m(3,3,nbcut_max),vec(3), ! basis vector for each intersectionp plane
140 . posj, pos, val, !position on the segment: unnormalized dot product
141 . interp(2), wet_ratio(0:4), z(3), stria(4), dsum, numgrav, pgrav(3),
142 . cumul_val1, cumul_val2
143
144 my_real,DIMENSION(3) :: pt1,pt2,pt3,pt4
145
146 INTEGER :: N_INTP_LAG(1:4,NBCUT_MAX) ! number of intersection point on lagrandian edge id=index1 !index3:ICUT
147! TYPE(LAGRANGIAN_EDGE) :: EDG_LAG(1:4, 1:6 , NBCUT_MAX) !index1: id !index2: up to 6 intersection points (exact value is N_INTP_LAG(index1)) is Scut edge id !index3:ICUT
148
149
150 LOGICAL BOOL, ATYP1, ATYP2, ATYP3
151
152 INTEGER :: ITYP
153
154 my_real :: xp1(3), xp2(3), xp3(3), xp4(3)
155
156
157
158 INTEGER,EXTERNAL :: ISONPOLYG
159
160! INTEGER IDX(7,6)
161 INTEGER IAD1,ILAG0,ISCUT0,ILAGk,ISCUTk,IIB,NTRIA,
162 . JPOS(6) !index of nonzero positions on iadd(id_lagrangian_edge,:)
163
164 TYPE(POLYGON) :: PolySCUT, PolyTria(4), POLYresult
165
166 INTERFACE
167 FUNCTION i22aera(NPTS,P, C)
168 INTEGER :: NPTS
169 my_real :: p(3,npts), i22aera(3), c(3)
170 END FUNCTION
171 END INTERFACE
172
173 integer,external :: GetProjectedFaceType
174C-----------------------------------------------
175C C o m m e n t s
176C-----------------------------------------------
177C Aim : compute Wet surface of a lagrangian face
178C method : lagrangian segment and Cut Surface are projected on Cut surface associated plane. Polygonal clipping also enables
179C to estimate wet ratio VAL1/VAL2 where VAL1 is clipped polygon and VAL2 is full 2D segment
180C
181C +--------------+
182C / \
183C / oooooooooooooo
184C + o \ o
185C | o \ o
186C | SCUT (2D) o + o
187C | o| o
188C | | o o
189C +------------------------+ o
190C o o
191C o o
192C o o Projected 3D-Segment (2D)
193C o o
194C oooooooooooo
195C
196C
197C
198C
199C
200 ! CE_LOC(1:MVSIZ) = CAND_E(NFT+(1:MVSIZ)) (local)
201 ! CB_LOC(1:MVSIZ) = CAND_B(NFT+(1:MVSIZ)) (local)
202
203 ! 3 3
204 ! Proj_ICUT : R -----> R
205 ! P |-----> Ph = Proj_ICUT(P)
206
207 ! 3 3 1
208 ! DIST_ICUT : ( R ,R ) -----> R
209 ! P,Ph |-----> <P-Ph, P-Ph>
210
211 ! PROJECTION CRITERIA
212 ! based on shortest distances.
213 ! Each point of lagrangian face are projected to all cut planes.
214 ! If one plane has the shortest distance, then the full face will be projected to the Cut Plane
215 ! (NBCUT = 2 but could be extended later):
216 ! +---+---+---+---+---+---+---+---+
217 ! | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | -> ICUT : cut plane id
218 ! +--------+---+---+---+---+---+---+---+---+
219 ! | Node_1 | | T | | | | | | |
220 ! +--------+---+---+---+---+---+---+---+---+
221 ! | Node_2 | | T | | | | | | |
222 ! +--------+---+---+---+---+---+---+---+---+
223 ! | Node_3 | | | | | | | | T |
224 ! +--------+---+---+---+---+---+---+---+---+
225 ! | Node_4 | | | | | | | | T |
226 ! +--------+---+---+---+---+---+---+---+---+
227 ! | |
228 ! IE going to be |
229 ! projected on Pcut_2 |
230 ! |
231 ! IE also going to be
232 ! projected on Pcut_8
233
234 ! EXAMPLE
235 !
236 ! Pcut_1
237 ! |
238 ! |
239 ! +----0--------------------+ Pcut_2
240 ! | | | /
241 ! | | |/
242 ! | | 0
243 ! | | /|
244 ! | | / |
245 ! | | / |
246 ! | | / |
247 ! | | / |
248 ! | | / |
249 ! Node_1 | / |
250 ! | + | Pcut>8 / |
251 ! | \| | / |
252 ! +----0----------0---------+
253 ! \ /
254 ! IE \ /
255 ! \
256 ! \
257 ! +
258 ! Node_2
259C--------------------------------------------------------
260C S o u r c e L i n e s
261C--------------------------------------------------------
262 nin = 1
263 ce_loc(1:jlt) = cand_e(1:jlt)
264 cb_loc(1:jlt) = cand_b(1:jlt)
265
266C--------------------------------------------------------
267C REFERENCE POINT : QUADRAGLE CENTROID OR TRIA THIRD NODE
268C--------------------------------------------------------
269
270 DO i=1,jlt
271 IF(ix3(i)/=ix4(i))THEN
272 np_rect(i) = 4
273 x0(i) = fourth*(x1(i)+x2(i)+x3(i)+x4(i))
274 y0(i) = fourth*(y1(i)+y2(i)+y3(i)+y4(i))
275 z0(i) = fourth*(z1(i)+z2(i)+z3(i)+z4(i))
276 ELSE
277 np_rect(i) = 3
278 x0(i) = x3(i)
279 y0(i) = y3(i)
280 z0(i) = z3(i)
281 ENDIF
282 ENDDO
283C--------------------------------------------------------
284C calculation of normals on lagrangian facets
285C--------------------------------------------------------
286C
287 DO i=1,jlt
288
289 IF(np_rect(i)==4)THEN
290
291 sel(i) = zero
292
293 x01(i) = x1(i) - x0(i)
294 y01(i) = y1(i) - y0(i)
295 z01(i) = z1(i) - z0(i)
296
297 x02(i) = x2(i) - x0(i)
298 y02(i) = y2(i) - y0(i)
299 z02(i) = z2(i) - z0(i)
300
301 x03(i) = x3(i) - x0(i)
302 y03(i) = y3(i) - y0(i)
303 z03(i) = z3(i) - z0(i)
304
305 x04(i) = x4(i) - x0(i)
306 y04(i) = y4(i) - y0(i)
307 z04(i) = z4(i) - z0(i)
308
309 sx0 = y01(i)*z02(i) - z01(i)*y02(i)
310 sy0 = z01(i)*x02(i) - x01(i)*z02(i)
311 sz0 = x01(i)*y02(i) - y01(i)*x02(i)
312 ss2(1) = sqrt(max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)) !2S_tria1
313 sel(i) = sel(i) + ss2(1)
314 ss2(1) = one/ss2(1)
315
316 nnx1(i) = sx0 * ss2(1)
317 nny1(i) = sy0 * ss2(1)
318 nnz1(i) = sz0 * ss2(1)
319
320 sx0 = y02(i)*z03(i) - z02(i)*y03(i)
321 sy0 = z02(i)*x03(i) - x02(i)*z03(i)
322 sz0 = x02(i)*y03(i) - y02(i)*x03(i)
323 ss2(2) = sqrt(max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)) !2S_tria2
324 sel(i) = sel(i) + ss2(2)
325 ss2(2) = one/ss2(2)
326
327 nnx2(i) = sx0 * ss2(2)
328 nny2(i) = sy0 * ss2(2)
329 nnz2(i) = sz0 * ss2(2)
330
331 sx0 = y03(i)*z04(i) - z03(i)*y04(i)
332 sy0 = z03(i)*x04(i) - x03(i)*z04(i)
333 sz0 = x03(i)*y04(i) - y03(i)*x04(i)
334 ss2(3) = sqrt(max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)) !2S_tria3
335 sel(i) = sel(i) + ss2(3)
336 ss2(3) = one/ss2(3)
337
338 nnx3(i) = sx0 * ss2(3)
339 nny3(i) = sy0 * ss2(3)
340 nnz3(i) = sz0 * ss2(3)
341
342 sx0 = y04(i)*z01(i) - z04(i)*y01(i)
343 sy0 = z04(i)*x01(i) - x04(i)*z01(i)
344 sz0 = x04(i)*y01(i) - y04(i)*x01(i)
345 ss2(4) = sqrt(max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)) !2S_tria4
346 sel(i) = sel(i) + ss2(4)
347 ss2(4) = one/ss2(4)
348
349 nnx4(i) = sx0 * ss2(4)
350 nny4(i) = sy0 * ss2(4)
351 nnz4(i) = sz0 * ss2(4)
352
353 sel(i) = half * sel(i) !(2*1/4 = 1/2)
354
355 ELSE ! => NP_RECT(I)==3
356
357 x01(i) = x1(i) - x0(i)
358 y01(i) = y1(i) - y0(i)
359 z01(i) = z1(i) - z0(i)
360
361 x02(i) = x2(i) - x0(i)
362 y02(i) = y2(i) - y0(i)
363 z02(i) = z2(i) - z0(i)
364
365 sx0 = y01(i)*z02(i) - z01(i)*y02(i)
366 sy0 = z01(i)*x02(i) - x01(i)*z02(i)
367 sz0 = x01(i)*y02(i) - y01(i)*x02(i)
368 ss2(1) = sqrt(max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0))
369 sel(i) = ss2(1)
370 ss2(1) = one/ss2(1)
371
372 nnx1(i) = sx0 * ss2(1)
373 nny1(i) = sy0 * ss2(1)
374 nnz1(i) = sz0 * ss2(1)
375
376 sel(i) = half*sel(i)
377
378 ENDIF
379
380 ENDDO
381
382
383
384C--------------------------------------------------------
385 !NORMAL BACKUP
386 ! for a quadrangle, normal vector is
387 ! averaged from its 4 trias.
388C--------------------------------------------------------
389#include "lockon.inc"
390 DO i=1,jlt
391 IF(np_rect(i)==4)THEN
392 surf(1,iabs(cand_e(i))) = fourth*(nnx1(i) + nnx2(i) + nnx3(i) + nnx4(i))
393 surf(2,iabs(cand_e(i))) = fourth*(nny1(i) + nny2(i) + nny3(i) + nny4(i))
394 surf(3,iabs(cand_e(i))) = fourth*(nnz1(i) + nnz2(i) + nnz3(i) + nnz4(i))
395 ELSE
396 surf(1,iabs(cand_e(i))) = nnx1(i)
397 surf(2,iabs(cand_e(i))) = nny1(i)
398 surf(3,iabs(cand_e(i))) = nnz1(i)
399 ENDIF
400 ENDDO
401#include "lockoff.inc"
402
403C--------------------------------------------------------
404 !SUBROUTINE OUTPUT
405C--------------------------------------------------------
406c DO I=1,JLT
407c IF(IX3(I)/=IX4(I))THEN
408c ! Average from the 4 triangles
409c NX1(I) = FOURTH*(NNX1(I)+NNX2(I)+NNX3(I)+NNX4(I))
410c NY1(I) = FOURTH*(NNY1(I)+NNY2(I)+NNY3(I)+NNY4(I))
411c NZ1(I) = FOURTH*(NNZ1(I)+NNZ2(I)+NNZ3(I)+NNZ4(I))
412c ELSE
413c NX1(I) = (NNX1(I)+NNX2(I)+NNX3(I)+NNX4(I))
414c NY1(I) = (NNY1(I)+NNY2(I)+NNY3(I)+NNY4(I))
415c NZ1(I) = (NNZ1(I)+NNZ2(I)+NNZ3(I)+NNZ4(I))
416c ENDIF
417c ENDDO
418
419! if(IBUG22_Swet==-11 .AND. INT22>0 )then
420! print *, " |----------i22wetsurf.F---------|"
421! print *, " | LAGRANGIAN PROJECTIONS |"
422! print *, " |-------------------------------|"
423! write (*,FMT='(A,I10)') " N_CAND_B =",NCANDB
424! write (*,FMT='(A,I10)') " N_CAND_E =",NCANDE
425! write (*,FMT='(A,I10)') " #COUPLES =",I_STOK
426! write (*,FMT='(A,I10)') " NB =",NB
427! Do I = 1, JLT!By curling on NB you must use global tables (less efficient, but it is for post/debug)
428! IB = CB_LOC(I)
429! IFT = BRICK_LIST(NIN,IB)%Seg_add_LFT
430! ILT = BRICK_LIST(NIN,IB)%Seg_add_LLT
431! ILT = MIN(ILT, JLT) !boucle MVSIZ
432!if(ilt<ift)cycle !looping here on bricks not on candidates, should not happen
433!
434! !IFT : [1,NBRIC] |-> 1st :si =I_stok alors pas candidate
435! !ILT : [1,NBRIC] |-> last :si =0 alors idem
436!
437! write (*,FMT='(A,I10)') " IB =",IB
438! write (*,FMT='(A,I10,I10)') " IFT,ILT =",IFT,ILT
439! write (*,FMT='(A,I5,A,100I10)')" CE_LOC(",I,")=", CE_LOC(I)
440! write (*,FMT='(A,100I10)') " N1=",ITAB(IRECT(1,CE_LOC(I)))
441! write (*,FMT='(A,100I10)') " N2=",ITAB(IRECT(2,CE_LOC(I)))
442! write (*,FMT='(A,100I10)') " N3=",ITAB(IRECT(3,CE_LOC(I)))
443! write (*,FMT='(A,100I10)') " N4=",ITAB(IRECT(4,CE_LOC(I)))
444! print *, ""
445! enddo
446! endif
447
448 dist_pcut(1:nbcut_max , 1:4 , 1:ncande) = ep30
449
450
451! +------------------------------------+
452! + +
453! + +
454! +\ +
455! + \ +
456! + \ +
457! + \ +
458! + \ + 8 real cuts from polyedra
459! + \ ICUT=1 + +6 fictitious cuts from closed surface hypothesis
460! + \ +
461! + \ +
462! + \ +
463! + \ +
464! + \ ICUT=2 /+
465! + \ / +
466! + \ / +
467! + \ ICUT=8+1 / +
468! +--------------+------------+--------+
469! ^
470! |_
471! brick_list()%Pol9woNodes(j) = 1 (cell 9 does not have nodes on this face = closed surface hypothesis)
472! These closed surfaces must be taken into account for contact forces (and also internal forces)
473
474
475 DO i=1,jlt
476C--------------------------------------------------------
477 !projection of nodes onto the cutting plane
478C--------------------------------------------------------
479 ib = cb_loc(i)
480 ift = brick_list(nin,ib)%Seg_add_LFT
481 ilt = brick_list(nin,ib)%Seg_add_LLT
482 nbcut = brick_list(nin,ib)%NBCUT
483 ie = i
484
485 !face coordinates
486 p(1:3,1) = (/ x1(ie), y1(ie), z1(ie) /)
487 p(1:3,2) = (/ x2(ie), y2(ie), z2(ie) /)
488 p(1:3,3) = (/ x3(ie), y3(ie), z3(ie) /)
489 p(1:3,4) = (/ x4(ie), y4(ie), z4(ie) /)
490
491 DO icut = 1,nbcut
492 delta(1:4,icut,i) = zero
493 normal(1:3) = brick_list(nin,ib)%PCUT(icut)%N(1:3)
494 basis(1:3) = brick_list(nin,ib)%PCUT(icut)%B(1:3)
495 n_scut(1:3,icut,i) = normal(1:3)
496 !(/a, b, c, d/) = (/normal (1: 3), -normal (1: 3)*Basis (1: 3)/)! Equation of the section plan (PCUT)
497 a = normal(1)
498 b = normal(2)
499 c = normal(3)
500 d = -a*basis(1) -b*basis(2) -c*basis(3)
501 a2 = a*a
502 b2 = b*b
503 c2 = c*c
504 l2 = one / max(em30,a2+b2+c2)
505 !projected face coordinates
506 DO j=1,4
507 ph( 1,icut,j) = (b2+c2)*p(1,j) - a*b*p(2,j) - a*c*p(3,j) - d*a
508 ph( 2,icut,j) = -a*b*p(1,j) + (a2+c2)*p(2,j) - b*c*p(3,j) - d*b
509 ph( 3,icut,j) = -a*c*p(1,j) - b*c*p(2,j) + (a2+b2)*p(3,j) - d*c
510 ph(1:3,icut,j) = l2 * ph(1:3,icut,j)
511 pph(1:3,j) = p(1:3,j) - ph(1:3,icut,j)
512 dist(j) = sum( pph(1:3,j) * pph(1:3,j) )
513 dist_pcut(icut,j,i) = dist(j) !L2 Norm ** 2
514 ENDDO
515 !projected cut surface nodes (coplanarity is assumption but not reality)
516 np = brick_list(nin,ib)%PCUT(icut)%NP
517 DO j=1,np
518 s(1:3,j) = brick_list(nin,ib)%PCUT(icut)%P(1:3,j)
519 ENDDO
520 DO j=1,np
521 sh( 1,icut,j) = (b2+c2)*s(1,j) - a*b*s(2,j) - a*c*s(3,j) - d*a
522 sh( 2,icut,j) = -a*b*s(1,j) + (a2+c2)*s(2,j) - b*c*s(3,j) - d*b
523 sh( 3,icut,j) = -a*c*s(1,j) - b*c*s(2,j) + (a2+b2)*s(3,j) - d*c
524 sh(1:3,icut,j) = l2 * sh(1:3,icut,j)
525 ENDDO
526 ENDDO !next ICUT
527
528 !Additional cut planes for closed surface hypothesis
529 nbcut_add = 0
530 DO icut = 8+1,8+6
531 j = icut - 8
532 np = brick_list(nin,ib)%PCUT(icut)%NP
533 IF(np<=0)cycle
534 nbcut_add = nbcut_add + 1
535 delta(1:4,icut,i) = zero
536 normal(1:3) = brick_list(nin,ib)%PCUT(icut)%N(1:3)
537 basis(1:3) = brick_list(nin,ib)%PCUT(icut)%B(1:3)
538 n_scut(1:3,icut,i) = normal(1:3)
539 !(/a, b, c, d/) = (/normal (1: 3), -normal (1: 3)*Basis (1: 3)/)! Equation of the section plan (PCUT)
540 a = normal(1)
541 b = normal(2)
542 c = normal(3)
543 d = -a*basis(1) -b*basis(2) -c*basis(3)
544 a2 = a*a
545 b2 = b*b
546 c2 = c*c
547 l2 = one / max(em30,a2+b2+c2)
548 !projected face coordinates
549 DO j=1,4
550 ph( 1,icut,j) = (b2+c2)*p(1,j) - a*b*p(2,j) - a*c*p(3,j) - d*a
551 ph( 2,icut,j) = -a*b*p(1,j) + (a2+c2)*p(2,j) - b*c*p(3,j) - d*b
552 ph( 3,icut,j) = -a*c*p(1,j) - b*c*p(2,j) + (a2+b2)*p(3,j) - d*c
553 ph(1:3,icut,j) = l2 * ph(1:3,icut,j)
554 pph(1:3,j) = p(1:3,j) - ph(1:3,icut,j)
555 dist(j) = sum( pph(1:3,j) * pph(1:3,j) )
556 dist_pcut(icut,j,i) = dist(j) !L2 Norm ** 2
557 ENDDO
558 !projected cut surface nodes (coplanarity is assumption but not reality)
559 np = brick_list(nin,ib)%PCUT(icut)%NP
560 DO j=1,np
561 s(1:3,j) = brick_list(nin,ib)%PCUT(icut)%P(1:3,j)
562 ENDDO
563 DO j=1,np
564 sh( 1,icut,j) = (b2+c2)*s(1,j) - a*b*s(2,j) - a*c*s(3,j) - d*a
565 sh( 2,icut,j) = -a*b*s(1,j) + (a2+c2)*s(2,j) - b*c*s(3,j) - d*b
566 sh( 3,icut,j) = -a*c*s(1,j) - b*c*s(2,j) + (a2+b2)*s(3,j) - d*c
567 sh(1:3,icut,j) = l2 * sh(1:3,icut,j)
568 ENDDO
569 ENDDO !next ICUT
570
571
572C TAG_NOD(1:4,1:NBCUT) = 0
573C DO J=1,4
574C !init before entering loop
575C ICUT = 1
576C VAL = DIST_PCUT(ICUT,J,I)
577C !loop on cut planes in current cut cell IB = CB_LOC(I)
578C DO ICUT=1,NBCUT
579C IF(DIST_PCUT(ICUT,J,I) <= VAL)THEN
580C TAG_NOD(J, ICUT) = 1
581C VAL = DIST_PCUT(ICUT,J,I)
582C ENDIF
583C ENDDO!next ICUT
584C ENDDO !next J
585C !tag for additional cuts (closed surface hypothesis)
586C IF(NBCUT_ADD > 0)THEN
587C DO J=1,4
588C !init before entering loop
589C ICUT = 8+1
590C VAL = DIST_PCUT(ICUT,J,I)
591C !loop on cut planes in current cut cell IB = CB_LOC(I)
592C DO ICUT=8+1,8+NBCUT_ADD
593C IF(DIST_PCUT(ICUT,J,I) <= VAL)THEN
594C TAG_NOD(J, ICUT)= 1
595C VAL = DIST_PCUT(ICUT,J,I)
596C ENDIF
597C ENDDO!next ICUT
598C ENDDO !next J
599C ENDIF!(NBCUT_ADD > 0)
600
601 !tag of lagrangian elements
602C TAG_EL(1:NBCUT_MAX) = 0
603C IF(NBCUT>0)THEN
604C DO ICUT=1,NBCUT
605C TAG_EL(ICUT) = SUM(TAG_NOD(:,ICUT))
606C if(IBUG22_Swet==-1 .AND. INT22>0 )write (*,FMT='(A,I10,A)') " |---COUPLE:",I,"-----------|"
607C if(IBUG22_Swet==-1 .AND. INT22>0 .AND. NBCUT<=1)print *, " TAG_EL(ICUT=1)=" ,TAG_EL(1)
608C if(IBUG22_Swet==-1 .AND. INT22>0 .AND. NBCUT<=2)print *, " TAG_EL(ICUT=2)=" ,TAG_EL(2)
609C if(IBUG22_Swet==-1 .AND. INT22>0 .AND. NBCUT<=3)print *, " TAG_EL(ICUT=3)=" ,TAG_EL(3)
610C if(IBUG22_Swet==-1 .AND. INT22>0 .AND. NBCUT<=4)print *, " TAG_EL(ICUT=4)=" ,TAG_EL(4)
611C if(IBUG22_Swet==-1 .AND. INT22>0 .AND. NBCUT<=5)print *, " TAG_EL(ICUT=5)=" ,TAG_EL(5)
612C if(IBUG22_Swet==-1 .AND. INT22>0 .AND. NBCUT<=6)print *, " TAG_EL(ICUT=6)=" ,TAG_EL(6)
613C if(IBUG22_Swet==-1 .AND. INT22>0 .AND. NBCUT<=7)print *, " TAG_EL(ICUT=7)=" ,TAG_EL(7)
614C if(IBUG22_Swet==-1 .AND. INT22>0 .AND. NBCUT<=8)print *, " TAG_EL(ICUT=8)=" ,TAG_EL(8)
615C ENDDO!next ICUT
616C ENDIF
617
618
619
620C--------------------------------------------------------
621 !injection into the proper space of the projection plane (i.e. intersection plane)
622C--------------------------------------------------------
623
624 !DIM = 2
625 !Base = Projection matrix Eingenvectors (eigenval = 1)
626 ! B1=(-b,a,0) B2=(-c,0,a)
627
628 DO icut=1,nbcut
629C if(TAG_EL(ICUT)==0)CYCLE
630 ib = cb_loc(i)
631 normal(1:3) = brick_list(nin,ib)%PCUT(icut)%N(1:3)
632 a = normal(1)
633 b = normal(2)
634 c = normal(3)
635 IF (a==zero . and. b==zero)THEN
636 v1 = (/ one,zero,zero /)
637 v2 = (/ zero,one,zero /)
638 ELSEIF(a==zero . and. c==zero)THEN
639 v1 = (/ one,zero,zero /)
640 v2 = (/ zero,zero,one /)
641 ELSEIF(b==zero . and. c==zero)THEN
642 v1 = (/ zero,one,zero /)
643 v2 = (/ zero,zero,one /)
644 ELSE
645 IF(a/=zero)THEN
646 v1 = (/ -b, a , zero /)
647 v2 = (/ -c, zero, a /)
648 ELSE
649 v1 = (/ one, zero, zero /)
650 v2 = (/ zero, -c, b /)
651 ENDIF
652 !ORTHONORMAL
653 !hyopothesis : V1, V2 are non colinear, ||V1||=1
654 !V2~ = labmda1 * V1 + lambda2 * V2 (V2~ is generated by (V1,V2) )
655 !let set lambda1 = <V1,V2>, then lambda2=-1
656 v1 = v1 / sqrt(sum(v1*v1))
657 ps = v1(1)*v2(1) + v1(2)*v2(2) +v1(3)*v2(3)
658 v2 = ps*v1 -v2
659 v2 = v2 / sqrt(sum(v2*v2))
660 ENDIF
661 v3 = normal / sqrt(sum(normal*normal))
662 !BASIS SWITCH
663 !third direction is Pcut normal
664 !change of basis matrix
665 m(1:3,1,icut) = v1
666 m(1:3,2,icut) = v2
667 m(1:3,3,icut) = v3
668 if(ibug22_swet==-1 .AND. int22>0 )then
669 print *," ORTHONORMAL BASIS :"
670 write (*,fmt='(a,f20.10,a1,f20.10,a1,f20.10)')" V1=", V1(1),",",V1(2),",",V1(3)
671 write (*,FMT='(a,f20.10,a1,f20.10,a1,f20.10)')" V2=", V2(1),",",V2(2),",",V2(3)
672 write (*,FMT='(a,f20.10,a1,f20.10,a1,f20.10)')" V3=", V3(1),",",V3(2),",",V3(3)
673 endif
674 !new coordinates in the projection plane (dim=2)
675 !no need to invert the matrix because it is orthogonal
676 DO J=1,4
677 VAL1 = M(1,1,ICUT)*Ph(1,ICUT,J) + M(2,1,ICUT)*Ph(2,ICUT,J) + M(3,1,ICUT)*Ph(3,ICUT,J)
678 VAL2 = M(1,2,ICUT)*Ph(1,ICUT,J) + M(2,2,ICUT)*Ph(2,ICUT,J) + M(3,2,ICUT)*Ph(3,ICUT,J)
679 VAL3 = ZERO
680 Ph(1,ICUT,J) = VAL1
681 Ph(2,ICUT,J) = VAL2
682 Ph(3,ICUT,J) = VAL3
683 ENDDO
684 NP = BRICK_LIST(NIN,IB)%PCUT(ICUT)%NP
685 DO J=1,NP
686 VAL1 = M(1,1,ICUT)*Sh(1,ICUT,J) + M(2,1,ICUT)*Sh(2,ICUT,J) + M(3,1,ICUT)*Sh(3,ICUT,J)
687 VAL2 = M(1,2,ICUT)*Sh(1,ICUT,J) + M(2,2,ICUT)*Sh(2,ICUT,J) + M(3,2,ICUT)*Sh(3,ICUT,J)
688 VAL3 = ZERO
689 SH(1,ICUT,J) = VAL1
690 SH(2,ICUT,J) = VAL2
691 Sh(3,ICUT,J) = VAL3
692 ENDDO!next NP
693 ENDDO!next ICUT
694
695
696 IF(NBCUT_ADD > 0)THEN
697 DO ICUT=8+1,8+6
698 J = ICUT - 8
699 NP = BRICK_LIST(NIN,IB)%PCUT(ICUT)%NP
700 IF(NP<=0)CYCLE
701C if(TAG_EL(ICUT)==0)CYCLE
702 IB = CB_LOC(I)
703 NORMAL(1:3) = BRICK_LIST(NIN,IB)%PCUT(ICUT)%N(1:3)
704 a = NORMAL(1)
705 b = NORMAL(2)
706 c = NORMAL(3)
707 IF (a==ZERO . AND. b==ZERO)THEN
708 V1 = (/ ONE,ZERO,ZERO /)
709 V2 = (/ ZERO,ONE,ZERO /)
710 ELSEIF(a==ZERO . AND. c==ZERO)THEN
711 V1 = (/ ONE,ZERO,ZERO /)
712 V2 = (/ ZERO,ZERO,ONE /)
713 ELSEIF(b==ZERO . AND. c==ZERO)THEN
714 V1 = (/ ZERO,ONE,ZERO /)
715 V2 = (/ ZERO,ZERO,ONE /)
716 ELSE
717 IF(a/=ZERO)THEN
718 V1 = (/ -b, a , ZERO /)
719 V2 = (/ -c, ZERO, a /)
720 ELSE
721 V1 = (/ ONE, ZERO, ZERO /)
722 V2 = (/ ZERO, -c, b /)
723 ENDIF
724 !ORTHONORMAL
725 !hyopothesis : V1, V2 are non colinear, ||V1||=1
726 !V2~ = labmda1 * V1 + lambda2 * V2 (V2~ is generated by (V1,V2) )
727 !let set lambda1 = <V1,V2>, then lambda2=-1
728 V1 = V1 / SQRT(SUM(V1*V1))
729 PS = V1(1)*V2(1) + V1(2)*V2(2) +V1(3)*V2(3)
730 V2 = PS*V1 -V2
731 V2 = V2 / SQRT(SUM(V2*V2))
732 ENDIF
733 V3 = NORMAL / SQRT(SUM(NORMAL*NORMAL))
734 !BASIS SWITCH
735 !third direction is Pcut normal
736 !change of basis matrix
737 M(1:3,1,ICUT) = V1
738 M(1:3,2,ICUT) = V2
739 M(1:3,3,ICUT) = V3
740.AND. if(IBUG22_Swet==-1 INT22>0 )then
741 print *," ORTHONORMAL BASIS _additional Scut (closed surf hypothesis) :"
742 write (*,FMT='(a,f20.10,a1,f20.10,a1,f20.10)')" V1=", V1(1),",",V1(2),",",V1(3)
743 write (*,FMT='(a,f20.10,a1,f20.10,a1,f20.10)')" V2=", V2(1),",",V2(2),",",V2(3)
744 write (*,FMT='(a,f20.10,a1,f20.10,a1,f20.10)')" V3=", V3(1),",",V3(2),",",V3(3)
745 endif
746 !new coordinates in the projection plane (dim=2)
747 !no need to invert the matrix because it is orthogonal
748 DO J=1,4
749 VAL1 = M(1,1,ICUT)*Ph(1,ICUT,J) + M(2,1,ICUT)*Ph(2,ICUT,J) + M(3,1,ICUT)*Ph(3,ICUT,J)
750 VAL2 = M(1,2,ICUT)*Ph(1,ICUT,J) + M(2,2,ICUT)*Ph(2,ICUT,J) + M(3,2,ICUT)*Ph(3,ICUT,J)
751 VAL3 = ZERO
752 Ph(1:3,ICUT,J) = (/VAL1,VAL2,VAL3/)
753 ENDDO
754 NP = BRICK_LIST(NIN,IB)%PCUT(ICUT)%NP
755 DO J=1,NP
756 VAL1 = M(1,1,ICUT)*Sh(1,ICUT,J) + M(2,1,ICUT)*Sh(2,ICUT,J) + M(3,1,ICUT)*Sh(3,ICUT,J)
757 VAL2 = M(1,2,ICUT)*Sh(1,ICUT,J) + M(2,2,ICUT)*Sh(2,ICUT,J) + M(3,2,ICUT)*Sh(3,ICUT,J)
758 VAL3 = ZERO
759 Sh(1:3,ICUT,J) = (/VAL1,VAL2,VAL3/)
760 ENDDO!next NP
761 ENDDO!next ICUT
762 ENDIF !(NBCUT_ADD > 0)
763
764
765
766C--------------------------------------------------------
767 !POLYGONAL CLIPPING (SCUT with projected segment)
768C--------------------------------------------------------
769! Projected Segment might be non convex. In this case it is decomposed into 2 triangles (hourglass shape)
770! If projected segment is convex then it is splitted into 4 triangles (required for parith on)
771!
772 DO ICUT = 1, NBCUT
773 SEFF(ICUT,I) = ZERO
774 P_Grav(1:3,ICUT,I) = ZERO
775 NumGrav = ZERO
776 NP = BRICK_LIST(NIN,IB)%PCUT(ICUT)%NP !BUG line was missing !
777C if(TAG_EL(ICUT)==0)CYCLE
778 PolySCUT%NumNodes = NP+1 !
779 DO J=1,NP
780 PolySCUT%node(J,1:2) = Sh(1:2,ICUT,J)
781 ENDDO
782 PolySCUT%node(NP+1,:) = PolySCUT%node(1,:)
783 CALL SetCounterClockWisePolyg( PolySCUT, VAL2 )
784 IF(VAL2==ZERO)CYCLE
785 ITYP = GetProjectedFaceType(
786 . Ph(1:2,ICUT,1), Ph(1:2,ICUT,2), Ph(1:2,ICUT,3), Ph(1:2,ICUT,4), InterP, ListP )
787 NTRIA = 4
788 IF(ITYP/=0)NTRIA = 2
789 IF(NP_RECT(I)==3) THEN
790 ITYP = 0
791 NTRIA = 1
792 ENDIF
793
794 Stria(1:4) = ZERO
795
796 !first triangle
797 PolyTria(1)%NumNodes = 3+1
798 PolyTria(1)%node(1,1:2) = Ph(1:2,ICUT,ListP(1))
799 PolyTria(1)%node(2,1:2) = Ph(1:2,ICUT,ListP(2))
800 PolyTria(1)%node(3,1:2) = InterP(1:2)
801 PolyTria(1)%node(4,:) = PolyTria(1)%node(1,:)
802 CALL SetCounterClockWisePolyg( PolyTria(1), Stria(1) )
803
804 !second triangle
805 PolyTria(2)%NumNodes = 3+1
806 PolyTria(2)%node(1,1:2) = Ph(1:2,ICUT,ListP(3))
807 PolyTria(2)%node(2,1:2) = Ph(1:2,ICUT,ListP(4))
808 PolyTria(2)%node(3,1:2) = InterP(1:2)
809 PolyTria(2)%node(4,:) = PolyTria(2)%node(1,:)
810 CALL SetCounterClockWisePolyg( PolyTria(2), Stria(2) )
811
812 IF(NTRIA == 4)THEN
813 !third triangle
814 PolyTria(3)%NumNodes = 3+1
815 PolyTria(3)%node(1,1:2) = Ph(1:2,ICUT,ListP(2))
816 PolyTria(3)%node(2,1:2) = Ph(1:2,ICUT,ListP(3))
817 PolyTria(3)%node(3,1:2) = InterP(1:2)
818 PolyTria(3)%node(4,:) = PolyTria(3)%node(1,:)
819 CALL SetCounterClockWisePolyg( PolyTria(3), Stria(3) )
820
821 !fourth triangle
822 PolyTria(4)%NumNodes = 3+1
823 PolyTria(4)%node(1,1:2) = Ph(1:2,ICUT,ListP(1))
824 PolyTria(4)%node(2,1:2) = Ph(1:2,ICUT,ListP(4))
825 PolyTria(4)%node(3,1:2) = InterP(1:2)
826 PolyTria(4)%node(4,:) = PolyTria(4)%node(1,:)
827 CALL SetCounterClockWisePolyg( PolyTria(4), Stria(4) )
828 ENDIF
829
830 Cumul_VAL2 = SUM(Stria(1:4))
831
832 !-------OUTPUT------
833.AND. if(IBUG22_Swet==-1 INT22>0 )then
834 print *, " |-----------i22wetsurf.F--------|"
835 print *, " | LAGRANGIAN PROJECTIONS |"
836 print *, " |-------------------------------|"
837 IB = CB_LOC(I)
838 brickID = BRICK_LISt(NIN,IB)%ID
839 write (*,FMT='(a,1i10,a,i2)') " --------brickID=",ixs(11,brickID), " ICUT=", ICUT
840 K = ABS(CE_LOC(I))
841 write (*,FMT='(a,1i10)') " N1=",ITAB(IRECT(1,K))
842 write (*,FMT='(a,1i10)') " N2=",ITAB(IRECT(2,K))
843 write (*,FMT='(a,1i10)') " N3=",ITAB(IRECT(3,K))
844 write (*,FMT='(a,1i10)') " N4=",ITAB(IRECT(4,K))
845 write (*,FMT='(a)' ) " Projected Tria"
846 write (*,FMT='(a,3f30.16)') " *createnode ",Ph(1:2,ICUT,1),ZERO
847 write (*,FMT='(a,3f30.16)') " *createnode ",Ph(1:2,ICUT,2),ZERO
848 write (*,FMT='(a,3f30.16)') " *createnode ",Ph(1:2,ICUT,3),ZERO
849 write (*,FMT='(a,3f30.16)') " *createnode ",Ph(1:2,ICUT,4),ZERO
850 write (*,FMT='(a,1i1)' ) " ityp=",ITYP
851 write (*,FMT='(a)' ) " PolySCUT"
852 write (*,FMT='(a,3f30.16)') " *createnode ",PolySCUT%node(1,1:2),ZERO
853 write (*,FMT='(a,3f30.16)') " *createnode ",PolySCUT%node(2,1:2),ZERO
854 write (*,FMT='(a,3f30.16)') " *createnode ",PolySCUT%node(3,1:2),ZERO
855 write (*,FMT='(a,3f30.16)') " *createnode ",PolySCUT%node(4,1:2),ZERO
856 write (*,FMT='(a,1i10)') " NTRIA=",NTRIA
857 DO J=1,NTRIA
858
859 ENDDO
860 print *, ""
861 endif
862
863 POLYresult%NumNodes = 0
864
865 WET_RATIO(:) = ZERO
866 VAL1 = ZERO
867 VAL2 = ZERO
868 Cumul_VAL1 = ZERO
869
870 DO J=1,NTRIA
871 CALL PolygonalClipping( PolyTria(J),PolySCUT, POLYresult)
872 IF(POLYresult%NumNodes==0)CYCLE
873 !full area
874 pt1(1:2) = PolyTria(J)%node(1,1:2) ; pt1(3)=ZERO
875 pt2(1:2) = PolyTria(J)%node(2,1:2) ; pt2(3)=ZERO
876 pt3(1:2) = PolyTria(J)%node(3,1:2) ; pt3(3)=ZERO
877 !wet area
878.AND. if(IBUG22_Swet==-1 INT22>0 )then
879 write (*,FMT='(a,1i10)') " --TRIA #", J
880 write (*,FMT='(a,3f30.16)') " *createnode ",PolyTria(J)%node(1,1:2),ZERO
881 write (*,FMT='(a,3f30.16)') " *createnode ",PolyTria(J)%node(2,1:2),ZERO
882 write (*,FMT='(a,3f30.16)') " *createnode ",PolyTria(J)%node(3,1:2),ZERO
883 write (*,FMT='(a,1i10)') " ->clip nodes", POLYresult%NumNodes
884 do k=1,POLYresult%NumNodes-1
885 write (*,FMT='(a,3f45.35)') " *createnode ",POLYresult%node(k,1:2),ZERO
886 enddo
887 endif
888
889 PGrav(1:3) = ZERO
890 !Gravity center estimation / To be updated later / must also be Parith ON
891 IF(POLYresult%NumNodes-1>0)THEN
892 DO k=1,POLYresult%NumNodes-1
893c P_Grav(1,ICUT,I) = P_Grav(1,ICUT,I) + POLYresult%node(k,1)
894c P_Grav(2,ICUT,I) = P_Grav(2,ICUT,I) + POLYresult%node(k,2)
895 PGrav(1) = PGrav(1) + POLYresult%node(k,1)
896 PGrav(2) = PGrav(2) + POLYresult%node(k,2)
897 ENDDO
898 PGrav(1) = PGrav(1) / (POLYresult%NumNodes-1)
899 PGrav(2) = PGrav(2) / (POLYresult%NumNodes-1)
900 ENDIF
901
902 P_Grav(1,ICUT,I) = P_Grav(1,ICUT,I) + PGrav(1)
903 P_Grav(2,ICUT,I) = P_Grav(2,ICUT,I) + PGrav(2)
904 NumGrav = Numgrav + ONE
905
906.AND. IF(POLYresult%NumNodes > 0 STria(J)/=ZERO)THEN
907 pt1(1:2) = PolyTria(J)%node(1,1:2) ; pt1(3)=ZERO
908 pt2(1:2) = PolyTria(J)%node(2,1:2) ; pt2(3)=ZERO
909 pt3(1:2) = PolyTria(J)%node(3,1:2) ; pt3(3)=ZERO
910 CALL SetCounterClockWisePolyg(POLYresult ,VAL)
911 VAL = MAX(ZERO,MIN(VAL,Stria(J)))
912 Cumul_VAL1 = Cumul_VAL1 + VAL
913.AND. if(IBUG22_Swet==-1 INT22>0 )write (*,FMT='(a,1f30.16)') " -> Tria Wet Surf=", VAL
914.AND. if(IBUG22_Swet==-1 INT22>0 )write (*,FMT='(a,1f30.16)') " -> Tria Surf =", Stria(J)
915.AND. if(IBUG22_Swet==-1 INT22>0 )write (*,FMT='(a,1f30.16)') " -> cumulated Wet Surf=", Cumul_VAL1
916.AND. if(IBUG22_Swet==-1 INT22>0 )write (*,FMT='(a,1f30.16)') " -> Sum of tria Surf=", Cumul_VAL2
917
918 !post-condition testing
919 WET_RATIO(J) = VAL/Stria(J)
920.OR. IF(WET_RATIO(J)<-EM06 WET_RATIO(J)-ONE>EM06)THEN
921 print *, "***error: inter22 topology issue while computing wet surface",WET_RATIO(J),VAL1,VAL2
922 ENDIF
923 ENDIF
924 ENDDO !next J=1,NTRIA
925
926 WET_RATIO(0) = Cumul_VAL1 !sum of wet surfaces for all triangles composing the lagrangian face.
927
928.OR. IF(WET_RATIO(0) == ZERO Cumul_VAL2==ZERO)THEN
929 SEFF(ICUT,I) = ZERO
930 CYCLE
931 ELSE
932 WET_RATIO(0) = WET_RATIO(0) / Cumul_VAL2
933 ENDIF
934
935 WET_RATIO(0) = MAX(ZERO,MIN(WET_RATIO(0),ONE))
936 SEFF(ICUT,I) = Cumul_VAL1
937
938 SEFF(ICUT,I) = WET_RATIO(0) * Sel(I)
939
940 !-------OUTPUT------
941.AND. if(IBUG22_Swet==-1 INT22>0 )then
942 write (*,FMT='(a,1f30.16)') " Sum Wet Surf =",Cumul_VAL1
943 write (*,FMT='(a,1f30.16)') " Total Surf =",Cumul_VAL2
944 write (*,FMT='(a,1f30.16)') " RATIO =",WET_RATIO(0)
945 print *, ""
946 endif
947
948 IF(NumGrav==0)CYCLE
949
950 P_Grav(1,ICUT,I) = P_Grav(1,ICUT,I) / NumGrav
951 P_Grav(2,ICUT,I) = P_Grav(2,ICUT,I) / NumGrav
952
953 !---sWeighting factors using pressure load center
954 !InterPolation using gravity centers of clipped polygons
955 !NP_RECT(I) = 3 or 4
956
957 !2d projected segment
958 pt1(1:2) = Ph(1:2,ICUT,ListP(1))
959 pt2(1:2) = Ph(1:2,ICUT,ListP(2))
960 pt3(1:2) = Ph(1:2,ICUT,ListP(3))
961 pt4(1:2) = Ph(1:2,ICUT,ListP(4))
962
963C write (*,FMT='(A )') "Segment ="
964C write (*,FMT='(A,2E30.16,A)') "*createnode ", Ph(1:2,ICUT,ListP(1)) ," 0.0"
965C write (*,FMT='(A,2E30.16,A)') "*createnode ", Ph(1:2,ICUT,ListP(2)) ," 0.0"
966C write (*,FMT='(A,2E30.16,A)') "*createnode ", Ph(1:2,ICUT,ListP(3)) ," 0.0"
967C write (*,FMT='(A,2E30.16,A)') "*createnode ", Ph(1:2,ICUT,ListP(4)) ," 0.0"
968
969C write (*,FMT='(A )') "Scut ="
970C do k=1,PolyScut%NumNodes-1
971C write (*,FMT='(A,2E30.16,A)')"*createnode ", PolyScut%node(k,1:2) ," 0.0"
972C enddo
973
974
975 !distance of each point from wet polygon gravity center
976 distance(1) = SQRT( (P_grav(1,ICUT,I)-Pt1(1))**2 + (P_grav(2,ICUT,I)-Pt1(2))**2 )
977 distance(2) = SQRT( (P_grav(1,ICUT,I)-Pt2(1))**2 + (P_grav(2,ICUT,I)-Pt2(2))**2 )
978 distance(3) = SQRT( (P_grav(1,ICUT,I)-Pt3(1))**2 + (P_grav(2,ICUT,I)-Pt3(2))**2 )
979 distance(4) = SQRT( (P_grav(1,ICUT,I)-Pt4(1))**2 + (P_grav(2,ICUT,I)-Pt4(2))**2 )
980
981C write (*,FMT='(A )') "Pgrav ="
982C write (*,FMT='(A,2E30.16,A)') "*createnode ", P_grav(1:2,ICUT,I) ," 0.0"
983
984
985 !segment edges
986 edge(1) = SQRT ( (Pt1(1)-Pt2(1))**2 + (Pt1(2)-Pt2(2))**2 )
987 edge(2) = SQRT ( (Pt2(1)-Pt3(1))**2 + (Pt2(2)-Pt3(2))**2 )
988 edge(3) = SQRT ( (Pt3(1)-Pt4(1))**2 + (Pt3(2)-Pt4(2))**2 )
989 edge(4) = SQRT ( (Pt4(1)-Pt1(1))**2 + (Pt4(2)-Pt1(2))**2 )
990 IF(NP_RECT(I)==3)THEN
991 edge(3) = SQRT ( (Pt3(1)-Pt1(1))**2 + (Pt3(2)-Pt1(2))**2 )
992 edge(4) = ZERO
993 ENDIF
994
995 a = edge(1)
996 b = distance(1)
997 c = distance(2)
998 s2 = HALF*(a+b+c)
999 Area_tria(1) = SQRT( ABS(s2*(s2-a)*(s2-b)*(s2-c)) )
1000
1001 a = edge(2)
1002 b = distance(2)
1003 c = distance(3)
1004 s2 = HALF*(a+b+c)
1005 Area_tria(2) = SQRT( ABS(s2*(s2-a)*(s2-b)*(s2-c)) )
1006
1007 IF(NP_RECT(I)==4)THEN
1008 a = edge(3)
1009 b = distance(3)
1010 c = distance(4)
1011 s2 = HALF*(a+b+c)
1012 Area_tria(3) = SQRT( ABS(s2*(s2-a)*(s2-b)*(s2-c)) )
1013
1014 a = edge(4)
1015 b = distance(4)
1016 c = distance(1)
1017 s2 = HALF*(a+b+c)
1018 Area_tria(4) = SQRT( ABS(s2*(s2-a)*(s2-b)*(s2-c)) )
1019 ELSE
1020 a = edge(3)
1021 b = distance(3)
1022 c = distance(1)
1023 s2 = HALF*(a+b+c)
1024 Area_tria(3) = SQRT( ABS(s2*(s2-a)*(s2-b)*(s2-c)) )
1025
1026 Area_tria(4) = ZERO
1027 ENDIF
1028
1029 Dsum = SUM(Area_tria(1:NP_RECT(I)))
1030
1031c22 IF(Dsum==ZERO)then
1032c22 print *, "Warning: inter22 - nodal forces computation"
1033c22 endif
1034
1035 IF(NP_RECT(I)==4)THEN
1036 !delta(1,ICUT,I) = HALF*(Area_tria(2)+Area_tria(3))/Dsum
1037 !delta(2,ICUT,I) = HALF*(Area_tria(3)+Area_tria(4))/Dsum
1038 !delta(3,ICUT,I) = HALF*(Area_tria(1)+Area_tria(4))/Dsum
1039 !delta(4,ICUT,I) = HALF*(Area_tria(1)+Area_tria(2))/Dsum
1040 !uniform loading: same weight
1041 delta(1:4,ICUT,I) = FOURTH
1042 ELSE
1043 !delta(1,ICUT,I) = (Area_tria(2))/Dsum
1044 !delta(2,ICUT,I) = (Area_tria(3))/Dsum
1045 !delta(3,ICUT,I) = (Area_tria(1))/Dsum
1046 !delta(4,ICUT,I) = ZERO
1047 !uniform loading : same weight
1048 delta(1:4,ICUT,I) = THIRD
1049 delta(4 ,ICUT,I) = ZERO
1050 ENDIF
1051
1052 !write (*,FMT='(a,3e30.16)') "CoG = ", P_grav(1:2,ICUT,I)
1053 !write (*,FMT='(a,4e30.16)') "dist = ", distance(1:4)
1054
1055c write (*,FMT='(A,4E30.16,A,1E30.16)') "A_tria= ", Area_tria(1:4)," sum= ",SUM(Area_tria)
1056c write (*,FMT='(A,4E30.16)') "delta= ", delta(1:4,ICUT,I)
1057
1058
1059 ENDDO!next ICUT
1060
1061
1062
1063
1064
1065
1066 DO ICUT = 8+1,8+6
1067 J = ICUT - 8
1068 SEFF(ICUT,I) = ZERO
1069C print *, "I=", I, " ICUT= ",ICUT
1070 NP = BRICK_LIST(NIN,IB)%PCUT(ICUT)%NP
1071 IF(NP<=0)CYCLE
1072 P_Grav(1:3,ICUT,I) = ZERO
1073 NumGrav = ZERO
1074
1075C if(TAG_EL(ICUT)==0)CYCLE
1076 PolySCUT%NumNodes = NP+1
1077 DO J=1,NP
1078 PolySCUT%node(J,1:2) = Sh(1:2,ICUT,J)
1079 ENDDO
1080 PolySCUT%node(NP+1,:) = PolySCUT%node(1,:)
1081 CALL SetCounterClockWisePolyg( PolySCUT, VAL2 )
1082
1083 IF(VAL2==ZERO)CYCLE
1084
1085 ITYP = GetProjectedFaceType(Ph(1:2,ICUT,1), Ph(1:2,ICUT,2), Ph(1:2,ICUT,3), Ph(1:2,ICUT,4), InterP, ListP)
1086
1087 NTRIA = 4
1088 IF(ITYP/=0)NTRIA = 2
1089 IF(NP_RECT(I)==3) THEN
1090 ITYP = 0
1091 NTRIA = 1
1092 ENDIF
1093
1094 Stria(1:4) = ZERO
1095
1096 !first triangle
1097 PolyTria(1)%NumNodes = 3+1
1098 PolyTria(1)%node(1,1:2) = Ph(1:2,ICUT,ListP(1))
1099 PolyTria(1)%node(2,1:2) = Ph(1:2,ICUT,ListP(2))
1100 PolyTria(1)%node(3,1:2) = InterP(1:2)
1101 PolyTria(1)%node(4,:) = PolyTria(1)%node(1,:)
1102 CALL SetCounterClockWisePolyg( PolyTria(1), Stria(1) )
1103
1104 !second triangle
1105 PolyTria(2)%NumNodes = 3+1
1106 PolyTria(2)%node(1,1:2) = Ph(1:2,ICUT,ListP(3))
1107 PolyTria(2)%node(2,1:2) = Ph(1:2,ICUT,ListP(4))
1108 PolyTria(2)%node(3,1:2) = InterP(1:2)
1109 PolyTria(2)%node(4,:) = PolyTria(2)%node(1,:)
1110 CALL SetCounterClockWisePolyg( PolyTria(2), Stria(2) )
1111
1112 IF(NTRIA == 4)THEN
1113 !third triangle
1114 PolyTria(3)%NumNodes = 3+1
1115 PolyTria(3)%node(1,1:2) = Ph(1:2,ICUT,ListP(2))
1116 PolyTria(3)%node(2,1:2) = Ph(1:2,ICUT,ListP(3))
1117 PolyTria(3)%node(3,1:2) = InterP(1:2)
1118 PolyTria(3)%node(4,:) = PolyTria(3)%node(1,:)
1119 CALL SetCounterClockWisePolyg( PolyTria(3), Stria(3) )
1120
1121 !fourth triangle
1122 PolyTria(4)%NumNodes = 3+1
1123 PolyTria(4)%node(1,1:2) = Ph(1:2,ICUT,ListP(1))
1124 PolyTria(4)%node(2,1:2) = Ph(1:2,ICUT,ListP(4))
1125 PolyTria(4)%node(3,1:2) = InterP(1:2)
1126 PolyTria(4)%node(4,:) = PolyTria(4)%node(1,:)
1127 CALL SetCounterClockWisePolyg( PolyTria(4), Stria(4) )
1128 ENDIF
1129
1130 Cumul_VAL2 = SUM(Stria(1:4))
1131
1132 !-------OUTPUT------
1133.AND. if(IBUG22_Swet==-1 INT22>0 )then
1134 print *, " |-----------i22wetsurf.F--------|"
1135 print *, " | LAGRANGIAN PROJECTIONS |"
1136 print *, " |-------------------------------|"
1137 IB = CB_LOC(I)
1138 brickID = BRICK_LISt(NIN,IB)%ID
1139 write (*,FMT='(a,1i10,a,i2)') " --------brickID=",ixs(11,brickID), " ICUT=", ICUT
1140 write (*,FMT='(a )') " additional Scut (closed surf. hypothesis."
1141 write (*,FMT='(a,1i10)') " N1=",ITAB(IRECT(1,CE_LOC(I)))
1142 write (*,FMT='(a,1i10)') " N2=",ITAB(IRECT(2,CE_LOC(I)))
1143 write (*,FMT='(a,1i10)') " N3=",ITAB(IRECT(3,CE_LOC(I)))
1144 write (*,FMT='(a,1i10)') " N4=",ITAB(IRECT(4,CE_LOC(I)))
1145 write (*,FMT='(a)' ) " Projected Tria"
1146 write (*,FMT='(a,3f30.16)') " *createnode ",Ph(1:2,ICUT,1),ZERO
1147 write (*,FMT='(a,3f30.16)') " *createnode ",Ph(1:2,ICUT,2),ZERO
1148 write (*,FMT='(a,3f30.16)') " *createnode ",Ph(1:2,ICUT,3),ZERO
1149 write (*,FMT='(a,3f30.16)') " *createnode ",Ph(1:2,ICUT,4),ZERO
1150 write (*,FMT='(a,1i1)' ) " ityp=",ITYP
1151 write (*,FMT='(a)' ) " PolySCUT"
1152 write (*,FMT='(a,3f30.16)') " *createnode ",PolySCUT%node(1,1:2),ZERO
1153 write (*,FMT='(a,3f30.16)') " *createnode ",PolySCUT%node(2,1:2),ZERO
1154 write (*,FMT='(a,3f30.16)') " *createnode ",PolySCUT%node(3,1:2),ZERO
1155 write (*,FMT='(a,3f30.16)') " *createnode ",PolySCUT%node(4,1:2),ZERO
1156 write (*,FMT='(a,1i10)') " NTRIA=",NTRIA
1157 DO J=1,NTRIA
1158
1159 ENDDO
1160 print *, ""
1161 endif
1162
1163 POLYresult%NumNodes = 0
1164
1165 WET_RATIO(:) = ZERO
1166 VAL1 = ZERO
1167 VAL2 = ZERO
1168 Cumul_VAL1 = ZERO
1169
1170 DO J=1,NTRIA
1171 CALL PolygonalClipping( PolyTria(J),PolySCUT, POLYresult)
1172 IF(POLYresult%NumNodes==0)CYCLE
1173 !full area
1174 pt1(1:2) = PolyTria(J)%node(1,1:2) ; pt1(3)=ZERO
1175 pt2(1:2) = PolyTria(J)%node(2,1:2) ; pt2(3)=ZERO
1176 pt3(1:2) = PolyTria(J)%node(3,1:2) ; pt3(3)=ZERO
1177 !wet area
1178.AND. if(IBUG22_Swet==-1 INT22>0 )then
1179 write (*,FMT='(a,1i10)') " --TRIA #", J
1180 write (*,FMT='(a,3f30.16)') " *createnode ",PolyTria(J)%node(1,1:2),ZERO
1181 write (*,FMT='(a,3f30.16)') " *createnode ",PolyTria(J)%node(2,1:2),ZERO
1182 write (*,FMT='(a,3f30.16)') " *createnode ",PolyTria(J)%node(3,1:2),ZERO
1183 write (*,FMT='(a,1i10)') " ->clip nodes", POLYresult%NumNodes
1184 do k=1,POLYresult%NumNodes-1
1185 write (*,FMT='(a,3f45.35)') " *createnode ",POLYresult%node(k,1:2),ZERO
1186 enddo
1187 endif
1188
1189 PGrav(1:3) = ZERO
1190 !Gravity center estimation / To be updated later / must also be Parith ON
1191 IF(POLYresult%NumNodes-1>0)THEN
1192 DO k=1,POLYresult%NumNodes-1
1193 P_Grav(1,ICUT,I) = P_Grav(1,ICUT,I) + POLYresult%node(k,1)
1194 P_Grav(2,ICUT,I) = P_Grav(2,ICUT,I) + POLYresult%node(k,2)
1195 PGrav(1) = PGrav(1) + POLYresult%node(k,1)
1196 PGrav(2) = PGrav(2) + POLYresult%node(k,2)
1197 ENDDO
1198 PGrav(1) = PGrav(1) / (POLYresult%NumNodes-1)
1199 PGrav(2) = PGrav(2) / (POLYresult%NumNodes-1)
1200 ENDIF
1201
1202 P_Grav(1,ICUT,I) = P_Grav(1,ICUT,I) + PGrav(1)
1203 P_Grav(2,ICUT,I) = P_Grav(2,ICUT,I) + PGrav(2)
1204 NumGrav = Numgrav + ONE
1205
1206.AND. IF(POLYresult%NumNodes > 0 STria(J)/=ZERO)THEN
1207 pt1(1:2) = PolyTria(J)%node(1,1:2) ; pt1(3)=ZERO
1208 pt2(1:2) = PolyTria(J)%node(2,1:2) ; pt2(3)=ZERO
1209 pt3(1:2) = PolyTria(J)%node(3,1:2) ; pt3(3)=ZERO
1210 CALL SetCounterClockWisePolyg(POLYresult ,VAL)
1211 VAL = MAX(ZERO,MIN(VAL,Stria(J)))
1212 Cumul_VAL1 = Cumul_VAL1 + VAL
1213.AND. if(IBUG22_Swet==-1 INT22>0 )write (*,FMT='(a,1f30.16)') " -> Tria Wet Surf=", VAL
1214.AND. if(IBUG22_Swet==-1 INT22>0 )write (*,FMT='(a,1f30.16)') " -> Tria Surf =", Stria(J)
1215.AND. if(IBUG22_Swet==-1 INT22>0 )write (*,FMT='(a,1f30.16)') " -> cumulated Wet Surf=", Cumul_VAL1
1216.AND. if(IBUG22_Swet==-1 INT22>0 )write (*,FMT='(a,1f30.16)') " -> Sum of tria Surf=", Cumul_VAL2
1217
1218 !post-condition testing
1219 WET_RATIO(J) = VAL/Stria(J)
1220.OR. IF(WET_RATIO(J)<-EM06 WET_RATIO(J)-ONE>EM06)THEN
1221 print *, "***error: inter22 topology issue while computing wet surface",WET_RATIO(J),VAL1,VAL2
1222 ENDIF
1223 ENDIF
1224 ENDDO !next J=1,NTRIA
1225
1226 WET_RATIO(0) = Cumul_VAL1 !sum of wet surfaces for all triangles composing the lagrangian face.
1227
1228.OR. IF(WET_RATIO(0) == ZERO Cumul_VAL2==ZERO)THEN
1229 SEFF(ICUT,I) = ZERO
1230 CYCLE
1231 ELSE
1232 WET_RATIO(0) = WET_RATIO(0) / Cumul_VAL2
1233 ENDIF
1234
1235 WET_RATIO(0) = MAX(ZERO,MIN(WET_RATIO(0),ONE))
1236 SEFF(ICUT,I) = Cumul_VAL1
1237
1238 SEFF(ICUT,I) = WET_RATIO(0) * Sel(I)
1239
1240 !-------OUTPUT------
1241.AND. if(IBUG22_Swet==-1 INT22>0 )then
1242 write (*,FMT='(a,1f30.16)') " Sum Wet Surf =",Cumul_VAL1
1243 write (*,FMT='(a,1f30.16)') " Total Surf =",Cumul_VAL2
1244 write (*,FMT='(a,1f30.16)') " RATIO =",WET_RATIO(0)
1245 print *, ""
1246 endif
1247
1248 IF(NumGrav==0)CYCLE
1249
1250 P_Grav(1,ICUT,I) = P_Grav(1,ICUT,I) /NumGrav
1251 P_Grav(2,ICUT,I) = P_Grav(2,ICUT,I) /NumGrav
1252
1253 !---sWeighting factors using pressure load center
1254 !InterPolation using gravity centers of clipped polygons
1255 !NP_RECT(I) = 3 or 4
1256
1257 !2d projected segment
1258 pt1(1:2) = Ph(1:2,ICUT,ListP(1))
1259 pt2(1:2) = Ph(1:2,ICUT,ListP(2))
1260 pt3(1:2) = Ph(1:2,ICUT,ListP(3))
1261 pt4(1:2) = Ph(1:2,ICUT,ListP(4))
1262
1263C write (*,FMT='(A )') "Segment ="
1264C write (*,FMT='(A,2E30.16,A)') "*createnode ", Ph(1:2,ICUT,ListP(1)) ," 0.0"
1265C write (*,FMT='(A,2E30.16,A)') "*createnode ", Ph(1:2,ICUT,ListP(2)) ," 0.0"
1266C write (*,FMT='(A,2E30.16,A)') "*createnode ", Ph(1:2,ICUT,ListP(3)) ," 0.0"
1267C write (*,FMT='(A,2E30.16,A)') "*createnode ", Ph(1:2,ICUT,ListP(4)) ," 0.0"
1268
1269C write (*,FMT='(A )') "Scut ="
1270C do k=1,PolyScut%NumNodes-1
1271C write (*,FMT='(A,2E30.16,A)')"*createnode ", PolyScut%node(k,1:2) ," 0.0"
1272C enddo
1273
1274
1275 !distance of each point from wet polygon gravity center
1276 distance(1) = SQRT( (P_grav(1,ICUT,I)-Pt1(1))**2 + (P_grav(2,ICUT,I)-Pt1(2))**2 )
1277 distance(2) = SQRT( (P_grav(1,ICUT,I)-Pt2(1))**2 + (P_grav(2,ICUT,I)-Pt2(2))**2 )
1278 distance(3) = SQRT( (P_grav(1,ICUT,I)-Pt3(1))**2 + (P_grav(2,ICUT,I)-Pt3(2))**2 )
1279 distance(4) = SQRT( (P_grav(1,ICUT,I)-Pt4(1))**2 + (P_grav(2,ICUT,I)-Pt4(2))**2 )
1280
1281C write (*,FMT='(A )') "Pgrav ="
1282C write (*,FMT='(A,2E30.16,A)') "*createnode ", P_grav(1:2,ICUT,I) ," 0.0"
1283
1284
1285 !segment edges
1286 edge(1) = SQRT ( (Pt1(1)-Pt2(1))**2 + (Pt1(2)-Pt2(2))**2 )
1287 edge(2) = SQRT ( (Pt2(1)-Pt3(1))**2 + (Pt2(2)-Pt3(2))**2 )
1288 edge(3) = SQRT ( (Pt3(1)-Pt4(1))**2 + (Pt3(2)-Pt4(2))**2 )
1289 edge(4) = SQRT ( (Pt4(1)-Pt1(1))**2 + (Pt4(2)-Pt1(2))**2 )
1290 IF(NP_RECT(I)==3)THEN
1291 edge(3) = SQRT ( (Pt3(1)-Pt1(1))**2 + (Pt3(2)-Pt1(2))**2 )
1292 edge(4) = ZERO
1293 ENDIF
1294
1295 a = edge(1)
1296 b = distance(1)
1297 c = distance(2)
1298 s2 = HALF*(a+b+c)
1299 Area_tria(1) = SQRT( ABS(s2*(s2-a)*(s2-b)*(s2-c)) )
1300
1301 a = edge(2)
1302 b = distance(2)
1303 c = distance(3)
1304 s2 = HALF*(a+b+c)
1305 Area_tria(2) = SQRT( ABS(s2*(s2-a)*(s2-b)*(s2-c)) )
1306
1307 IF(NP_RECT(I)==4)THEN
1308 a = edge(3)
1309 b = distance(3)
1310 c = distance(4)
1311 s2 = HALF*(a+b+c)
1312 Area_tria(3) = SQRT( ABS(s2*(s2-a)*(s2-b)*(s2-c)) )
1313
1314 a = edge(4)
1315 b = distance(4)
1316 c = distance(1)
1317 s2 = HALF*(a+b+c)
1318 Area_tria(4) = SQRT( ABS(s2*(s2-a)*(s2-b)*(s2-c)) )
1319 ELSE
1320 a = edge(3)
1321 b = distance(3)
1322 c = distance(1)
1323 s2 = HALF*(a+b+c)
1324 Area_tria(3) = SQRT( ABS(s2*(s2-a)*(s2-b)*(s2-c)) )
1325
1326 Area_tria(4) = ZERO
1327 ENDIF
1328
1329 Dsum = SUM(Area_tria(1:NP_RECT(I)))
1330
1331 IF(Dsum==ZERO)then
1332 print *, "Warning: inter22 - nodal forces computation"
1333 endif
1334
1335 IF(NP_RECT(I)==4)THEN
1336 delta(1,ICUT,I) = HALF*(Area_tria(2)+Area_tria(3))/Dsum
1337 delta(2,ICUT,I) = HALF*(Area_tria(3)+Area_tria(4))/Dsum
1338 delta(3,ICUT,I) = HALF*(Area_tria(1)+Area_tria(4))/Dsum
1339 delta(4,ICUT,I) = HALF*(Area_tria(1)+Area_tria(2))/Dsum
1340 !uniform loading : same weight
1341 delta(1:4,ICUT,I) = FOURTH
1342 ELSE
1343 delta(1,ICUT,I) = (Area_tria(2))/Dsum
1344 delta(2,ICUT,I) = (Area_tria(3))/Dsum
1345 delta(3,ICUT,I) = (Area_tria(1))/Dsum
1346 delta(4,ICUT,I) = ZERO
1347 !uniform loading : same weight
1348 delta(1:4,ICUT,I) = THIRD
1349 delta(4,ICUT,I) = ZERO
1350 ENDIF
1351
1352 !write (*,FMT='(a,3e30.16)') "CoG = ", P_grav(1:2,ICUT,I)
1353 !write (*,FMT='(a,4e30.16)') "dist = ", distance(1:4)
1354
1355C write (*,FMT='(A,4E30.16,A,1E30.16)') "A_tria= ", Area_tria(1:4)," sum= ",SUM(Area_tria)
1356C write (*,FMT='(A,4E30.16,A,I2,A,I6)') "delta= ", delta(1:4,ICUT,I), "ICUT=", ICUT, " I=", I
1357
1358 ENDDO!next ICUT
1359
1360 ENDDO!next I (couple)
1361
1362C--------------------------------------------------------
1363C COMMENTAIRES
1364C -2- in the case where nbcut > 1 (scut not connected) need to handle elements that do not have a directed projection but are in the continuity of the surface (generally there is a connected part in the neighboring element)
1365C solution: extend scut by continuity by closing the facet
1366C (Useless if similar mesh size between Lagrane and Euler: Cicular Smaes)
1367C--------------------------------------------------------
1368
1369
1370C
1371 RETURN
1372 END
1373C===============================================================================
1374
1375
1376!||====================================================================
1377!|| getprojectedfacetype ../engine/source/interfaces/int22/i22wetsurf.F
1378!||--- called by ------------------------------------------------------
1379!|| i22wetsurf ../engine/source/interfaces/int22/i22wetsurf.F
1380!||--- calls -----------------------------------------------------
1381!|| intersectp ../engine/source/interfaces/int22/i22clip_tools.F
1382!||====================================================================
1383 FUNCTION GetProjectedFaceType(P1,P2,P3,P4,I,ListP)
1384C-----------------------------------------------
1385C I m p l i c i t T y p e s
1386C-----------------------------------------------
1387#include "implicit_f.inc"
1388C-----------------------------------------------
1389C D u m m y A r g u m e n t s
1390C-----------------------------------------------
1391 my_real :: P1(2),P2(2),P3(2),P4(2)
1392 my_real,intent(inout) :: I(2)
1393 INTEGER,intent(inout) :: ListP(4)
1394 INTEGER :: GetProjectedFaceType
1395C-----------------------------------------------
1396C L o c a l A r g u m e n t s
1397C-----------------------------------------------
1398 my_real :: PA(2), PB(2)
1399C-----------------------------------------------
1400C I n t e r f a c e B l o c k s
1401C-----------------------------------------------
1402 INTERFACE
1403 FUNCTION intersectP(a,b,c,d,iflg)
1404 my_real, intent(inout) :: a(2),b(2),c(2),d(2)
1405 my_real :: intersectP(2)
1406 integer,intent(in) :: iflg
1407 END FUNCTION
1408 END INTERFACE
1409C-----------------------------------------------
1410C S o u r c e L i n e s
1411C-----------------------------------------------
1412C
1413C
1414C 3 2 3 4 4 3
1415C +---------+ +---------+ +----------+
1416C \ / \ / | |
1417C \ / \ / | |
1418C \ / \ / | |
1419C \ / \ / | |
1420C + I + I | + I |
1421C / \ / \ | |
1422C / \ / \ | |
1423C / \ / \ | |
1424C / \ / \ | |
1425C +---------+ +---------+ +----------+
1426C 1 4 1 2 1 2
1427C
1428C GetProjectedFaceType = 2 GetProjectedFaceType = 1 GetProjectedFaceType = 0
1429C PA /= EP30 PA = EP30 PA = EP30
1430C PB = EP30 PB /= EP30 PB = EP30
1431C
1432C-----------------------------------------------
1433
1434 PA(1:2) = intersectP( P1, P2, P3, P4 ,0) ! [P1P2] & (P3P4)
1435 PB(1:2) = intersectP( P1, P4, P2, P3 ,0) ! [P1P4] & (P2P3)
1436
1437.AND. IF (PA(1)/=EP30 PB(1)==EP30)THEN
1438 GetProjectedFaceType = 2
1439 I(:) = PA(:)
1440 ListP(1:2) = (/1,4/) !defining tiangles
1441 ListP(3:4) = (/2,3/)
1442.AND. ELSEIF(PA(1)==EP30 PB(1)/=EP30)THEN
1443 GetProjectedFaceType = 1
1444 I(:) = PB(:)
1445 ListP(1:2) = (/1,2/) !defining tiangles
1446 ListP(3:4) = (/3,4/)
1447 ELSE
1448 GetProjectedFaceType = 0
1449 I(1) = FOURTH*(P1(1)+P2(1)+P3(1)+P4(1))
1450 I(2) = FOURTH*(P1(2)+P2(2)+P3(2)+P4(2))
1451 ListP(1:4) = (/1,2,3,4/)
1452 ENDIF
1453
1454
1455 END FUNCTION
1456
i22aera
function i22aera(npts, p, c)
Definition i22subvol.F:2382
i22wetsurf
subroutine i22wetsurf(jlt, cand_b, cand_e, cb_loc, ce_loc, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, xi, yi, zi, nx1, nx2, nx3, nx4, ny1, ny2, ny3, ny4, nz1, nz2, nz3, nz4, lb1, lb2, lb3, lb4, lc1, lc2, lc3, lc4, p1, p2, p3, p4, ix1, ix2, ix3, ix4, nsvg, stif, jlt_new, gapv, inacti, cand_p, n_scut, index, vxi, vyi, vzi, msi, kini, surf, ibag, itab, irect, i_stok, ixs, nft, cog, seff, delta, x)
Definition i22wetsurf.F:52
interp
subroutine interp(tf, tt, npoint, f, tg)
Definition interp.F:35
max
#define max(a, b)
Definition macros.h:21
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::nbcut_max
integer nbcut_max
Definition cut-cell-search_mod.F:517
i22tri_mod::ibug22_swet
integer ibug22_swet
Definition cut-cell-search_mod.F:342
i22tri_mod::ncande
integer ncande
Definition cut-cell-search_mod.F:499
area_tria
subroutine area_tria(x, n1, n2, n3, a2)
Definition s4mass3.F:410