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!|| i22bufbric_mod ../common_source/modules/interfaces/cut-cell-search_mod.F
33!|| i22edge_mod ../common_source/modules/interfaces/cut-cell-buffer_mod.F
34!|| i22tri_mod ../common_source/modules/interfaces/cut-cell-search_mod.F
35!|| realloc_mod ../engine/share/modules/realloc_mod.F
36!||====================================================================
37 SUBROUTINE i22intersect(
38 1 X ,II_STOK, CAND_B ,CAND_E ,ITASK,
39 2 NBRIC ,ITAB , BUFBRIC ,NCAND,
40 3 IXS ,NIN )
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.AND. !if (ixs(11,brick_list(nin,i)%id)==5882 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.AND. IF((ABS(POW(1))<=TOL)(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.AND..OR..AND. IF( ((POW(1)<-TOL)(POW(2)>TOL)) ((POW(1)>TOL)(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.AND. IF((CUTCOOR<=ONE+TOL)(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.AND. ELSEIF((ABS(POW(1))<=TOL)(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.OR..OR..OR. !IF(J==1 J==5 J==8 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.OR..OR..OR. !IF(J==1 J==5 J==8 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.OR..OR..OR. !IF(J==1 J==5 J==8 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.AND. if(itask==0 debug_outp)then
548.or. if(ibug22_intersect==-1 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.AND. IF(NBCUT2>0 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.AND. if(itask==0 debug_outp)then
563.or. if(ibug22_intersect==-1 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.OR. IF(CUTCOOR==EM06 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.AND..AND. IF(NBCUT==1 CUTCOOR>EM06 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.AND. if(ITASK==0 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.or..or. if(ibug22_intersect==ixs(11,brick_list(nin,i)%id) ibug22_intersect==-1 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.or..or. if(ibug22_intersect==ixs(11,brick_list(nin,i)%id) ibug22_intersect==-1 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
761 END
762
763
764
765
766
767
768!||====================================================================
769!|| isonsh3n ../engine/source/interfaces/int22/i22intersect.F
770!||--- called by ------------------------------------------------------
771!|| i22intersect ../engine/source/interfaces/int22/i22intersect.F
772!|| i22trivox ../engine/source/interfaces/intsort/i22trivox.F
773!||====================================================================
774 INTEGER FUNCTION ISONSH3N(Z, A, ZA, ZB, P, IFLG_DB)
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
866 END FUNCTION
867
my_real
#define my_real
Definition cppsort.cpp:32
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
i22intersect
subroutine i22intersect(x, ii_stok, cand_b, cand_e, itask, nbric, itab, bufbric, ncand, ixs, nin)
Definition i22intersect.F:41
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
id
initmumps id
Definition multiplerhs_example.m:4
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::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
second
real function second()
SECOND Using ETIME
Definition second_EXT_ETIME.f:35
my_barrier
subroutine my_barrier
Definition machine.F:31
i22edge_mod::list_secnd
Definition cut-cell-buffer_mod.F:99