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