41
42
43
45
46
47
48#include "implicit_f.inc"
49
50
51
52#include "mvsiz_p.inc"
53
54
55
56#include "param_c.inc"
57
58
59
60
61
62
63
64
65
66
67
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
106
107
108
109
110
111
112
113
114 INTEGER NRTM,NRTSR,I_ADD,MAXSIZ,I_MEM,ESHIFT,NRTS,
115 . NSN4,NB_N_B,NOINT,I_ADD_MAX,IAUTO ,NIN,
116 . ADD(2,*),IRECTS(2,*),IRECTM(2,*),
117 . CAND_S(*),CAND_M(*),ADDCM(*),CHAINE(2,*),ITAB(*),
118 . NB_OLD(2,*),IFPEN(*),IFORM,II_STOK
119
121 . x(3,*),xyzm(6,*),stfs(*),stfm(*),
122 . tzinf,maxbox,minbox
123
124
125
126 INTEGER NB_NCN,NB_NCN1,NB_ECN,ADDNN,ADDNE,I,J,DIR,NN1,NN2,
127 . N1,N2,N3,N4,NN,NE,K_STOK,K,L,NCAND_PROV,J_STOK,NI,
128 . ISTOP,NB_ECN1,PROV_S(2*MVSIZ),PROV_M(2*MVSIZ),
129 . NB_NC_OLD, NB_EC_OLD, NB_NC, NB_EC,JJ,KK,
130
131
132 . BPE(NRTM+100),PE(MAXSIZ),BPN(NRTS+NRTSR+100),PN(MAXSIZ)
133
135 . dx,dy,dz,dsup,seuil,seuils,seuili, xx1, xx2, xx3, xx4,
136 . xmin, xmax,ymin,
ymax,zmin, zmax, xx,yy,zz,
137 . xmins,ymins,zmins,xmaxs,ymaxs,zmaxs,
138 . yy1,yy2,zz1,zz2,dmx,dmy,dmz,
139 . xy1,xy2,xz1,xz2,ximin,ximax,xjmin,xjmax,xkmin,xkmax,
140 . timin,timax,tjmin,tjmax,tkmin,tkmax,tsmin,tsmax,
141 . txmin, txmax,tymin, tymax,tzmin, tzmax
143 LOGICAL I11INSID
144
145
146
147
148
149
150 xmin = xyzm(1,i_add)
151 ymin = xyzm(2,i_add)
152 zmin = xyzm(3,i_add)
153 xmax = xyzm(4,i_add)
155 zmax = xyzm(6,i_add)
156
157
158 nb_ec = 0
159 DO i=1,nrtm
160
161 IF(stfm(i)/=zero)THEN
162 nb_ec = nb_ec + 1
163 bpe(nb_ec) = i
164 END IF
165 ENDDO
166
167
168
169
170 nb_nc = 0
171 DO i=1,nrts
172
173 IF(stfs(i)/=zero)THEN
174 n1=irects(1,i)
175 n2=irects(2,i)
176 xmins =
min(x(1,n1),x(1,n2))
177 ymins =
min(x(2,n1),x(2,n2))
178 zmins =
min(x(3,n1),x(3,n2))
179 xmaxs =
max(x(1,n1),x(1,n2))
180 ymaxs =
max(x(2,n1),x(2,n2))
181 zmaxs =
max(x(3,n1),x(3,n2))
182 IF(xmaxs>=xmin.AND.xmins<=xmax.AND.
183 . ymaxs>=ymin.AND.ymins<=
ymax.AND.
184 . zmaxs>=zmin.AND.zmins<=zmax)THEN
185 nb_nc = nb_nc + 1
186 bpn(nb_nc) = i
187 ENDIF
188 END IF
189 ENDDO
190
191
192
193 DO i = nrts+1, nrts+nrtsr
194 nb_nc = nb_nc + 1
195 bpn(nb_nc) = i
196 ENDDO
197
198
199
200
201 j_stok = 0
202 istop = 0
203 nb_nc_old = 0
204 nb_ec_old = 0
205
206 nb_old(1,i_add) = 0
207 nb_old(2,i_add) = 0
208
209 dx = xyzm(4,i_add) - xyzm(1,i_add)
210 dy = xyzm(5,i_add) - xyzm(2,i_add)
211 dz = xyzm(6,i_add) - xyzm(3,i_add)
213 GOTO 200
214
215 100 CONTINUE
216
217
218
219
220
221
222
223
224
225
226
227 xmin = 1.e30
228 xmax = -1.e30
229
230 ymin = 1.e30
232
233 zmin = 1.e30
234 zmax = -1.e30
235
236 DO i=1,nb_ec
237 ne = bpe(i)
238 xx1=x(1, irectm(1,ne))
239 xx2=x(1, irectm(2,ne))
240 xmin=
min(xmin,xx1,xx2)
241 xmax=
max(xmax,xx1,xx2)
242
243 yy1=x(2, irectm(1,ne))
244 yy2=x(2, irectm(2,ne))
245 ymin=
min(ymin,yy1,yy2)
247
248 zz1=x(3, irectm(1,ne))
249 zz2=x(3, irectm(2,ne))
250 zmin=
min(zmin,zz1,zz2)
251 zmax=
max(zmax,zz1,zz2)
252 ENDDO
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
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 xmin =
max(xmin - tzinf , xyzm(1,i_add))
302 ymin =
max(ymin - tzinf , xyzm(2,i_add))
303 zmin =
max(zmin - tzinf , xyzm(3,i_add))
304 xmax =
min(xmax + tzinf , xyzm(4,i_add))
306 zmax =
min(zmax + tzinf , xyzm(6,i_add))
307
308 txmin = xmin - tzinf
309 tymin = ymin - tzinf
310 tzmin = zmin - tzinf
311 txmax = xmax + tzinf
313 tzmax = zmax + tzinf
314
315 dmx = xmax-xmin
317 dmz = zmax-zmin
318
319 dsup =
max(dmx,dmy,dmz)
320
321 IF(dmy==dsup) THEN
322 dir = 2
323 jj = 3
324 kk = 1
325 seuil = (ymin+
ymax)*0.5
326 ximin = ymin
327 xjmin = zmin
328 xkmin = xmin
330 xjmax = zmax
331 xkmax = xmax
332 timin = tymin
333 tjmin = tzmin
334 tkmin = txmin
335 timax = tymax
336 tjmax = tzmax
337 tkmax = txmax
338 ELSE IF(dmz==dsup) THEN
339 dir = 3
340 jj = 1
341 kk = 2
342 seuil = (zmin+zmax)*0.5
343 ximin = zmin
344 xjmin = xmin
345 xkmin = ymin
346 ximax = zmax
347 xjmax = xmax
349 timin = tzmin
350 tjmin = txmin
351 tkmin = tymin
352 timax = tzmax
353 tjmax = txmax
354 tkmax = tymax
355 ELSE
356 dir = 1
357 jj = 2
358 kk = 3
359 seuil = (xmin+xmax)*0.5
360 ximin = xmin
361 xjmin = ymin
362 xkmin = zmin
363 ximax = xmax
365 xkmax = zmax
366 timin = txmin
367 tjmin = tymin
368 tkmin = tzmin
369 timax = txmax
370 tjmax = tymax
371 tkmax = tzmax
372 ENDIF
373
374 tsmin = seuil - tzinf
375 tsmax = seuil + tzinf
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393 nb_ncn= 0
394 nb_ncn1= 0
395 addnn= add(1,i_add)
396 DO i=1,nb_nc
397 nn = bpn(i)
398 IF(nn<=nrts) THEN
399 xx1=x(dir,irects(1,nn))
400 xx2=x(dir,irects(2,nn))
401 xy1=x(jj, irects(1,nn))
402 xy2=x(jj, irects(2,nn))
403 xz1=x(kk, irects(1,nn))
404 xz2=x(kk, irects(2,nn))
405 ELSE
406 ni = nn-nrts
407 xx1=xrem(dir,ni)
408 xx2=xrem(dir+7,ni)
409 xy1=xrem(jj ,ni)
410 xy2=xrem(jj+7 ,ni)
411 xz1=xrem(kk ,ni)
412 xz2=xrem(kk+7 ,ni)
413 END IF
416 IF(xmin<seuil.AND.xmax>=ximin) THEN
417 IF(
i11insid(xx1,xx2,xy1,xy2,xz1,xz2,
418 . ximin,seuil,xjmin,xjmax,xkmin,xkmax)) THEN
419
420 nb_ncn1 = nb_ncn1 + 1
421 addnn = addnn + 1
422 pn(addnn) = nn
423 END IF
424 END IF
425 ENDDO
426 DO i=1,nb_nc
427 nn = bpn(i)
428 IF(nn<=nrts) THEN
429 xx1=x(dir,irects(1,nn))
430 xx2=x(dir,irects(2,nn))
431 xy1=x(jj, irects(1,nn))
432 xy2=x(jj, irects(2,nn))
433 xz1=x(kk, irects(1,nn))
434 xz2=x(kk, irects(2,nn))
435 ELSE
436 ni = nn-nrts
437 xx1=xrem(dir,ni)
438 xx2=xrem(dir+7,ni)
439 xy1=xrem(jj ,ni)
440 xy2=xrem(jj+7 ,ni)
441 xz1=xrem(kk ,ni)
442 xz2=xrem(kk+7 ,ni)
443 END IF
446 IF(xmax>=seuil.AND.xmin<=ximax) THEN
447 IF(
i11insid(xx1,xx2,xy1,xy2,xz1,xz2,
448 . seuil,ximax,xjmin,xjmax,xkmin,xkmax)) THEN
449
450 nb_ncn = nb_ncn + 1
451 bpn(nb_ncn) = nn
452 ENDIF
453 ENDIF
454 ENDDO
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481 nb_ecn= 0
482 nb_ecn1= 0
483 addne= add(2,i_add)
484 IF(nb_ncn1==0) THEN
485 DO i=1,nb_ec
486 ne = bpe(i)
487 xx1=x(dir, irectm(1,ne))
488 xx2=x(dir, irectm(2,ne))
489 IF(
max(xx1,xx2)>=tsmin)
THEN
490 xy1=x(jj, irectm(1,ne))
491 xy2=x(jj, irectm(2,ne))
492 xz1=x(kk, irectm(1,ne))
493 xz2=x(kk, irectm(2,ne))
494 IF(
i11insid(xx1,xx2,xy1,xy2,xz1,xz2,
495 . tsmin,timax,tjmin,tjmax,tkmin,tkmax)) THEN
496
497 nb_ecn = nb_ecn + 1
498 bpe(nb_ecn) = ne
499 ENDIF
500 ENDIF
501 ENDDO
502 ELSEIF(nb_ncn==0) THEN
503 DO i=1,nb_ec
504 ne = bpe(i)
505 xx1=x(dir, irectm(1,ne))
506 xx2=x(dir, irectm(2,ne))
507 IF(
min(xx1,xx2)<tsmax)
THEN
508 xy1=x(jj, irectm(1,ne))
509 xy2=x(jj, irectm(2,ne))
510 xz1=x(kk, irectm(1,ne))
511 xz2=x(kk, irectm(2,ne))
512 IF(
i11insid(xx1,xx2,xy1,xy2,xz1,xz2,
513 . timin,tsmax,tjmin,tjmax,tkmin,tkmax)) THEN
514
515 addne = addne + 1
516 nb_ecn1= nb_ecn1 + 1
517 pe(addne) = ne
518 ENDIF
519 ENDIF
520 ENDDO
521 ELSE
522 DO i=1,nb_ec
523 ne = bpe(i)
524 xx1=x(dir, irectm(1,ne))
525 xx2=x(dir, irectm(2,ne))
526 IF(
min(xx1,xx2)<tsmax)
THEN
527 xy1=x(jj, irectm(1,ne))
528 xy2=x(jj, irectm(2,ne))
529 xz1=x(kk, irectm(1,ne))
530 xz2=x(kk, irectm(2,ne))
531 IF(
i11insid(xx1,xx2,xy1,xy2,xz1,xz2,
532 . timin,tsmax,tjmin,tjmax,tkmin,tkmax)) THEN
533
534 addne = addne + 1
535 nb_ecn1= nb_ecn1 + 1
536 pe(addne) = ne
537 ENDIF
538 ENDIF
539 ENDDO
540 DO i=1,nb_ec
541 ne = bpe(i)
542 xx1=x(dir, irectm(1,ne))
543 xx2=x(dir, irectm(2,ne))
544 IF(
max(xx1,xx2)>=tsmin)
THEN
545 xy1=x(jj, irectm(1,ne))
546 xy2=x(jj, irectm(2,ne))
547 xz1=x(kk, irectm(1,ne))
548 xz2=x(kk, irectm(2,ne))
549 IF(
i11insid(xx1,xx2,xy1,xy2,xz1,xz2,
550 . tsmin,timax,tjmin,tjmax,tkmin,tkmax)) THEN
551
552 nb_ecn = nb_ecn + 1
553 bpe(nb_ecn) = ne
554 ENDIF
555 ENDIF
556 ENDDO
557 ENDIF
558
559
560
561 add(1,i_add+1) = addnn
562 add(2,i_add+1) = addne
563
564
565
566
567
568
569 xyzm(1,i_add+1) = xyzm(1,i_add)
570 xyzm(2,i_add+1) = xyzm(2,i_add)
571 xyzm(3,i_add+1) = xyzm(3,i_add)
572 xyzm(4,i_add+1) = xyzm(4,i_add)
573 xyzm(5,i_add+1) = xyzm(5,i_add)
574 xyzm(6,i_add+1) = xyzm(6,i_add)
575
576 xyzm(dir ,i_add) = ximin
577 xyzm(dir+3,i_add) = seuil
578 xyzm(dir ,i_add+1) = seuil
579 xyzm(dir+3,i_add+1) = ximax
580
581 nb_old(1,i_add)=nb_nc
582 nb_old(2,i_add)=nb_ec
583 nb_old(1,i_add+1)=nb_nc
584 nb_old(2,i_add+1)=nb_ec
585
586 nb_nc = nb_ncn
587 nb_ec = nb_ecn
588
589 i_add = i_add + 1
590 IF(i_add+1>=i_add_max) THEN
591 i_mem = 3
592 RETURN
593 ENDIF
594
595 200 CONTINUE
596
597
598
599
600
601
602
603
604
605
606
607
608 IF(add(1,i_add)+nb_nc>maxsiz) THEN
609
610 i_mem = 1
611 RETURN
612 ENDIF
613 IF(add(2,i_add)+nb_ec>maxsiz) THEN
614
615 i_mem = 1
616 RETURN
617 ENDIF
618
619
620
621 IF(nb_ec/=0.AND.nb_nc/=0) THEN
622
623 dx = xyzm(4,i_add) - xyzm(1,i_add)
624 dy = xyzm(5,i_add) - xyzm(2,i_add)
625 dz = xyzm(6,i_add) - xyzm(3,i_add)
627
628
629
630
631
632 IF(nb_ec+nb_nc<=128) THEN
633 ncand_prov = nb_ec*nb_nc
634 ELSE
635 ncand_prov = 129
636 ENDIF
637
638 nb_nc_old = nb_old(1,i_add)
639 nb_ec_old = nb_old(2,i_add)
640
641 IF(dsup<minbox.OR.
642 . nb_nc<=nb_n_b.OR.nb_ec<=nb_n_b.OR.
643 . ncand_prov<=128.OR.(nb_ec==nb_ec_old
644 . .AND.nb_nc==nb_nc_old)) THEN
645
646 ncand_prov = nb_ec*nb_nc
647 DO k=1,ncand_prov,nvsiz
648 DO l=k,
min(k-1+nvsiz,ncand_prov)
649 i = 1+(l-1)/nb_nc
650 j = l-(i-1)*nb_nc
651 ne = bpe(i)
652 nn = bpn(j)
653 n1=irectm(1,ne)
654 n2=irectm(2,ne)
655 IF(nn<=nrts) THEN
656 nn1=irects(1,nn)
657 nn2=irects(2,nn)
658 IF(iauto==0 .OR. itab(n1)>itab(nn1) )THEN
659 IF(nn1/=n1.AND.nn1/=n2.AND.
660 . nn2/=n1.AND.nn2/=n2) THEN
661 j_stok = j_stok + 1
662 prov_s(j_stok) = nn
663 prov_m(j_stok) = ne
664 ENDIF
665 ENDIF
666 ELSE
667 ni = nn-nrts
670 n1 = itab(n1)
671 n2 = itab(n2)
672 IF(iauto==0 .OR. n1>nn1 )THEN
673 IF(nn1/=n1.AND.nn1/=n2.AND.
674 . nn2/=n1.AND.nn2/=n2) THEN
675 j_stok = j_stok + 1
676 prov_s(j_stok) = nn
677 prov_m(j_stok) = ne
678 ENDIF
679 ENDIF
680 END IF
681 ENDDO
682 IF(j_stok>=nvsiz)THEN
684 1 nvsiz,irects,irectm,x ,ii_stok,
685 2 cand_s,cand_m,nsn4 ,noint ,tzinf ,
686 3 i_mem ,prov_s,prov_m,eshift,addcm ,
687 4 chaine,nrts ,itab ,ifpen ,iform )
688 IF(i_mem==2)RETURN
689 j_stok = j_stok-nvsiz
690#include "vectorize.inc"
691 DO j=1,j_stok
692 prov_s(j) = prov_s(j+nvsiz)
693 prov_m(j) = prov_m(j+nvsiz)
694 ENDDO
695 ENDIF
696 ENDDO
697 ELSE
698
699 GOTO 100
700
701 ENDIF
702 ENDIF
703
704
705
706
707
708 i_add = i_add - 1
709 IF (i_add/=0) THEN
710
711
712
713
714 CALL i7dstk(nb_nc,nb_ec,add(1,i_add),bpn,pn,bpe,pe)
715
716 GOTO 200
717
718 ENDIF
719
720
721
723 1 j_stok,irects,irectm,x ,ii_stok,
724 2 cand_s,cand_m,nsn4 ,noint ,tzinf ,
725 3 i_mem ,prov_s,prov_m,eshift,addcm ,
726 4 chaine,nrts ,itab ,ifpen ,iform )
727
728 RETURN
integer, dimension(:,:), allocatable irem
subroutine i11sto(j_stok, irects, irectm, x, ii_stok, cand_n, cand_e, nsn, noint, tzinf, i_mem, prov_n, prov_e, multimp, addcm, chaine, iadfin)
subroutine i7dstk(i_add, nb_nc, nb_ec, add, bpn, pn, bpe, pe)