45
46
47
51
52
53
54#include "implicit_f.inc"
55
56
57
58#include "mvsiz_p.inc"
59
60
61
62#include "param_c.inc"
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 INTEGER ::
113 . NRTM,NRTS,MULTIMP,IADFIN,IGAP,
114 . NSN,NOINT,ITAB(*),NBX,NBY,NBZ,IAUTO,
115 . IRECTS(2,NRTS),IRECTM(2,NRTM),FLAGREMNODE
116 INTEGER ITASK
117 INTEGER, INTENT(INOUT) ::
118 . CAND_S(*),CAND_M(*),ADDCM(*),CHAINE(2,*),
119 . VOXEL(1:NBX+2,1:NBY+2,1:NBZ+2), I_MEM,II_STOK,
120 . KREMNODE(*),REMNODE(*)
122 . ,INTENT(IN) ::
123 . x(3,*),xyzm(6,*),
124 . gapmin, drad, marge, tzinf, dgapload,
125 . gap_s(*), gap_m(*), gap_s_l(*), gap_m_l(*)
126
127
128
129 INTEGER
130 . I,J,SS1,SS2,IBUG,
131 . N1,N2,MM1,MM2, iN1, iN2, iM1, iM2, K,L,
132 . PROV_S(2*MVSIZ),PROV_M(2*MVSIZ),
133 . IX1,IY1,IZ1,IX2,IY2,IZ2,
134 . IX,IY,IZ, FIRST_ADD,
135 . I_STOK, I_STOK_BAK, IEDG,
136 . PREV_ADD, CHAIN_ADD, CURRENT_ADD,
137 . NEDG, DEJA , MAX_ADD ,II_STOK0, M
139 . xx1, xx2,
140 . xmin, xmax,ymin,
ymax,zmin, zmax,
141 . yy1,yy2,zz1,zz2,
142 . aaa, dd,
143 . xminb,yminb,zminb,xmaxb,ymaxb,zmaxb
144 my_real,
dimension(:),
ALLOCATABLE :: xmax_edgs, xmin_edgs, ymax_edgs, ymin_edgs, zmax_edgs, zmin_edgs
145 my_real,
dimension(:),
ALLOCATABLE :: xmax_edgm, xmin_edgm, ymax_edgm, ymin_edgm, zmax_edgm, zmin_edgm
146 INTEGER, DIMENSION(:), ALLOCATABLE :: TAGREMLINE
147
148 ALLOCATE(xmax_edgs(nrts), xmin_edgs(nrts), ymax_edgs(nrts))
149 ALLOCATE(ymin_edgs(nrts), zmax_edgs(nrts), zmin_edgs(nrts))
150 ALLOCATE(xmax_edgm(nrtm), xmin_edgm(nrtm), ymax_edgm(nrtm))
151 ALLOCATE(ymin_edgm(nrtm), zmax_edgm(nrtm), zmin_edgm(nrtm))
152
153 IF(flagremnode==2)THEN
154 ALLOCATE(tagremline(nrts))
155 tagremline(1:nrts) = 0
156 ENDIF
157
158 aaa = zero
159
166
167
168
169
170 max_add =
max(1,4*(nrts))
174
175
176 IF(nrtm==0.OR.nrts==0)THEN
177
181 END IF
182
183
184
185
186 xmin = xyzm(1,1)
187 ymin = xyzm(2,1)
188 zmin = xyzm(3,1)
189 xmax = xyzm(4,1)
191 zmax = xyzm(6,1)
192
193 xminb = xmin
194 yminb = ymin
195 zminb = zmin
196 xmaxb = xmax
198 zmaxb = zmax
199
200
201
202
203
204
205
206
207
208
209 current_add=1
210
211 DO i = 1,nrts
212
213
214
215
216
217
218 n1=irects(1,i)
219 n2=irects(2,i)
220
221
222
223
224
225 xx1=x(1,n1)
226 xx2=x(1,n2)
227 xmax_edgs(i)=
max(xx1,xx2);
IF(xmax_edgs(i) < xmin) cycle
228 xmin_edgs(i)=
min(xx1,xx2);
IF(xmin_edgs(i) > xmax) cycle
229 yy1=x(2,n1)
230 yy2=x(2,n2)
231 ymax_edgs(i)=
max(yy1,yy2);
IF(ymax_edgs(i) < ymin) cycle
232 ymin_edgs(i)=
min(yy1,yy2);
IF(ymin_edgs(i) >
ymax) cycle
233 zz1=x(3,n1)
234 zz2=x(3,n2)
235 zmax_edgs(i)=
max(zz1,zz2);
IF(zmax_edgs(i) < zmin) cycle
236 zmin_edgs(i)=
min(zz1,zz2);
IF(zmin_edgs(i) > zmax) cycle
237
238
239
240
241
242 ix1=int(nbx*(xmin_edgs(i)-xminb)/(xmaxb-xminb))
243 iy1=int(nby*(ymin_edgs(i)-yminb)/(ymaxb-yminb))
244 iz1=int(nbz*(zmin_edgs(i)-zminb)/(zmaxb-zminb))
248
249 ix2=int(nbx*(xmax_edgs(i)-xminb)/(xmaxb-xminb))
250 iy2=int(nby*(ymax_edgs(i)-yminb)/(ymaxb-yminb))
251 iz2=int(nbz*(zmax_edgs(i)-zminb)/(zmaxb-zminb))
255
256
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 DO iz = iz1,iz2
294 DO iy = iy1,iy2
295 DO ix = ix1,ix2
296
297 first_add = voxel(ix,iy,iz)
298
299 IF(first_add == 0)THEN
300
301 voxel(ix,iy,iz) = current_add
305 ELSE
306
312 ENDIF
313
314 current_add = current_add+1
315
316 IF( current_add>=max_add)THEN
317
318
319 max_add = 2 * max_add
320
324 ENDIF
325
326 ENDDO
327 ENDDO
328 ENDDO
329
330 ENDDO
331
332
333
334
335
336
337
338
339 nedg = 0
340 i_stok = 0
341
342 DO iedg=1,nrtm
343
344 aaa = zero
345
346
347
348
349 n1 = irectm(1,iedg)
350 n2 = irectm(2,iedg)
351 mm1 = itab(n1)
352 mm2 = itab(n2)
353
354
355
356
357 xx1=x(1,n1)
358 xx2=x(1,n2)
359 xmax_edgm(iedg)=
max(xx1,xx2)+tzinf
360 xmin_edgm(iedg)=
min(xx1,xx2)-tzinf
361 yy1=x(2,n1)
362 yy2=x(2,n2)
363 ymax_edgm(iedg)=
max(yy1,yy2)+tzinf
364 ymin_edgm(iedg)=
min(yy1,yy2)-tzinf
365 zz1=x(3,n1)
366 zz2=x(3,n2)
367 zmax_edgm(iedg)=
max(zz1,zz2)+tzinf
368 zmin_edgm(iedg)=
min(zz1,zz2)-tzinf
369
370
371
372
373
374 ix1=int(nbx*(xmin_edgm(iedg)-aaa-xminb)/(xmaxb-xminb))
375 iy1=int(nby*(ymin_edgm(iedg)-aaa-yminb)/(ymaxb-yminb))
376 iz1=int(nbz*(zmin_edgm(iedg)-aaa-zminb)/(zmaxb-zminb))
380
381 ix2=int(nbx*(xmax_edgm(iedg)+aaa-xminb)/(xmaxb-xminb))
382 iy2=int(nby*(ymax_edgm(iedg)+aaa-yminb)/(ymaxb-yminb))
383 iz2=int(nbz*(zmax_edgm(iedg)+aaa-zminb)/(zmaxb-zminb))
387
388 deja = 0
389 i_stok_bak = i_stok
390
391
392 IF(flagremnode==2)THEN
393 k = kremnode(iedg)
394 l = kremnode(iedg+1)-1
395 DO m=k,l
396 tagremline(remnode(m)) = 1
397 ENDDO
398 ENDIF
399
400
401
402
403 DO iz = iz1,iz2
404 DO iy = iy1,iy2
405 DO ix = ix1,ix2
406
407 chain_add = voxel(ix,iy,iz)
408 DO WHILE(chain_add /= 0)
410
411
412 ss1=itab(irects(1,i))
413 ss2=itab(irects(2,i))
414
415 IF( (ss1==mm1).OR.(ss1==mm2).OR.
416 . (ss2==mm1).OR.(ss2==mm2) )THEN
418 cycle
419 END IF
420
421
422 IF(iauto==1 .AND. mm1<ss1 )THEN
424 cycle
425 END IF
426
427
428 IF(flagremnode==2)THEN
429 IF(tagremline(i)==1) THEN
431 cycle
432 ENDIF
433 ENDIF
434
435 i_stok = i_stok + 1
436 prov_s(i_stok) = i
437 prov_m(i_stok) = iedg
438
439
440 IF(deja==0) nedg = nedg + 1
441 deja=1
443
444 IF(i_stok>=nvsiz)THEN
446 1 nvsiz,irects,irectm,x ,ii_stok ,
447 2 cand_s,cand_m,nsn ,noint ,marge ,
448 3 i_mem ,prov_s,prov_m,multimp,addcm ,
449 4 chaine,iadfin,gapmin,drad ,igap ,
450 5 gap_s ,gap_m ,gap_s_l,gap_m_l,dgapload)
451
452 IF(i_mem==2) THEN
453 ii_stok=zero
454 GOTO 1000
456 i_stok = i_stok-nvsiz
457 DO j=1,i_stok
458 prov_s(j) = prov_s(j+nvsiz)
459 prov_m(j) = prov_m(j+nvsiz)
460 ENDDO
461 ENDIF
462
463
464 ENDDO
465 ENDDO
466 ENDDO
467 ENDDO
468
469
470 IF(flagremnode==2)THEN
471 k = kremnode(iedg)
472 l = kremnode(iedg+1)-1
473 DO m=k,l
474 tagremline(remnode(m)) = 0
475 ENDDO
476 ENDIF
477
478 ENDDO
479
480
481
482
483
484
486 1 i_stok,irects,irectm,x ,ii_stok,
487 2 cand_s,cand_m,nsn ,noint ,marge ,
488 3 i_mem ,prov_s,prov_m,multimp,addcm ,
489 4 chaine,iadfin,gapmin,drad ,igap ,
490 5 gap_s ,gap_m ,gap_s_l,gap_m_l,dgapload)
491
492
493
494
495
496
497
498 1000 CONTINUE
499
500
501
502
503
504
508 2 nbx, nby, nbz, voxel )
509
513 IF(flagremnode==2) DEALLOCATE(tagremline)
514
515 DEALLOCATE(xmax_edgs, xmin_edgs, ymax_edgs)
516 DEALLOCATE(ymin_edgs, zmax_edgs, zmin_edgs)
517 DEALLOCATE(xmax_edgm, xmin_edgm, ymax_edgm)
518 DEALLOCATE(ymin_edgm, zmax_edgm, zmin_edgm)
519
520 RETURN
if(complex_arithmetic) id
subroutine i11resetvoxel1(tmin, tmax, nbx, nby, nbz, voxel)
subroutine ymax(idn, fac, npc, pld, stiffmin, stiffmax, stiffini, stiffavg)
integer function, dimension(:), pointer ireallocate(ptr, new_size)
integer, dimension(:), pointer lchain_elem
integer, dimension(:), pointer lchain_last
integer, dimension(:), pointer lchain_next
subroutine i11sto_vox1(j_stok, irects, irectm, x, ii_stok, cand_n, cand_e, nsn, noint, marge, i_mem, prov_n, prov_e, multimp, addcm, chaine, iadfin, gapmin, drad, igap, gap_s, gap_m, gap_s_l, gap_m_l, dgapload)