48
49
50
51 USE my_alloc_mod
54 use inter_save_candidate_mod , only : inter_save_candidate
56 use intbufdef_mod
57
58
59
60#include "implicit_f.inc"
61
62
63
64#include "mvsiz_p.inc"
65
66 INTEGER NVECSZ
67 parameter(nvecsz = mvsiz)
68
69
70
71#include "com04_c.inc"
72#include "param_c.inc"
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 integer, intent(in) :: nrtm
114 integer, intent(in) :: nrtm_l
115 INTEGER I_MEM,NSN,
116 . MULNSN,NOINT,IGAP,NBX,NBY,NBZ,IREMNODE,FLAGREMNODE,
117 . NSV(NSN),
118 . IRECT(4,NRTM), VOXEL(NBX+2,NBY+2,NBZ+2),ISTF,
119 . I_STOK ,J_STOK,
120 . INDEX(*),KREMNODE(*),REMNODE(*)
122 . x(3,*),xyzm(6,2),stf(*),stfn(*),gap_s(*),gap_m(*),
123 . tzinf,marge,gap,gapmin,gapmax,bgapsmx,drad,
124 . curv_max(*),gap_s_l(*),gap_m_l(*)
125 INTEGER ID
126 CHARACTER(LEN=NCHARTITLE)::TITR
127 integer, intent(in) :: nin
128 integer, dimension(npari), intent(inout) :: ipari
129 type(intbuf_struct_), intent(inout) :: intbuf_tab
130 INTEGER, dimension(nsn), intent(inout) :: iix,iiy,iiz,local_next_nod
131 LOGICAL,INTENT(IN) :: IS_USED_WITH_LAW151
132
133
134
135 INTEGER NB_NCN,NB_NCN1,NB_ECN,I,J,DIR,NB_NC,NB_EC,
136 . N1,N2,N3,N4,NN,NE,K,L,NCAND_PROV,II,JJ,KK,
137 . NSNF, NSNL, ,DELNOD
139 . dx,dy,dz,xs,ys,zs,xx,sx,sy,sz,s2,
140 . xmin, xmax,ymin,
ymax,zmin, zmax, tz, gapsmx, gapl,
141 . xx1,xx2,xx3,xx4,yy1,yy2,yy3,yy4,zz1,zz2,zz3,zz4,
142 . d1x,d1y,d1z,d2x,d2y,d2z,dd1,dd2,d2,a2,gs,gapv(mvsiz)
143 INTEGER LAST_NOD(NSN)
144 INTEGER*8 meanjj8
145 INTEGER IX,IY,IZ,NEXT,M1,M2,M3,M4,
146 . IX1,IY1,IZ1,IX2,IY2,IZ2
147 integer, dimension(mvsiz) :: prov_n,prov_e
149 . xminb,yminb,zminb,xmaxb,ymaxb,zmaxb,
150 . xmine,ymine,zmine,xmaxe,ymaxe,zmaxe,aaa,tstart,tstop
151 my_real ,
INTENT(IN) :: dgapload
152 INTEGER FIRST,NEW,LAST,M
153 INTEGER, DIMENSION(MVSIZ) :: IX11,IX12,IX13,IX14,NSVG
154 my_real,
DIMENSION(MVSIZ) :: x1,x2,x3,x4
155 my_real,
DIMENSION(MVSIZ) :: y1,y2,y3,y4
156 my_real,
DIMENSION(MVSIZ) :: z1,z2,z3,z4
157 my_real,
DIMENSION(MVSIZ) :: xi,yi,zi
158 my_real,
DIMENSION(MVSIZ) :: x0,y0,z0
159 my_real,
DIMENSION(MVSIZ) :: nx1,ny1,nz1
160 my_real,
DIMENSION(MVSIZ) :: nx2,ny2,nz2
161 my_real,
DIMENSION(MVSIZ) :: nx3,ny3,nz3
162 my_real,
DIMENSION(MVSIZ) :: nx4,ny4,nz4
163 my_real,
DIMENSION(MVSIZ) :: p1,p2,p3,p4
164 my_real,
DIMENSION(MVSIZ) :: lb1,lb2,lb3,lb4
165 my_real,
DIMENSION(MVSIZ) :: lc1,lc2,lc3,lc4
166 my_real,
DIMENSION(MVSIZ) :: n11,n21,n31
167 my_real,
DIMENSION(MVSIZ) :: stif,pene
168
169 INTEGER, DIMENSION(:), ALLOCATABLE :: TAGNOD
170
171 integer , external :: omp_get_thread_num,omp_get_num_threads
172 integer :: itask,nthreads
173 integer :: my_old_size,my_address
174 integer :: local_i_stok,multimp
175 integer :: local_cand_n_size,
176
177 integer, save :: i_stok_old
178 integer, dimension(:), allocatable, save :: cand_n_size,cand_e_size
179 integer, dimension(:), allocatable, save :: address_cand_n,address_cand_e
180 type(array_type_int_1d) :: local_cand_n
181 type(array_type_int_1d) :: local_cand_e
182
183 integer :: my_size,mode
184 integer, dimension(:), allocatable :: my_index
185 integer, dimension(:,:), allocatable :: sort_array,save_array
186 integer, dimension(70000) :: work
187 integer stagnod
188
189
190 itask = omp_get_thread_num()
191 nthreads = omp_get_num_threads()
192 local_cand_n_size = size(intbuf_tab%cand_n) / nthreads
193 local_cand_e_size = size(intbuf_tab%cand_e) / nthreads
194 local_i_stok = 0
195 local_cand_n%size_int_array_1d = local_cand_n_size
196 local_cand_e%size_int_array_1d = local_cand_e_size
199
200 xmin = xyzm(1,1)
201 ymin = xyzm(2,1)
202 zmin = xyzm(3,1)
203 xmax = xyzm(4,1)
205 zmax = xyzm(6,1)
206
207 xminb = xyzm(1,2)
208 yminb = xyzm(2,2)
209 zminb = xyzm(3,2)
210 xmaxb = xyzm(4,2)
211 ymaxb = xyzm(5,2)
212 zmaxb = xyzm(6,2)
213
214
215
216
217
218
219
220 allocate( cand_n_size(nthreads+1),cand_e_size(nthreads+1) )
221 allocate( address_cand_n(nthreads+1),address_cand_e(nthreads+1) )
222 cand_n_size(1:nthreads+1) = 0
223 cand_e_size(1:nthreads+1) = 0
224 address_cand_n(1:nthreads+1) = 0
225 address_cand_e(1:nthreads+1) = 0
226
227 DO i=1,nsn
228 iix(i)=0
229 iiy(i)=0
230 iiz(i)=0
231 IF(stfn(i) == zero)cycle
232 j=nsv(i)
233
234
235 ix=int(
lrvoxel*(x(1,j)-xmin)/(xmax-xmin))
236 IF(ix < 0 .OR. ix >
lrvoxel) cycle
238 IF(iy < 0 .OR. iy >
lrvoxel) cycle
239 iz=int(
lrvoxel*(x(3,j)-zmin)/(zmax-zmin))
240 IF(iz < 0 .OR. iz >
lrvoxel) cycle
241 IF(.NOT.(btest(
crvoxel(iy,iz),ix))) cycle
242
243
244 IF( (x(1,j)-xminb)/(xmaxb-xminb) > one ) THEN
245 iix(i) = nbx
246 ELSE
247 iix(i)=int(
max(nbx*(x(1,j)-xminb)/(xmaxb-xminb),-one))
248 ENDIF
249 IF( (x(2,j)-yminb)/(ymaxb-yminb) > one ) THEN
250 iiy(i) = nby
251 ELSE
252 iiy(i)=int(
max(nby*(x(2,j)-yminb)/(ymaxb-yminb),-one))
253 ENDIF
254 IF( (x(3,j)-zminb)/(zmaxb-zminb) > one ) THEN
255 iiz(i) = nbz
256 ELSE
257 iiz(i)=int(
max(nbz*(x(3,j)-zminb)/(zmaxb-zminb),-one))
258 ENDIF
259
260
261 iix(i)=
max(1,2+iix(i))
262 iiy(i)=
max(1,2+iiy(i))
263 iiz(i)=
max(1,2+iiz(i))
264
265 first = voxel(iix(i),iiy(i),iiz(i))
266 IF(first == 0)THEN
267
268 voxel(iix(i),iiy(i),iiz(i)) = i
269 local_next_nod(i) = 0
270 last_nod(i) = 0
271 ELSEIF(last_nod(first) == 0)THEN
272
273
274 local_next_nod(first) = i
275 last_nod(first) = i
276 local_next_nod(i) = 0
277 ELSE
278
279
280 last = last_nod(first)
281 local_next_nod(last) = i
282 last_nod(first) = i
283 local_next_nod(i) = 0
284 ENDIF
285 ENDDO
286
287
288
289
290
291
292
293
294 stagnod = numnod+numfakenodigeo
295 IF(is_used_with_law151) stagnod = stagnod + numels
297
298 local_i_stok = 0
299 j_stok = 0
300
301 DO kk=1,nrtm_l
302 ne=index(kk)
303 IF(stf(ne) == zero)cycle
304
305 IF(flagremnode==2.AND.iremnode==2)THEN
306 k = kremnode(ne)+1
307 l = kremnode(ne+1)
308 DO m=k,l
310 ENDDO
311 ENDIF
312
313 IF(igap == 0)THEN
314 aaa = tzinf+curv_max(ne)
315 ELSE
316 aaa = marge+curv_max(ne)+
max(
min(gapmax,
max(gapmin,bgapsmx+gap_m(ne)))+dgapload,drad)
317 ENDIF
318
319 m1 = irect(1,ne)
320 m2 = irect(2,ne)
321 m3 = irect(3,ne)
322 m4 = irect(4,ne)
323
324 xx1=x(1,m1)
325 xx2=x(1,m2)
326 xx3=x(1,m3)
327 xx4=x(1,m4)
328 xmaxe=
max(xx1,xx2,xx3,xx4)
329 xmine=
min(xx1,xx2,xx3,xx4)
330
331 yy1=x(2,m1)
332 yy2=x(2,m2)
333 yy3=x(2,m3)
334 yy4=x(2,m4)
335 ymaxe=
max(yy1,yy2,yy3,yy4)
336 ymine=
min(yy1,yy2,yy3,yy4)
337
338 zz1=x(3,m1)
339 zz2=x(3,m2)
340 zz3=x(3,m3)
341 zz4=x(3,m4)
342 zmaxe=
max(zz1,zz2,zz3,zz4)
343 zmine=
min(zz1,zz2,zz3,zz4)
344
345
346
347
348 sx = (yy3-yy1)*(zz4-zz2) - (zz3-zz1)*(yy4-yy2)
349 sy = (zz3-zz1)*(xx4-xx2) - (xx3-xx1)*(zz4-zz2)
350 sz = (xx3-xx1)*(yy4-yy2) - (yy3-yy1)*(xx4-xx2)
351 s2 = sx*sx + sy*sy + sz*sz
352
353
354 IF( (xmine - aaa - xminb)/(xmaxb-xminb) > one ) THEN
355 ix1 = nbx
356 ELSE
357 ix1=int(
max(nbx*(xmine-aaa-xminb)/(xmaxb-xminb),-one))
358 ENDIF
359 IF( (ymine - aaa - yminb)/(ymaxb-yminb) > one ) THEN
360 iy1 = nby
361 ELSE
362 iy1=int(
max(nby*(ymine-aaa-yminb)/(ymaxb-yminb),-one))
363 ENDIF
364 IF( (zmine - aaa - zminb)/(zmaxb-zminb) > one ) THEN
365 iz1 = nbz
366 ELSE
367 iz1=int(
max(nbz*(zmine-aaa-zminb)/(zmaxb-zminb),-one))
368 ENDIF
369
373
374 IF( (xmaxe + aaa - xminb)/(xmaxb-xminb) > one ) THEN
375 ix2 = nbx
376 ELSE
377 ix2=int(
max(nbx*(xmaxe+aaa-xminb)/(xmaxb-xminb),-one))
378 ENDIF
379 IF( (ymaxe + aaa - yminb)/(ymaxb-yminb) > one ) THEN
380 iy2 = nby
381 ELSE
382 iy2=int(
max(nby*(ymaxe+aaa-yminb)/(ymaxb-yminb),-one))
383 ENDIF
384 IF( (zmaxe + aaa - zminb)/(zmaxb-zminb) > one ) THEN
385 iz2 = nbz
386 ELSE
387 iz2=int(
max(nbz*(zmaxe+aaa-zminb)/(zmaxb-zminb),-one))
388 ENDIF
389
393
394 DO iz = iz1,iz2
395 DO iy = iy1,iy2
396 DO ix = ix1,ix2
397
398 jj = voxel(ix,iy,iz)
399
400 DO WHILE(jj /= 0)
401
402 nn=nsv(jj)
403 IF(
tagnod(nn) == 1 )
GOTO 300
404 IF(nn == m1)GOTO 300
405 IF(nn == m2)GOTO 300
406 IF(nn == m3)GOTO 300
407 IF(nn == m4)GOTO 300
408
409 xs = x(1,nn)
410 ys = x(2,nn)
411 zs = x(3,nn)
412 IF(igap /= 0)THEN
413 aaa = marge+curv_max(ne)+
max(
min(gapmax,
max(gapmin,gap_s(jj)+gap_m(ne)))+dgapload,drad)
414 ENDIF
415 IF(xs<=xmine-aaa)GOTO 300
416 IF(xs>=xmaxe+aaa)GOTO 300
417 IF(ys<=ymine-aaa)GOTO 300
418 IF(ys>=ymaxe+aaa)GOTO 300
419 IF(zs<=zmine-aaa)GOTO 300
420 IF(zs>=zmaxe+aaa)GOTO 300
421
422
423
424 d1x = xs - xx1
425 d1y = ys - yy1
426 d1z = zs - zz1
427 d2x = xs - xx2
428 d2y = ys - yy2
429 d2z = zs - zz2
430 dd1 = d1x*sx+d1y*sy+d1z*sz
431 dd2 = d2x*sx+d2y*sy+d2z*sz
432 IF(dd1*dd2 > zero)THEN
433 d2 =
min(dd1*dd1,dd2*dd2)
434 a2 = aaa*aaa*s2
435 IF(d2 > a2)GOTO 300
436 ENDIF
437 j_stok = j_stok + 1
438 prov_n(j_stok) = jj
439 prov_e(j_stok) = ne
440 IF(j_stok == nvsiz) THEN
441 CALL i7cor3(x ,irect,nsv ,prov_e ,prov_n,
442 . stf ,stfn ,gapv ,igap ,gap ,
443 . gap_s,gap_m,istf ,gapmin ,gapmax,
444 . gap_s_l,gap_m_l ,drad ,ix11,ix12,
445 5 ix13 ,ix14 ,nsvg,x1 ,x2 ,
446 6 x3 ,x4 ,y1 ,y2 ,y3 ,
447 7 y4 ,z1 ,z2 ,z3 ,z4 ,
448 8 xi ,yi ,zi ,stif ,dgapload,
449 9 j_stok)
450 CALL i7dst3(ix13,ix14,x1 ,x2 ,x3 ,
451 1 x4 ,y1 ,y2 ,y3 ,y4 ,
452 2 z1 ,z2 ,z3 ,z4 ,xi ,
453 3 yi ,zi ,x0 ,y0 ,z0 ,
454 4 nx1,ny1,nz1,nx2,ny2,
455 5 nz2,nx3,ny3,nz3,nx4,
456 6 ny4,nz4,p1 ,p2 ,p3 ,
457 7 p4 ,lb1,lb2,lb3,lb4,
458 8 lc1,lc2,lc3,lc4,j_stok)
459 CALL i7pen3(marge,gapv,n11,n21,n31,
460 1 pene ,nx1 ,ny1,nz1,nx2,
461 2 ny2 ,nz2 ,nx3,ny3,nz3,
462 3 nx4 ,ny4 ,nz4,p1 ,p2 ,
463 4 p3 ,p4,j_stok)
464 call inter_save_candidate( local_i_stok,j_stok,prov_n,prov_e,pene,local_cand_n,local_cand_e )
465 j_stok = 0
466 ENDIF
467
468 300 CONTINUE
469 jj = local_next_nod(jj)
470 ENDDO
471 ENDDO
472 ENDDO
473 ENDDO
474
475 IF(flagremnode==2.AND.iremnode==2)THEN
476 k = kremnode(ne)+1
477 l = kremnode(ne+1)
478 DO m=k,l
480 ENDDO
481 ENDIF
482
483 ENDDO
484
486
487 IF(j_stok/=0) THEN
488 CALL i7cor3(x ,irect,nsv ,prov_e ,prov_n,
489 . stf ,stfn ,gapv ,igap ,gap ,
490 . gap_s,gap_m,istf ,gapmin ,gapmax,
491 . gap_s_l,gap_m_l ,drad ,ix11,ix12,
492 5 ix13 ,ix14 ,nsvg,x1 ,x2 ,
493 6 x3 ,x4 ,y1 ,y2 ,y3 ,
494 7 y4 ,z1 ,z2 ,z3 ,z4 ,
495 8 xi ,yi ,zi ,stif ,dgapload,
496 9 j_stok)
497 CALL i7dst3(ix13,ix14,x1 ,x2 ,x3 ,
498 1 x4 ,y1 ,y2 ,y3 ,y4 ,
499 2 z1 ,z2 ,z3 ,z4 ,xi ,
500 3 yi ,zi ,x0 ,y0 ,z0 ,
501 4 nx1,ny1,nz1,nx2,ny2,
502 5 nz2,nx3,ny3,nz3,nx4,
503 6 ny4,nz4,p1 ,p2 ,p3 ,
504 7 p4 ,lb1,lb2,lb3,lb4,
505 8 lc1,lc2,lc3,lc4,j_stok)
506 CALL i7pen3(marge,gapv,n11,n21,n31,
507 1 pene ,nx1 ,ny1,nz1,nx2,
508 2 ny2 ,nz2 ,nx3,ny3,nz3,
509 3 nx4 ,ny4 ,nz4,p1 ,p2 ,
510 4 p3 ,p4,j_stok)
511 call inter_save_candidate( local_i_stok,j_stok,prov_n,prov_e,pene,local_cand_n,local_cand_e )
512 j_stok = 0
513 ENDIF
514
515
516 cand_n_size(itask+1) = local_i_stok
517 cand_e_size(itask+1) = local_i_stok
518
519
520
521
522 ! compute
the address
for each threads &
the total number of candidates
523 address_cand_n(1) = 0
524 address_cand_e(1) = 0
525
526 do i=1,nthreads
527 address_cand_n(i+1) = cand_n_size(i) + address_cand_n(i)
528 address_cand_e(i+1) = cand_e_size(i) + address_cand_e(i)
529
530 cand_n_size(nthreads+1) = cand_n_size(nthreads+1) + cand_n_size(i)
531 cand_e_size(nthreads+1) = cand_e_size(nthreads+1) + cand_e_size(i)
532 enddo
533
534
535
536
537 my_old_size = ipari(18)*ipari(23)
538 i_stok_old = i_stok
539 i_stok = i_stok + cand_n_size(nthreads+1)
540 if(i_stok > my_old_size) then
541 multimp = i_stok/ipari(18) + 1
543 endif
544
545
546
547
548
549
550 my_address = i_stok_old + address_cand_n(itask+1)
551 intbuf_tab%cand_n(my_address+1:my_address+local_i_stok) = local_cand_n%int_array_1d(1:local_i_stok)
552 my_address = i_stok_old + address_cand_e(itask+1)
553 intbuf_tab%cand_e(my_address+1:my_address+local_i_stok) = local_cand_e%int_array_1d(1:local_i_stok)
554
555
556
559
560
561
562
563
564
565 my_size = cand_n_size(nthreads+1)
566 allocate(my_index(2*my_size))
567 allocate(sort_array(2,my_size))
568 allocate(save_array(2,my_size))
569
570 my_address = i_stok_old + address_cand_n(1)
571 sort_array(1,1:my_size) = intbuf_tab%cand_n(my_address+1:my_address+my_size)
572 my_address = i_stok_old + address_cand_e(1)
573 sort_array(2,1:my_size) = intbuf_tab%cand_e(my_address+1:my_address+my_size)
574 do i=1,my_size
575 my_index(i) = i
576 enddo
577 save_array(1:2,1:my_size) = sort_array(1:2,1:my_size)
578 mode = 0
579
580 call my_orders( mode,work,sort_array,my_index,my_size,2)
581 my_address = i_stok_old + address_cand_n(1)
582 do i=1,my_size
583 intbuf_tab%cand_n(my_address+i) = save_array(1,my_index(i))
584 enddo
585 my_address = i_stok_old + address_cand_e(1)
586 do i=1,my_size
587 intbuf_tab%cand_e(my_address+i) = save_array(2,my_index(i))
588 enddo
589 deallocate(my_index)
590 deallocate(sort_array)
591 deallocate(save_array)
592
593
594
595
596 DO i=1,nsn
597 IF(iix(i)/=0)THEN
598 voxel(iix(i),iiy(i),iiz(i))=0
599 ENDIF
600 ENDDO
601
602
603
604 deallocate( cand_n_size,cand_e_size )
605 deallocate( address_cand_n,address_cand_e )
606
607
608 RETURN
end diagonal values have been computed in the(sparse) matrix id.SOL
subroutine ymax(idn, fac, npc, pld, stiffmin, stiffmax, stiffini, stiffavg)
for(i8=*sizetab-1;i8 >=0;i8--)
void my_orders(int *mode, int *iwork, int *data, int *index, int *n, int *irecl)
integer, parameter nchartitle
integer, dimension(0:lrvoxel, 0:lrvoxel) crvoxel
subroutine i7cor3(x, irect, nsv, cand_e, cand_n, stf, stfn, gapv, igap, gap, gap_s, gap_m, istf, gapmin, gapmax, gap_s_l, gap_m_l, drad, ix1, ix2, ix3, ix4, nsvg, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, xi, yi, zi, stif, dgapload, last)
subroutine i7dst3(ix3, ix4, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, xi, yi, zi, x0, y0, z0, nx1, ny1, nz1, nx2, ny2, nz2, nx3, ny3, nz3, nx4, ny4, nz4, p1, p2, p3, p4, lb1, lb2, lb3, lb4, lc1, lc2, lc3, lc4, last)
subroutine i7pen3(marge, gapv, n1, n2, n3, pene, nx1, ny1, nz1, nx2, ny2, nz2, nx3, ny3, nz3, nx4, ny4, nz4, p1, p2, p3, p4, last)
subroutine tagnod(ix, nix, nix1, nix2, numel, iparte, tagbuf, npart)
subroutine upgrade_multimp(ni, multimp_parameter, intbuf_tab)