46
47
48
50 use my_alloc_mod
51
52
53
54#include "implicit_f.inc"
55
56
57
58#include "mvsiz_p.inc"
59
60 INTEGER NVECSZ
61 parameter(nvecsz = mvsiz)
62
63
64
65#include "com04_c.inc"
66#include "param_c.inc"
67#include "task_c.inc"
68
69
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 INTEGER I_MEM,ESHIFT,NSN,ISZNSNR,NRTM,NIN,ITASK,IGAP,
105 . MULNSN,NOINT,NSNR,NBX,NBY,NBZ,
106 . NSV(*),CAND_N(*),CAND_E(*),
107 . IRECT(4,*), VOXEL(NBX+2,NBY+2,NBZ+2),II_STOK,ITAB(*),
108 . NBINFLG(*),MBINFLG(*),ILEV,MSEGTYP(*),
109 . FLAGREMNODE,KREMNOD(*),REMNOD(*), ICODT(*), ISKEW(*)
110
112 . x(3,*),v(3,*),xyzm(6),stf(*),stfn(*),gap_s(*),gap_m(*),
113 . curv_max(*),pene_old(5,nsn),gap_s_l(*),gap_m_l(*),
114 . marge,bgapsmx,pmax_gap,vmaxdt
115 my_real ,
INTENT(IN) :: dgapload ,drad
116
117
118
119 INTEGER I,J,
120 . NN,NE,K,L,J_STOK,JJ,
121 . PROV_N(MVSIZ),PROV_E(MVSIZ),
122 . NSNF, NSNL,M,DELNOD
123
125 . xs,ys,zs,sx,sy,sz,s2,
126 . xmin, xmax,ymin,
ymax,zmin, zmax,
127 . xx1,xx2,xx3,xx4,yy1,yy2,yy3,yy4,zz1,zz2,zz3,zz4,
128 . d1x,d1y,d1z,d2x,d2y,d2z,dd1,dd2,d2,a2
129
130 INTEGER IX,IY,IZ,M1,M2,M3,M4,
131 . IX1,IY1,IZ1,IX2,IY2,IZ2
132 INTEGER, DIMENSION(:),ALLOCATABLE :: LAST_NOD
133 INTEGER, DIMENSION(:),ALLOCATABLE :: IIX,IIY,IIZ
135 . xminb,yminb,zminb,xmaxb,ymaxb,zmaxb,
136 . xmine,ymine,zmine,xmaxe,ymaxe,zmaxe,aaa
137 INTEGER FIRST,LAST
138 SAVE iix,iiy,iiz
139 INTEGER, DIMENSION(:),ALLOCATABLE :: TAG
140
141 CALL my_alloc(tag,numnod)
142 CALL my_alloc(last_nod,nsn+nsnr)
143 IF(itask == 0)THEN
145 ALLOCATE(iix(nsn+nsnr))
146 ALLOCATE(iiy(nsn+nsnr))
147 ALLOCATE(iiz(nsn+nsnr))
148 END IF
149
151
152
153 xmin = xyzm(1)
154 ymin = xyzm(2)
155 zmin = xyzm(3)
156 xmax = xyzm(4)
158 zmax = xyzm(6)
159
160
161 xminb = xmin
162 yminb = ymin
163 zminb = zmin
164 xmaxb = xmax
166 zmaxb = zmax
167
168
169
170
171 IF(itask == 0)THEN
172 DO i=1,nsn
173 iix(i)=0
174 iiy(i)=0
175 iiz(i)=0
176 IF(stfn(i) <= zero)cycle
177 j=nsv(i)
178
179
180 IF(x(1,j) < xmin) cycle
181 IF(x(1,j) > xmax) cycle
182 IF(x(2,j) < ymin) cycle
183 IF(x(2,j) >
ymax) cycle
184 IF(x(3,j) < zmin) cycle
185 IF(x(3,j) > zmax) cycle
186
187 iix(i)=int(nbx*(x(1,j)-xminb)/(xmaxb-xminb))
188 iiy(i)=int(nby*(x(2,j)-yminb)/(ymaxb-yminb))
189 iiz(i)=int(nbz*(x(3,j)-zminb)/(zmaxb-zminb))
190
191 iix(i)=
max(1,2+
min(nbx,iix(i)))
192 iiy(i)=
max(1,2+
min(nby,iiy(i)))
193 iiz(i)=
max(1,2+
min(nbz,iiz(i)))
194
195 first = voxel(iix(i),iiy(i),iiz(i))
196 IF(first == 0)THEN
197
198 voxel(iix(i),iiy(i),iiz(i)) = i
200 last_nod(i) = 0
201 ELSEIF(last_nod(first) == 0)THEN
202
203
205 last_nod(first) = i
207 ELSE
208
209
210 last = last_nod(first)
212 last_nod(first) = i
214 ENDIF
215 ENDDO
216
217
218
219
220 DO j = 1, nsnr
221 iix(nsn+j)=int(nbx*(xrem(1,j)-xminb)/(xmaxb-xminb))
222 iiy(nsn+j)=int(nby*(xrem(2,j)-yminb)/(ymaxb-yminb))
223 iiz(nsn+j)=int(nbz*(xrem(3,j)-zminb)/(zmaxb-zminb))
224 iix(nsn+j)=
max(1,2+
min(nbx,iix(nsn+j)))
225 iiy(nsn+j)=
max(1,2+
min(nby,iiy(nsn+j)))
226 iiz(nsn+j)=
max(1,2+
min(nbz,iiz(nsn+j)))
227
228 first = voxel(iix(nsn+j),iiy(nsn+j),iiz(nsn+j))
229 IF(first == 0)THEN
230
231 voxel(iix(nsn+j),iiy(nsn+j),iiz(nsn+j)) = nsn+j
233 last_nod(nsn+j) = 0
234 ELSEIF(last_nod(first) == 0)THEN
235
236
238 last_nod(first) = nsn+j
240 ELSE
241
242
243 last = last_nod(first)
245 last_nod(first) = nsn+j
247 ENDIF
248 ENDDO
249 END IF
250
252
253
254
255
256 j_stok = 0
257
258 IF(flagremnode == 2)THEN
259 DO i=1,numnod
260 tag(i) = 0
261 ENDDO
262
263 DO ne=1,nrtm
264
265 IF(stf(ne) <= zero)cycle
266 k = kremnod(2*(ne-1)+1)+1
267 l = kremnod(2*(ne-1)+2)
268 DO i=k,l
269 tag(remnod(i)) = 1
270 ENDDO
271
272 aaa = marge+curv_max(ne)+
max(
max(bgapsmx+gap_m(ne),pmax_gap)+dgapload,drad)+vmaxdt
273 +
274
275
276
277
278
279 m1 = irect(1,ne)
280 m2 = irect(2,ne)
281 m3 = irect(3,ne)
282 m4 = irect(4,ne)
283
284 xx1=x(1,m1)
285 xx2=x(1,m2)
286 xx3=x(1,m3)
287 xx4=x(1,m4)
288 xmaxe=
max(xx1,xx2,xx3,xx4)
289 xmine=
min(xx1,xx2,xx3,xx4)
290
291 yy1=x(2,m1)
292 yy2=x(2,m2)
293 yy3=x(2,m3)
294 yy4=x(2,m4)
295 ymaxe=
max(yy1,yy2,yy3,yy4)
296 ymine=
min(yy1,yy2,yy3,yy4)
297
298 zz1=x(3,m1)
299 zz2=x(3,m2)
300 zz3=x(3,m3)
301 zz4=x(3,m4)
302 zmaxe=
max(zz1,zz2,zz3,zz4)
303 zmine=
min(zz1,zz2,zz3,zz4)
304
305
306
307
308 sx = (yy3-yy1)*(zz4-zz2) - (zz3-zz1)*(yy4-yy2)
309 sy = (zz3-zz1)*(xx4-xx2) - (xx3-xx1)*(zz4-zz2)
310 sz = (xx3-xx1)*(yy4-yy2) - (yy3-yy1)*(xx4-xx2)
311 s2 = sx*sx + sy*sy + sz*sz
312
313
314
315 ix1=int(nbx*(xmine-aaa-xminb)/(xmaxb-xminb))
316 iy1=int(nby*(ymine-aaa-yminb)/(ymaxb-yminb))
317 iz1=int(nbz*(zmine-aaa-zminb)/(zmaxb-zminb))
318
322
323 ix2=int(nbx*(xmaxe+aaa-xminb)/(xmaxb-xminb))
324 iy2=int(nby*(ymaxe+aaa-yminb)/(ymaxb-yminb))
325 iz2=int(nbz*(zmaxe+aaa-zminb)/(zmaxb-zminb))
326
330
331
332
333
334
335
336 DO iz = iz1,iz2
337 DO iy = iy1,iy2
338 DO ix = ix1,ix2
339
340
341
342 jj = voxel(ix,iy,iz)
343
344 DO WHILE(jj /= 0)
345 delnod = 0
346
347
348 IF(jj<=nsn)THEN
349 nn=nsv(jj)
350 IF(nn == m1)GOTO 200
351 IF(nn == m2)GOTO 200
352 IF(nn == m3)GOTO 200
353 IF(nn == m4)GOTO 200
354 IF(tag(nn) == 1)GOTO 200
355 xs = x(1,nn)
356 ys = x(2,nn)
357 zs = x(3,nn)
358
359
360
361
362
363 aaa = marge + curv_max(ne)
364 + +
max(gap_s(jj)+gap_m(ne)+dgapload,drad)
365 + +vmaxdt
366 ELSE
367 j=jj-nsn
368
369 k = kremnod(2*(ne-1)+2) + 1
370 l = kremnod(2*(ne-1)+3)
371
372 DO m=k,l
373 IF(remnod(m) == -
irem(2,j) )
THEN
374 delnod = delnod + 1
375 EXIT
376 ENDIF
377 ENDDO
378
379 IF(delnod /= 0)GOTO 200
380
381 xs = xrem(1,j)
382 ys = xrem(2,j)
383 zs = xrem(3,j)
384 aaa = marge+curv_max(ne)
385
386
387
389 + + vmaxdt
390 ENDIF
391
392 IF(xs<=xmine-aaa)GOTO 200
393 IF(xs>=xmaxe+aaa)GOTO 200
394 IF(ys<=ymine-aaa)GOTO 200
395 IF(ys>=ymaxe+aaa)GOTO 200
396 IF(zs<=zmine-aaa)GOTO 200
397 IF(zs>=zmaxe+aaa)GOTO 200
398
399
400
401
402
403 d1x = xs - xx1
404 d1y = ys - yy1
405 d1z = zs - zz1
406 d2x = xs - xx2
407 d2y = ys - yy2
408 d2z = zs - zz2
409 dd1 = d1x*sx+d1y*sy+d1z*sz
410 dd2 = d2x*sx+d2y*sy+d2z*sz
411 IF(dd1*dd2 > zero)THEN
412 d2 =
min(dd1*dd1,dd2*dd2)
413 a2 = aaa*aaa*s2
414 IF(d2 > a2)GOTO 200
415 ENDIF
416
417
418
419 j_stok = j_stok + 1
420 prov_n(j_stok) = jj
421 prov_e(j_stok) = ne
422 IF(j_stok == nvsiz)THEN
423
425 1 nvsiz ,irect ,x ,nsv ,ii_stok,
426 2 cand_n,cand_e ,mulnsn,noint ,marge ,
427 3 i_mem ,prov_n ,prov_e,eshift,v ,
428 4 nsn ,nrtm ,gap_s ,gap_m ,curv_max,nin ,
429 5 pene_old,nbinflg ,mbinflg,ilev,msegtyp,
430 6 itab ,igap ,gap_s_l,gap_m_l,icodt,iskew,
431 7 drad ,dgapload)
432 IF(i_mem==2)GOTO 100
433 j_stok = 0
434 ENDIF
435
436 200 CONTINUE
437
439
440 ENDDO
441
442 ENDDO
443 ENDDO
444 ENDDO
445 k = kremnod(2*(ne-1)+1)+1
446 l = kremnod(2*(ne-1)+2)
447 DO i=k,l
448 tag(remnod(i)) = 0
449 ENDDO
450
451
452
453
454
455 ENDDO
456
457 ELSE
458 DO ne=1,nrtm
459
460 IF(stf(ne) <= zero)cycle
461
462 aaa = marge+curv_max(ne)+
max(
max(bgapsmx+gap_m(ne),pmax_gap)+dgapload,drad)+vmaxdt
463 +
464
465
466
467
468
469 m1 = irect(1,ne)
470 m2 = irect(2,ne)
471 m3 = irect(3,ne)
472 m4 = irect(4,ne)
473
474 xx1=x(1,m1)
475 xx2=x(1,m2)
476 xx3=x(1,m3)
477 xx4=x(1,m4)
478 xmaxe=
max(xx1,xx2,xx3,xx4)
479 xmine=
min(xx1,xx2,xx3,xx4)
480
481 yy1=x(2,m1)
482 yy2=x(2,m2)
483 yy3=x(2,m3)
484 yy4=x(2,m4)
485 ymaxe=
max(yy1,yy2,yy3,yy4)
486 ymine=
min(yy1,yy2,yy3,yy4)
487
488 zz1=x(3,m1)
489 zz2=x(3,m2)
490 zz3=x(3,m3)
491 zz4=x(3,m4)
492 zmaxe=
max(zz1,zz2,zz3,zz4)
493 zmine=
min(zz1,zz2,zz3,zz4)
494
495
496
497
498 sx = (yy3-yy1)*(zz4-zz2) - (zz3-zz1)*(yy4-yy2)
499 sy = (zz3-zz1)*(xx4-xx2) - (xx3-xx1)*(zz4-zz2)
500 sz = (xx3-xx1)*(yy4-yy2) - (yy3-yy1)*(xx4-xx2)
501 s2 = sx*sx + sy*sy + sz*sz
502
503
504
505 ix1=int(nbx*(xmine-aaa-xminb)/(xmaxb-xminb))
506 iy1=int(nby*(ymine-aaa-yminb)/(ymaxb-yminb))
507 iz1=int(nbz*(zmine-aaa-zminb)/(zmaxb-zminb))
508
512
513 ix2=int(nbx*(xmaxe+aaa-xminb)/(xmaxb-xminb))
514 iy2=int(nby*(ymaxe+aaa-yminb)/(ymaxb-yminb))
515 iz2=int(nbz*(zmaxe+aaa-zminb)/(zmaxb-zminb))
516
520
521
522
523
524
525
526 DO iz = iz1,iz2
527 DO iy = iy1,iy2
528 DO ix = ix1,ix2
529
530
531
532 jj = voxel(ix,iy,iz)
533
534 DO WHILE(jj /= 0)
535
536
537
538 IF(jj<=nsn)THEN
539 nn=nsv(jj)
540 IF(nn == m1)GOTO 300
541 IF(nn == m2)GOTO 300
542 IF(nn == m3)GOTO 300
543 IF(nn == m4)GOTO 300
544 xs = x(1,nn)
545 ys = x(2,nn)
546 zs = x(3,nn)
547
548
549
550
551
552 aaa = marge + curv_max(ne)
553 + +
max(gap_s(jj)+gap_m(ne)+dgapload,drad)
554 + +vmaxdt
555 ELSE
556 j=jj-nsn
557
558 xs = xrem(1,j)
559 ys = xrem(2,j)
560 zs = xrem(3,j)
561 aaa = marge+curv_max(ne)
562
563
564
566 + + vmaxdt
567 ENDIF
568
569 IF(xs<=xmine-aaa)GOTO 300
570 IF(xs>=xmaxe+aaa)GOTO 300
571 IF(ys<=ymine-aaa)GOTO 300
572 IF(ys>=ymaxe+aaa)GOTO 300
573 IF(zs<=zmine-aaa)GOTO 300
574 IF(zs>=zmaxe+aaa)GOTO 300
575
576
577
578
579
580 d1x = xs - xx1
581 d1y = ys - yy1
582 d1z = zs - zz1
583 d2x = xs - xx2
584 d2y = ys - yy2
585 d2z = zs - zz2
586 dd1 = d1x*sx+d1y*sy+d1z*sz
587 dd2 = d2x*sx+d2y*sy+d2z*sz
588 IF(dd1*dd2 > zero)THEN
589 d2 =
min(dd1*dd1,dd2*dd2)
590 a2 = aaa*aaa*s2
591 IF(d2 > a2)GOTO 300
592 ENDIF
593
594
595
596 j_stok = j_stok + 1
597 prov_n(j_stok) = jj
598 prov_e(j_stok) = ne
599 IF(j_stok == nvsiz)THEN
600
602 1 nvsiz ,irect ,x ,nsv ,ii_stok,
603 2 cand_n,cand_e ,mulnsn,noint ,marge ,
604 3 i_mem ,prov_n ,prov_e,eshift,v ,
605 4 nsn ,nrtm ,gap_s ,gap_m ,curv_max,nin ,
606 5 pene_old,nbinflg ,mbinflg,ilev,msegtyp,
607 6 itab ,igap ,gap_s_l,gap_m_l,icodt,iskew,
608 7 drad ,dgapload)
609 IF(i_mem==2)GOTO 100
610 j_stok = 0
611 ENDIF
612
613 300 CONTINUE
614
616
617 ENDDO
618
619 ENDDO
620 ENDDO
621 ENDDO
622
623
624
625
626
627 ENDDO
628 END IF
629
630
631
633 1 j_stok,irect ,x ,nsv ,ii_stok,
634 2 cand_n,cand_e ,mulnsn,noint ,marge ,
635 3 i_mem ,prov_n ,prov_e,eshift,v ,
636 4 nsn ,nrtm ,gap_s ,gap_m ,curv_max,nin ,
637 5 pene_old,nbinflg,mbinflg,ilev ,msegtyp,
638 6 itab ,igap ,gap_s_l,gap_m_l,icodt,iskew,
639 7 drad ,dgapload)
640
641
642
643
644 100 CONTINUE
645
646
648 nsnf = 1 + itask*nsn / nthread
649 nsnl = (itask+1)*nsn / nthread
650
651 DO i=nsnf,nsnl
652 IF(iix(i)/=0)THEN
653 voxel(iix(i),iiy(i),iiz(i))=0
654 ENDIF
655 ENDDO
656
657
658
659
660 nsnf = 1 + itask*nsnr / nthread
661 nsnl = (itask+1)*nsnr / nthread
662 DO j = nsnf, nsnl
663 voxel(iix(nsn+j),iiy(nsn+j),iiz(nsn+j))=0
664 ENDDO
665
666
668 IF(itask == 0)THEN
670 DEALLOCATE(iix)
671 DEALLOCATE(iiy)
672 DEALLOCATE(iiz)
673 ENDIF
674
675 DEALLOCATE(tag)
676 DEALLOCATE(last_nod)
677 RETURN
subroutine ymax(idn, fac, npc, pld, stiffmin, stiffmax, stiffini, stiffavg)
integer, dimension(:), allocatable next_nod
integer, dimension(:,:), allocatable irem
subroutine i25sto(j_stok, irect, x, nsv, local_i_stok, local_cand_n, local_cand_e, marge, prov_n, prov_e, eshift, nsn, nrtm, gap_s, gap_m, nbinflg, mbinflg, ilev, msegtyp, igap, gap_s_l, gap_m_l, edge_l2, icode, iskew, drad, dgapload, nrtmt)