OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i22intersect.F File Reference
#include "implicit_f.inc"
#include "comlock.inc"
#include "mvsiz_p.inc"
#include "task_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine i22intersect (x, ii_stok, cand_b, cand_e, itask, nbric, itab, bufbric, ncand, ixs, nin)
integer function isonsh3n (z, a, za, zb, p, iflg_db)

Function/Subroutine Documentation

◆ i22intersect()

subroutine i22intersect ( x,
integer ii_stok,
integer, dimension(ncand) cand_b,
integer, dimension(ncand) cand_e,
integer itask,
integer nbric,
integer, dimension(*) itab,
integer, dimension(nbric) bufbric,
integer ncand,
integer, dimension(nixs,*) ixs,
integer nin )

Definition at line 38 of file i22intersect.F.

42C============================================================================
43C-----------------------------------------------
44C D e s c r i p t i o n
45C-----------------------------------------------
46C Interface Type22 (/INTER/TYPE22) is an FSI coupling method based on cut cell method.
47C This experimental cut cell method is not completed, abandoned, and is not an official option.
48C
49C X(3,*) is used for local nodes (cells) for surface nodes use IRECTL
50C-----------------------------------------------
51C M o d u l e s
52C-----------------------------------------------
53 USE i22bufbric_mod
54 USE i22tri_mod
55 USE i22edge_mod
56 USE realloc_mod
57 use element_mod , only : nixs
58C-----------------------------------------------
59C I m p l i c i t T y p e s
60C-----------------------------------------------
61#include "implicit_f.inc"
62#include "comlock.inc"
63C-----------------------------------------------
64C G l o b a l P a r a m e t e r s
65C-----------------------------------------------
66#include "mvsiz_p.inc"
67C-----------------------------------------------
68C C o m m o n B l o c k s
69C-----------------------------------------------
70#include "task_c.inc"
71C-----------------------------------------------
72C=======================================================================
73C 6 creation of edge list and addresses for each candidate cell
74C================================================================
75C-----------------------------------------------
76C D u m m y A r g u m e n t s
77C-----------------------------------------------
78 INTEGER CAND_B(NCAND),CAND_E(NCAND), NCAND, NIN,
79 . ITASK, NBRIC, ITAB(*),
80 . BUFBRIC(NBRIC), IXS(NIXS,*), II_STOK
81 my_real :: x(3,*)
82C-----------------------------------------------
83C L o c a l V a r i a b l e s
84C-----------------------------------------------
85 INTEGER NB_NCN,NB_NCN1,NB_ECN,I,J,K,L,DIR,NB_NC,NB_EC,
86 . N1,N2,N3,N4,NN,NE,NCAND_PROV,J_STOK,II,JJ,TT,
87 . NSNF, NSNL, TANGENT(12),
88 . PROV_B(2*MVSIZ), PROV_E(2*MVSIZ), LAST_NE,
89 . VOXBND(2*MVSIZ,0:1,1:3) !voxel bounds storage for shell: comp1=id, comp2=lbound/ubound, comp3=direction.
90
91 my_real
92 . dx,dy,dz,xs,ys,zs,xx,sx,sy,sz,s2,
93 . xmin, xmax,ymin, ymax,zmin, zmax, tz, gapsmx, gapl,
94 . d1x,d1y,d1z,d2x,d2y,d2z,dd1,dd2,d2,a2,gs, point(3),point2(3),d_1,d_2,d_3,
95 . on1(3),n1n2(3)
96
97 INTEGER IX,IY,IZ,NEXT,M(8),
98 . IX1,IY1,IZ1,IX2,IY2,IZ2,IBUG,IBUG2,I_LOC,
99 . BIX1(NBRIC),BIY1(NBRIC),BIZ1(NBRIC),
100 . BIX2(NBRIC),BIY2(NBRIC),BIZ2(NBRIC),
101 . FIRST_ADD, PREV_ADD, LCHAIN_ADD, I_STOK , idb_ID
102
103 INTEGER, DIMENSION(1) :: SHELL_ADD !address of surface face which is candidate in CAND_E (temp)
104 !possible to decrease allocated size later
105 INTEGER :: NC, I_STOK_BAK, IPA,IPB
106 my_real
107 . xminb,yminb,zminb,xmaxb,ymaxb,zmaxb,
108 . dxb,dyb,dzb,
109 . aaa, basisconst2,ns,
110 . power(8), cutcoor,cutcoor2,cut(2),
111 . pow(2), old_cutcoor, old_cutshell, cutnode(2)
112
113 LOGICAL :: IsSH3N
114
115 LOGICAL, DIMENSION(NBRIC) :: COUNT
116! Logical, dimension (12*nbric) :: ledge
117 LOGICAL :: BOOL(NIRECT_L)
118 INTEGER NBCUT, NBCUT2,DEJA, ISONSHELL, ISONSH3N
119 INTEGER :: COUNTER, NEDGE, NFACE, NODES8(8), COUNTER_BRICK(NBRIC)
120
121 INTEGER :: iN(2), iN1a, iN2a, iN1b, iN2b , iN3, iN4
122 INTEGER :: POS, IAD, IADE, IB ,IBG , NBF, NBL
123 INTEGER :: I_12bits, nbits, npqts, pqts(4), SOM, SECTION
124 INTEGER :: I_bits(12), MAX_ADD, IMIN_LOC, IMAX_LOC
125
126 my_real :: aeradiag, debugtab(24*ncand,3)
127 LOGICAL db_FLAG, TAGnode(8), debug_outp
128
129 CHARACTER*12 :: sectype
130 CHARACTER*12 ::filename
131 LOGICAL :: IsSecDouble
132
133 CHARACTER(LEN=1) filenum
134
135 INTEGER ::
136 . MIN_IX_LOC, MIN_IY_LOC, MIN_IZ_LOC, !indice voxel min used
137 . MAX_IX_LOC, MAX_IY_LOC, MAX_IZ_LOC !indice voxel max used
138
139 INTEGER :: ISHIFT, IDX
140
141 my_real, dimension(:), allocatable :: powb
142
143 INTEGER :: A(5), IE, N_CUT_EDGE
144
145 INTEGER :: TAG_INDEX(NBRIC), I8(8,NBRIC),IFLG_DB
146 my_real :: r8(8,nbric), denom, norm, tolcrit,tol
147
148 TYPE(LIST_SECND) :: OLD_SECND_LIST
149
150C----------------------------------------------------------------
151 a(1:5)=((/1,2,3,4,1/))
152
153 idb_id = 0
154
155
156C--------------
157C
158C
159C 60 61 62
160C +-----------+-----------+
161C /| /| /|
162C / | / | / |
163C / | / | / |
164C 50+---------51+---------52+ |
165C | |-------|---+-------|---+42
166C | /| | /| | /|
167C | / | | / | | / |
168C |/ | |/ | |/ |
169C 30+---------31+---------32+ |
170C | |-------|---+-------|---+
171C | / 20 | / 21 | / 22
172C | / | / | /
173C |/ |/ |/
174C +-----------+-----------+
175C 10 11 12
176C .
177C
178C-----------------------------------------------
179
180 nb = ncandb
181 nbf = 1+itask*nb/nthread
182 nbl = (itask+1)*nb/nthread
183 nin = 1
184
185 cutcoor = -ep30
186
187C=======================================================================
188C 1 Allocation/INIT list of edges
189C================================================================
190
191 !---------------------------------------------------------!
192 ! DEBUG OUTPUT !
193 !---------------------------------------------------------!
194 !INTERFACE 22 ONLY - OUTPUT---------------!
195 debug_outp = .false.
196 if(ibug22_intersect/=0)then
197 if(ibug22_intersect>0)then
198 do ib=nbf,nbl
199 ie=bufbric(list_b(ib))
200 if(ixs(11,ie)==ibug22_intersect)then
201 debug_outp=.true.
202 exit
203 endif
204 enddo
205 elseif(ibug22_intersect==-1)then
206 debug_outp = .true.
207 endif
208 endif
209 if(itask==0.AND.debug_outp)then
210 print *, ""
211 print *, "================================="
212 print *, "==== BRICK INTERSECTIONS ====="
213 print *, "================================="
214 endif
215
216 nb = ncandb ! decreasing list of cells which are really intersected. Above we needed also uncut cells
217 nbf = 1+itask*nb/nthread
218 nbl = (itask+1)*nb/nthread
219 nin = 1
220
221 IF (itask==0)THEN
222 ALLOCATE(vne(ncande,3,4))
223 ALLOCATE(basisconst(ncande,4))
224 ALLOCATE(nbsubtriangles(ncande))
225 ALLOCATE(diag22(nb))
226 ALLOCATE(ptz(ncande,3))
227 ALLOCATE(vza(ncande,3,4))
228 ALLOCATE(vzb(ncande,3,4))
229 ALLOCATE(pta(ncande,3,4))
230 END IF
231 CALL my_barrier
232
233!internal cell number : BRICK_LIST(NIN,I)%ID
234!local number in group : LIST_B(I)
235
236
237 !temp array (init)
238 DO i=nbf, nbl ! 1,NB (threads)
239 brick_list(nin,i)%ID = bufbric(list_b(i)) ! internal cell id
240 d_1 = zero
241 d_2 = zero
242 d_3 = zero
243 DO j=1,12
244 k=(i-1)*12+j
245 edge_list(nin,k)%NODE(1:2) = ixs(int22_buf%iEDGE(1:2,j)+1,brick_list(nin,i)%ID) ! ids cell nodes
246 edge_list(nin,k)%NBCUT = 0
247 edge_list(nin,k)%CUTSHELL = 0
248 edge_list(nin,k)%CUTCOOR(1:2) = zero
249 edge_list(nin,k)%VECTOR(1:3) = x(1:3,edge_list(nin,k)%NODE(2))-x(1:3,edge_list(nin,k)%NODE(1)) ! vector (iN1->iN2) to be store to compute intersections
250 d_1 = max(d_1, abs(edge_list(nin,k)%VECTOR(1)))
251 d_2 = max(d_2, abs(edge_list(nin,k)%VECTOR(2)))
252 d_3 = max(d_3, abs(edge_list(nin,k)%VECTOR(3)))
253 END DO
254 diag22(i) = max( d_1,d_2,d_3 )
255 END DO
256 CALL my_barrier !waiting for the list to be created before going on
257
258 if(itask==0 .AND. debug_outp)print *, ""
259
260
261
262C=======================================================================
263C 2 a.Normal & basis point
264C Storage in vNE(1:NCANDE,1:3,TT)
265C --> 1st arg : candidate Face (structure)
266C --> 2nd arg : components x,y,z
267C --> 3rd arg : normal number (1shell=4tria)
268C B.PLANE EQUATION (Constant Value)
269C BasisCONST(TT)
270C --> 1st arg : normal number (1shell=4tria)
271C================================================================
272 nbf = 1+itask*ncande/nthread
273 nbl = (itask+1)*ncande/nthread
274
275 DO i=nbf,nbl!1,NCANDE
276 ne = list_e(i) !shell number IRECT(1:4,*)
277 m(3:4)=irect_l(3:4,ne)
278 IF(m(3)==m(4))THEN !type : sh4n or sh3n
279 issh3n=.true.
280 nbsubtriangles(i)=1
281 ELSE
282 issh3n=.false.
283 nbsubtriangles(i)=4 ;
284 END IF
285 IF(.NOT.issh3n)THEN !centre
286 !ptZ(I,1:3)=(/(FOURTH*SUM( IRECT_L( 4*J+(/1:4/),NE) ),J=1,3)/) ! Z = mean X(5:8), mean Y(9:12), mean Z(13:16)
287 ptz(i,1)=fourth*sum( irect_l( 05:08,ne) )
288 ptz(i,2)=fourth*sum( irect_l( 09:12,ne) )
289 ptz(i,3)=fourth*sum( irect_l( 13:16,ne) )
290 ELSE
291 !ptZ(I,1:3)=IRECT_L((/8,12,16/),NE) !M(4)
292 ptz(i,1)=irect_l(08,ne) !M(4)
293 ptz(i,2)=irect_l(12,ne) !M(4)
294 ptz(i,3)=irect_l(16,ne) !M(4)
295 END IF
296 DO tt=1,nbsubtriangles(i) !Triangle comping (Z, IPA, IPB) = (Z, 1,2) or (Z, 2.3) or (Z, 3.4) or (Z, 4.1)
297 ipa=a(tt) ! Z is centroid of 4node-shell.
298 ipb=a(tt+1) !; IF(IPB==5) IPB=1,for cyclic permutation
299 !IPA=IRECT_L(IPA,NE)
300 !sIPB=IRECT_L(IPB,NE)
301 pta(i,1:3,tt) = irect_l((/4,8,12/)+ipa,ne) !A coordinates (for intersection computation ISONSH3N) <=> X(1:3,IPA)
302 vza(i,1:3,tt) = irect_l((/4,8,12/)+ipa,ne)-ptz(i,1:3) !vector ZA et ZB for each triangle.
303 vzb(i,1:3,tt) = irect_l((/4,8,12/)+ipb,ne)-ptz(i,1:3) !vector ZA et ZB for each triangle.
304 vne(i,1:3,tt) = vza(i,(/2,3,1/),tt)*vzb(i,(/3,1,2/),tt) - ! vNE = vZA (x) vZB
305 . vza(i,(/3,1,2/),tt)*vzb(i,(/2,3,1/),tt)
306 !normalizing the normal vector to normalize the cartesian plane equation
307 norm = vne(i,1,tt)*vne(i,1,tt)+vne(i,2,tt)*vne(i,2,tt)+vne(i,3,tt)*vne(i,3,tt)
308 norm = sqrt(norm)
309 IF(norm/=zero)THEN
310 vne(i,1,tt) = vne(i,1,tt) / norm
311 vne(i,2,tt) = vne(i,2,tt) / norm
312 vne(i,3,tt) = vne(i,3,tt) / norm
313 ENDIF
314 basisconst(i,tt) = sum(ptz(i,1:3)*vne(i,1:3,tt)) !Constant basis (Equation surface)
315 !basis point is Zentrum, it is common for each triangle
316 !PLANE EQUATION : nx.X + nY.Y + nZ.Z + k = 0
317 ! where n is vNE
318 ! k such as ptZ satisfying equation => k = -nx.Zx - ny.Zy - nz.Zz is BasisCONST
319 END DO !next TT
320 END DO !next I
321
322 CALL my_barrier !wait for complete initialization
323
324C=======================================================================
325C 3 Calculation Intersections
326C================================================================
327 if(itask==0 .AND. debug_outp)then
328 print *," Calcul des Intersections sur Proc=", itask+1
329 endif
330
331 !to be optimized : loop over cells from list_b (need to get corresponding addresses)
332 nbf = 1+itask*ncandb/nthread
333 nbl = (itask+1)*ncandb/nthread
334
335 !NCANDB bricks are present in candidates. These NCANDB bricks are multi-threaded. Each thread will handle ( (/CAND_B(IADF(I):IADL(I)) , CAND_B(IADF(I):IADL(I))/), I=NBF,NBL)
336
337 DO i=nbf,nbl !each thread handle its bricks (no duplicate work)
338 brick_list(nin,i)%Seg_add_LFT = iadf(i) ! corresponding window in CAND_B(i) : CAND_B(IADFL(I):IADL(I)) = I
339 brick_list(nin,i)%Seg_add_LLT = iadl(i)
340 cand_b(iadf(i):iadl(i)) = i !TRANSFORM CAND_B(I) into ADD into brick_list()
341c print *, "IB =", I
342c print *, "....IADF,IADL=", IADF(I), IADL(I)
343 ib = list_b(i) !local ID in group IBUFSSG(IAD+(/1:NBRIC/))
344 ibg = bufbric(ib) !global ID
345 tagnode=.false.
346 tangent(1:12) = 0 !detect tangent edge : keep them uncut, its adjacent connected edges will be cut
347
348 if (ixs(11,brick_list(nin,i)%id)==idb_id)then
349 print *, "idb_ID====="
350 print *, "CAND_E =", cand_e(iadf(i):iadl(i))
351 endif
352 ! if (ixs(11,brick_list(nin,i)%id)==0)then
353 ! print *, "481====="
354 ! print *, "CAND_E =", CAND_E(IADF(I):IADL(I))
355 ! endif
356
357 DO iad=iadf(i), iadl(i)
358 ie = cand_e(iad)
359 IF(ie<=0)cycle
360 iade = get_list_e_pos_from_cand_e_pos(iad)
361C !non intersection criterion (couple is retained anyway in CAND_B/CAND_E because we may not search/sort on each cycle)
362C IF((XMAXS(IB)<IRECT_L(17,IE)).OR. ! ignore if intersection with cell (without margin)is null, IRECT_L <=> xmine,xmaxe,...zmine,zmaxe
363C . (XMINS(IB)>IRECT_L(20,IE)).OR.
364C . (YMAXS(IB)<IRECT_L(18,IE)).OR.
365C . (YMINS(IB)>IRECT_L(21,IE)).OR.
366C . (ZMAXS(IB)<IRECT_L(19,IE)).OR.
367C . (ZMINS(IB)>IRECT_L(22,IE)) ) THEN
368C! CAND_E(IAD) = 0
369C CYCLE !next couple candidat
370C END IF
371 DO tt=1,nbsubtriangles(iade) !shell decomposition to get rid of warping errors
372 power(1:8)=(/(sum(vne(iade,1:3,tt) * x(1:3,ixs(ii,ibg)))- basisconst(iade,tt),ii=2,9)/) !power of nodes for plane generated by current triangle : POW = distance = origin ordinate in working system
373 n_cut_edge = 0
374 DO j=1,12 !12 edges !opposite sign for nodes power ?
375 k = (i-1)*12+j !k : address number in EDGE_LIST
376 pow(1:2)= power(int22_buf%iEDGE(1:2, j))
377 in(1:2) = edge_list(nin,k)%NODE(1:2) !IXS(1+iEDGE(1, J),NS)
378 deja = edge_list(nin,k)%NBCUT !is edge already cut by another shell ?
379 nbcut = -1
380 nbcut2 = -1
381
382 if(ixs(11,brick_list(nin,i)%id)==idb_id)then
383 print *, "J=", j, itab(in(1:2))
384 write(*,fmt='(A,4I20)') "shell N1-N2-N3 :",int(irect_l(01:04,ie))
385 write(*,fmt='(A,3F20.12)') " shell N1 :",irect_l( (/05,09,13/),ie)
386 write(*,fmt='(A,3F20.12)') " shell N2 :",irect_l( (/06,10,14/),ie)
387 write(*,fmt='(A,3F20.12)') " shell N3 :",irect_l( (/07,11,15/),ie)
388 write(*,fmt='(A,2F40.20)')" POW(1:2)=", pow(1:2)
389 print *,""
390 endif
391 !(2 intersection points : if edge is on the plane)
392 if (ixs(11,brick_list(nin,i)%id)==idb_id )then
393 print *, "idb_ID====="
394 write(*,fmt='(A,4I20)') "shell N1-N2-N3-N4 :",int(irect_l(01:04,ie))
395 print *, "IE =", ie
396 print *, "TT =", tt
397 print *, "J =", j
398 print *, "POW1,POW2", pow(1:2)
399 endif
400 !if (ixs(11,brick_list(nin,i)%id)==5882 .AND. J==12)then
401 ! print *, "5882====="
402 ! write(*,FMT='(A,4I20)') "shell N1-N2-N3-N4 :",INT(IRECT_L(01:04,IE))
403 ! print *, "IE =", IE
404 ! print *, "TT =", TT
405 ! print *, "J =", J
406 ! print *, "POW1,POW2", POW(1:2)
407 !endif
408
409 tolcrit = em06
410 tol = (one+em04)*tolcrit*diag22(i) !EM04 in ISONSH3N
411 ! POWER is in mm corresponds to the shifted distance of the intercept (ordonnee a lorigine) it must be normalized
412 IF((abs(pow(1))<=tol).AND.(abs(pow(2))<=tol))THEN !edge is on the shell then no intersection with the face
413 tangent(j) = 1
414 !print *, "TANGENT=1",ixs(11,brick_list(nin,i)%id)
415 cycle
416 ENDIF
417 IF(deja==2)cycle !hypothesis: 2 intersections max. Warning possible dpedency of elem numbering
418 IF( ((pow(1)<-tol).AND.(pow(2)>tol)) .OR.((pow(1)>tol).AND.(pow(2)<-tol)) )THEN !intersection with current shell
419 on1(1:3) = x(1:3,in(1))
420 n1n2(1:3)= edge_list(nin,k)%VECTOR(1:3)
421 denom = sum( vne(iade,1:3,tt) * n1n2(1:3) )
422 IF(abs(denom)>em12)THEN
423 cutcoor = ( basisconst(iade,tt) - sum( vne(iade,1:3,tt) * on1(1:3) ) ) / denom !N.N1 change all 3 cycles (possible optimization, cf definition of iedge)
424 ELSE
425 !in this case edge is intersected with a plane which is quasi tangent. possibly infinite solutions. We choose cutcoor=0.5, error related to this description is neglictible
426 cutcoor = half
427 CALL arret(2)
428 ENDIF
429
430 IF((cutcoor<=one+tol).AND.(cutcoor>=-tol))THEN !check if intersection point is one on the shell
431 cutcoor = min(one-em06,cutcoor)
432 cutcoor = max(em06 ,cutcoor)
433 point(1:3)=on1(1:3) + cutcoor * n1n2(1:3)
434 ELSE
435 print *, " CUTCOOR =", cutcoor
436 CALL arret(2)
437 END IF
438
439 nbcut2 = 0
440 ELSEIF((abs(pow(1))<=tol).AND.(abs(pow(2))<=tol))THEN !edge is on the shell then no intersection with the face
441 !DON'T CUT EDGE WHICH IS LAYING ON THE CUT PLANE. ITS ADJACENT EDGES WILL BE CUT.
442 ! specific case : edge is on the structural face
443 !IF(J==1 .OR. J==5 .OR. J==8 .OR. J==12)THEN
444 ! ON1(1:3) = X(1:3,iN(1))
445 ! N1N2(1:3)= EDGE_LIST(NIN,K)%VECTOR(1:3)
446 ! CUTCOOR = ONE-EM06 !intserection on edge nodes1,5,8 ou 12 (necessary case to get a corresponding topology without risking to loose the lagrangian surface while it moves)
447 ! CUTCOOR2 = EM06
448
449
450 !ELSE
451 !-!ON1(1:3) = X(1:3,iN(1))
452 !-!N1N2(1:3)= EDGE_LIST(NIN,K)%VECTOR(1:3)
453 !-!CUTCOOR = HALF
454 !-!POINT(1:3) = ON1(1:3) + CUTCOOR * N1N2(1:3)
455 !-!NBCUT2 = 0
456 !ENDIF
457 tangent(j) = 1
458 ELSEIF (abs(pow(1))<=tol)THEN !intersection with summit #1
459 !SET CUTCOOR = EM06 TO HAVE A EPSILON CUT SURFACE. OTHERWISE CRITERIA FOR AMBIGUS COMBINATION WILL GET A DIVISION BY 0 (double penta check in i22subvol)
460 !IF(J==1 .OR. J==5 .OR. J==8 .OR. J==12)THEN
461 on1(1:3) = x(1:3,in(1))
462 n1n2(1:3)= edge_list(nin,k)%VECTOR(1:3)
463 cutcoor = em06 !intserection on edge nodes 1,5,8 ou 12 (necessary case to get a corresponding topology without risking to loose the lagrangian surface while it moves)
464 point(1:3) = on1(1:3) + zero * n1n2(1:3)
465 nbcut2 = 0
466 !ELSE
467 ! NBCUT = 0
468 ! NBCUT2 = 0
469 !ENDIF
470! IF(TAGnode(iEDGE(1,J))==.FALSE.)THEN
471! POINT(1:3)=X(1:3,iN(1))
472! CUTCOOR=ZERO
473! TAGnode(iEDGE(1,J))=.TRUE.
474! ELSE
475! NBCUT = 0
476! ENDIF
477! NBCUT2 = 0
478 ELSEIF (abs(pow(2))<=tol)THEN !intersection with summit #2
479 !SET CUTCOOR = ONE-EM06 TO HAVE A EPSILON CUT SURFACE. OTHERWISE CRITERIA FOR AMBIGUS COMBINATION WILL GET A DIVISION BY 0 (double penta check in i22subvol)
480 !IF(J==1 .OR. J==5 .OR. J==8 .OR. J==12)THEN
481 on1(1:3) = x(1:3,in(1))
482 n1n2(1:3)= edge_list(nin,k)%VECTOR(1:3)
483 cutcoor = one-em06 !intserection with edge nodes 1,5,8 ou 12 (necessary case to get a corresponding topology without risking to loose the lagrangian surface while it moves)
484 point(1:3) = on1(1:3) + one * n1n2(1:3)
485 nbcut2 = 0
486 !ELSE
487 ! NBCUT = 0
488 ! NBCUT2 = 0
489 !ENDIF
490! IF(TAGnode(iEDGE(2,J))==.FALSE.)THEN
491! POINT(1:3)=X(1:3,iN(2))
492! CUTCOOR=ONE
493! TAGnode(iEDGE(2,J))=.TRUE.
494! ELSE
495! NBCUT = 0
496! ENDIF
497! NBCUT2 = 0
498 ELSE !2 nodes of the edges are on the same side of the intersection plane
499 nbcut = 0
500 nbcut2 = 0
501 END IF
502
503 if (ixs(11,brick_list(nin,i)%id)==idb_id )then
504 print *, "cutcoor=", cutcoor
505 iflg_db=1
506 else
507 iflg_db=0
508 endif
509
510 IF(nbcut==-1) nbcut =isonsh3n( ptz(iade,1:3),pta(iade,1:3,tt),vza(iade,1:3,tt),vzb(iade,1:3,tt),point(1:3) ,iflg_db) ! Check if point (:) belongs to the triangle
511 IF(nbcut2==-1)nbcut2=isonsh3n( ptz(iade,1:3),pta(iade,1:3,tt),vza(iade,1:3,tt),vzb(iade,1:3,tt),point2(1:3),iflg_db) ! Check if point (:) belongs to the triangle
512
513 if (ixs(11,brick_list(nin,i)%id)==idb_id )then
514 print *, "NBCUT, NBCUT2=", nbcut,nbcut2
515 endif
516
517 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
518 !If intersection with infinite plane, we calculate coordinates on the edge : COOR = N1+ t N1N2
519 ! if t is not within [0,1] we are not on the segment but on the open entity (N1N2)\‍[N1N2]
520 IF(nbcut>0) THEN !if we detect intersection on current shell
521 IF(deja==0)THEN !edge is not yet cut : first intersection detected
522 !print *, "nouveau point", CUTCOOR
523 nbcut=1
524 edge_list(nin,k)%CUTCOOR(1) = cutcoor ! storage
525 edge_list(nin,k)%CUTSHELL(1) = ie
526 ELSEIF (deja==1)THEN ! already cut => 2nd intersection
527 old_cutcoor = edge_list(nin,k)%CUTCOOR(1) ! first store, the ordering the 2 intersections
528 old_cutshell = edge_list(nin,k)%CUTSHELL(1)
529 IF (abs(cutcoor-old_cutcoor)>em6) THEN
530 nbcut=2
531 IF(cutcoor>old_cutcoor)THEN
532 edge_list(nin,k)%CUTCOOR(1) = old_cutcoor
533 edge_list(nin,k)%CUTCOOR(2) = cutcoor
534 edge_list(nin,k)%CUTSHELL(1) = old_cutshell
535 edge_list(nin,k)%CUTSHELL(2) = ie
536 !print *, "second points point",OLD_CUTCOOR,CUTCOOR
537 ELSEIF(cutcoor<old_cutcoor)THEN
538 edge_list(nin,k)%CUTCOOR(1) = cutcoor
539 edge_list(nin,k)%CUTCOOR(2) = old_cutcoor
540 edge_list(nin,k)%CUTSHELL(1) = ie
541 edge_list(nin,k)%CUTSHELL(2) = old_cutshell
542 !print *, "second points point",CUTCOOR,OLD_CUTCOOR
543 ELSE
544 !print *, "point deja trouve", CUTCOOR
545 nbcut=1 ! the point is the same, no need to record it
546 END IF
547 END IF
548 ELSEIF(deja==2)THEN ! already 2 intersections, no more expected (hypothesis)
549 if(itask==0 .AND. debug_outp)then
550 if(ibug22_intersect==-1 .or. ibug22_intersect==ixs(11,brick_list(nin,i)%id))then
551 print *, "THREE INTERSECTION SUR UNE ARRETE - STOP"
552 CALL arret(2)
553 endif
554 endif
555 END IF
556 edge_list(nin,k)%NBCUT = nbcut
557 edge_list(nin,k)%LEN = sqrt(n1n2(1)*n1n2(1)+n1n2(2)*n1n2(2)+n1n2(3)*n1n2(3))
558 END IF !NBCUT>0
559
560 IF(nbcut2>0 .AND. deja==0)THEN !if already intersected
561 nbcut=2
562 edge_list(nin,k)%CUTCOOR(1)=zero ! storage
563 edge_list(nin,k)%CUTCOOR(2)=one ! storage
564 if(itask==0 .AND. debug_outp)then
565 if(ibug22_intersect==-1 .or. ibug22_intersect==ixs(11,brick_list(nin,i)%id))then
566 print *, "edge fully on intersection plane",j
567 endif
568 endif
569 edge_list(nin,k)%CUTSHELL(1) = ie
570 edge_list(nin,k)%CUTSHELL(2) = ie
571 edge_list(nin,k)%NBCUT = 2
572 ENDIF
573
574 END DO !next J=1,12 // next edge
575 END DO !next TT=1,NbSubTriangles // next sub triangles
576 END DO !next IAD=IADF(I),IADL(I)
577
578 if (ixs(11,brick_list(nin,i)%id)==idb_id )print *, "TANGENT 1-12=", tangent(1:12)
579
580 DO j=1,12
581 IF(tangent(j)==1)THEN
582 k =(i-1)*12+j
583 nbcut = edge_list(nin,k)%NBCUT
584 !IF(NBCUT==0)THEN
585 ! EDGE_LIST(NIN,K)%CUTSHELL(1) = 0
586 ! EDGE_LIST(NIN,K)%CUTSHELL(2) = 0
587 ! EDGE_LIST(NIN,K)%NBCUT = 0
588 ! EDGE_LIST(NIN,K)%CUTCOOR(1) = 0
589 ! EDGE_LIST(NIN,K)%CUTCOOR(2) = 0
590 !ELSEIF(NBCUT==1)THEN
591 !
592 ! +-------------+ +-------------+
593 ! | | | |
594 ! | | | |
595 ! | o-----------------1----o dble -> O O
596 ! | | | | |
597 ! | | | | |
598 ! +-|-----------+ +-O-----------+
599 ! |
600 ! o
601 ! face1 cuts edge within tolerance, intersection point must be doubled to get PENTA HEXA, otherwise, max(secID) = HEXA
602 !
603 !cleaning if cutcoor = EM06 ou CUTCOOR = ONE-EM06 (CASE : lagrangian face == euler face => HEXA cutting)
604 IF(nbcut>=1)THEN
605 cutcoor = edge_list(nin,k)%CUTCOOR(1)
606 IF(cutcoor==em06 .OR. cutcoor==one-em06)THEN
607 edge_list(nin,k)%CUTSHELL(1) = edge_list(nin,k)%CUTSHELL(2)
608 edge_list(nin,k)%NBCUT = max(0,edge_list(nin,k)%NBCUT-1)
609 edge_list(nin,k)%CUTCOOR(1) = edge_list(nin,k)%CUTCOOR(2)
610 nbcut = edge_list(nin,k)%NBCUT
611 IF(nbcut==1)THEN
612 edge_list(nin,k)%CUTSHELL(2) = 0
613 edge_list(nin,k)%CUTCOOR(2) = zero
614 ENDIF
615 ENDIF
616 ENDIF
617 !duplicated nodes
618 nbcut = edge_list(nin,k)%NBCUT
619 cutcoor = edge_list(nin,k)%CUTCOOR(1)
620 IF(nbcut==1 .AND. cutcoor>em06 .AND. cutcoor<one-em06)THEN
621 edge_list(nin,k)%CUTSHELL(2) = edge_list(nin,k)%CUTSHELL(1)
622 edge_list(nin,k)%NBCUT = 2
623 edge_list(nin,k)%CUTCOOR(2) = edge_list(nin,k)%CUTCOOR(1)
624 ENDIF
625
626 !ENDIF
627 ENDIF
628 ENDDO
629
630 END DO !next I=NBF,NBL
631
632C=======================================================================
633C 4 Calcalculation of intersection points
634C + HM script (*createnode) to visualize them with HM
635C================================================================
636
637 CALL my_barrier !wait for all nodes to be computed
638
639 if(itask==0 .AND. ibug22_outp_intpoint==1)THEN
640
641 if(debug_outp)then
642 print *, " ===== intersection_nodes.txt ======="
643 endif
644 ipa = 100+ispmd
645 filename(1:12) = "cut_nod0.txt"
646 write(filename(8:8),'(i1.1)')ispmd+1
647 !===script HM=====!
648 open( unit=ipa, file = filename(1:12) )
649
650 !===script HM=====!
651 som=0
652 DO i=1,nb
653 !===script HM=====!
654 if(ibug22_intersect==ixs(11,brick_list(nin,i)%id) .or. ibug22_intersect==-1 .or. ibug22_outp_intpoint == 1)then
655 write (unit=ipa,fmt='(A,I10)') "cell ID = ", ixs(11,brick_list(nin,i)%id)
656 write (* ,fmt='(A,I10)') "cell ID = ", ixs(11,brick_list(nin,i)%id)
657 endif
658 DO j=1,12
659 iad = (i-1)*12+j
660 nbcut = edge_list(nin,iad)%NBCUT
661 cut(1:2) = edge_list(nin,iad)%CUTCOOR(1:2)
662 in(1:2) = edge_list(nin,iad)%NODE(1:2)
663 n1n2(1:3) = edge_list(nin,iad)%VECTOR(1:3)
664 on1(1:3) = x(1:3,in(1))
665 DO k=1,nbcut
666 !intersection coordinates
667 point(1:3)= on1(1:3) + cut(k) * n1n2(1:3)
668 ! EDGE_LIST(NIN,IAD)%CUTPOINT(1:3) = POINT(1:3)
669 !===script HM=====!
670 if(ibug22_intersect==ixs(11,brick_list(nin,i)%id) .or. ibug22_intersect==-1 .or. ibug22_outp_intpoint == 1)then
671 write(unit=ipa,
672 . fmt='(A12,F20.14,A1,F20.14,A1,F20.14,A13)')"*createnode ",point(1) ," ", point(2)," ",point(3)," 0 0 0 "
673 write(*,
674 . fmt='(A12,F20.14,A1,F20.14,A1,F20.14,A13)')"*createnode ",point(1) ," ", point(2)," ",point(3)," 0 0 0 "
675 endif
676 !===script HM=====!
677
678 db_flag=.true.
679 if (som>0)then
680 do l=1,som
681 if(sum(abs(point(1:3)-debugtab(l,1:3)))<em06)
682 . db_flag=.false.
683 end do
684 end if
685 if(db_flag)then
686 som=som+1 !debug
687 debugtab(som,1:3) =point(1:3)
688 end if
689
690 END DO ! (DO K=1,NBCUT <=> NBCUT>0)
691 END DO !next J=1,12
692 END do! next I=1,NB
693 !===script HM=====!
694 CLOSE(ipa)
695 !===script HM=====!
696 end if !ITASK==0
697
698! if(debug_outp)then
699! print *, ""
700! print *, " |--------i22intersect.F---------|"
701! print *, " | WRITE : SCRIPT HM NODES |"
702! print *, " |-------------------------------|"
703! print *, SOM , "nodes written"
704! do L=1,SOM
705! print *, debugTAB(L,1:3)
706! end do
707! print *, ""
708! endif
709
710C=======================================================================
711C 5 ***
712C================================================================
713
714 if(debug_outp)then
715 if(itask==0)then
716 print *, ""
717 print *, " |--------i22intersect.F---------|"
718 print *, " | EDGES |"
719 print *, " |-------------------------------|"
720 print *, 12*nb , " edges (12*NBRIC)"
721 DO i=1, nb ! 1,nb split across different threads
722 if(ibug22_intersect/=-1 .and. ibug22_intersect/=ixs(11,brick_list(nin,i)%id))cycle
723 print *, " ** CELL **", ixs(11,brick_list(nin,i)%ID)
724 DO j=1,12
725 k=(i-1)*12+j
726 IF( edge_list(nin,k)%NBCUT==0)THEN
727 WRITE(*,fmt='(A10,I10,A1,I12,I12,A8)') " edge ",k,":",
728 . itab(edge_list(nin,k)%NODE(1)),
729 . itab(edge_list(nin,k)%NODE(2))," "
730 ELSEIF( edge_list(nin,k)%NBCUT==1)THEN
731 WRITE(*,fmt='(A10,I10,A1,I12,I12,A8,1F30.16)') " edge ",k,":",
732 . itab(edge_list(nin,k)%NODE(1)),
733 . itab(edge_list(nin,k)%NODE(2))," CUTTED :" ,edge_list(nin,k)%CUTCOOR(1)
734 ELSE
735 WRITE(*,fmt='(A10,I10,A1,I12,I12,A8,2F30.16)') " edge ",k,":",
736 . itab(edge_list(nin,k)%NODE(1)),
737 . itab(edge_list(nin,k)%NODE(2))," 2CUTTED :" ,edge_list(nin,k)%CUTCOOR(1:2)
738 END IF
739 END DO
740 END DO
741 end if
742 end if
743
744
745C=======================================================================
746C 5 DEALLOCATE
747C================================================================
748
749 IF (itask==0)THEN
750 DEALLOCATE(vne)
751 DEALLOCATE(basisconst)
752 DEALLOCATE(nbsubtriangles)
753 DEALLOCATE(ptz)
754 DEALLOCATE(vza)
755 DEALLOCATE(vzb)
756 DEALLOCATE(pta)
757 DEALLOCATE(diag22)
758 END IF
759
760
761
762 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
isonsh3n
integer function isonsh3n(z, a, za, zb, p, iflg_db)
Definition i22intersect.F:777
ymax
subroutine ymax(idn, fac, npc, pld, stiffmin, stiffmax, stiffini, stiffavg)
Definition law100_upd.F:274
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::edge_list
type(edge_entity), dimension(:,:), allocatable, target edge_list
Definition cut-cell-search_mod.F:106
i22bufbric_mod::brick_list
type(brick_entity), dimension(:,:), allocatable, target brick_list
Definition cut-cell-search_mod.F:105
i22bufbric_mod::int22_buf
type(int22_buf_t) int22_buf
Definition cut-cell-search_mod.F:130
i22edge_mod
Definition cut-cell-buffer_mod.F:48
i22tri_mod
Definition cut-cell-search_mod.F:212
i22tri_mod::nb
integer nb
Definition cut-cell-search_mod.F:487
i22tri_mod::list_e
integer, dimension(:), allocatable list_e
Definition cut-cell-search_mod.F:495
i22tri_mod::iadf
integer, dimension(:), allocatable iadf
Definition cut-cell-search_mod.F:495
i22tri_mod::get_list_e_pos_from_cand_e_pos
integer, dimension(:), allocatable get_list_e_pos_from_cand_e_pos
Definition cut-cell-search_mod.F:495
i22tri_mod::ibug22_outp_intpoint
integer ibug22_outp_intpoint
Definition cut-cell-search_mod.F:367
i22tri_mod::iadl
integer, dimension(:), allocatable iadl
Definition cut-cell-search_mod.F:495
i22tri_mod::ibug22_intersect
integer ibug22_intersect
Definition cut-cell-search_mod.F:241
i22tri_mod::diag22
integer, dimension(:), allocatable diag22
Definition cut-cell-search_mod.F:502
i22tri_mod::nbsubtriangles
integer, dimension(:), allocatable nbsubtriangles
Definition cut-cell-search_mod.F:510
i22tri_mod::list_b
integer, dimension(:), allocatable list_b
Definition cut-cell-search_mod.F:495
i22tri_mod::ncande
integer ncande
Definition cut-cell-search_mod.F:499
i22tri_mod::ncandb
integer ncandb
Definition cut-cell-search_mod.F:499
realloc_mod
Definition realloc_mod.F:29
arret
subroutine arret(nn)
Definition arret.F:86
my_barrier
subroutine my_barrier
Definition machine.F:31

