OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i22wetsurf.F File Reference
#include "implicit_f.inc"
#include "comlock.inc"
#include "inter22.inc"
#include "com04_c.inc"
#include "mvsiz_p.inc"
#include "lockon.inc"
#include "lockoff.inc"

Go to the source code of this file.

Functions/Subroutines

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)
integer function getprojectedfacetype (p1, p2, p3, p4, i, listp)

Function/Subroutine Documentation

◆ getprojectedfacetype()

integer function getprojectedfacetype ( p1,
p2,
p3,
p4,
dimension(2), intent(inout) i,
integer, dimension(4), intent(inout) listp )

Definition at line 1383 of file i22wetsurf.F.

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 IF (pa(1)/=ep30 .AND. 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 ELSEIF(pa(1)==ep30 .AND. 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
my_real
#define my_real
Definition cppsort.cpp:32
intersectp
function intersectp(seg_p1, seg_p2, line_p1, line_p2, iflg)
Definition i22clip_tools.F:189
getprojectedfacetype
integer function getprojectedfacetype(p1, p2, p3, p4, i, listp)
Definition i22wetsurf.F:1384

◆ i22wetsurf()

subroutine i22wetsurf ( integer jlt,
integer, dimension(*) cand_b,
integer, dimension(*) cand_e,
integer, dimension(mvsiz) cb_loc,
integer, dimension(mvsiz) 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,
integer, dimension(mvsiz) ix1,
integer, dimension(mvsiz) ix2,
integer, dimension(mvsiz) ix3,
integer, dimension(mvsiz) ix4,
integer, dimension(mvsiz) nsvg,
stif,
integer jlt_new,
gapv,
integer inacti,
cand_p,
n_scut,
integer, dimension(mvsiz) index,
vxi,
vyi,
vzi,
msi,
integer, dimension(*) kini,
surf,
integer ibag,
integer, dimension(*) itab,
integer, dimension(4,*) irect,
integer i_stok,
integer, dimension(nixs,*) ixs,
integer nft,
cog,
seff,
delta,
x )

Definition at line 36 of file i22wetsurf.F.

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 IF(polyresult%NumNodes > 0 .AND. 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 if(ibug22_swet==-1 .AND. int22>0 )write (*,fmt='(A,1F30.16)') " -> Tria Wet Surf=", val
914 if(ibug22_swet==-1 .AND. int22>0 )write (*,fmt='(A,1F30.16)') " -> Tria Surf =", stria(j)
915 if(ibug22_swet==-1 .AND. int22>0 )write (*,fmt='(A,1F30.16)') " -> cumulated Wet Surf=", cumul_val1
916 if(ibug22_swet==-1 .AND. 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 IF(wet_ratio(j)<-em06 .OR. 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 IF(wet_ratio(0) == zero .OR. 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 if(ibug22_swet==-1 .AND. 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 if(ibug22_swet==-1 .AND. 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 if(ibug22_swet==-1 .AND. 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 IF(polyresult%NumNodes > 0 .AND. 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 if(ibug22_swet==-1 .AND. int22>0 )write (*,fmt='(A,1F30.16)') " -> Tria Wet Surf=", val
1214 if(ibug22_swet==-1 .AND. int22>0 )write (*,fmt='(A,1F30.16)') " -> Tria Surf =", stria(j)
1215 if(ibug22_swet==-1 .AND. int22>0 )write (*,fmt='(A,1F30.16)') " -> cumulated Wet Surf=", cumul_val1
1216 if(ibug22_swet==-1 .AND. 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 IF(wet_ratio(j)<-em06 .OR. 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 IF(wet_ratio(0) == zero .OR. 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 if(ibug22_swet==-1 .AND. 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
setcounterclockwisepolyg
subroutine setcounterclockwisepolyg(polyg, area)
Definition i22clip_tools.F:350
polygonalclipping
subroutine polygonalclipping(basepolygon, clippolygon, resultpolygon)
Definition i22clip_tools.F:60
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
min
#define min(a, b)
Definition macros.h:20
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