40
41
42
43#include "implicit_f.inc"
44
45
46
47 INTEGER NN, NHOL, IPOLY(6+2*NN+1+NHOL,*), IPOLY_OLD(*),
48 . INEDGE(6,*), NBNEDGE, INS(2,*), INZ, IZ(3,*), NNS3,
49 . NPOLY, NS, IELNOD(2*NN,*), IELNOD_OLD(*)
51 . rpoly(4+3*2*nn+3*nhol,*), rpoly_old(*), rnedge(6,*), nx,
52 . ny, nz, x0, y0, z0, rns(4,*)
53
54
55
56 INTEGER NN1, I, IADHOL(NHOL+1), NSEG, ITAG(2*NN+1), II,
57 . TSEG(3,2*NN), NHOL_OLD, ICUT, J, JJ, NSEG_INI, NSEG1,
58 . REDIR(NN), N1, N2, ISTOP, NNP, ICLOSE, ISEG,
59 . POLY(2*NN,NN), IN, IN1, IN2, LENPOLY(NN), K, KK,
60 . THOL(NHOL), JN, NHOLP, IADHOLP(NHOL), KN,
61 . ITAG2(2*NN), ITAGN(NN), ITAGS(NN), ISEG_OLD
63 . x1, y1, z1, x2, y2, z2, zl1, zl2,
alpha, npx, npy, npz,
64 . xp0, yp0, zp0, vx, vy, vz, xx, yy, zz, xl(nn), xlmin,
65 . xlc, xx1, yy1, zz1, xx2, yy2, zz2, vx1, vy1,
66 . vz1, vx2, vy2, vz2, nr1, nr2, ss, vvx, vvy, vvz,
67 . ss1, zl, xns(3,nn), tole, ll, zlm
68
69 tole=epsilon(zero)*0.5
70
71
72 IF (nhol==0) THEN
73 nn1=nn
74 ELSE
75 nn1=ipoly_old(6+nn+1+1)-1
76 DO i=1,nhol
77 iadhol(i)=ipoly_old(6+nn+1+i)
78 ENDDO
79 iadhol(nhol+1)=nn+1
80 ENDIF
81 ns=0
82 nseg=0
83
84 DO i=1,2*nn
85 itag(i)=0
86 itag2(i)=0
87 ENDDO
88 itag(2*nn+1)=0
89 DO i=1,nn
90 itagn(i)=0
91 itags(i)=0
92 ENDDO
93
94 DO i=1,nn1
95 ii=i+1
96 IF (i==nn1) ii=1
97 x1=rpoly_old(4+3*(i-1)+1)
98 y1=rpoly_old(4+3*(i-1)+2)
99 z1=rpoly_old(4+3*(i-1)+3)
100 x2=rpoly_old(4+3*(ii-1)+1)
101 y2=rpoly_old(4+3*(ii-1)+2)
102 z2=rpoly_old(4+3*(ii-1)+3)
103 zl1=(x1-x0)*nx+(y1-y0)*ny+(z1-z0)*nz
104 zl2=(x2-x0)*nx+(y2-y0)*ny+(z2-z0)*nz
105 IF (zl1-zl2/=zero) THEN
108 ns=ns+1
109 ins(1,ns)=ipoly_old(6+i)
110 ins(2,ns)=ipoly_old(6+ii)
115 nseg=nseg+1
116 tseg(1,nseg)=i
117 tseg(2,nseg)=-ns
118 tseg(3,nseg)=nseg+1
119 nseg=nseg+1
120 tseg(1,nseg)=-ns
121 tseg(2,nseg)=ii
122 tseg(3,nseg)=nseg+1
123 ELSEIF (
alpha==zero)
THEN
124 IF (itagn(ii)==0) THEN
125 ns=ns+1
126 itagn(ii)=ns
127 ins(1,ns)=ipoly_old(6+i)
128 ins(2,ns)=ipoly_old(6+ii)
129 rns(1,ns)=zero
130 xns(1,ns)=x2
131 xns(2,ns)=y2
132 xns(3,ns)=z2
133 itags(ns)=1
134 nseg=nseg+1
135 tseg(1,nseg)=i
136 tseg(2,nseg)=-ns
137 tseg(3,nseg)=nseg+1
138 nseg=nseg+1
139 tseg(1,nseg)=-ns
140 tseg(2,nseg)=ii
141 tseg(3,nseg)=nseg+1
142 ELSE
143 nseg=nseg+1
144 tseg(1,nseg)=i
145 tseg(2,nseg)=ii
146 tseg(3,nseg)=nseg+1
147 ENDIF
148 ELSEIF (
alpha==one)
THEN
149 IF (itagn(i)==0) THEN
150 ns=ns+1
151 itagn(i)=ns
152 ins(1,ns)=ipoly_old(6+i)
153 ins(2,ns)=ipoly_old(6+ii)
154 rns(1,ns)=one
155 xns(1,ns)=x1
156 xns(2,ns)=y1
157 xns(3,ns)=z1
158 itags(ns)=1
159 nseg=nseg+1
160 tseg(1,nseg)=i
161 tseg(2,nseg)=-ns
162 tseg(3,nseg)=nseg+1
163 nseg=nseg+1
164 tseg(1,nseg)=-ns
165 tseg(2,nseg)=ii
166 tseg(3,nseg)=nseg+1
167 ELSE
168 nseg=nseg+1
169 tseg(1,nseg)=i
170 tseg(2,nseg)=ii
171 tseg(3,nseg)=nseg+1
172 ENDIF
173 ELSE
174 nseg=nseg+1
175 tseg(1,nseg)=i
176 tseg(2,nseg)=ii
177 tseg(3,nseg)=nseg+1
178 ENDIF
179 ELSE
180 nseg=nseg+1
181 tseg(1,nseg)=i
182 tseg(2,nseg)=ii
183 tseg(3,nseg)=nseg+1
184 ENDIF
185 ENDDO
186 tseg(3,nseg)=1
187
188 nhol_old=nhol
189 nhol=0
190 DO i=1,nhol_old
191 ii = i+1
192 icut=0
193 nseg_ini=nseg
194 DO j=iadhol(i),iadhol(i+1)-1
195 jj=j+1
196 IF (j==(iadhol(i+1)-1)) jj=iadhol(i)
197 x1=rpoly_old(4+3*(j-1)+1)
198 y1=rpoly_old(4+3*(j-1)+2)
199 z1=rpoly_old(4+3*(j-1)+3)
200 x2=rpoly_old(4+3*(jj-1)+1)
201 y2=rpoly_old(4+3*(jj-1)+2)
202 z2=rpoly_old(4+3*(jj-1)+3)
203 zl1=(x1-x0)*nx+(y1-y0)*ny+(z1-z0)*nz
204 zl2=(x2-x0)*nx+(y2-y0)*ny+(z2-z0)*nz
205 IF (zl1-zl2/=zero) THEN
208 icut=1
209 ns=ns+1
210 ins(1,ns)=ipoly_old(6+j)
211 ins(2,ns)=ipoly_old(6+jj)
216 nseg=nseg+1
217 tseg(1,nseg)=j
218 tseg(2,nseg)=-ns
219 tseg(3,nseg)=nseg+1
220 nseg=nseg+1
221 tseg(1,nseg)=-ns
222 tseg(2,nseg)=jj
223 tseg(3,nseg)=nseg+1
224 ELSEIF (
alpha==zero)
THEN
225 icut=1
226 IF (itagn(ii)==0) THEN
227 ns=ns+1
228 itagn(jj)=ns
229 ins(1,ns)=ipoly_old(6+j)
230 ins(2,ns)=ipoly_old(6+jj)
231 rns(1,ns)=zero
232 xns(1,ns)=x2
233 xns(2,ns)=y2
234 xns(3,ns)=z2
235 itags(ns)=1
236 nseg=nseg+1
237 tseg(1,nseg)=j
238 tseg(2,nseg)=-ns
239 tseg(3,nseg)=nseg+1
240 nseg=nseg+1
241 tseg(1,nseg)=-ns
242 tseg(2,nseg)=jj
243 tseg(3,nseg)=nseg+1
244 ELSE
245 nseg=nseg+1
246 tseg(1,nseg)=j
247 tseg(2,nseg)=jj
248 tseg(3,nseg)=nseg+1
249 ENDIF
250 ELSEIF (
alpha==one)
THEN
251 icut=1
252 IF (itagn(i)==0) THEN
253 ns=ns+1
254 itagn(j)=ns
255 ins(1,ns)=ipoly_old(6+j)
256 ins(2,ns)=ipoly_old(6+jj)
257 rns(1,ns)=one
258 xns(1,ns)=x1
259 xns(2,ns)=y1
260 xns(3,ns)=z1
261 itags(ns)=1
262 nseg=nseg+1
263 tseg(1,nseg)=j
264 tseg(2,nseg)=-ns
265 tseg(3,nseg)=nseg+1
266 nseg=nseg+1
267 tseg(1,nseg)=-ns
268 tseg(2,nseg)=jj
269 tseg(3,nseg)=nseg+1
270 ELSE
271 nseg=nseg+1
272 tseg(1,nseg)=j
273 tseg(2,nseg)=jj
274 tseg(3,nseg)=nseg+1
275 ENDIF
276 ELSE
277 nseg=nseg+1
278 tseg(1,nseg)=j
279 tseg(2,nseg)=jj
280 tseg(3,nseg)=nseg+1
281 ENDIF
282 ELSE
283 nseg=nseg+1
284 tseg(1,nseg)=j
285 tseg(2,nseg)=jj
286 tseg(3,nseg)=nseg+1
287 ENDIF
288 ENDDO
289 tseg(3,nseg)=nseg_ini+1
290
291 IF (icut==0) THEN
292
293 nhol=nhol+1
294 DO j=nseg_ini+1,nseg
295 itag(j)=-nhol
296 ENDDO
297 ENDIF
298 ENDDO
299
300 IF(ns <= 1) THEN
301 npoly=0
302 RETURN
303 ENDIF
304 nseg1=nseg
305
306 npx=rpoly_old(2)
307 npy=rpoly_old(3)
308 npz=rpoly_old(4)
309 xp0=rpoly_old(5)
310 yp0=rpoly_old(6)
311 zp0=rpoly_old(7)
312 vx=ny*npz-nz*npy
313 vy=nz*npx-nx*npz
314 vz=nx*npy-ny*npx
315 DO i=1,ns
316 xx=xns(1,i)
317 yy=xns(2,i)
318 zz=xns(3,i)
319 xl(i)=(xx-xp0)*vx+(yy-yp0)*vy+(zz-zp0)*vz
320 ENDDO
321 xlmin=ep30
322 DO i=1,ns
323 redir(i)=i
324 ENDDO
325 DO i=1,ns
326 xlmin=xl(redir(i))
327 DO j=i+1,ns
328 xlc=xl(redir(j))
329 IF (xlc<xlmin) THEN
330 jj=redir(j)
331 redir(j)=redir(i)
332 redir(i)=jj
333 xlmin=xlc
334 ENDIF
335 ENDDO
336 ENDDO
337 ii=0
338 DO i=1,ns
339 ii=ii+1
340 IF (ii==1) THEN
341 n1=redir(i)
342 n2=redir(i+1)
343 x1=xns(1,n1)
344 y1=xns(2,n1)
345 z1=xns(3,n1)
346 x2=xns(1,n2)
347 y2=xns(2,n2)
348 z2=xns(3,n2)
349 ll=sqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)
350
351 IF (ll>tole) THEN
352 nseg=nseg+1
353 tseg(1,nseg)=-n1
354 tseg(2,nseg)=-n2
355 tseg(3,nseg)=-1
356 ENDIF
357 ELSE
358 ii=0
359 ENDIF
360 ENDDO
361
362 istop=0
363 npoly=1
364 nnp=0
365 DO WHILE (istop==0)
366 i=1
367 DO WHILE (itag(i)/=0.AND.i<=nseg)
368 i=i+1
369 ENDDO
370 IF (i==nseg+1) THEN
371 istop=1
372 cycle
373 ENDIF
374
375 iclose=0
376 iseg=i
377 DO WHILE (iclose==0)
378 nnp=nnp+1
379 poly(nnp,npoly)=tseg(1,iseg)
380 in=tseg(2,iseg)
381 IF (tseg(3,iseg)>0) THEN
382 itag(iseg)=1
383 IF (in<0) THEN
384 iseg_old=iseg
385 iseg=0
386 DO j=nseg1+1,nseg
387 IF (itag(j)/=0) cycle
388 in1=tseg(1,j)
389 in2=tseg(2,j)
390 IF (in1==in) THEN
391 iseg=j
392 ELSEIF (in2==in) THEN
393 iseg=j
394 tseg(1,j)=in2
395 tseg(2,j)=in1
396 ENDIF
397 ENDDO
398 IF (iseg==0) iseg=tseg(3,iseg_old)
399 ELSE
400 iseg=tseg(3,iseg)
401 ENDIF
402 ELSE
403 IF (itag2(iseg)==1) itag(iseg)=1
404 itag2(iseg)=1
405 iseg=0
406 DO j=1,nseg1
407 in1=tseg(1,j)
408 IF (in1==in) iseg=j
409 ENDDO
410 ENDIF
411
412 IF (itag(iseg)/=0) THEN
413 iclose=1
414 lenpoly(npoly)=nnp
415 npoly=npoly+1
416 nnp=0
417 ENDIF
418 ENDDO
419 ENDDO
420 npoly=npoly-1
421
422 DO i=1,nhol
423 j=1
424 DO WHILE (itag(j)/=-i)
425 j=j+1
426 ENDDO
427 n1=tseg(1,j)
428 xx=rpoly_old(4+3*(n1-1)+1)
429 yy=rpoly_old(4+3*(n1-1)+2)
430 zz=rpoly_old(4+3*(n1-1)+3)
431 DO j=1,npoly
433 DO k=1,lenpoly(j)
434 kk=k+1
435 IF (k==lenpoly(j)) kk=1
436 n1=poly(k,i)
437 n2=poly(kk,i)
438 IF (n1>0) THEN
439 xx1=rpoly_old(4+3*(n1-1)+1)
440 yy1=rpoly_old(4+3*(n1-1)+2)
441 zz1=rpoly_old(4+3*(n1-1)+3)
442 ELSE
443 xx1=xns(1,-n1)
444 yy1=xns(2,-n1)
445 zz1=xns(3,-n1)
446 ENDIF
447 IF (n2>0) THEN
448 xx2=rpoly_old(4+3*(n2-1)+1)
449 yy2=rpoly_old(4+3*(n2-1)+2)
450 zz2=rpoly_old(4+3*(n2-1)+3)
451 ELSE
452 xx2=xns(1,-n2)
453 yy2=xns(2,-n2)
454 zz2=xns(3,-n2)
455 ENDIF
456 vx1=xx1-xx
457 vy1=yy1-yy
458 vz1=zz1-zz
459 vx2=xx2-xx
460 vy2=yy2-yy
461 vz2=zz2-zz
462 nr1=sqrt(vx1**2+vy1**2+vz1**2)
463 nr2=sqrt(vx2**2+vy2**2+vz2**2)
464 vx1=vx1/nr1
465 vy1=vy1/nr1
466 vz1=vz1/nr1
467 vx2=vx2/nr2
468 vy2=vy2/nr2
469 vz2=vz2/nr2
470 ss=vx1*vx2+vy1*vy2+vz1*vz2
471 vvx=vy1*vz2-vz1*vy2
472 vvy=vz1*vx2-vx1*vz2
473 vvz=vx1*vy2-vy1*vx2
474 ss1=npx*vvx+npy*vvy+npz*vvz
475 IF (ss1>=zero) THEN
477 ELSE
479 ENDIF
480 ENDDO
481
482 IF (abs(
alpha)>=one) thol(i)=j
483 ENDDO
484 ENDDO
485
486 DO i=1,ns
487 rns(2,i)=xns(1,i)
488 rns(3,i)=xns(2,i)
489 rns(4,i)=xns(3,i)
490 ENDDO
491
492 DO i=1,npoly
493 ipoly(1,i)=ipoly_old(1)
494 ipoly(3,i)=ipoly_old(3)
495 ipoly(4,i)=ipoly_old(4)
496 ipoly(5,i)=ipoly_old(5)
497 ipoly(6,i)=ipoly_old(6)
498 rpoly(1,i)=zero
499 rpoly(2,i)=npx
500 rpoly(3,i)=npy
501 rpoly(4,i)=npz
502 nnp=0
503 DO j=1,lenpoly(i)
504 jn=poly(j,i)
505 IF (jn>0) THEN
506 nnp=nnp+1
507 ipoly(6+nnp,i)=ipoly_old(6+jn)
508 rpoly(4+3*(nnp-1)+1,i)=rpoly_old(4+3*(jn-1)+1)
509 rpoly(4+3*(nnp-1)+2,i)=rpoly_old(4+3*(jn-1)+2)
510 rpoly(4+3*(nnp-1)+3,i)=rpoly_old(4+3*(jn-1)+3)
511 ielnod(nnp,i)=ielnod_old(jn)
512 ELSE
513 IF (itags(-jn)==0) THEN
514 nnp=nnp+1
515 ipoly(6+nnp,i)=-nns3+jn
516 rpoly(4+3*(nnp-1)+1,i)=xns(1,-jn)
517 rpoly(4+3*(nnp-1)+2,i)=xns(2,-jn)
518 rpoly(4+3*(nnp-1)+3,i)=xns(3,-jn)
519 ielnod(nnp,i)=-1
520 ELSE
521 jj=poly(j+1,i)
522 IF (j==lenpoly(i)) jj=poly(1,i)
523 x1=xns(1,-jn)
524 y1=xns(2,-jn)
525 z1=xns(3,-jn)
526 IF (jj>0) THEN
527 x2=rpoly_old(4+3*(jj-1)+1)
528 y2=rpoly_old(4+3*(jj-1)+2)
529 z2=rpoly_old(4+3*(jj-1)+3)
530 ELSE
531 x2=xns(1,-jj)
532 y2=xns(2,-jj)
533 z2=xns(3,-jj)
534 ENDIF
535 ll=sqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)
536 IF (ll>tole) THEN
537 nnp=nnp+1
538 ipoly(6+nnp,i)=-nns3+jn
539 rpoly(4+3*(nnp-1)+1,i)=xns(1,-jn)
540 rpoly(4+3*(nnp-1)+2,i)=xns(2,-jn)
541 rpoly(4+3*(nnp-1)+3,i)=xns(3,-jn)
542 ielnod(nnp,i)=-1
543 ENDIF
544 ENDIF
545 ENDIF
546 ENDDO
547
548 nholp=0
549
550 DO j=1,nhol
551 IF (thol(j)==i) THEN
552 nholp=nholp+1
553 iadholp(nholp)=nnp+1
554 DO k=1,nseg
555 IF (tseg(3,k)==-j) THEN
556 nnp=nnp+1
557 kn=tseg(1,k)
558 ipoly(6+nnp,i)=kn
559 rpoly(4+3*(nnp-1)+1,i)=rpoly_old(4+3*(kn-1)+1)
560 rpoly(4+3*(nnp-1)+2,i)=rpoly_old(4+3*(kn-1)+2)
561 rpoly(4+3*(nnp-1)+3,i)=rpoly_old(4+3*(kn-1)+3)
562 ENDIF
563 ENDDO
564 ENDIF
565 ENDDO
566 ipoly(2,i)=nnp
567 ipoly(6+nnp+1,i)=nholp
568 DO j=1,nholp
569 ipoly(6+nnp+1+j,i)=iadholp(j)
570 ENDDO
571 ENDDO
572
573 DO i=nseg1+1,nseg
574 nbnedge=nbnedge+1
575 n1=-tseg(1,i)
576 n2=-tseg(2,i)
577 inedge(1,nbnedge)=ipoly_old(1)
578 inedge(2,nbnedge)=nns3+n1
579 inedge(3,nbnedge)=nns3+n2
580 inedge(4,nbnedge)=ipoly_old(3)
581 inedge(5,nbnedge)=ipoly_old(4)
582 inedge(6,nbnedge)=inz
583
584 xx1=xns(1,n1)
585 yy1=xns(2,n1)
586 zz1=xns(3,n1)
587
588 xx2=xns(1,n2)
589 yy2=xns(2,n2)
590 zz2=xns(3,n2)
591
592 rnedge(1,nbnedge)=xx1
593 rnedge(2,nbnedge)=yy1
594 rnedge(3,nbnedge)=zz1
595 rnedge(4,nbnedge)=xx2
596 rnedge(5,nbnedge)=yy2
597 rnedge(6,nbnedge)=zz2
598 ENDDO
599
600 DO i=1,npoly
601 zl=zero
602 zlm=zero
603 DO j=1,ipoly(2,i)
604 xx=rpoly(4+3*(j-1)+1,i)
605 yy=rpoly(4+3*(j-1)+2,i)
606 zz=rpoly(4+3*(j-1)+3,i)
607 zl=(xx-x0)*nx+(yy-y0)*ny+(zz-z0)*nz
608 IF (abs(zl)>abs(zlm)) zlm=zl
609 ENDDO
610 iz(1,i)=1
611 IF (zlm>zero) THEN
612 iz(2,i)=inz+1
613 ELSEIF (zlm<zero) THEN
614 iz(2,i)=inz
615 ENDIF
616 ENDDO
617
618 RETURN