46
47
48
52
53
54
55#include "implicit_f.inc"
56#include "comlock.inc"
57
58
59
60#include "mvsiz_p.inc"
61
62
63
64#include "param_c.inc"
65
66
67
68#ifdef MPI
69#endif
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114 INTEGER ::
115 . NRTM,NRTSR,ESHIFT,NRTS,IGAP,
116 . NSN4,NOINT,ITAB(*),NBX,NBY,NBZ,IAUTO,
117 . IRECTS(2,NRTS),IRECTM(2,NRTM)
118 INTEGER ITASK,IFORM
119 INTEGER, INTENT(INOUT) ::
120 . CAND_S(*),CAND_M(*),ADDCM(*),CHAINE(2,*),
121 . VOXEL(1:NBX+2,1:NBY+2,1:NBZ+2), I_MEM,IFPEN(*),II_STOK,
122 . FLAGREMNODE,KREMNODE(*),REMNODE(*)
124 . ,INTENT(IN) ::
125 . x(3,*),xyzm(6,*),
126 . stfs(nrts),stfm(nrtm), tzinf, gap
127 my_real ,
INTENT(IN) :: dgapload,drad
129 . gapmin,marge,bgapsmx,
130 . gap_s(*),gap_m(*), gap_s_l(*), gap_m_l(*)
131
132
133
134 INTEGER
135 . I,J,SS1,SS2,
136 . ,N2,MM1,MM2, K,,
137 . PROV_S(2*MVSIZ),PROV_M(2*MVSIZ),
138 . IX1,IY1,IZ1,IX2,IY2,IZ2,
139 . IX,IY,IZ, FIRST_ADD,
140 . I_STOK, I_STOK_BAK, IEDG,
141 . PREV_ADD, CHAIN_ADD, CURRENT_ADD,
142 . NEDG, DEJA , MAX_ADD , M,REMOVE_REMOTE
143 INTEGER, DIMENSION(3) :: TMIN,TMAX
145 . xx1, xx2,
146 . xmin, xmax,ymin,
ymax,zmin, zmax,
147 . yy1,yy2,zz1,zz2,
148 . aaa,
149 . xmax_edgs(nrts+nrtsr), xmin_edgs(nrts+nrtsr),
150 . ymax_edgs(nrts+nrtsr), ymin_edgs(nrts+nrtsr),
151 . zmax_edgs(nrts+nrtsr), zmin_edgs(nrts+nrtsr),
152 . xmax_edgm(nrtm), xmin_edgm(nrtm),
153 . ymax_edgm(nrtm), ymin_edgm(nrtm),
154 . zmax_edgm(nrtm), zmin_edgm(nrtm),
155 . xminb,yminb,zminb,xmaxb,ymaxb,zmaxb
156
157 INTEGER, DIMENSION(:), ALLOCATABLE :: TAGREMLINE
158
159
160 IF(flagremnode==2) THEN
161 ALLOCATE(tagremline(nrts))
162 tagremline(1:nrts) = 0
163 ENDIF
164
165 aaa = zero
166
173
174
175
176
177 IF(itask == 0)THEN
178 max_add =
max(1,4*(nrts+nrtsr))
182 END IF
183
185
186 IF(nrtm==0.OR.nrts==0)THEN
187
191 END IF
192
193
194
195
196 xmin = xyzm(1,1)
197 ymin = xyzm(2,1)
198 zmin = xyzm(3,1)
199 xmax = xyzm(4,1)
201 zmax = xyzm(6,1)
202
203 xminb = xmin
204 yminb = ymin
205 zminb = zmin
206 xmaxb = xmax
208 zmaxb = zmax
209
210
211
212
213
214
215
216
217 IF(itask == 0)THEN
218
219 current_add=1
220
221 DO i = 1,nrts
222
223 IF(stfs(i)==zero)cycle
224
225
226
227
228 n1=irects(1,i)
229 n2=irects(2,i)
230
231
232
233
234
235 xx1=x(1,n1)
236 xx2=x(1,n2)
237 xmax_edgs(i)=
max(xx1,xx2);
IF(xmax_edgs(i) < xmin) cycle
238 xmin_edgs(i)=
min(xx1,xx2);
IF(xmin_edgs(i) > xmax) cycle
239 yy1=x(2,n1)
240 yy2=x(2,n2)
241 ymax_edgs(i)=
max(yy1,yy2);
IF(ymax_edgs(i) < ymin) cycle
242 ymin_edgs(i)=
min(yy1,yy2);
IF(ymin_edgs(i) >
ymax) cycle
243 zz1=x(3,n1)
244 zz2=x(3,n2)
245 zmax_edgs(i)=
max(zz1,zz2);
IF(zmax_edgs(i) < zmin) cycle
246 zmin_edgs(i)=
min(zz1,zz2);
IF(zmin_edgs(i) > zmax) cycle
247
248
249
250
251
252 ix1=int(nbx*(xmin_edgs(i)-xminb)/(xmaxb-xminb))
253 iy1=int(nby*(ymin_edgs(i)-yminb)/(ymaxb-yminb))
254 iz1=int(nbz*(zmin_edgs(i)-zminb)/(zmaxb-zminb))
258
259 ix2=int(nbx*(xmax_edgs(i)-xminb)/(xmaxb-xminb))
260 iy2=int(nby*(ymax_edgs(i)-yminb)/(ymaxb-yminb))
261 iz2=int(nbz*(zmax_edgs(i)-zminb)/(zmaxb-zminb))
265
266
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303 DO iz = iz1,iz2
304 DO iy = iy1,iy2
305 DO ix = ix1,ix2
306
307 first_add = voxel(ix,iy,iz)
308
309 IF(first_add == 0)THEN
310
311 voxel(ix,iy,iz) = current_add
315 ELSE
316
322 ENDIF
323
324 current_add = current_add+1
325
326 IF( current_add>=max_add)THEN
327
328
329 max_add = 2 * max_add
330
334 ENDIF
335
336 ENDDO
337 ENDDO
338 ENDDO
339
340 ENDDO
341
342
343
344
345
346
347
348 DO i = nrts+1,nrts+nrtsr
349
350 j=i-nrts
351
352
353
354
355
356 xx1=xrem(1,j)
357 xx2=xrem(8,j)
358 xmax_edgs(i)=
max(xx1,xx2) ;
IF(xmax_edgs(i) < xmin) cycle
359 xmin_edgs(i)=
min(xx1,xx2) ;
IF(xmin_edgs(i) > xmax) cycle
360 yy1=xrem(2,j)
361 yy2=xrem(9,j)
362 ymax_edgs(i)=
max(yy1,yy2) ;
IF(ymax_edgs(i) < ymin) cycle
363 ymin_edgs(i)=
min(yy1,yy2) ;
IF(ymin_edgs(i) >
ymax) cycle
364 zz1=xrem(3,j)
365 zz2=xrem(10,j)
366 zmax_edgs(i)=
max(zz1,zz2) ;
IF(zmax_edgs(i) < zmin) cycle
367 zmin_edgs(i)=
min(zz1,zz2) ;
IF(zmin_edgs(i) > zmax) cycle
368
369
370
371
372
373 ix1=int(nbx*(xmin_edgs(i)-xminb)/(xmaxb-xminb))
374 iy1=int(nby*(ymin_edgs(i)-yminb)/(ymaxb-yminb))
375 iz1=int(nbz*(zmin_edgs(i)-zminb)/(zmaxb-zminb))
379
380 ix2=int(nbx*(xmax_edgs(i)-xminb)/(xmaxb-xminb))
381 iy2=int(nby*(ymax_edgs(i)-yminb)/(ymaxb-yminb))
382 iz2=int(nbz*(zmax_edgs(i)-zminb)/(zmaxb-zminb))
386
387
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424 DO iz = iz1,iz2
425 DO iy = iy1,iy2
426 DO ix = ix1,ix2
427 first_add = voxel(ix,iy,iz)
428 IF(first_add == 0)THEN
429
430 voxel(ix,iy,iz) = current_add
434 ELSE
435 !voxel contenant deja une edge
441 ENDIF
442 current_add = current_add+1
443 IF( current_add>=max_add)THEN
444
445
446 max_add = 2 * max_add
447
451 ENDIF
452 ENDDO
453 ENDDO
454 ENDDO
455
456 ENDDO
457
458
459 END IF
460
462
463
464
465
466
467
468
469 nedg = 0
470 i_stok = 0
471 marge = tzinf -
max(gap+dgapload,drad)
472
473 DO iedg=1,nrtm
474
475 IF(stfm(iedg) == zero)cycle
476
477
478 aaa = tzinf
479 IF(igap == 0)THEN
480 aaa = tzinf
481 ELSE
482 aaa = marge+
483 .
max(
max(gapmin,bgapsmx+gap_m(iedg))+dgapload,drad)
484 ENDIF
485
486
487
488
489
490 n1 = irectm(1,iedg)
491 n2 = irectm(2,iedg)
492 mm1 = itab(n1)
493 mm2 = itab(n2)
494
495
496
497
498 xx1=x(1,n1)
499 xx2=x(1,n2)
500 yy1=x(2,n1)
501 yy2=x(2,n2)
502 zz1=x(3,n1)
503 zz2=x(3,n2)
504
505 xmax_edgm(iedg)=
max(xx1,xx2)
506 xmin_edgm(iedg)=
min(xx1,xx2)
507 ymax_edgm(iedg)=
max(yy1,yy2)
508 ymin_edgm(iedg)=
min(yy1,yy2)
509 zmax_edgm(iedg)=
max(zz1,zz2)
510 zmin_edgm(iedg)=
min(zz1,zz2)
511
512
513
514
515
516 ix1=int(nbx*(xmin_edgm(iedg)-aaa-xminb)/(xmaxb-xminb))
517 iy1=int(nby*(ymin_edgm(iedg)-aaa-yminb)/(ymaxb-yminb))
518 iz1=int(nbz*(zmin_edgm(iedg)-aaa-zminb)/(zmaxb-zminb))
522
523 ix2=int(nbx*(xmax_edgm(iedg)+aaa-xminb)/(xmaxb-xminb))
524 iy2=int(nby*(ymax_edgm(iedg)+aaa-yminb)/(ymaxb-yminb))
525 iz2=int(nbz*(zmax_edgm(iedg)+aaa-zminb)/(zmaxb-zminb))
529
530 deja = 0
531 i_stok_bak = i_stok
532
533
534 IF(flagremnode==2)THEN
535 k = kremnode(2*(iedg-1)+1)
536 l = kremnode(2*(iedg-1)+2)-1
537 DO m=k,l
538 tagremline(remnode(m)) = 1
539 ENDDO
540 ENDIF
541
542
543
544
545 DO iz = iz1,iz2
546 DO iy = iy1,iy2
547 DO ix = ix1,ix2
548
549 chain_add = voxel(ix,iy,iz)
550 DO WHILE(chain_add /= 0)
552
553
554 IF (i<=nrts)THEN
555 ss1=itab(irects(1,i))
556 ss2=itab(irects(2,i))
557 ELSE
560 END IF
561
562 IF( (ss1==mm1).OR.(ss1==mm2).OR.
563 . (ss2==mm1).OR.(ss2==mm2) )THEN
565 cycle
566 END IF
567
568
569 IF(iauto==1 .AND. mm1<ss1 )THEN
571 cycle
572 END IF
573
574
575 IF (flagremnode == 2) THEN
576 IF (i <= nrts) THEN
577
578 IF(tagremline(i)==1) THEN
580 cycle
581 ENDIF
582 ELSE
583
584 k = kremnode(2*(iedg-1)+2)
585 l = kremnode(2*(iedg-1)+3)-1
586 remove_remote = 0
587 DO m=k,l,2
588 IF ((ss1==remnode(m)).AND.(ss2==remnode(m+1))) remove_remote = 1
589 ENDDO
590 IF (remove_remote==1) THEN
592 cycle
593 ENDIF
594 ENDIF
595 ENDIF
596
597 i_stok = i_stok + 1
598 prov_s(i_stok) = i
599 prov_m(i_stok) = iedg
600
601
602 IF(deja==0) nedg = nedg + 1
603 deja=1
605
606 IF(i_stok>=nvsiz)THEN
608 1 nvsiz ,irects,irectm,x ,ii_stok,
609 2 cand_s,cand_m,nsn4 ,noint ,marge,
610 3 i_mem ,prov_s,prov_m,eshift,addcm ,
611 4 chaine,nrts ,itab ,ifpen ,iform,
612 5 gapmin,drad ,igap, gap_s, gap_m,
613 7 gap_s_l, gap_m_l ,dgapload)
614
615 IF(i_mem==2) THEN
616
617
618 GOTO 1000
620 i_stok = i_stok-nvsiz
621
622 DO j=1,i_stok
623 prov_s(j) = prov_s(j+nvsiz)
624 prov_m(j) = prov_m(j+nvsiz)
625 ENDDO
626 ENDIF
627
628
629 ENDDO
630 ENDDO
631 ENDDO
632 ENDDO
633
634
635 IF(flagremnode==2)THEN
636 k = kremnode(2*(iedg-1)+1)
637 l = kremnode(2*(iedg-1)+2)-1
638 DO m=k,l
639 tagremline(remnode(m)) = 0
640 ENDDO
641 ENDIF
642
643 ENDDO
644
645
646
647
648
650 1 i_stok,irects,irectm,x ,ii_stok,
651 2 cand_s,cand_m,nsn4 ,noint ,marge ,
652 3 i_mem ,prov_s,prov_m,eshift,addcm ,
653 4 chaine,nrts ,itab ,ifpen ,iform ,
654 5 gapmin,drad ,igap, gap_s ,gap_m ,
655 7 gap_s_l, gap_m_l ,dgapload)
656
657
658
659
660
661
662 1000 CONTINUE
663
665
666
667
668
672
676
677 IF (itask==0)THEN
678
679 DO k= tmin(3),tmax(3)
680 DO j= tmin(2),tmax(2)
681 DO i= tmin(1),tmax(1)
682 voxel(i,j,k) = 0
683 END DO
684 END DO
685 END DO
686
690 IF(flagremnode==2) DEALLOCATE(tagremline)
691 END IF
692
693
694
695 RETURN
if(complex_arithmetic) id
subroutine i11sto_vox(j_stok, irects, irectm, x, ii_stok, cand_s, cand_m, nsn4, noint, marge, i_mem, prov_s, prov_m, eshift, addcm, chaine, nrts, itab, ifpen, iform, gapmin, drad, igap, gap_s, gap_m, gap_s_l, gap_m_l, dgapload)
subroutine ymax(idn, fac, npc, pld, stiffmin, stiffmax, stiffini, stiffavg)
integer function, dimension(:), pointer ireallocate(ptr, new_size)
integer, dimension(:), pointer lchain_elem
integer, dimension(:), pointer lchain_last
integer, dimension(:), pointer lchain_next
integer, dimension(:,:), allocatable irem