◆ isonsh3n()

integer function isonsh3n ( dimension(3), intent(in) z,
dimension(3), intent(in) a,
dimension(3), intent(in) za,
dimension(3), intent(in) zb,
dimension(3), intent(in) p,
integer, intent(in) iflg_db )

Definition at line 776 of file i22intersect.F.

777C-----------------------------------------------
778C P r e c o n d i t i o n s
779C-----------------------------------------------
780! Precondition : Point P is on the infinite plane generated by 3-node-shell (Z,A,B)
781C-----------------------------------------------
782C D e s c r i p t i o n
783C-----------------------------------------------
784!Is point in interior or exterior ?
785! 1 :interior
786! 0 :exterior
787!
788! Z+-----+B CRITERIA :
789! | P / ----------
790! | + / (ZA^ZP).(ZP^ZB) & (AZ^AP).(AP^AB) > 0
791! | / <=> / (ZA^U ).( U^ZB) & (AZ^V ).( V^AB) > 0
792! | / | U=ZP
793! |/ \ V=AP
794! A+
795!
796! CRITERION
797! coordinates of P in frame (ZA,ZB) must be
798! -TOL < x < 1+ TOL
799! -TOL < y < 1+ TOL
800C-----------------------------------------------
801C I m p l i c i t T y p e s
802C-----------------------------------------------
803#include "implicit_f.inc"
804C-----------------------------------------------
805C D u m m y A r g u m e n t s
806C-----------------------------------------------
807 my_real , intent(in) ::
808 . z(3), a(3), p(3) !points
809 my_real , intent(in) ::
810 . za(3),zb(3) !vectors
811 INTEGER , intent(in) :: IFLG_DB
812C-----------------------------------------------
813C L o c a l V a r i a b l e s
814C-----------------------------------------------
815 my_real
816 . u(3),v(3),ab(3),pv1(3),pv2(3),ps1,ps2,ps3,u2(3),zp(3),ap(3),
817 . pz(3),pa(3),pb(3),
818 . norm, k , norm1, norm2, eps1,eps2, lref1,lref2 , tol,
819 . norm_pz,norm_pa,norm_pb,norm_za_2,norm_zb_2,sin1,sin2,tol_sin,norm_ab,
820 . pab, pbz, paz, abz,
821 . coor1, coor2, coor3,coef,
822 . t, s, tolcrit
823C-----------------------------------------------
824C S o u r c e L i n e s
825C-----------------------------------------------
826 isonsh3n = 0
827
828 tolcrit = em06
829 tol = tolcrit
830
831 norm_za_2 = (za(1)*za(1)+za(2)*za(2)+za(3)*za(3))
832 norm_zb_2 = (zb(1)*zb(1)+zb(2)*zb(2)+zb(3)*zb(3))
833
834 zp(1) = p(1) - z(1)
835 zp(2) = p(2) - z(2)
836 zp(3) = p(3) - z(3)
837
838 ps1 = za(1)*zp(1)+za(2)*zp(2)+za(3)*zp(3)
839
840 ps2 = zb(1)*zp(1)+zb(2)*zp(2)+zb(3)*zp(3)
841
842 ps3 = zb(1)*za(1)+zb(2)*za(2)+zb(3)*za(3)
843
844 coef = one-ps3*ps3/norm_za_2/norm_zb_2
845 coef = one / coef
846
847 t = coef * ( ps2/norm_zb_2 - ps3*ps1/norm_za_2/norm_zb_2)
848 s = coef * ( -ps3*ps2/norm_za_2/norm_zb_2 + ps1/norm_za_2)
849
850 IF(iflg_db == 1)THEN
851 print *, "coor ZA,ZB =", t,s
852 write(*,fmt='(A12,3F30.16)')"*createnode ",z(1:3)
853 write(*,fmt='(A12,3F30.16)')"*createnode ",a(1:3)
854 write(*,fmt='(A12,3F30.16)')"*createnode ",z(1:3)+zb(1:3)
855 write(*,fmt='(A12,3F30.16)')"*createnode ",p(1:3)
856 ENDIF
857
858 IF(t>=-tol)THEN
859 IF(s>=-tol )THEN
860 IF(s+t<=one+tol)THEN
861 isonsh3n = 1
862 ENDIF
863 ENDIF
864 ENDIF
865
866 RETURN
867