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
115
116 INTEGER ::
117 . NRTM,NRTSR,ESHIFT,NRTS,IGAP,
118 . NSN4,NOINT,ITAB(*),NBX,NBY,NBZ,IAUTO,
119 . IRECTS(2,NRTS),IRECTM(2,NRTM)
120 INTEGER ITASK,IFORM
121 INTEGER, INTENT(INOUT) ::
122 . CAND_S(*),CAND_M(*),ADDCM(*),CHAINE(2,*),
123 . VOXEL(1:NBX+2,1:NBY+2,1:NBZ+2), I_MEM,IFPEN(*),II_STOK,
124 . FLAGREMNODE,KREMNODE(*),REMNODE(*)
126 . ,INTENT(IN) ::
127 . x(3,*),xyzm(6,*),
128 . stfs(nrts),stfm(nrtm), tzinf, gap
129 my_real ,
INTENT(IN) :: dgapload,drad
131 . gapmin,marge,bgapsmx,
132 . gap_s(*),gap_m(*), gap_s_l(*), gap_m_l(*)
133
134
135
136 INTEGER
137 . I,J,SS1,SS2,IBUG,
138 . N1,N2,MM1,MM2, iN1, iN2, iM1, , K,L,
139 . PROV_S(2*MVSIZ),PROV_M(2*MVSIZ),
140 . IX1,IY1,IZ1,IX2,IY2,IZ2,
141 . IX,IY,IZ, FIRST_ADD,
142 . I_STOK, I_STOK_BAK, IEDG,
143 . PREV_ADD, CHAIN_ADD, CURRENT_ADD,
144 . NEDG, DEJA , MAX_ADD ,II_STOK0, M,REMOVE_REMOTE
145 INTEGER, DIMENSION(3) :: TMIN,
147 . xx1, xx2,
148 . xmin, xmax,ymin,
ymax,zmin, zmax,
149 . yy1,yy2,zz1,zz2,
150 . aaa, dd,
151 . xmax_edgs(nrts+nrtsr), xmin_edgs(nrts+nrtsr),
152 . ymax_edgs(nrts+nrtsr), ymin_edgs(nrts+nrtsr),
153 . zmax_edgs(nrts+nrtsr), zmin_edgs(nrts+nrtsr),
154 . xmax_edgm(nrtm), xmin_edgm(nrtm),
155 . ymax_edgm(nrtm), ymin_edgm(nrtm),
156 . zmax_edgm(nrtm), zmin_edgm(nrtm),
157 . xminb,yminb,zminb,xmaxb,ymaxb,zmaxb
158
159 INTEGER, DIMENSION(:), ALLOCATABLE :: TAGREMLINE
160
161
162 IF(flagremnode==2) THEN
163 ALLOCATE(tagremline(nrts))
164 tagremline(1:nrts) = 0
165 ENDIF
166
167 aaa = zero
168
175
176
177
178
179 IF(itask == 0)THEN
180 max_add =
max(1,4*(nrts+nrtsr))
184 END IF
185
187
188 IF(nrtm==0.OR.nrts==0)THEN
189
193 END IF
194
195
196
197
198 xmin = xyzm(1,1)
199 ymin = xyzm(2,1)
200 zmin = xyzm(3,1)
201 xmax = xyzm(4,1)
203 zmax = xyzm(6,1)
204
205 xminb = xmin
206 yminb = ymin
207 zminb = zmin
208 xmaxb = xmax
210 zmaxb = zmax
211
212
213
214
215
216
217
218
219 IF(itask == 0)THEN
220
221 current_add=1
222
223 DO i = 1,nrts
224
225 IF(stfs(i)==zero)cycle
226
227
228
229
230 n1=irects(1,i)
231 n2=irects(2,i)
232
233
234
235
236
237 xx1=x(1,n1)
238 xx2=x(1,n2)
239 xmax_edgs(i)=
max(xx1,xx2);
IF(xmax_edgs(i) < xmin) cycle
240 xmin_edgs(i)=
min(xx1,xx2);
IF(xmin_edgs(i) > xmax) cycle
241 yy1=x(2,n1)
242 yy2=x(2,n2)
243 ymax_edgs(i)=
max(yy1,yy2);
IF(ymax_edgs(i) < ymin) cycle
244 ymin_edgs(i)=
min(yy1,yy2);
IF(ymin_edgs(i) >
ymax) cycle
245 zz1=x(3,n1)
246 zz2=x(3,n2)
247 zmax_edgs(i)=
max(zz1,zz2);
IF(zmax_edgs(i) < zmin) cycle
248 zmin_edgs(i)=
min(zz1,zz2);
IF(zmin_edgs(i) > zmax) cycle
249
250
251
252
253
254 ix1=int(nbx*(xmin_edgs(i)-xminb)/(xmaxb-xminb))
255 iy1=int(nby*(ymin_edgs(i)-yminb)/(ymaxb-yminb))
256 iz1=int(nbz*(zmin_edgs(i)-zminb)/(zmaxb-zminb))
260
261 ix2=int(nbx*(xmax_edgs(i)-xminb)/(xmaxb-xminb))
262 iy2=int(nby*(ymax_edgs(i)-yminb)/(ymaxb-yminb))
263 iz2=int(nbz*(zmax_edgs(i)-zminb)/(zmaxb-zminb))
267
268
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
304
305 DO iz = iz1,iz2
306 DO iy = iy1,iy2
307 DO ix = ix1,ix2
308
309 first_add = voxel(ix,iy,iz)
310
311 IF(first_add == 0)THEN
312
313 voxel(ix,iy,iz) = current_add
317 ELSE
318
324 ENDIF
325
326 current_add = current_add+1
327
328 IF( current_add>=max_add)THEN
329
330
331 max_add = 2 * max_add
332
336 ENDIF
337
338 ENDDO
339 ENDDO
340 ENDDO
341
342 ENDDO
343
344
345
346
347
348
349
350
351 DO i = nrts+1,nrts+nrtsr
352
353 j=i-nrts
354
355
356
357
358
359 xx1=xrem(1,j)
360 xx2=xrem(8,j)
361 xmax_edgs(i)=
max(xx1,xx2) ;
IF(xmax_edgs(i) < xmin) cycle
362 xmin_edgs(i)=
min(xx1,xx2) ;
IF(xmin_edgs(i) > xmax) cycle
363 yy1=xrem(2,j)
364 yy2=xrem(9,j)
365 ymax_edgs(i)=
max(yy1,yy2) ;
IF(ymax_edgs(i) < ymin) cycle
366 ymin_edgs(i)=
min(yy1,yy2) ;
IF(ymin_edgs(i) >
ymax) cycle
367 zz1=xrem(3,j)
368 zz2=xrem(10,j)
369 zmax_edgs(i)=
max(zz1,zz2) ;
IF(zmax_edgs(i) < zmin) cycle
370 zmin_edgs(i)=
min(zz1,zz2) ;
IF(zmin_edgs(i) > zmax) cycle
371
372
373
374
375
376 ix1=int(nbx*(xmin_edgs(i)-xminb)/(xmaxb-xminb))
377 iy1=int(nby*(ymin_edgs(i)-yminb)/(ymaxb-yminb))
378 iz1=int(nbz*(zmin_edgs(i)-zminb)/(zmaxb-zminb))
382
383 ix2=int(nbx*(xmax_edgs(i)-xminb)/(xmaxb-xminb))
384 iy2=int(nby*(ymax_edgs(i)-yminb)/(ymaxb-yminb
385 iz2=int(nbz*(zmax_edgs(i)-zminb)/(zmaxb-zminb))
389
390
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
425
426
427 DO iz = iz1,iz2
428 DO iy = iy1,iy2
429 DO ix = ix1,ix2
430 first_add = voxel(ix,iy,iz)
431 IF(first_add == 0)THEN
432
433 voxel(ix,iy,iz) = current_add
437 ELSE
438
444 ENDIF
445 current_add = current_add+1
446 IF( current_add>=max_add)THEN
447
448
449 max_add = 2 * max_add
450
454 ENDIF
455 ENDDO
456 ENDDO
457 ENDDO
458
459 ENDDO
460
461
462 END IF
463
465
466
467
468
469
470
471
472 nedg = 0
473 i_stok = 0
474 marge = tzinf -
max(gap+dgapload,drad)
475
476 DO iedg=1,nrtm
477
478 IF(stfm(iedg) == zero)cycle
479
480
481 aaa = tzinf
482 IF(igap == 0)THEN
483 aaa = tzinf
484 ELSE
485 aaa = marge+
486 .
max(
max(gapmin,bgapsmx+gap_m(iedg))+dgapload,drad)
487 ENDIF
488
489
490
491
492
493 n1 = irectm(1,iedg)
494 n2 = irectm(2,iedg)
495 mm1 = itab(n1)
496 mm2 = itab(n2)
497
498
499
500
501 xx1=x(1,n1)
502 xx2=x(1,n2)
503 yy1=x(2,n1)
504 yy2=x(2,n2)
505 zz1=x(3,n1)
506 zz2=x(3,n2)
507
508 xmax_edgm(iedg)=
max(xx1,xx2)
509 xmin_edgm(iedg)=
min(xx1,xx2)
510 ymax_edgm(iedg)=
max(yy1,yy2)
511 ymin_edgm(iedg)=
min(yy1,yy2)
512 zmax_edgm(iedg)=
max(zz1,zz2)
513 zmin_edgm(iedg)=
min(zz1,zz2)
514
515
516
517
518
519 ix1=int(nbx*(xmin_edgm(iedg)-aaa-xminb)/(xmaxb-xminb))
520 iy1=int(nby*(ymin_edgm(iedg)-aaa-yminb)/(ymaxb-yminb))
521 iz1=int(nbz*(zmin_edgm(iedg)-aaa-zminb)/(zmaxb-zminb))
525
526 ix2=int(nbx*(xmax_edgm(iedg)+aaa-xminb)/(xmaxb-xminb))
527 iy2=int(nby*(ymax_edgm(iedg)+aaa-yminb)/(ymaxb-yminb))
528 iz2=int(nbz*(zmax_edgm(iedg)+aaa-zminb)/(zmaxb-zminb))
532
533 deja = 0
534 i_stok_bak = i_stok
535
536
537 IF(flagremnode==2)THEN
538 k = kremnode(2*(iedg-1)+1)
539 l = kremnode(2*(iedg-1)+2)-1
540 DO m=k,l
541 tagremline(remnode(m)) = 1
542 ENDDO
543 ENDIF
544
545
546
547
548 DO iz = iz1,iz2
549 DO iy = iy1,iy2
550 DO ix = ix1,ix2
551
552 chain_add = voxel(ix,iy,iz)
553 DO WHILE(chain_add /= 0)
555
556
557 IF (i<=nrts)THEN
558 ss1=itab(irects(1,i))
559 ss2=itab(irects(2,i))
560 ELSE
563 END IF
564
565 IF( (ss1==mm1).OR.(ss1==mm2).OR.
566 . (ss2==mm1).OR.(ss2==mm2) )THEN
568 cycle
569 END IF
570
571
572 IF(iauto==1 .AND. mm1<ss1 )THEN
574 cycle
575 END IF
576
577
578 IF (flagremnode == 2) THEN
579 IF (i <= nrts) THEN
580
581 IF(tagremline(i)==1) THEN
583 cycle
584 ENDIF
585 ELSE
586
587 k = kremnode(2*(iedg-1)+2)
588 l = kremnode(2*(iedg-1)+3)-1
589 remove_remote = 0
590 DO m=k,l,2
591 IF ((ss1==remnode(m)).AND.(ss2==remnode(m+1))) remove_remote = 1
592 ENDDO
593 IF (remove_remote==1) THEN
595 cycle
596 ENDIF
597 ENDIF
598 ENDIF
599
600 i_stok = i_stok + 1
601 prov_s(i_stok) = i
602 prov_m(i_stok) = iedg
603
604
605 IF(deja==0) nedg = nedg + 1
606 deja=1
608
609 IF(i_stok>=nvsiz)THEN
611 1 nvsiz ,irects,irectm,x ,ii_stok,
612 2 cand_s,cand_m,nsn4 ,noint ,marge,
613 3 i_mem ,prov_s,prov_m,eshift,addcm ,
614 4 chaine,nrts ,itab ,ifpen ,iform,
615 5 gapmin,drad ,igap, gap_s, gap_m,
616 7 gap_s_l, gap_m_l ,dgapload)
617
618 IF(i_mem==2) THEN
619
620
621 GOTO 1000
623 i_stok = i_stok-nvsiz
624
625 DO j=1,i_stok
626 prov_s(j) = prov_s(j+nvsiz)
627 prov_m(j) = prov_m(j+nvsiz)
628 ENDDO
629 ENDIF
630
631
632 ENDDO
633 ENDDO
634 ENDDO
635 ENDDO
636
637
638 IF(flagremnode==2)THEN
639 k = kremnode(2*(iedg-1)+1)
640 l = kremnode(2*(iedg-1)+2)-1
641 DO m=k,l
642 tagremline(remnode(m)) = 0
643 ENDDO
644 ENDIF
645
646 ENDDO
647
648
649
650
651
653 1 i_stok,irects,irectm,x ,ii_stok,
654 2 cand_s,cand_m,nsn4 ,noint ,marge ,
655 3 i_mem ,prov_s,prov_m,eshift,addcm ,
656 4 chaine,nrts ,itab ,ifpen ,iform ,
657 5 gapmin,drad ,igap, gap_s ,gap_m ,
658 7 gap_s_l, gap_m_l ,dgapload)
659
660
661
662
663
664
665 1000 CONTINUE
666
668
669
670
671
675
679
680 IF (itask==0)THEN
681
682 DO k= tmin(3),tmax(3)
683 DO j= tmin(2),tmax(2)
684 DO i= tmin(1),tmax(1)
685 voxel(i,j,k) = 0
686 END DO
687 END DO
688 END DO
689
693 IF(flagremnode==2) DEALLOCATE(tagremline)
694 END IF
695
696
697
698 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