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 1381 of file i22wetsurf.F.

1382C-----------------------------------------------
1383C I m p l i c i t T y p e s
1384C-----------------------------------------------
1385#include "implicit_f.inc"
1386C-----------------------------------------------
1387C D u m m y A r g u m e n t s
1388C-----------------------------------------------
1389 my_real :: p1(2),p2(2),p3(2),p4(2)
1390 my_real,intent(inout) :: i(2)
1391 INTEGER,intent(inout) :: ListP(4)
1392 INTEGER :: GetProjectedFaceType
1393C-----------------------------------------------
1394C L o c a l A r g u m e n t s
1395C-----------------------------------------------
1396 my_real :: pa(2), pb(2)
1397C-----------------------------------------------
1398C I n t e r f a c e B l o c k s
1399C-----------------------------------------------
1400 INTERFACE
1401 FUNCTION intersectp(a,b,c,d,iflg)
1402 my_real, intent(inout) :: a(2),b(2),c(2),d(2)
1403 my_real :: intersectp(2)
1404 integer,intent(in) :: iflg
1405 END FUNCTION
1406 END INTERFACE
1407C-----------------------------------------------
1408C S o u r c e L i n e s
1409C-----------------------------------------------
1410C
1411C
1412C 3 2 3 4 4 3
1413C +---------+ +---------+ +----------+
1414C \ / \ / | |
1415C \ / \ / | |
1416C \ / \ / | |
1417C \ / \ / | |
1418C + I + I | + I |
1419C / \ / \ | |
1420C / \ / \ | |
1421C / \ / \ | |
1422C / \ / \ | |
1423C +---------+ +---------+ +----------+
1424C 1 4 1 2 1 2
1425C
1426C GetProjectedFaceType = 2 GetProjectedFaceType = 1 GetProjectedFaceType = 0
1427C PA /= EP30 PA = EP30 PA = EP30
1428C PB = EP30 PB /= EP30 PB = EP30
1429C
1430C-----------------------------------------------
1431
1432 pa(1:2) = intersectp( p1, p2, p3, p4 ,0) ! [P1P2] & (P3P4)
1433 pb(1:2) = intersectp( p1, p4, p2, p3 ,0) ! [P1P4] & (P2P3)
1434
1435 IF (pa(1)/=ep30 .AND. pb(1)==ep30)THEN
1436 getprojectedfacetype = 2
1437 i(:) = pa(:)
1438 listp(1:2) = (/1,4/) !defining tiangles
1439 listp(3:4) = (/2,3/)
1440 ELSEIF(pa(1)==ep30 .AND. pb(1)/=ep30)THEN
1441 getprojectedfacetype = 1
1442 i(:) = pb(:)
1443 listp(1:2) = (/1,2/) !defining tiangles
1444 listp(3:4) = (/3,4/)
1445 ELSE
1446 getprojectedfacetype = 0
1447 i(1) = fourth*(p1(1)+p2(1)+p3(1)+p4(1))
1448 i(2) = fourth*(p1(2)+p2(2)+p3(2)+p4(2))
1449 listp(1:4) = (/1,2,3,4/)
1450 ENDIF
1451
1452
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:1382

◆ 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 35 of file i22wetsurf.F.

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