31 1 JLT ,CAND_N ,CAND_E ,
35 6 IX2 ,IX3 ,IX4 ,NSVG ,STIF ,
36 7 INACTI ,MSEGLO ,GAPS ,GAPM ,GAPMXL ,
37 8 IRECT ,IRTLM ,TIME_S ,GAP_NM ,ITAB ,
38 9 ICONT_I,NNX ,NNY ,NNZ ,
39 A FAR ,PENT ,DIST ,LB ,LC ,
40 B LBP ,LCP ,SUBTRIA ,MVOISN,IBOUND,
41 C VTX_BISECTOR ,DRAD ,DGAPLOAD)
49#include "implicit_f.inc"
58 INTEGER JLT, NIN, NSN, INACTI,
60 . CAND_E(*),NSVG(MVSIZ)
61 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
62 . MSEGLO(*),IRTLM(4,NSN), SUBTRIA(MVSIZ),
63 . MVOISN(MVSIZ,4), IBOUND(4,*)
65 . TIME_S(2,NSN), GAPS(*), GAPM(*), GAPMXL(MVSIZ)
66 my_real ,
INTENT(IN) :: DGAPLOAD ,DRAD
68 . XX(MVSIZ,5),YY(MVSIZ,5),ZZ(MVSIZ,5),
69 . xi(mvsiz), yi(mvsiz), zi(mvsiz), stif(mvsiz),
70 . nnx(mvsiz,5), nny(mvsiz,5), nnz(mvsiz,5),
71 . pent(mvsiz,4), dist(mvsiz), lb(mvsiz,4), lc(mvsiz,4),
72 . lbp(mvsiz,4), lcp(mvsiz,4), gap_nm(4,mvsiz), dd(mvsiz,4)
73 INTEGER IRECT(4,*),ITAB(
74REAL*4 VTX_BISECTOR(3,2,*)
78 INTEGER I, J, K, JJ, IT, ITOLD
79 INTEGER I4N(MVSIZ), I1, I2,
80 . ITRIA(2,4), INEIGH(3,4), MARQUE(4,MVSIZ),
81 . SUBTRIA_N(MVSIZ), IB1, IB2, IB3, IBX, IX, IY, IZ
83 . XI1(MVSIZ), XI0V(MVSIZ),
84 . yi1(mvsiz), yi0v(mvsiz),
85 . zi1(mvsiz), zi0v(mvsiz),
90 . gap_mm(mvsiz), gap, gap2,
93 . lap, hla, hlb(mvsiz,4), hlc(mvsiz,4),
100 . x0n(mvsiz,4), y0n(mvsiz,4), z0n(mvsiz,4),
101 . xn(mvsiz,4),yn(mvsiz,4),zn(mvsiz,4),
103 . x12(mvsiz), y12(mvsiz), z12(mvsiz),
104 . px, py, pz, pp, p1, p2, xh, yh, zh, d1, d2, d3, vx, vy, vz,
106 . lx, lax, lbx, lcx, ld(mvsiz), dmin
108 DATA itria/1,2,2,3,3,4,4,1/,
109 . ineigh/4,3,2,1,4,3,2,1,4,3,2,1/
114 marque(1:4,1:jlt) = 0
115 pent(1:jlt,1:4) = ep20
122 IF(stif(i) <= zero)cycle
124 x0n(i,1) = xx(i,1) - xx(i,5)
125 y0n(i,1) = yy(i,1) - yy(i,5)
126 z0n(i,1) = zz(i,1) - zz(i,5)
128 x0n(i,2) = xx(i,2) - xx(i,5)
129 y0n(i,2) = yy(i,2) - yy(i,5)
130 z0n(i,2) = zz(i,2) - zz(i,5)
132 x0n(i,3) = xx(i,3) - xx(i,5)
133 y0n(i,3) = yy(i,3) - yy(i,5)
134 z0n(i,3) = zz(i,3) - zz(i,5)
136 x0n(i,4) = xx(i,4) - xx(i,5)
137 y0n(i,4) = yy(i,4) - yy(i,5)
138 z0n(i,4) = zz(i,4) - zz(i,5)
140 IF(ix3(i)/=ix4(i))
THEN
141 gap_mm(i)=fourth*(gap_nm(1,i)+gap_nm(2,i)+gap_nm(3,i)+gap_nm(4,i))
143 gap_mm(i)=gap_nm(3,i)
150#include "vectorize.inc"
153 IF(stif(i) <= zero)cycle
155 xi0v(i) = xx(i,5) - xi(i)
156 yi0v(i) = yy(i,5) - yi(i)
157 zi0v(i) = zz(i,5) - zi(i)
159 xij(i,1) = xx(i,1) - xi(i)
160 yij(i,1) = yy(i,1) - yi(i)
161 zij(i,1) = zz(i,1) - zi(i)
163 xij(i,2) = xx(i,2) - xi(i)
164 yij(i,2) = yy(i,2) - yi(i)
165 zij(i,2) = zz(i,2) - zi(i)
167 xij(i,3) = xx(i,3) - xi(i)
168 yij(i,3) = yy(i,3) - yi(i)
169 zij(i,3) = zz(i,3) - zi(i)
171 xij(i,4) = xx(i,4) - xi(i)
172 yij(i,4) = yy(i,4) - yi(i)
173 zij(i,4) = zz(i,4) - zi(i)
175 sx1 = yi0v(i)*zij(i,1) - zi0v(i)*yij(i,1)
176 sy1 = zi0v(i)*xij(i,1) - xi0v(i)*zij(i,1)
177 sz1 = xi0v(i)*yij(i,1) - yi0v(i)*xij(i,1)
179 sx2 = yi0v(i)*zij(i,2) - zi0v(i)*yij(i,2)
180 sy2 = zi0v(i)*xij(i,2) - xi0v(i)*zij(i,2)
181 sz2 = xi0v(i)*yij(i,2) - yi0v(i)*xij(i,2)
183 xn(i,1) = y0n(i,1)*z0n(i,2) - z0n(i,1)*y0n(i,2)
184 yn(i,1) = z0n(i,1)*x0n(i,2) - x0n(i,1)*z0n(i,2)
185 zn(i,1) = x0n(i,1)*y0n(i,2) - y0n(i,1)*x0n(i,2)
187 .
max(em30,sqrt(xn(i,1)*xn(i,1)+ yn(i,1)*yn(i,1)+ zn(i,1)*zn(i,1)))
191 lb(i,1) = -(xn(i,1)*sx2 + yn(i,1)*sy2 + zn(i,1)*sz2)*nn
192 lc(i,1) = (xn(i,1)*sx1 + yn(i,1)*sy1 + zn(i,1)*sz1)*nn
194 sx3 = yi0v(i)*zij(i,3) - zi0v(i)*yij(i,3)
195 sy3 = zi0v(i)*xij(i,3) - xi0v(i)*zij(i,3)
196 sz3 = xi0v(i)*yij(i,3) - yi0v(i)*xij(i,3)
198 xn(i,2) = y0n(i,2)*z0n(i,3) - z0n(i,2)*y0n(i,3)
199 yn(i,2) = z0n(i,2)*x0n(i,3) - x0n(i,2)*z0n(i,3)
200 zn(i,2) = x0n(i,2)*y0n(i,3) - y0n(i,2)*x0n(i,3)
202 .
max(em30,sqrt(xn(i,2)*xn(i,2)+ yn(i,2)*yn(i
206 lb(i,2) = -(xn(i,2)*sx3 + yn(i,2)*sy3 + zn(i,2)*sz3)*nn
207 lc(i,2) = (xn(i,2)*sx2 + yn(i,2)*sy2 + zn(i,2)*sz2)*nn
209 sx4 = yi0v(i)*zij(i,4) - zi0v(i)*yij(i,4)
210 sy4 = zi0v(i)*xij(i,4) - xi0v(i
211 sz4 = xi0v(i)*yij(i,4) - yi0v(i)*xij(i,4)
213 xn(i,3) = y0n(i,3)*z0n(i,4) - z0n(i,3)*y0n(i,4)
214 yn(i,3) = z0n(i,3)*x0n(i,4) - x0n(i,3)*z0n(i,4)
215 zn(i,3) = x0n(i,3)*y0n(i,4) - y0n(i,3)*x0n(i,4)
217 .
max(em30,sqrt(xn(i,3)*xn(i,3)+ yn(i,3)*yn(i,3)+ zn(i,3)*zn(i,3)))
221 lb(i,3) = -(xn(i,3)*sx4 + yn(i,3)*sy4 + zn(i,3)*sz4)*nn
222 lc(i,3) = (xn(i,3)*sx3 + yn(i,3)*sy3 + zn(i,3)*sz3)*nn
224 xn(i,4) = y0n(i,4)*z0n(i,1) - z0n(i,4)*y0n(i,1)
225 yn(i,4) = z0n(i,4)*x0n(i,1) - x0n(i,4)*z0n(i,1)
226 zn(i,4) = x0n(i,4)*y0n(i,1) - y0n(i,4)*x0n(i,1)
228 .
max(em30,sqrt(xn(i,4)*xn(i,4)+ yn(i,4)*yn(i,4)+ zn(i,4)*zn
232 lb(i,4) = -(xn(i,4)*sx1 + yn(i,4)*sy1 + zn(i,4)*sz1)*nn
233 lc(i,4) = (xn(i,4)*sx4 + yn(i,4)*sy4 + zn(i,4)*sz4)*nn
237#include "vectorize.inc"
240 IF(stif(i) <= zero)cycle
242 aaa = one/
max(em30,x0n(i,1)*x0n(i,1)+y0n(i,1)*y0n(i,1)+z0n(i,1)*z0n(i,1))
243 hlc(i,1) = lc(i,1)*abs(lc(i,1))*aaa
244 hlb(i,4) = lb(i,4)*abs(lb(i,4))*aaa
245 al(i,1) = -(xi0v(i)*x0n(i,1)+yi0v(i)*y0n(i,1)+zi0v(i)*z0n(i,1))*aaa
246 al(i,1) =
max(zero,
min(one,al(i,1)))
247 aaa = one/
max(em30,x0n(i,2)*x0n(i,2)+y0n(i,2)*y0n(i,2)+z0n(i,2)*z0n(i,2))
248 hlc(i,2) = lc(i,2)*abs(lc(i,2))*aaa
249 hlb(i,1) = lb(i,1)*abs(lb(i,1))*aaa
250 al(i,2) = -(xi0v(i)*x0n(i,2)+yi0v(i)*y0n(i,2)+zi0v(i)*z0n(i,2))*aaa
251 al(i,2) =
max(zero,
min(one,al(i,2)))
252 aaa = one/
max(em30,x0n(i,3)*x0n(i,3)+y0n(i,3)*y0n(i,3)+z0n(i,3)*z0n(i,3))
253 hlc(i,3) = lc(i,3)*abs(lc(i,3))*aaa
254 hlb(i,2) = lb(i,2)*abs(lb(i,2))*aaa
255 al(i,3) = -(xi0v(i)*x0n(i,3)+yi0v(i)*y0n(i,3)+zi0v(i)*z0n(i,3))*aaa
256 al(i,3) =
max(zero,
min(one,al(i,3)))
257 aaa = one/
max(em30,x0n(i,4)*x0n(i,4)+y0n(i,4)*y0n(i,4)+z0n(i,4)*z0n(i,4))
258 hlc(i,4) = lc(i,4)*abs(lc(i,4))*aaa
259 hlb(i,3) = lb(i,3)*abs(lb(i,3))*aaa
260 al(i,4) = -(xi0v(i)*x0n(i,4)+yi0v(i)*y0n(i,4)+zi0v(i)*z0n(i,4))*aaa
261 al(i,4) =
max(zero,
min(one,al(i,4)))
268#include "vectorize.inc"
271 IF(stif(i) <= zero)cycle
273 x12(i) = xx(i,i2) - xx(i,i1)
274 y12(i) = yy(i,i2) - yy(i,i1)
275 z12(i) = zz(i,i2) - zz(i,i1)
279 lap = one-lbp(i,it)-lcp(i,it)
281 aaa = one /
max(em20,x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i))
282 hla= lap*abs(lap) * aaa
284 + hla<=hlb(i,it).AND.hla<=hlc(i,it))
THEN
285 lbp(i,it) = (xij(i,i2)*x12(i)+yij(i,i2)*y12(i)+zij(i,i2)*z12(i)) * aaa
286 lbp(i,it) =
max(zero,
min(one,lbp(i,it)))
287 lcp(i,it) = one - lbp(i,it)
288 ELSEIF(lbp(i,it)<zero.AND.
289 + hlb(i,it)<=hlc(i,it).AND.hlb(i,it)<=hla)
THEN
292 ELSEIF(lcp(i,it)<zero.AND.
293 + hlc(i,it)<=hla.AND.hlc(i,it)<=hlb(i,it))
THEN
298 lap = one-lbp(i,it)-lcp(i,it)
300 yp =lap*yy(i,5)+lbp(i,it)*yy(i,i1) + lcp(i,it)*yy(i,i2)
301 zp =lap*zz(i,5)+lbp(i,it)*zz(i,i1) + lcp(i,it)*zz(i,i2)
305 dd(i,it)=dx*dx+dy*dy+dz*dz
311#include "vectorize.inc"
314 IF (stif(i) <= zero)cycle
316 lap = one-lbp(i,it)-lcp(i,it)
318 gap =
min(
max(drad,gaps(i)+lap*gap_mm(i)+lbp(i,it)*gap_nm(i1,i)+lcp(i,it)*gap_nm(i2,i)+dgapload),
319 .
max(drad,gapmxl(i)+dgapload))
322 bb(i,it) =((xx(i,5)-xi(i))*xn(i,it)+(yy(i,5)-yi(i))*yn
324 IF(bb(i,it) > zero)
THEN
327 pent(i,it)=
max(zero,gap+bb(i,it))
329 pent(i,it)=
max(zero,gap-sqrt(dd(i,it)))
334#include "vectorize.inc"
337 IF (stif(i) <= zero)cycle
340 IF(lb(i,it) < -em03 .OR. lc(i,it) < -em03 .OR.
341 . lb(i,it)+lc(i,it) > one+em03) far(i,it)=1
343 xh=xi(i)+bb(i,it)*xn(i,it)
344 yh=yi(i)+bb(i,it)*yn(i,it)
345 zh=zi(i)+bb(i,it)*zn(i,it)
347 IF(ix3(i) /= ix4(i))
THEN
351 IF(mvoisn(i,it)==0)
THEN
353 IF( (xh-xx(i,i1))* nnx(i,it)+
354 . (yh-yy(i,i1))* nny(i,it)+
355 . (zh-zz(i,i1))* nnz(i,it) >= gaps(i)) far(i,it) =2
357 ELSEIF((ib1 /= 0 .AND. ib2 == 0).OR.
358 . (ib2 /= 0 .AND. ib1 == 0))
THEN
368 IF(vtx_bisector(1,1,ibx)/=zero.OR.
369 . vtx_bisector(2,1,ibx)/=zero.OR.
370 . vtx_bisector(3,1,ibx)/=zero.OR.
371 . vtx_bisector(1,2,ibx)/=zero.OR.
372 . vtx_bisector(2,2,ibx)/=zero.OR.
373 . vtx_bisector(3,2,ibx)/=zero)
THEN
374 p1 = (xh-xx(i,ix))* vtx_bisector(1,1,ibx)+
375 . (yh-yy(i,ix))* vtx_bisector(2,1,ibx)+
376 . (zh-zz(i,ix))* vtx_bisector(3,1,ibx)
377 p2 = (xh-xx(i,ix))* vtx_bisector(1,2,ibx)+
378 . (yh-yy(i,ix))* vtx_bisector(2,2,ibx)+
379 . (zh-zz(i,ix))* vtx_bisector(3,2,ibx)
380 IF(p1 >= gaps(i) .AND. p2 >= gaps(i)) far(i,it) =2
385 nn = one/
max(em20,sqrt(vx*vx+vy*vy+vz*vz))
387 IF(pn >= gaps(i)) far(i,it) =2
393 IF(far(i,it)==1 .OR. bb(i,it) <= zero)
THEN
395 x12(i)=xx(i,i2)-xx(i,i1)
396 y12(i)=yy(i,i2)-yy(i,i1)
397 z12(i)=zz(i,i2)-zz(i,i1)
400 px = z12(i)*nny(i,it)-y12(i)*nnz(i,it)
401 py = x12(i)*nnz(i,it)-z12(i)*nnx(i,it)
402 pz = y12(i)*nnx(i,it)-x12(i)*nny(i,it)
403 pp = px*px+py*py+pz*pz
405 ll = xij(i,i1)*xij(i,i1)+yij(i,i1)*yij(i,i1)+zij(i,i1)*zij(i,i1)
407 side(i,it)=-(xij(i,i1)*px+yij(i,i1)*py+zij(i,i1)*pz)*sqrt(one/
max(em30,ll*pp))
410 IF(side(i,it) < -zep01) far(i,it) =2
413 IF(far(i,it)==2) pent(i,it)=zero
421 d1=(xh-xx(i,1))* nnx(i,1)+
422 . (yh-yy(i,1))* nny(i,1)+
423 . (zh-zz(i,1))* nnz(i,1)
424 d2=(xh-xx(i,2))* nnx(i,2)+
425 . (yh-yy(i,2))* nny(i,2)+
426 . (zh-zz(i,2))* nnz(i,2)
427 d3=(xh-xx(i,3))* nnx(i,4)+
428 . (yh-yy(i,3))* nny(i,4)+
429 . (zh-zz(i,3))* nnz(i,4)
431 IF ( mvoisn(i,1)==0 .AND. d1 >= gaps(i))
THEN
433 ELSEIF( mvoisn(i,2)==0 .AND. d2 >= gaps(i))
THEN
435 ELSEIF( mvoisn(i,4)==0 .AND. d3 >= gaps(i))
THEN
437 ELSEIF((ib1/=0 .AND. ib2==0 .AND. ib3==0).OR.
438 . (ib2/=0 .AND. ib3==0 .AND. ib1==0).OR.
439 . (ib3/=0 .AND. ib1==0 .AND. ib2==0))
THEN
456 IF(vtx_bisector(1,1,ibx)/=zero.OR.
457 . vtx_bisector(2,1,ibx)/=zero.OR.
458 . vtx_bisector(3,1,ibx)/=zero.OR.
459 . vtx_bisector(1,2,ibx)/=zero.OR.
460 . vtx_bisector(2,2,ibx)/=zero.OR.
461 . vtx_bisector(3,2,ibx)/=zero)
THEN
462 p1 = (xh-xx(i,ix))* vtx_bisector(1,1,ibx)+
463 . (yh-yy(i,ix))* vtx_bisector(2,1,ibx)+
464 . (zh-zz(i,ix))* vtx_bisector(3,1,ibx)
465 p2 = (xh-xx(i,ix))* vtx_bisector(1,2,ibx)+
466 . (yh-yy(i,ix))* vtx_bisector(2,2,ibx)+
467 . (zh-zz(i,ix))* vtx_bisector(3,2,ibx)
468 IF(p1 >= gaps(i) .AND. p2 >= gaps(i)) far(i,it) =2
470 vx = two*xx(i,ix)-(xx(i,iy)+xx(i,iz))
471 vy = two*yy(i,ix)-(yy(i,iy)+yy(i,iz))
472 vz = two*zz(i,ix)-(zz(i,iy)+zz(i,iz))
473 nn = one/
max(em20,sqrt(vx*vx+vy*vy+vz*vz))
474 pn = ((xh-xx(i,ix))*vx+(yh-yy(i,ix))*vy+(zh-zz(i,ix))*vz)*nn
475 IF(pn >= gaps(i)) far(i,it) =2
480 IF(far(i,it)==1 .OR. bb(i,it) <= zero)
THEN
482 IF(mvoisn(i,1)/=0)
THEN
484 x12(i)=xx(i,2)-xx(i,1)
485 y12(i)=yy(i,2)-yy(i,1)
486 z12(i)=zz(i,2)-zz(i,1)
489 px = z12(i)*nny(i,1)-y12(i)*nnz(i,1)
490 py = x12(i)*nnz(i,1)-z12(i)*nnx(i,1)
491 pz = y12(i)*nnx(i,1)-x12(i)*nny(i,1)
492 pp = px*px+py*py+pz*pz
494 xi1(i) = xx(i,1) - xi(i)
495 yi1(i) = yy(i,1) - yi(i)
496 zi1(i) = zz(i,1) - zi(i)
497 ll = xi1(i)*xi1(i)+yi1(i)*yi1(i)+zi1(i)*zi1(i)
499 side(i,1)=-(xi1(i)*px+yi1(i)*py+zi1(i)*pz)*sqrt(one/
max(em30,ll*pp))
500 IF(side(i,1) < -zep01) far(i,1)=2
504 IF(mvoisn(i,2)/=0)
THEN
506 x12(i)=xx(i,3)-xx(i,2)
507 y12(i)=yy(i,3)-yy(i,2)
508 z12(i)=zz(i,3)-zz(i,2)
511 px = z12(i)*nny(i,2)-y12(i)*nnz(i,2)
512 py = x12(i)*nnz(i,2)-z12(i)*nnx(i,2)
513 pz = y12(i)*nnx(i,2)-x12(i)*nny(i,2)
514 pp = px*px+py*py+pz*pz
516 xi1(i) = xx(i,2) - xi(i)
517 yi1(i) = yy(i,2) - yi(i)
518 zi1(i) = zz(i,2) - zi(i)
519 ll = xi1(i)*xi1(i)+yi1(i)*yi1(i)+zi1(i)*zi1(i)
521 side(i,2)=-(xi1(i)*px+yi1(i)*py+zi1(i)*pz)*sqrt(one/
max(em30,ll*pp))
522 IF(side(i,2) < -zep01) far(i,2)=2
526 IF(mvoisn(i,4)/=0)
THEN
528 x12(i)=xx(i,1)-xx(i,3)
529 y12(i)=yy(i,1)-yy(i,3)
530 z12(i)=zz(i,1)-zz(i,3)
533 px = z12(i)*nny(i,4)-y12(i)*nnz(i,4)
534 py = x12(i)*nnz(i,4)-z12(i)*nnx(i,4)
535 pz = y12(i)*nnx(i,4)-x12(i)*nny(i,4)
536 pp = px*px+py*py+pz*pz
538 xi1(i) = xx(i,3) - xi(i)
539 yi1(i) = yy(i,3) - yi(i)
540 zi1(i) = zz(i,3) - zi(i)
541 ll = xi1(i)*xi1(i)+yi1(i)*yi1(i)+zi1(i)*zi1(i)
543 side(i,4)=-(xi1(i)*px+yi1(i)*py+zi1(i)*pz)*sqrt(one/
max(em30,ll*pp))
544 IF(side(i,4) < -zep01) far(i,3)=2
550 IF(far(i,1)==2.OR.far(i,2)==2.OR.far(i,3)==2) pent(i,it)=zero
558 IF(stif(i) <= zero)cycle
560 IF(ix3(i)/=ix4(i))
THEN
571#include "vectorize.inc"
575 x12(i) = xx(i,i2) - xx(i,i1)
576 y12(i) = yy(i,i2) - yy(i,i1)
577 z12(i) = zz(i,i2) - zz(i,i1)
581 lap = one-lbp(i,it)-lcp(i,it)
583 aaa = one /
max(em20,x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i))
584 hla= lap*abs(lap) * aaa
586 + hla<=hlb(i,it).AND.hla<=hlc(i,it))
THEN
587 lbp(i,it) = (xij(i,i2)*x12(i)+yij(i,i2)*y12(i)+zij(i,i2)*z12(i)) * aaa
588 lbp(i,it) =
max(zero,
min(one,lbp(i,it)))
589 lcp(i,it) = one - lbp(i,it)
590 ELSEIF(lbp(i,it)<zero.AND.
591 + hlb(i,it)<=hlc(i,it).AND.hlb(i,it)<=hla)
THEN
594 ELSEIF(lcp(i,it)<zero.AND.
595 + hlc(i,it)<=hla.AND.hlc(i,it)<=hlb(i,it))
THEN
600 lap = one-lbp(i,it)-lcp(i,it)
601 xp =lap*xx(i,5)+lbp(i,it)*xx(i,i1) + lcp(i,it)*xx(i,i2)
602 yp =lap*yy(i,5)+lbp(i,it)*yy(i,i1) + lcp(i,it)*yy(i,i2)
603 zp =lap*zz(i,5)+lbp(i,it)*zz(i,i1) + lcp(i,it)*zz(i,i2)
607 dd(i,it)=dx*dx+dy*dy+dz*dz
613#include "vectorize.inc"
617 lap = one-lbp(i,it)-lcp(i,it)
619 gap =
min(
max(gaps(i)+lap*gap_mm(i)+lbp(i,it)*gap_nm(i1,i)+lcp(i,it)*gap_nm(i2,i)+dgapload,drad),
620 .
max(drad,gapmxl(i)+dgapload))
623 bb(i,it) =((xx(i,5)-xi(i))*xn(i,it)+(yy(i,5)-yi(i))*yn(i,it)+(zz(i,5)-zi(i))*zn(i,it))
625 IF(bb(i,it) > zero)
THEN
628 pent(i,it)=
max(zero,gap+bb(i,it))
630 pent(i,it)=
max(zero,gap-sqrt(dd(i,it)))
635#include "vectorize.inc"
640 IF(lb(i,it) < -em03 .OR. lc(i,it) < -em03 .OR.
641 . lb(i,it)+lc(i,it) > one+em03) far(i,it)=1
643 xh=xi(i)+bb(i,it)*xn(i,it)
644 yh=yi(i)+bb(i,it)*yn(i,it)
645 zh=zi(i)+bb(i,it)*zn(i,it)
649 IF(mvoisn(i,it)==0)
THEN
651 IF( (xh-xx(i,i1))* nnx(i,it)+
652 . (yh-yy(i,i1))* nny(i,it)+
653 . (zh-zz(i,i1))* nnz(i,it) >= gaps(i)) far(i,it) =2
655 ELSEIF((ib1 /= 0 .AND. ib2 == 0).OR.
656 . (ib2 /= 0 .AND. ib1 == 0))
THEN
665 IF(vtx_bisector(1,1,ibx)/=zero.OR.
666 . vtx_bisector(2,1,ibx)/=zero.OR.
667 . vtx_bisector(3,1,ibx)/=zero.OR.
668 . vtx_bisector(1,2,ibx)/=zero.OR.
669 . vtx_bisector(2,2,ibx)/=zero.OR.
670 . vtx_bisector(3,2,ibx)/=zero)
THEN
671 p1 = (xh-xx(i,ix))* vtx_bisector(1,1,ibx)+
672 . (yh-yy(i,ix))* vtx_bisector(2,1,ibx)+
673 . (zh-zz(i,ix))* vtx_bisector(3,1,ibx)
674 p2 = (xh-xx(i,ix))* vtx_bisector(1,2,ibx)+
675 . (yh-yy(i,ix))* vtx_bisector(2,2,ibx)+
676 . (zh-zz(i,ix))* vtx_bisector(3,2,ibx)
677 IF(p1 >= gaps(i) .AND. p2 >= gaps(i)) far(i,it) =2
682 nn = one/
max(em20,sqrt(vx*vx+vy*vy+vz*vz))
683 pn = ((xh-xx(i,ix))*vx+(yh-yy(i,ix))*vy+(zh-zz(i,ix))*vz)*nn
684 IF(pn >= gaps(i)) far(i,it) =2
690 IF(far(i,it)==1 .OR. bb(i,it) <= zero)
THEN
694 x12(i)=xx(i,i2)-xx(i,i1)
695 y12(i)=yy(i,i2)-yy(i,i1)
696 z12(i)=zz(i,i2)-zz(i,i1)
699 px = z12(i)*nny(i,it)-y12(i)*nnz(i,it)
700 py = x12(i)*nnz(i,it)-z12(i)*nnx(i,it)
701 pz = y12(i)*nnx(i,it)-x12(i)*nny(i,it)
702 pp = px*px+py*py+pz*pz
704 ll = xij(i,i1)*xij(i,i1)+yij(i,i1)*yij(i,i1)+zij(i,i1)*zij(i,i1)
706 side(i,it)=-(xij(i,i1)*px+yij(i,i1)*py+zij(i,i1)*pz)*sqrt(one/
max(em30,ll*pp))
707 IF(side(i,it) < -zep01) far(i,it) =2
712 IF(far(i,it)==2) pent(i,it)=zero
718#include "vectorize.inc"
721 IF (stif(i) <= zero)cycle
723 subtria_n(i)= subtria(i)
725 IF(ix3(i)/=ix4(i))
THEN
726 dmin=
min(dd(i,1),dd(i,2),dd(i,3),dd(i,4))
728 IF(dd(i,jj) <= onep03*dmin)
THEN
731 lax = one-lb(i,jj)-lc(i,jj)
734 IF(lbx >= zero .AND. lcx >= zero)
THEN
737 lx =
max(zero,dd(i,jj)-bb(i,jj)*bb(i,jj))
751 IF(dd(i,1) <= dist(i))
THEN
759#include "vectorize.inc"
762 IF(stif(i) <= zero)cycle
767 IF(it/=0.AND.pent(i,it)==zero) it=0
770 IF(j /= it) pent(i,j)=zero
777 irtlm(2,n) = it+5*subtria(i)
780 time_s(1,n) = dist(i)
783 IF(ix3(i)/=ix4(i))
THEN
784 IF(pent(i,itold)==zero.OR.far(i,itold)>=2)
THEN
785 irtlm(2,n) = -irtlm(2,n)
788 IF(pent(i,itold)==zero.OR.far(i,1)>=2.OR.far(i,2)>=2.OR.far(i,3)>=2)
THEN
789 irtlm(2,n) = -irtlm(2,n)
804 irtlm_fi(nin)%P(2,n-nsn) = it+5*subtria(i)
807 time_sfi(nin)%P(2*(n-nsn-1)+1) = dist(i)
810 IF(ix3(i)/=ix4(i))
THEN
811 IF(pent(i,itold)==zero.OR.far(i,itold)>=2)
THEN
815 IF(pent(i,itold)==zero.OR.far(i,1)>=2.OR.far(i,2)>=2.OR.far(i,3)>=2)
THEN