OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i25trivox_edg.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!|| i25trivox_edg ../engine/source/interfaces/intsort/i25trivox_edg.F
25!||--- called by ------------------------------------------------------
26!|| i25buce_edg ../engine/source/interfaces/intsort/i25buce_edg.F
27!||--- calls -----------------------------------------------------
28!|| bitget ../engine/source/interfaces/intsort/i20sto.F
29!|| i25sto_e2s ../engine/source/interfaces/intsort/i25sto_e2s.F
30!|| i25sto_edg ../engine/source/interfaces/intsort/i25sto_edg.F
31!|| int_checksum ../engine/share/modules/debug_mod.F
32!|| ireallocate ../engine/share/modules/realloc_mod.F
33!|| my_barrier ../engine/source/system/machine.F
34!||--- uses -----------------------------------------------------
35!|| debug_mod ../engine/share/modules/debug_mod.F
36!|| realloc_mod ../engine/share/modules/realloc_mod.F
37!|| tri11 ../engine/share/modules/tri11_mod.F
38!|| tri25ebox ../engine/share/modules/tri25ebox.F
39!|| tri7box ../engine/share/modules/tri7box.F
40!||====================================================================
41 SUBROUTINE i25trivox_edg(
42 1 I_MEM ,VMAXDT ,INACTI ,IRECT ,
43 2 X ,V ,STF ,STFE ,XYZM ,
44 3 II_STOK ,CANDS_E2E ,ESHIFT ,NEDGE_T ,CANDM_E2E ,
45 4 MULNSNE,NOINT ,BGAPEMX ,SSHIFT ,NRTM_T ,
46 5 VOXEL ,NBX ,NBY ,NBZ ,
47 6 IGAP ,GAP_M ,GAP_M_L ,DRAD ,MARGE ,
48 7 ITASK ,ITAB ,LL_STOK ,MULNSNS ,
49 8 MBINFLG ,EBINFLG,ILEV ,CAND_A ,CAND_P ,
50 9 FLAGREMNODE,KREMNODE_EDG,REMNODE_EDG,KREMNODE_E2S,
51 . REMNODE_E2S,
52 A IEDGE ,NEDGE ,LEDGE ,MSEGTYP ,IGAP0 ,
53 B ADMSR,EDG_BISECTOR,VTX_BISECTOR,
54 C CANDM_E2S,CANDS_E2S,CAND_B,CAND_PS,GAPE ,
55 D GAP_E_L,NEDGE_LOCAL,IFQ,CANDE2E_FX ,CANDE2E_FY,
56 E CANDE2E_FZ,CANDE2S_FX ,CANDE2S_FY,CANDE2S_FZ,IFPEN_E,IFPEN_E2S,
57 F KREMNODE_EDG_SIZ,REMNODE_EDG_SIZ,KREMNODE_E2S_SIZ,REMNODE_E2S_SIZ,
58 G DGAPLOAD )
59C============================================================================
60C M o d u l e s
61C-----------------------------------------------
62 USE realloc_mod
63 USE tri25ebox
64 USE tri7box
65 USE tri11
66#ifdef WITH_ASSERT
67 USE debug_mod
68#endif
69C-----------------------------------------------
70C I m p l i c i t T y p e s
71C-----------------------------------------------
72#include "implicit_f.inc"
73#include "comlock.inc"
74C-----------------------------------------------
75C G l o b a l P a r a m e t e r s
76C-----------------------------------------------
77#include "mvsiz_p.inc"
78c parameter setting the size for the vector (orig version is 128)
79 INTEGER NVECSZ
80 parameter(nvecsz = mvsiz)
81C-----------------------------------------------
82C C o m m o n B l o c k s
83C-----------------------------------------------
84#include "param_c.inc"
85#include "assert.inc"
86#include "i25edge_c.inc"
87C-----------------------------------------------
88C ROLE DE LA ROUTINE:
89C ===================
90C CLASSE LES EDGES DANS DES BOITES
91C RECHERCHE POUR CHAQUE FACETTE DES BOITES CONCERNES
92C RECHERCHE DES CANDIDATS
93C-----------------------------------------------
94C D u m m y A r g u m e n t s
95C-----------------------------------------------
96 INTEGER I_MEM(2),INACTI,ITASK,IGAP,IEDGE,NEDGE,ESHIFT,NEDGE_T,SSHIFT,NRTM_T,IGAP0,
97 . mulnsne,mulnsns,noint,nbx,nby,nbz,ifq,
98 . cands_e2e(*),candm_e2e(*),
99 . irect(4,*), voxel(nbx+2,nby+2,nbz+2),ii_stok,ll_stok,itab(*),
100 . mbinflg(*),ebinflg(*),ilev,cand_a(*),ledge(nledge,*),admsr(4,*),msegtyp(*),
101 . candm_e2s(*),cands_e2s(*),cand_b(*),ifpen_e(*),ifpen_e2s(*)
102C INTEGER :: NEDGE_REMOTE_OLD, RENUM(*)
103 INTEGER , INTENT(IN) :: KREMNODE_EDG_SIZ,REMNODE_EDG_SIZ,KREMNODE_E2S_SIZ,REMNODE_E2S_SIZ,
104 . FLAGREMNODE, KREMNODE_EDG(KREMNODE_EDG_SIZ), REMNODE_EDG(REMNODE_EDG_SIZ),
105 . KREMNODE_E2S(KREMNODE_E2S_SIZ), REMNODE_E2S(REMNODE_E2S_SIZ)
106C REAL
107 my_real , INTENT(IN) :: DGAPLOAD ,DRAD
108 my_real
109 . X(3,*),V(3,*),XYZM(6),STF(*), STFE(NEDGE), GAP_M(*), GAP_M_L(*), GAPE(*), GAP_E_L(*),
110 . CAND_P(*),CAND_PS(*),MARGE,BGAPEMX,VMAXDT,
111 . CANDE2E_FX(*) ,CANDE2E_FY(*),CANDE2E_FZ(*),
112 . CANDE2S_FX(4,*) ,CANDE2S_FY(4,*),CANDE2S_FZ(4,*)
113 real*4 edg_bisector(3,4,*), vtx_bisector(3,2,*)
114 INTEGER, INTENT(IN) :: NEDGE_LOCAL
115C-----------------------------------------------
116C L o c a l V a r i a b l e s
117C-----------------------------------------------
118 INTEGER I,J,I_STOK, SOL_EDGE, SH_EDGE,
119 . N1,N2,NN,NE,K,L,J_STOK,II,JJ,NA,NB,
120 . prov_s(mvsiz),prov_m(mvsiz),
121 . m,ns1,ns2,nse,ns,SIZE,z_first,z_last
122C REAL
123 my_real
124 . dx,dy,dz,xs,ys,zs,xx,sx,sy,sz,s2,
125 . xmin, xmax,ymin, ymax,zmin, zmax, tz, gapsmx, gapl,
126 . xx1,xx2,xx3,xx4,yy1,yy2,yy3,yy4,zz1,zz2,zz3,zz4,
127 . d1x,d1y,d1z,d2x,d2y,d2z,dd1,dd2,d2,a2,gs,drad2
128c provisoire
129 INTEGER IX,IY,IZ,IEDG,IE,
130 . m1, m2, m3, m4, mm1,mm2,mm3,mm4,ss1,ss2,
131 . ims1,ims2,iss1,iss2,
132 . am1,am2,as1,as2,
133 . ix1,iy1,iz1,ix2,iy2,iz2,remove_remote
134 INTEGER, DIMENSION(3) :: TMIN,TMAX
135 my_real
136 . xminb,yminb,zminb,xmaxb,ymaxb,zmaxb,aaa,
137 . xmax_edgs, xmin_edgs, !cotes min/max des aretes seconds et mains
138 . ymax_edgs, ymin_edgs,
139 . zmax_edgs, zmin_edgs,
140 . xmax_edgm, xmin_edgm,
141 . ymax_edgm, ymin_edgm,
142 . zmax_edgm, zmin_edgm
143 my_real :: g ! gap
144 INTEGER, DIMENSION(:), ALLOCATABLE :: TAGEDG
145 INTEGER :: EDGE_TYPE
146 INTEGER :: EID
147 INTEGER FIRST_ADD, PREV_ADD, CHAIN_ADD, CURRENT_ADD, MAX_ADD
148 INTEGER BITGET
149 EXTERNAL BITGET
150
151C-----------------------------------------------
152 INTEGER IDS(4), PROV_IDS(2,MVSIZ)
153
154C-----------------------------------------------
155 INTEGER, DIMENSION(:), ALLOCATABLE :: TAGREMLINE
156C-----------------------------------------------
157C
158 DRAD2 =zero
159
160 IF(flagremnode==2) THEN
161 ALLOCATE(tagremline(nedge))
162 tagremline(1:nedge) = 0
163 ENDIF
164
165 ids(1:4) = 0
166 prov_ids(1:2,1:mvsiz) = 0
167
168 sol_edge =iedge/10 ! solids
169 sh_edge =iedge-10*sol_edge ! shells
170
171 min_ix=nbx+2
172 min_iy=nby+2
173 min_iz=nbz+2
174 max_ix=1
175 max_iy=1
176 max_iz=1
177
178 !---------------------------------------------------------!
179 ! Allocation des tableaux chaines !
180 !---------------------------------------------------------!
181 max_add = max(1,4*(nedge+nedge_remote))
182 IF(itask==0)THEN
183 ALLOCATE(lchain_elem(1:max_add))
184 ALLOCATE(lchain_next(1:max_add))
185 ALLOCATE(lchain_last(1:max_add))
186 END IF
187
188C Barrier to wait init voxel and allocation
189 CALL my_barrier
190C
191 xmin = xyzm(1)
192 ymin = xyzm(2)
193 zmin = xyzm(3)
194 xmax = xyzm(4)
195 ymax = xyzm(5)
196 zmax = xyzm(6)
197
198c dev future: xminb plus grand que xmin...
199 xminb = xmin
200 yminb = ymin
201 zminb = zmin
202 xmaxb = xmax
203 ymaxb = ymax
204 zmaxb = zmax
205
206c IF( NSPMD > 1) THEN
207c CALL SPMD_OLDNUMCD(RENUM,OLDNUM,ISZNSNR,NSNROLD)
208c ENDIF
209C=======================================================================
210C 1 mise des edges dans les boites
211C=======================================================================
212 IF(itask == 0)THEN
213
214 current_add=1 ! premiere adresse
215 DO i=1,nedge + nedge_remote
216
217 IF(i <= nedge_local) THEN
218 ne = ledge(1,i)
219C PRINTIF((STFE(I) == 0 .AND.LEDGE(6,I) == D_ES),STFE(I))
220 IF(stfe(i)==zero) cycle ! on ne retient pas les facettes detruites
221
222 IF(ledge(7,i) < 0) cycle ! larete nest pas une arete second
223 n1 = ledge(5,i)
224 n2 = ledge(6,i)
225 eid = ledge(8,i)
226
227 xx1=x(1,n1)
228 xx2=x(1,n2)
229 yy1=x(2,n1)
230 yy2=x(2,n2)
231 zz1=x(3,n1)
232 zz2=x(3,n2)
233 debug_e2e(eid == d_es,eid)
234 ELSE IF(i > nedge) THEN
235 xx1=xrem_edge(e_x1,i-nedge)
236 xx2=xrem_edge(e_x2,i-nedge)
237 yy1=xrem_edge(e_y1,i-nedge)
238 yy2=xrem_edge(e_y2,i-nedge)
239 zz1=xrem_edge(e_z1,i-nedge)
240 zz2=xrem_edge(e_z2,i-nedge)
241 eid = irem_edge(e_global_id,i-nedge)
242 debug_e2e(eid == d_es,eid)
243 ELSE
244 ! Secondary edge is boundary between domains
245 ! ISPMD is not the owner of this edge
246 assert(nspmd > 1)
247 cycle
248 ENDIF
249 debug_e2e(eid==d_es,igap0)
250
251 IF(igap0 == 0)THEN
252 xmax_edgs=max(xx1,xx2);
253 xmin_edgs=min(xx1,xx2);
254 ymax_edgs=max(yy1,yy2);
255 ymin_edgs=min(yy1,yy2);
256 zmax_edgs=max(zz1,zz2);
257 zmin_edgs=min(zz1,zz2);
258 debug_e2e(eid==d_es,xmin_edgs)
259 debug_e2e(eid==d_es,ymin_edgs)
260 debug_e2e(eid==d_es,zmin_edgs)
261 debug_e2e(eid==d_es,xmax_edgs)
262 debug_e2e(eid==d_es,ymax_edgs)
263 debug_e2e(eid==d_es,zmax_edgs)
264 debug_e2e(eid==d_es,xmin)
265 debug_e2e(eid==d_es,ymin)
266 debug_e2e(eid==d_es,zmin)
267 debug_e2e(eid==d_es,xmax)
268 debug_e2e(eid==d_es,ymax)
269 debug_e2e(eid==d_es,zmax)
270 IF(xmax_edgs < xmin) cycle
271 IF(xmin_edgs > xmax) cycle
272 IF(ymax_edgs < ymin) cycle
273 IF(ymin_edgs > ymax) cycle
274 IF(zmax_edgs < zmin) cycle
275 IF(zmin_edgs > zmax) cycle
276
277 ELSE
278 IF(i <= nedge) THEN
279 g = gape(i)
280 ELSE
281 g = xrem_edge(e_gap,i-nedge)
282 END IF
283
284
285 xmax_edgs=max(xx1,xx2)+g;
286 xmin_edgs=min(xx1,xx2)-g;
287 ymax_edgs=max(yy1,yy2)+g;
288 ymin_edgs=min(yy1,yy2)-g;
289 zmax_edgs=max(zz1,zz2)+g;
290 zmin_edgs=min(zz1,zz2)-g;
291
292
293 debug_e2e(eid==d_es,xmin_edgs)
294 debug_e2e(eid==d_es,ymin_edgs)
295 debug_e2e(eid==d_es,zmin_edgs)
296 debug_e2e(eid==d_es,xmax_edgs)
297 debug_e2e(eid==d_es,ymax_edgs)
298 debug_e2e(eid==d_es,zmax_edgs)
299
300
301 END IF
302
303
304
305 !-------------------------------------------!
306 ! VOXEL OCCUPIED BY THE EDGE !
307 !-------------------------------------------!
308 !Voxel_lower_left_bound for this edge
309 ix1=int(nbx*(xmin_edgs-xminb)/(xmaxb-xminb))
310 iy1=int(nby*(ymin_edgs-yminb)/(ymaxb-yminb))
311 iz1=int(nbz*(zmin_edgs-zminb)/(zmaxb-zminb))
312 ix1=max(1,2+min(nbx,ix1))
313 iy1=max(1,2+min(nby,iy1))
314 iz1=max(1,2+min(nbz,iz1))
315 !Voxel_upper_right_bound for this edge
316 ix2=int(nbx*(xmax_edgs-xminb)/(xmaxb-xminb))
317 iy2=int(nby*(ymax_edgs-yminb)/(ymaxb-yminb))
318 iz2=int(nbz*(zmax_edgs-zminb)/(zmaxb-zminb))
319 ix2=max(1,2+min(nbx,ix2))
320 iy2=max(1,2+min(nby,iy2))
321 iz2=max(1,2+min(nbz,iz2))
322
323 !pour reset des voxel
324 min_ix = min(min_ix,ix1)
325 min_iy = min(min_iy,iy1)
326 min_iz = min(min_iz,iz1)
327 max_ix = max(max_ix,ix2)
328 max_iy = max(max_iy,iy2)
329 max_iz = max(max_iz,iz2)
330
331 !----------------------------------------------!
332 ! EDGE STORAGE FOR EACH VOXEL (CHAINED ARRAY) !
333 !----------------------------------------------!
334C
335C VOXEL(i,j,k) LCHAIN_LAST(FIRST)
336C +-----------+------------+
337C | =>FIRST | =>LAST |
338C +--+--------+--+---------+
339C | |
340C | |
341C | |
342C | | LCHAIN_ELEM(*) LCHAIN_NEXT(*)
343C | | +------------+-----------+
344C +-------------->| edge_id | iadd 3 | 1:FIRST --+
345C | +------------+-----------+ |
346C | | | | 2 |
347C | +------------+-----------+ |
348C | | edge_id | iadd 4 | 3 <-------+
349C | +------------+-----------+ |
350C | | edge_id | iadd 6 | 4 <-------+
351C | +------------+-----------+ |
352C | | | | 5 |
353C | +------------+-----------+ |
354C +-->| edge_id | 0 | 6:LAST <--+
355C +------------+-----------+
356C | | | MAX_ADD
357C +------------+-----------+
358 DO iz = iz1,iz2
359 DO iy = iy1,iy2
360 DO ix = ix1,ix2
361
362 first_add = voxel(ix,iy,iz)
363
364 IF(first_add == 0)THEN
365 !voxel encore vide
366 voxel(ix,iy,iz) = current_add ! adresse dans le tableau chaine de la premiere eddge trouvee occupant le voxel
367 lchain_last(current_add) = current_add ! dernier=adresse pour l edge courante
368 lchain_elem(current_add) = i ! edge ID
369 lchain_next(current_add) = 0 ! pas de suivant car dernier de la liste !
370 ELSE
371 !voxel contenant deja une edge
372 prev_add = lchain_last(first_add) ! devient l'avant-dernier
373 lchain_last(first_add) = current_add ! maj du dernier
374 lchain_elem(current_add) = i ! edge ID
375 lchain_next(prev_add) = current_add ! maj du suivant 0 -> CURRENT_ADD
376 lchain_next(current_add) = 0 ! pas de suivant car dernier de la liste
377 ENDIF
378
379 current_add = current_add+1
380
381 IF( current_add>=max_add)THEN
382 !OPTIMISATION : suprresion du deallocate/GOTO debut.
383 !REALLOCATE SI PAS ASSEZ DE PLACE : inutile de recommencer de 1 a MAX_ADD-1, on poursuit de MAX_ADD a 2*MAX_ADD
384 max_add = 2 * max_add
385 !print *, "reallocate"
386 lchain_next => ireallocate(lchain_next, max_add)
387 lchain_elem => ireallocate(lchain_elem, max_add)
388 lchain_last => ireallocate(lchain_last, max_add)
389 ENDIF
390
391 ENDDO !IX
392 ENDDO !IY
393 ENDDO !IZ
394
395 ENDDO
396
397 END IF
398C Barrier to wait task0 treatment
399 CALL my_barrier
400C
401! Attention: allocation en NTHREADS x (NEDGE+NEDGE_REMOTE)
402 ALLOCATE(tagedg(1:nedge+nedge_remote))
403 tagedg(1:nedge+nedge_remote)=0
404C=======================================================================
405C Sorting vs main shell edges
406C=======================================================================
407 IF(sh_edge==0) GOTO 300
408C=======================================================================
409C 3 A partir des voxels occupes par une edge main, on est en
410C mesure de connaitre toutes les edges escalves dans ce voisinage.
411C Ce qui permet de creer des couples cancidats pour le contact
412C Si la penetration est positive.
413C=======================================================================
414
415 j_stok = 0
416
417 DO i=1,nedge_t
418
419 iedg=eshift+i
420
421 IF(stfe(iedg)==zero) cycle ! on ne retient pas les facettes detruites
422 ne=ledge(1,iedg)
423
424 IF(iabs(ledge(7,iedg))==1) cycle ! Main solid edge
425
426 !-------------------------------------------!
427 ! (N1,N2) is the current main edge !
428 !-------------------------------------------!
429
430 aaa = marge+bgapemx+gape(iedg)+dgapload
431
432 n1 = ledge(5,iedg)
433 n2 = ledge(6,iedg)
434 mm1 = itab(n1)
435 mm2 = itab(n2)
436 am1 = min(mm1,mm2)
437 am2 = max(mm1,mm2)
438
439 IF(ilev==2)THEN
440 ims1 = bitget(ebinflg(iedg),0)
441 ims2 = bitget(ebinflg(iedg),1)
442 END IF
443
444 !-------------------------------------------!
445 ! X-coordinates of the four nodes !
446 !-------------------------------------------!
447
448 xx1=x(1,n1)
449 xx2=x(1,n2)
450 yy1=x(2,n1)
451 yy2=x(2,n2)
452 zz1=x(3,n1)
453 zz2=x(3,n2)
454 xmax_edgm=max(xx1,xx2)+gape(iedg) ! +TZINF
455 xmin_edgm=min(xx1,xx2)-gape(iedg) ! -TZINF
456 ymax_edgm=max(yy1,yy2)+gape(iedg) ! +TZINF
457 ymin_edgm=min(yy1,yy2)-gape(iedg) ! -TZINF
458 zmax_edgm=max(zz1,zz2)+gape(iedg) ! +TZINF
459 zmin_edgm=min(zz1,zz2)-gape(iedg) ! -TZINF
460 !-------------------------------------------!
461 ! VOXEL OCCUPIED BY THE BRICK !
462 !-------------------------------------------!
463 !Voxel_lower_left_bound for this element---+
464 ix1=int(nbx*(xmin_edgm-aaa-xminb)/(xmaxb-xminb))
465 iy1=int(nby*(ymin_edgm-aaa-yminb)/(ymaxb-yminb))
466 iz1=int(nbz*(zmin_edgm-aaa-zminb)/(zmaxb-zminb))
467 ix1=max(1,2+min(nbx,ix1))
468 iy1=max(1,2+min(nby,iy1))
469 iz1=max(1,2+min(nbz,iz1))
470 !Voxel_upper_right_bound for this element---+
471 ix2=int(nbx*(xmax_edgm+aaa-xminb)/(xmaxb-xminb))
472 iy2=int(nby*(ymax_edgm+aaa-yminb)/(ymaxb-yminb))
473 iz2=int(nbz*(zmax_edgm+aaa-zminb)/(zmaxb-zminb))
474 ix2=max(1,2+min(nbx,ix2))
475 iy2=max(1,2+min(nby,iy2))
476 iz2=max(1,2+min(nbz,iz2))
477
478C--- IREMGAP - tag of deactivated lines
479 IF(flagremnode==2)THEN
480 k = kremnode_edg(2*(iedg-1)+1)
481 l = kremnode_edg(2*(iedg-1)+2)-1
482 DO m=k,l
483 tagremline(remnode_edg(m)) = 1
484 ENDDO
485 ENDIF
486
487 DO iz = iz1,iz2
488 DO iy = iy1,iy2
489 DO ix = ix1,ix2
490
491 chain_add = voxel(ix,iy,iz) ! adresse dans le tableau chaine de la premiere edge stoquee dans le voxel
492 DO WHILE(chain_add /= 0) ! BOUCLE SUR LES EDGES DU VOXEL COURANT
493 jj = lchain_elem(chain_add) ! numeros des edge_id balayes dans le voxel courant
494
495 IF(tagedg(jj)/=0)THEN ! edge deja traitee vs cette arete main
496
497 chain_add = lchain_next(chain_add)
498 cycle
499 END IF
500 tagedg(jj)=1
501
502 !secnd edge nodes, exclure couples avec noeud commun
503 IF (jj<=nedge)THEN
504 ss1= itab(ledge(5,jj))
505 ss2= itab(ledge(6,jj))
506 eid = ledge(8,jj)
507 ELSE
508 ss1=irem_edge(e_node1_globid,jj-nedge)
509 ss2=irem_edge(e_node2_globid,jj-nedge)
510 eid = irem_edge(e_global_id,jj-nedge)
511 END IF
512
513 IF( (ss1==mm1).OR.(ss1==mm2).OR.
514 . (ss2==mm1).OR.(ss2==mm2) )THEN
515 chain_add = lchain_next(chain_add)
516 cycle
517 END IF
518
519 IF(ilev==2)THEN
520 IF(jj <= nedge) THEN
521 iss1=bitget(ebinflg(jj),0)
522 iss2=bitget(ebinflg(jj),1)
523 ELSE
524C double-check
525 iss1 = bitget(irem_edge(e_ebinflg,jj-nedge),0)
526 iss2 = bitget(irem_edge(e_ebinflg,jj-nedge),1)
527 ENDIF
528
529 IF(.NOT.((ims1 == 1 .and. iss2==1).or.
530 . (ims2 == 1 .and. iss1==1)))THEN
531 chain_add = lchain_next(chain_add)
532 cycle
533 ENDIF
534 ENDIF
535
536 IF( jj <= nedge) THEN
537 edge_type = ledge(7,jj)
538 ELSE
539 edge_type = irem_edge(e_type ,jj - nedge)
540 ENDIF
541
542 IF(iabs(ledge(7,iedg))/=1 .AND. edge_type /= 1 )THEN
543 ! attention les traitements de i25dst3e pour les solides
544 ! ne sont pas symetriques main second.
545 as1 = min(ss1,ss2)
546 as2 = max(ss1,ss2)
547 ! unicite des couples
548 IF(am1 < as1 .OR. (am1 == as1 .AND. am2 < as2))THEN
549 chain_add = lchain_next(chain_add)
550 cycle
551 ENDIF
552 ENDIF
553C IREMPGAP
554 IF (flagremnode == 2) THEN
555 IF (jj <= nedge) THEN
556C- Local Taged lines are removed
557 IF(tagremline(jj)==1) THEN
558 chain_add = lchain_next(chain_add)
559 cycle
560 ENDIF
561
562 IF(tagremline(jj)==0) THEN
563C- Even if it is not Remote lines have to be looked in remote list: Edge oin 2 procs
564 k = kremnode_edg(2*(iedg-1)+2)
565 l = kremnode_edg(2*(iedg-1)+3)-1
566 remove_remote = 0
567 DO m=k,l,2
568 IF ((ss1==remnode_edg(m)).AND.(ss2==remnode_edg(m+1))) remove_remote = 1
569 ENDDO
570 IF (remove_remote==1) THEN
571 chain_add = lchain_next(chain_add)
572 cycle
573 ENDIF
574 ENDIF
575 ELSE
576C- Remote lines are identified by nodes
577 k = kremnode_edg(2*(iedg-1)+2)
578 l = kremnode_edg(2*(iedg-1)+3)-1
579 remove_remote = 0
580 DO m=k,l,2
581 IF ((ss1==remnode_edg(m)).AND.(ss2==remnode_edg(m+1))) remove_remote = 1
582 ENDDO
583 IF (remove_remote==1) THEN
584 chain_add = lchain_next(chain_add)
585 cycle
586 ENDIF
587 ENDIF
588 ENDIF
589
590 j_stok = j_stok + 1 !on dispose d'un candidat
591 assert(jj > 0)
592 assert(jj <= nedge + nedge_remote)
593 prov_s(j_stok) = jj !edge secnd
594 prov_m(j_stok) = iedg !edge main
595
596 debug_e2e(ledge(8,iedg) == d_em .AND. eid == d_es,eid)
597
598c IF(DEJA==0) NEDG = NEDG + 1 !nombre d edges candidate au calcul de contact (debug)
599c DEJA=1 !l edge main IEDG fait l'objet d'une ecriture de candidat. On compte les edges main faisant l'objet dun couple candidate : on ne doit plus incrementer NEDG pour les autres edge secnd testees.
600 chain_add = lchain_next(chain_add)
601C-----------------------------------------------------
602 IF(j_stok==nvsiz)THEN
603 CALL i25sto_edg(
604 1 nvsiz ,irect ,x ,ii_stok,inacti,
605 2 cands_e2e,candm_e2e ,mulnsne,noint ,marge ,
606 3 i_mem(1) ,prov_s ,prov_m ,igap0,cand_a,
607 4 nedge ,ledge ,itab ,drad2 ,igap ,
608 5 gape ,gap_e_l,admsr ,edg_bisector,vtx_bisector ,
609 6 cand_p,ifq,cande2e_fx ,cande2e_fy,cande2e_fz,ifpen_e,
610 7 dgapload)
611 IF(i_mem(1)/=0) GOTO 300
612 j_stok = 0
613 ENDIF
614C-----------------------------------------------------
615
616 ENDDO ! WHILE(CHAIN_ADD /= 0)
617
618 ENDDO !NEXT IZ
619 ENDDO !NEXT IY
620 ENDDO !NEXT IZ
621
622C Reset TAGEDG
623 DO iz = iz1,iz2
624 DO iy = iy1,iy2
625 DO ix = ix1,ix2
626
627 chain_add = voxel(ix,iy,iz)
628 DO WHILE(chain_add /= 0) ! BOUCLE SUR LES EDGES DU VOXEL COURANT
629
630 jj = lchain_elem(chain_add) ! numeros des edge_id balayes dans le voxel courant
631 tagedg(jj)=0
632
633 chain_add = lchain_next(chain_add)
634
635 END DO
636
637 ENDDO !NEXT IZ
638 ENDDO !NEXT IY
639 ENDDO !NEXT IZ
640
641
642C--- IREMGAP - clean of tagremline
643 IF(flagremnode==2)THEN
644 k = kremnode_edg(2*(iedg-1)+1)
645 l = kremnode_edg(2*(iedg-1)+2)-1
646 DO m=k,l
647 tagremline(remnode_edg(m)) = 0
648 ENDDO
649 ENDIF
650
651 ENDDO !NEXT IEDG
652
653
654C-------------------------------------------------------------------------
655C FIN DU TRI vs Main shell edges
656C-------------------------------------------------------------------------
657 IF(j_stok/=0)CALL i25sto_edg(
658 1 j_stok ,irect ,x ,ii_stok,inacti,
659 2 cands_e2e,candm_e2e ,mulnsne,noint ,marge ,
660 3 i_mem(1) ,prov_s ,prov_m ,igap0,cand_a,
661 4 nedge ,ledge ,itab ,drad2 ,igap ,
662 5 gape ,gap_e_l,admsr ,edg_bisector,vtx_bisector ,
663 6 cand_p,ifq,cande2e_fx ,cande2e_fy,cande2e_fz,ifpen_e,
664 7 dgapload)
665
666 300 CONTINUE
667C=======================================================================
668C Sorting vs main solid edges
669C=======================================================================
670 IF(sol_edge==0) GOTO 400
671C=======================================================================
672C 3bis A partir des voxels occupes par une edge main, on est en
673C mesure de connaitre toutes les edges escalves dans ce voisinage.
674C Ce qui permet de creer des couples cancidats pour le contact
675C Si la penetration est positive.
676C=======================================================================
677
678 j_stok = 0
679
680 DO i=1,nrtm_t
681
682 ne =sshift+i
683
684 IF(msegtyp(ne)/=0) cycle ! not a solid edge
685 IF(stf(ne)==zero) cycle ! on ne retient pas les facettes detruites
686
687 m1 = irect(1,ne)
688 m2 = irect(2,ne)
689 m3 = irect(3,ne)
690 m4 = irect(4,ne)
691
692 mm1= itab(m1)
693 mm2= itab(m2)
694 mm3= itab(m3)
695 mm4= itab(m4)
696
697 xx1=x(1,m1)
698 yy1=x(2,m1)
699 zz1=x(3,m1)
700 xx2=x(1,m2)
701 yy2=x(2,m2)
702 zz2=x(3,m2)
703 xx3=x(1,m3)
704 yy3=x(2,m3)
705 zz3=x(3,m3)
706 xx4=x(1,m4)
707 yy4=x(2,m4)
708 zz4=x(3,m4)
709
710 xmax_edgm=max(xx1,xx2,xx3,xx4) ! +TZINF
711 xmin_edgm=min(xx1,xx2,xx3,xx4) ! -TZINF
712 ymax_edgm=max(yy1,yy2,yy3,yy4) ! +TZINF
713 ymin_edgm=min(yy1,yy2,yy3,yy4) ! -TZINF
714 zmax_edgm=max(zz1,zz2,zz3,zz4) ! +TZINF
715 zmin_edgm=min(zz1,zz2,zz3,zz4) ! -TZINF
716
717 dx=em02*(xmax_edgm-xmin_edgm)
718 dy=em02*(ymax_edgm-ymin_edgm)
719 dz=em02*(zmax_edgm-zmin_edgm)
720 xmax_edgm=xmax_edgm+dx
721 xmin_edgm=xmin_edgm-dx
722 ymax_edgm=ymax_edgm+dy
723 ymin_edgm=ymin_edgm-dy
724 zmax_edgm=zmax_edgm+dz
725 zmin_edgm=zmin_edgm-dz
726
727 aaa = marge+bgapemx+dgapload ! filtrer vs GAPE(JJ) dans i25pen3_edg !
728
729 !-------------------------------------------!
730 ! VOXEL OCCUPIED BY THE BRICK !
731 !-------------------------------------------!
732 !Voxel_lower_left_bound for this element---+
733 ix1=int(nbx*(xmin_edgm-aaa-xminb)/(xmaxb-xminb))
734 iy1=int(nby*(ymin_edgm-aaa-yminb)/(ymaxb-yminb))
735 iz1=int(nbz*(zmin_edgm-aaa-zminb)/(zmaxb-zminb))
736 ix1=max(1,2+min(nbx,ix1))
737 iy1=max(1,2+min(nby,iy1))
738 iz1=max(1,2+min(nbz,iz1))
739 !Voxel_upper_right_bound for this element---+
740 ix2=int(nbx*(xmax_edgm+aaa-xminb)/(xmaxb-xminb))
741 iy2=int(nby*(ymax_edgm+aaa-yminb)/(ymaxb-yminb))
742 iz2=int(nbz*(zmax_edgm+aaa-zminb)/(zmaxb-zminb))
743 ix2=max(1,2+min(nbx,ix2))
744 iy2=max(1,2+min(nby,iy2))
745 iz2=max(1,2+min(nbz,iz2))
746
747 IF(ilev==2)THEN
748 ims1 = bitget(mbinflg(ne),0)
749 ims2 = bitget(mbinflg(ne),1)
750 END IF
751
752#ifdef WITH_ASSERT
753C debug only
754 ids(1) = itab(irect(1,ne))
755 ids(2) = itab(irect(2,ne))
756 ids(3) = itab(irect(3,ne))
757 ids(4) = itab(irect(4,ne))
758 debug_e2e(int_checksum(ids,4,1)==d_em,xmin_edgm)
759 debug_e2e(int_checksum(ids,4,1)==d_em,ymin_edgm)
760 debug_e2e(int_checksum(ids,4,1)==d_em,zmin_edgm)
761 debug_e2e(int_checksum(ids,4,1)==d_em,xmax_edgm)
762 debug_e2e(int_checksum(ids,4,1)==d_em,ymax_edgm)
763 debug_e2e(int_checksum(ids,4,1)==d_em,zmax_edgm)
764#endif
765
766C--- IREMGAP - tag of deactivated lines
767 IF(flagremnode==2)THEN
768 k = kremnode_e2s(2*(ne-1)+1)
769 l = kremnode_e2s(2*(ne-1)+2)-1
770 DO m=k,l
771 tagremline(remnode_e2s(m)) = 1
772 ENDDO
773 ENDIF
774
775 DO iz = iz1,iz2
776 DO iy = iy1,iy2
777 DO ix = ix1,ix2
778
779 chain_add = voxel(ix,iy,iz) ! adresse dans le tableau chaine de la premiere edge stoquee dans le voxel
780 DO WHILE(chain_add /= 0) ! BOUCLE SUR LES EDGES DU VOXEL COURANT
781 jj = lchain_elem(chain_add) ! numeros des edge_id balayes dans le voxel courant
782
783
784 IF (jj<=nedge)THEN
785 eid = ledge(8,jj)
786 ELSE
787 eid = irem_edge(e_global_id,jj-nedge)
788 END IF
789
790 IF(tagedg(jj)/=0)THEN ! edge deja traitee vs cette arete main
791 chain_add = lchain_next(chain_add)
792 cycle
793 END IF
794 tagedg(jj)=1
795
796 !secnd edge nodes, exclure couples avec noeud commun
797 IF (jj<=nedge)THEN
798 ss1= itab(ledge(5,jj))
799 ss2= itab(ledge(6,jj))
800 ELSE
801 ss1=irem_edge(e_node1_globid,jj-nedge)
802 ss2=irem_edge(e_node2_globid,jj-nedge)
803 END IF
804
805 IF((ss1==mm1).OR.(ss1==mm2).OR.(ss1==mm3).OR.(ss1==mm4).OR.
806 . (ss2==mm1).OR.(ss2==mm2).OR.(ss2==mm3).OR.(ss2==mm4))THEN
807 chain_add = lchain_next(chain_add)
808 cycle
809 END IF
810
811 IF(ilev==2)THEN
812 IF(jj <= nedge) THEN
813 iss1=bitget(ebinflg(jj),0)
814 iss2=bitget(ebinflg(jj),1)
815 ELSE
816 iss1 = bitget(irem_edge(e_ebinflg,jj-nedge),0)
817 iss2 = bitget(irem_edge(e_ebinflg,jj-nedge),1)
818 ENDIF
819 IF(.NOT.((ims1 == 1 .and. iss2==1).or.
820 . (ims2 == 1 .and. iss1==1)))THEN
821 chain_add = lchain_next(chain_add)
822 cycle
823 ENDIF
824 ENDIF
825
826C IREMPGAP
827 IF (flagremnode == 2) THEN
828 IF (jj<=nedge)THEN
829C- Local Taged lines are removed
830 IF(tagremline(jj)==1) THEN
831 chain_add = lchain_next(chain_add)
832 cycle
833 ENDIF
834 ELSE
835C- Remote lines are identified by nodes
836 k = kremnode_e2s(2*(ne-1)+2)
837 l = kremnode_e2s(2*(ne-1)+3)-1
838 remove_remote = 0
839 DO m=k,l,2
840 IF ((ss1==remnode_e2s(m)).AND.(ss2==remnode_e2s(m+1))) remove_remote = 1
841 ENDDO
842 IF (remove_remote==1) THEN
843 chain_add = lchain_next(chain_add)
844 cycle
845 ENDIF
846 ENDIF
847 ENDIF
848
849CCC ================== DEBUG PRINT =====================
850C IF(JJ > NEDGE) THEN
851C WRITE(6,"(A,X,2I20)") "VOX REM",
852C . INT_CHECKSUM(IDS,4,1),IREM_EDGE(E_GLOBAL_ID,JJ-NEDGE)
853C ELSE
854C WRITE(6,"(A,X,2I20)") "VOX LOC",
855C . INT_CHECKSUM(IDS,4,1),LEDGE(8,JJ)
856C ENDIF
857CCC ================== DEBUG PRINT =====================
858C DEBUG_E2E(EID==D_ES.AND.INT_CHECKSUM(IDS,4,1)==D_EM,0)
859
860
861C ===================================================
862C-----------------------------------------------------
863 j_stok = j_stok + 1 !on dispose d'un candidat
864 prov_s(j_stok) = jj !edge secnd
865 prov_m(j_stok) = ne !segment main
866
867
868C DEBUG ONLY
869#ifdef WITH_ASSERT
870 prov_ids(2,j_stok) = eid
871 prov_ids(1,j_stok) = int_checksum(ids,4,1)
872#endif
873
874
875 assert(jj > 0)
876 assert(jj <= nedge + nedge_remote)
877C-----------------------------------------------------
878 IF(j_stok==nvsiz)THEN
879 CALL i25sto_e2s(
880 1 nvsiz ,irect ,x ,ll_stok,inacti,
881 2 cands_e2s,candm_e2s,mulnsns,noint ,marge ,
882 3 i_mem(2) ,prov_s ,prov_m ,igap0 ,cand_b,
883 4 nedge ,ledge ,itab ,drad2 ,igap ,
884 5 gap_m ,gap_m_l,gape ,gap_e_l,admsr ,
885 6 edg_bisector,vtx_bisector ,cand_ps,prov_ids,
886 7 ifq,cande2s_fx ,cande2s_fy,cande2s_fz,ifpen_e2s,
887 8 dgapload)
888
889 IF(i_mem(2)/=0) GOTO 400
890 j_stok = 0
891 ENDIF
892C-----------------------------------------------------
893
894 chain_add = lchain_next(chain_add) ! Next RTM
895
896 ENDDO ! WHILE(CHAIN_ADD /= 0)
897
898 ENDDO !NEXT IZ
899 ENDDO !NEXT IY
900 ENDDO !NEXT IZ
901
902C Reset TAGEDG
903 DO iz = iz1,iz2
904 DO iy = iy1,iy2
905 DO ix = ix1,ix2
906
907 chain_add = voxel(ix,iy,iz)
908 DO WHILE(chain_add /= 0) ! BOUCLE SUR LES EDGES DU VOXEL COURANT
909
910 jj = lchain_elem(chain_add) ! numeros des edge_id balayes dans le voxel courant
911 tagedg(jj)=0
912
913 chain_add = lchain_next(chain_add)
914
915 END DO
916
917 ENDDO !NEXT IZ
918 ENDDO !NEXT IY
919 ENDDO !NEXT IZ
920
921C--- IREMGAP - clean of tagremline
922 IF(flagremnode==2)THEN
923 k = kremnode_e2s(2*(ne-1)+1)
924 l = kremnode_e2s(2*(ne-1)+2)-1
925 DO m=k,l
926 tagremline(remnode_e2s(m)) = 0
927 ENDDO
928 ENDIF
929
930 ENDDO !NEXT IEDG
931C-------------------------------------------------------------------------
932C FIN DU TRI vs Main solid edges
933C-------------------------------------------------------------------------
934 IF(j_stok/=0)CALL i25sto_e2s(
935 1 j_stok ,irect ,x ,ll_stok,inacti,
936 2 cands_e2s,candm_e2s,mulnsns,noint ,marge ,
937 3 i_mem(2) ,prov_s ,prov_m ,igap0 ,cand_b,
938 4 nedge ,ledge ,itab ,drad2 ,igap ,
939 5 gap_m ,gap_m_l,gape ,gap_e_l,admsr ,
940 6 edg_bisector,vtx_bisector ,cand_ps,prov_ids,
941 7 ifq,cande2s_fx ,cande2s_fy,cande2s_fz,ifpen_e2s,
942 8 dgapload)
943
944C=======================================================================
945C 4 remise a zero des noeuds dans les boites
946C=======================================================================
947
948CC=============== DEBUG
949C DO I = 1, LL_STOK
950C JJ = CANDS_E2S(I)
951C NE = CANDM_E2S(I)
952C IDS(1) = ITAB(IRECT(1,NE))
953C IDS(2) = ITAB(IRECT(2,NE))
954C IDS(3) = ITAB(IRECT(3,NE))
955C IDS(4) = ITAB(IRECT(4,NE))
956C IF(JJ > NEDGE) THEN
957C CRITE(6,"(A,X,2I20)") "VOX REM",
958C . INT_CHECKSUM(IDS,4,1),IREM_EDGE(E_GLOBAL_ID,JJ-NEDGE)
959C ELSE
960C WRITE(6,"(A,X,2I20)") "VOX LOC",
961C . INT_CHECKSUM(IDS,4,1),LEDGE(8,JJ)
962C ENDIF
963C ENDDO
964C=======================================================================
965
966
967 400 CONTINUE
968
969
970
971C Barrier to avoid reinitialization before end of sorting
972 CALL my_barrier
973
974 tmin(1) = min_ix
975 tmin(2) = min_iy
976 tmin(3) = min_iz
977
978 tmax(1) = max_ix
979 tmax(2) = max_iy
980 tmax(3) = max_iz
981
982 IF (itask==0)THEN
983 !RESET VOXEL WITHIN USED RANGE ONLY
984 DO k= tmin(3),tmax(3)
985 DO j= tmin(2),tmax(2)
986 DO i= tmin(1),tmax(1)
987 voxel(i,j,k) = 0
988 END DO
989 END DO
990 END DO
991 !CHAINED LIST DEALLOCATION
992 DEALLOCATE(lchain_next)
993 DEALLOCATE(lchain_elem)
994 DEALLOCATE(lchain_last)
995 IF(flagremnode==2) DEALLOCATE(tagremline)
996 ENDIF
997
998 DEALLOCATE(tagedg)
999
1000C=======================================================================
1001
1002 RETURN
1003 END
i25trivox_edg
subroutine i25trivox_edg(i_mem, vmaxdt, inacti, irect, x, v, stf, stfe, xyzm, ii_stok, cands_e2e, eshift, nedge_t, candm_e2e, mulnsne, noint, bgapemx, sshift, nrtm_t, voxel, nbx, nby, nbz, igap, gap_m, gap_m_l, drad, marge, itask, itab, ll_stok, mulnsns, mbinflg, ebinflg, ilev, cand_a, cand_p, flagremnode, kremnode_edg, remnode_edg, kremnode_e2s, remnode_e2s, iedge, nedge, ledge, msegtyp, igap0, admsr, edg_bisector, vtx_bisector, candm_e2s, cands_e2s, cand_b, cand_ps, gape, gap_e_l, nedge_local, ifq, cande2e_fx, cande2e_fy, cande2e_fz, cande2s_fx, cande2s_fy, cande2s_fz, ifpen_e, ifpen_e2s, kremnode_edg_siz, remnode_edg_siz, kremnode_e2s_siz, remnode_e2s_siz, dgapload)
Definition i25trivox_edg.F:59
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
debug_mod
Definition debug_mod.F:45
debug_mod::int_checksum
pure integer function int_checksum(a, siz1, siz2)
Definition debug_mod.F:167
realloc_mod
Definition realloc_mod.F:29
realloc_mod::ireallocate
integer function, dimension(:), pointer ireallocate(ptr, new_size)
Definition realloc_mod.F:39
tri11::lchain_elem
integer, dimension(:), pointer lchain_elem
Definition tri11_mod.F:37
tri11::max_iz
integer max_iz
Definition tri11_mod.F:33
tri11::min_ix
integer min_ix
Definition tri11_mod.F:33
tri11::lchain_last
integer, dimension(:), pointer lchain_last
Definition tri11_mod.F:39
tri11::min_iz
integer min_iz
Definition tri11_mod.F:33
tri11::min_iy
integer min_iy
Definition tri11_mod.F:33
tri11::lchain_next
integer, dimension(:), pointer lchain_next
Definition tri11_mod.F:38
tri11::max_iy
integer max_iy
Definition tri11_mod.F:33
tri11::max_ix
integer max_ix
Definition tri11_mod.F:33
tri25ebox::nedge_remote
integer nedge_remote
Definition tri25ebox.F:73
tri25ebox::irem_edge
integer, dimension(:,:), allocatable irem_edge
Definition tri25ebox.F:64
i25sto_e2s
subroutine i25sto_e2s(j_stok, irect, x, ii_stok, inacti, cand_s, cand_m, mulnsn, noint, marge, i_mem, prov_s, prov_m, igap0, cand_a, nedge, ledge, itab, drad, igap, gap_m, gap_m_l, gape, gap_e_l, admsr, edg_bisector, vtx_bisector, cand_p, dgapload)
Definition i25sto_e2s.F:39
i25sto_edg
subroutine i25sto_edg(j_stok, irect, x, ii_stok, inacti, cand_s, cand_m, mulnsn, noint, marge, i_mem, prov_s, prov_m, igap0, cand_a, nedge, ledge, itab, drad, igap, gape, gap_e_l, admsr, edg_bisector, vtx_bisector, cand_p, dgapload)
Definition i25sto_edg.F:39
my_barrier
subroutine my_barrier
Definition machine.F:31