43
44
45
48
49
50
51#include "implicit_f.inc"
52
53
54
55#include "mvsiz_p.inc"
56
57 INTEGER NVECSZ ,NSP2SORTF, NSP2SORTL
58 parameter(nvecsz = mvsiz)
59
60
61
62#include "com01_c.inc"
63#include "com04_c.inc"
64#include "param_c.inc"
65#include "task_c.inc"
66#include "sphcom.inc"
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 INTEGER NSN,ITASK,NSNR,NBX,NBY,NBZ,
92 . NLIST(*),NOD2SP(*) ,
93 . VOXEL(NBX+2,NBY+2,NBZ+2),JVOIS(*) ,JSTOR(*), JPERM(*) ,
94 . IREDUCE,KXSP(NISP,*), IXSP(KVOISPH,*), KREDUCE(*),
95 . LGAUGE(3,*),KXSPR(*),(KVOISPH,*)
96
98 . x(3,*),bminma(12),
99 . marge ,spbuf(nspbuf,*), dvois(*), gauge(llgauge,*)
100
101
102
103 INTEGER NB_NCN,NB_NCN1,NB_ECN,I,J,DIR,NB_NC,NB_EC,
104 . N1,N2,N3,N4,NN,NE,K,L,II,JJ,JS,NS,N,
105 . NSNF, NSNL,NVOIS, IG, IL
106
108 . dx,dy,dz,xs,ys,zs,xx,sx,sy,sz,s2,xn,yn,zn,
109 . xmin, xmax,ymin,
ymax,zmin, zmax, tz,
110 . d1x,d1y,d1z,d2,a2,alpha_marge,distmax
111 my_real,
DIMENSION(:),
ALLOCATABLE :: tab_dk
112
113 INTEGER LAST_NOD(NSN+NSNR)
114 INTEGER IX,IY,IZ,NEXT,
115 . IX1,IY1,IZ1,IX2,IY2,IZ2
116 INTEGER, DIMENSION(:),ALLOCATABLE :: IIX,IIY,IIZ
118 . xminb,yminb,zminb,xmaxb,ymaxb,zmaxb,
119 . xmine,ymine,zmine,xmaxe,ymaxe,zmaxe,aaa,bbb,
120 . aaa2
121 INTEGER FIRST,NEW,LAST,REQ_RECV(NSPMD)
122 SAVE iix,iiy,iiz,distmax,tab_dk
123 INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: NEXT_NOD_LOCAL
124 INTEGER, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: VOXEL_LOCAL
125
126 IF(itask == 0)THEN
128 ALLOCATE(iix(nsn+nsnr))
129 ALLOCATE(iiy(nsn+nsnr))
130 ALLOCATE(iiz(nsn+nsnr))
131 ALLOCATE(tab_dk(numsph))
132 ALLOCATE(next_nod_local(nsn))
133 ALLOCATE( voxel_local(nbx+2,nby+2,nbz+2) )
134 END IF
135
137
138
139 alpha_marge = sqrt(one +spasort)
140
141 xmin = bminma(4)
142 ymin = bminma(5)
143 zmin = bminma(6)
144 xmax = bminma(1)
146 zmax = bminma(3)
147
148 xminb = bminma(10)
149 yminb = bminma(11)
150 zminb = bminma(12)
151 xmaxb = bminma(7)
152 ymaxb = bminma(8)
153 zmaxb = bminma(9)
154
155
156
157
158 IF(itask == 0)THEN
159
160 distmax = zero
161
162 DO i=1,nsn
163 iix(i)=0
164 iiy(i)=0
165 iiz(i)=0
166
167 j=nlist(i)
168 nn = kxsp(3,j)
169
170 distmax =
max(distmax,spbuf(1,j))
171
172
173 IF(x(1,nn) < xmin) cycle
174 IF(x(1,nn) > xmax) cycle
175 IF(x(2,nn) < ymin) cycle
176 IF(x(2,nn) >
ymax) cycle
177 IF(x(3,nn) < zmin) cycle
178 IF(x(3,nn) > zmax) cycle
179
180 iix(i)=int(nbx*(x(1,nn)-xminb)/(xmaxb-xminb))
181 iiy(i)=int(nby*(x(2,nn)-yminb)/(ymaxb-yminb))
182 iiz(i)=int(nbz*(x(3,nn)-zminb)/(zmaxb-zminb))
183
184 iix(i)=
max(1,2+
min(nbx,iix(i)))
185 iiy(i)=
max(1,2+
min(nby,iiy(i)))
186 iiz(i)=
max(1,2+
min(nbz,iiz(i)))
187
188 first = voxel(iix(i),iiy(i),iiz(i))
189 IF(first == 0)THEN
190
191 voxel(iix(i),iiy(i),iiz(i)) = i
193 last_nod(i) = 0
194 ELSEIF(last_nod(first) == 0)THEN
195
196
198 last_nod(first) = i
200 ELSE
201
202
203 last = last_nod(first)
205 last_nod(first) = i
207 ENDIF
208 ENDDO
209 voxel_local(1:nbx+2,1:nby+2,1:nbz+2) = voxel(1:nbx+2,1:nby+2,1:nbz+2)
210 next_nod_local(1:nsn) =
next_nod(1:nsn)
211
212
213
214
215 DO j = 1, nsnr
216 i = j+nsn
217 n = nlist(i)-numsph
218 distmax =
max(distmax,xsphr(2,n))
219 iix(i)=int(nbx*(xsphr(3,n)-xminb)/(xmaxb-xminb))
220 iiy(i)=int(nby*(xsphr(4,n)-yminb)/(ymaxb-yminb))
221 iiz(i)=int(nbz*(xsphr(5,n)-zminb)/(zmaxb-zminb))
222 iix(i)=
max(1,2+
min(nbx,iix(i)))
223 iiy(i)=
max(1,2+
min(nby,iiy(i)))
224 iiz(i)=
max(1,2+
min(nbz,iiz(i)))
225
226 first = voxel(iix(i),iiy(i),iiz(i))
227 IF(first == 0)THEN
228
229 voxel(iix(i),iiy(i),iiz(i)) = i
231 last_nod(i) = 0
232 ELSEIF(last_nod(first) == 0)THEN
233
234
236 last_nod(first) = i
238 ELSE
239
240
241 last = last_nod(first)
243 last_nod(first) = i
245 ENDIF
246 ENDDO
247 END IF
248
249
250
252 IF(itask == 0 ) THEN
254 ENDIF
256
257
258
259
260
261 DO ne = nsp2sortf,nsp2sortl
262
263
264
265 j=nlist(ne)
266 nn = kxsp(3,j)
267
268
269 aaa = two*spbuf(1,j)* alpha_marge
270
271
272
273 ix1=int(nbx*(x(1,nn)-aaa-xminb)/(xmaxb-xminb))
274 iy1=int(nby*(x(2,nn)-aaa-yminb)/(ymaxb-yminb))
275 iz1=int(nbz*(x(3,nn)-aaa-zminb)/(zmaxb-zminb))
276
280
281 ix2=int(nbx*(x(1,nn)+aaa-xminb)/(xmaxb-xminb))
282 iy2=int(nby*(x(2,nn)+aaa-yminb)/(ymaxb-yminb))
283 iz2=int(nbz*(x(3,nn)+aaa-zminb)/(zmaxb-zminb))
284
288
289
290
291
292
293
294 DO iz = iz1,iz2
295 DO iy = iy1,iy2
296 DO ix = ix1,ix2
297
298
299
300 jj = voxel(ix,iy,iz)
301
302 DO WHILE(jj /= 0)
303
304
305
306 js=nlist(jj)
307 IF(jj<=nsn)THEN
308 IF(spbuf(1,js) > spbuf(1,j) .OR.
309 . (spbuf(1,js) == spbuf(1,j) .AND. kxsp(8,js)>=kxsp(8,j)))GOTO 200
310
311 ns=kxsp(3,js)
312 xs = x(1,ns)
313 ys = x(2,ns)
314 zs = x(3,ns)
315
316 aaa = spbuf(1,j)+spbuf(1,js)
317 ELSE
318 n = nlist(jj)-numsph
319 IF(xsphr(2,n) > spbuf(1,j).OR.
320 . (xsphr(2,n) == spbuf(1,j) .AND.(nint(xsphr(6,n))>=kxsp(8,j) ) ) ) GOTO 200
321
322 xs = xsphr(3,n)
323 ys = xsphr(4,n)
324 zs = xsphr(5,n)
325 aaa = spbuf(1,j)+xsphr(2,n)
326 ENDIF
327
328 bbb = aaa * alpha_marge
329
330 IF(xs<=x(1,nn)-bbb)GOTO 200
331 IF(xs>=x(1,nn)+bbb)GOTO 200
332 IF(ys<=x(2,nn)-bbb)GOTO 200
333 IF(ys>=x(2,nn)+bbb)GOTO 200
334 IF(zs<=x(3,nn)-bbb)GOTO 200
335 IF(zs>=x(3,nn)+bbb)GOTO 200
336
337
338
339
340 d1x = xs - x(1,nn)
341 d1y = ys - x(2,nn)
342 d1z = zs - x(3,nn)
343 d2 = d1x*d1x+d1y*d1y+d1z*d1z
344 a2 = bbb*bbb
345 IF(d2 > a2)GOTO 200
346
347
348 aaa2 = aaa*aaa
349 nvois=kxsp(5,j)+1
350 jvois(nvois)=js
351 dvois(nvois)=d2/aaa2
352
353 kxsp(5,j) =nvois
354
355 200 CONTINUE
356
358
359 ENDDO
360
361 ENDDO
362 ENDDO
363 ENDDO
364
365
366
367
368
369 CALL sppro3(j ,kxsp ,ixsp ,nod2sp,jvois,
370 . jstor,jperm ,dvois,ireduce,kreduce,
371 . kxspr,ixspr,tab_dk)
372
373 ENDDO
374
375
376
377
379 IF(itask == 0 ) THEN
381 ENDIF
383
384
385
386
387
388 DO j = itask+1, nsnr, nthread
389
390 i = j+nsn
391 n = nlist(i)-numsph
392
393 aaa = two * xsphr(2,n) * alpha_marge
394
395
396
397 ix1=int(nbx*(xsphr(3,n)-aaa-xminb)/(xmaxb-xminb))
398 iy1=int(nby*(xsphr(4,n)-aaa-yminb)/(ymaxb-yminb))
399 iz1=int(nbz*(xsphr(5,n)-aaa-zminb)/(zmaxb-zminb))
400
404
405 ix2=int(nbx*(xsphr(3,n)+aaa-xminb)/(xmaxb-xminb))
406 iy2=int(nby*(xsphr(4,n)+aaa-yminb)/(ymaxb-yminb))
407 iz2=int(nbz*(xsphr(5,n)+aaa-zminb)/(zmaxb-zminb))
408
412
413
414
415
416
417
418 DO iz = iz1,iz2
419 DO iy = iy1,iy2
420 DO ix = ix1,ix2
421
422
423
424 jj = voxel_local(ix,iy,iz)
425
426 DO WHILE(jj /= 0)
427
428
429
430 js=nlist(jj)
431 IF(jj<=nsn)THEN
432 IF(xsphr(2,n) < spbuf(1,js).OR.
433 . (xsphr(2,n)==spbuf(1,js).AND.nint(xsphr(6,n))<kxsp(8,js)) )GOTO 250
434
435 ns=kxsp(3,js)
436 xs = x(1,ns)
437 ys = x(2,ns)
438 zs = x(3,ns)
439
440 aaa = xsphr(2,n)+spbuf(1,js)
441
442
443 bbb = aaa * alpha_marge
444
445 IF(xs<=xsphr(3,n)-bbb)GOTO 250
446 IF(xs>=xsphr(3,n)+bbb)GOTO 250
447 IF(ys<=xsphr(4,n)-bbb)GOTO 250
448 IF(ys>=xsphr(4,n)+bbb)GOTO 250
449 IF(zs<=xsphr(5,n)-bbb)GOTO 250
450 IF(zs>=xsphr(5,n)+bbb)GOTO 250
451
452
453
454
455 d1x = xsphr(3,n) - xs
456 d1y = xsphr(4,n) - ys
457 d1z = xsphr(5,n) - zs
458 d2 = d1x*d1x+d1y*d1y+d1z*d1z
459 a2 = bbb*bbb
460 IF(d2 > a2)GOTO 250
461
462
463 aaa2 = aaa*aaa
464 nvois=kxspr(j)+1
465 jvois(nvois)=js
466 dvois(nvois)=d2/aaa2
467
468 kxspr(j) =nvois
469
470 250 CONTINUE
471 END IF
472
473 jj = next_nod_local(jj)
474
475 ENDDO
476
477 ENDDO
478 ENDDO
479 ENDDO
480
481
482
483
484
485 CALL sppro3(j+numsph,kxsp ,ixsp ,nod2sp,jvois,
486 . jstor,jperm ,dvois,ireduce,kreduce,
487 . kxspr,ixspr,tab_dk)
488
489 ENDDO
490
491
492
493
494
495 DO ig = 1,nbgauge
496
497 IF(lgauge(1,ig) > -(numels+1))cycle
498
499 j = numsph+ig
500 xn=gauge(2,ig)
501 yn=gauge(3,ig)
502 zn=gauge(4,ig)
503
504 aaa = (two*distmax)* alpha_marge
505
506
507
508 ix1=int(nbx*(xn-aaa-xminb)/(xmaxb-xminb))
509 iy1=int(nby*(yn-aaa-yminb)/(ymaxb-yminb))
510 iz1=int(nbz*(zn-aaa-zminb)/(zmaxb-zminb))
511
515
516 ix2=int(nbx*(xn+aaa-xminb)/(xmaxb-xminb))
517 iy2=int(nby*(yn+aaa-yminb)/(ymaxb-yminb
518 iz2=int(nbz*(zn+aaa-zminb)/(zmaxb-zminb))
519
523
524
525
526
527
528
529 DO iz = iz1,iz2
530 DO iy = iy1,iy2
531 DO ix = ix1,ix2
532
533
534
535 jj = voxel_local(ix,iy,iz)
536
537 DO WHILE(jj /= 0)
538
539
540
541 js=nlist(jj)
542 IF(jj<=nsn)THEN
543 ns=kxsp(3,js)
544 xs = x(1,ns)
545 ys = x(2,ns)
546 zs = x(3,ns)
547
548 aaa = two*spbuf(1,js)
549 ELSE
550 GOTO 300
551 ENDIF
552
553 bbb = aaa * alpha_marge
554
555 IF(xs<=xn-bbb)GOTO 300
556 IF(xs>=xn+bbb)GOTO 300
557 IF(ys<=yn-bbb)GOTO 300
558 IF(ys>=yn+bbb)GOTO 300
559 IF(zs<=zn-bbb)GOTO 300
560 IF(zs>=zn+bbb)GOTO 300
561
562
563
564 d1x = xs - xn
565 d1y = ys - yn
566 d1z = zs - zn
567 d2 = d1x*d1x+d1y*d1y+d1z*d1z
568 a2 = bbb*bbb
569 IF(d2 > a2)GOTO 300
570
571
572 aaa2 = aaa*aaa
573 nvois=kxsp(5,j)+1
574 jvois(nvois)=js
575 dvois(nvois)=d2/aaa2
576
577 kxsp(5,j) =nvois
578
579 300 CONTINUE
580
581 jj = next_nod_local(jj)
582
583 ENDDO
584
585 ENDDO
586 ENDDO
587 ENDDO
588
589
590
591
592
593 il=-j
594 CALL sppro3(il ,kxsp ,ixsp ,nod2sp,jvois,
595 . jstor,jperm ,dvois,ireduce,kreduce,
596 . kxspr,ixspr, tab_dk)
597
598 ENDDO
599
600
601
602
603
604
605
606
607 100 CONTINUE
608
609
611
612 DO i=nsp2sortf,nsp2sortl
613 IF(iix(i)/=0)THEN
614 voxel(iix(i),iiy(i),iiz(i))=0
615 ENDIF
616 ENDDO
617
618
619
620
621 nsnf = 1 + itask*nsnr / nthread
622 nsnl = (itask+1)*nsnr / nthread
623 DO j = nsnf, nsnl
624 voxel(iix(nsn+j),iiy(nsn+j),iiz(nsn+j))=0
625 ENDDO
626
627
629 IF(itask == 0)THEN
631 DEALLOCATE(iix)
632 DEALLOCATE(iiy)
633 DEALLOCATE(iiz)
634 DEALLOCATE(tab_dk)
635 DEALLOCATE( next_nod_local )
636 DEALLOCATE( voxel_local )
637 ENDIF
638
639 RETURN
subroutine ymax(idn, fac, npc, pld, stiffmin, stiffmax, stiffini, stiffavg)
integer, dimension(:), allocatable next_nod
subroutine spmd_sphgetdk(tab_dk, act, req_recv)
subroutine sppro3(il, kxsp, ixsp, nod2sp, jvois, jstor, jperm, dvois, ireduce, kreduce, kxspr, ixspr, tab_dk)