OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i22intersect.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| i22intersect ../engine/source/interfaces/int22/i22intersect.F
25!||--- called by ------------------------------------------------------
26!|| i22buce ../engine/source/interfaces/intsort/i22buce.F
27!||--- calls -----------------------------------------------------
28!|| arret ../engine/source/system/arret.F
29!|| isonsh3n ../engine/source/interfaces/int22/i22intersect.F
30!|| my_barrier ../engine/source/system/machine.F
31!||--- uses -----------------------------------------------------
32!|| element_mod ../common_source/modules/elements/element_mod.F90
33!|| i22bufbric_mod ../common_source/modules/interfaces/cut-cell-search_mod.F
34!|| i22edge_mod ../common_source/modules/interfaces/cut-cell-buffer_mod.F
35!|| i22tri_mod ../common_source/modules/interfaces/cut-cell-search_mod.F
36!|| realloc_mod ../engine/share/modules/realloc_mod.F
37!||====================================================================
38 SUBROUTINE i22intersect(
39 1 X ,II_STOK, CAND_B ,CAND_E ,ITASK,
40 2 NBRIC ,ITAB , BUFBRIC ,NCAND,
41 3 IXS ,NIN )
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
763 END
764
765
766
767
768
769
770!||====================================================================
771!|| isonsh3n ../engine/source/interfaces/int22/i22intersect.F
772!||--- called by ------------------------------------------------------
773!|| i22intersect ../engine/source/interfaces/int22/i22intersect.F
774!|| i22trivox ../engine/source/interfaces/intsort/i22trivox.F
775!||====================================================================
776 INTEGER FUNCTION isonsh3n(Z, A, ZA, ZB, P, IFLG_DB)
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
868 END FUNCTION
869
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
i22intersect
subroutine i22intersect(x, ii_stok, cand_b, cand_e, itask, nbric, itab, bufbric, ncand, ixs, nin)
Definition i22intersect.F:42
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
i22edge_mod::list_secnd
Definition cut-cell-buffer_mod.F:99