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 37 of file i22intersect.F.

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

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