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 . INTTH,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(*),ICONT_I(*), FAR(MVSIZ,4)
74 REAL*4 VTX_BISECTOR(3,2,*)
78 INTEGER I, J, K, L, JJ, INTERSECT, N, ITQ, N4N, IT, IKEEP, MGLOB, NINDX, ITOLD
79 INTEGER I4N(MVSIZ), INDX(MVSIZ), I1, I2, II,
80 . ITRIA(2,4), INEIGH(3,4), MARQUE(4,MVSIZ),
81 . SUBTRIA_N(MVSIZ), IT1, IB1, IB2, IB3, IBX, IX, IY, IZ
83 . XI1(MVSIZ), XI2(MVSIZ), XI0V(MVSIZ), XI5,
84 . yi1(mvsiz), yi2(mvsiz), yi0v(mvsiz), yi5,
85 . zi1(mvsiz), zi2(mvsiz), zi0v(mvsiz), zi5,
92 . xo1, xo2, xo3, xo4, xo5, xoi,
93 . yo1, yo2, yo3, yo4, yo5, yoi,
94 . zo1, zo2, zo3, zo4, zo5, zoi,
95 . vo12, vo23, vo34, vo41, pene
97 . gap_mm(mvsiz), gap, gap2,
100 . lap, hla, hlb(mvsiz,4), hlc(mvsiz,4),
107 . x0n(mvsiz,4), y0n(mvsiz,4), z0n(mvsiz,4),
108 . xn(mvsiz,4),yn(mvsiz,4),zn(mvsiz,4),
110 . x12(mvsiz), y12(mvsiz), z12(mvsiz),
111 . px, py, pz, pp, p1, p2, xh, yh, zh, d1, d2, d3, vx, vy, vz,
115 . lx, lax, lbx, lcx, ld(mvsiz), dmin
117 DATA ITRIA/1,2,2,3,3,4,4,1/,
118 . INEIGH/4,3,2,1,4,3,2,1,4,3,2,1/
123 marque(1:4,1:jlt) = 0
124 pent(1:jlt,1:4) = ep20
131 IF(stif(i) <= zero)cycle
133 x0n(i,1) = xx(i,1) - xx(i,5)
134 y0n(i,1) = yy(i,1) - yy(i,5)
135 z0n(i,1) = zz(i,1) - zz(i,5)
137 x0n(i,2) = xx(i,2) - xx(i,5)
138 y0n(i,2) = yy(i,2) - yy(i,5)
139 z0n(i,2) = zz(i,2) - zz(i,5)
141 x0n(i,3) = xx(i,3) - xx(i,5)
142 y0n(i,3) = yy(i,3) - yy(i,5)
143 z0n(i,3) = zz(i,3) - zz(i,5)
145 x0n(i,4) = xx(i,4) - xx(i,5)
146 y0n(i,4) = yy(i,4) - yy(i,5)
147 z0n(i,4) = zz(i,4) - zz(i,5)
149 IF(ix3(i)/=ix4(i))
THEN
150 gap_mm(i)=fourth*(gap_nm(1,i)+gap_nm(2,i)+gap_nm(3,i)+gap_nm(4,i))
152 gap_mm(i)=gap_nm(3,i)
159#include "vectorize.inc"
162 IF(stif(i) <= zero)cycle
164 xi0v(i) = xx(i,5) - xi(i)
165 yi0v(i) = yy(i,5) - yi(i)
166 zi0v(i) = zz(i,5) - zi(i)
168 xij(i,1) = xx(i,1) - xi(i)
169 yij(i,1) = yy(i,1) - yi(i)
170 zij(i,1) = zz(i,1) - zi(i)
172 xij(i,2) = xx(i,2) - xi(i)
173 yij(i,2) = yy(i,2) - yi(i)
174 zij(i,2) = zz(i,2) - zi(i)
176 xij(i,3) = xx(i,3) - xi(i)
177 yij(i,3) = yy(i,3) - yi(i)
178 zij(i,3) = zz(i,3) - zi(i)
180 xij(i,4) = xx(i,4) - xi(i)
181 yij(i,4) = yy(i,4) - yi(i)
182 zij(i,4) = zz(i,4) - zi(i)
184 sx1 = yi0v(i)*zij(i,1) - zi0v(i)*yij(i,1)
185 sy1 = zi0v(i)*xij(i,1) - xi0v(i)*zij(i,1)
186 sz1 = xi0v(i)*yij(i,1) - yi0v(i)*xij(i,1)
188 sx2 = yi0v(i)*zij(i,2) - zi0v(i)*yij(i,2)
189 sy2 = zi0v(i)*xij(i,2) - xi0v(i)*zij(i,2)
190 sz2 = xi0v(i)*yij(i,2) - yi0v(i)*xij(i,2)
192 xn(i,1) = y0n(i,1)*z0n(i,2) - z0n(i,1)*y0n(i,2)
193 yn(i,1) = z0n(i,1)*x0n(i,2) - x0n(i,1)*z0n(i,2)
194 zn(i,1) = x0n(i,1)*y0n(i,2) - y0n(i,1)*x0n(i,2)
196 .
max(em30,sqrt(xn(i,1)*xn(i,1)+ yn(i,1)*yn(i,1)+ zn(i,1)*zn(i,1)))
200 lb(i,1) = -(xn(i,1)*sx2 + yn(i,1)*sy2 + zn(i,1)*sz2)*nn
201 lc(i,1) = (xn(i,1)*sx1 + yn(i,1)*sy1 + zn(i,1)*sz1)*nn
203 sx3 = yi0v(i)*zij(i,3) - zi0v(i)*yij(i,3)
204 sy3 = zi0v(i)*xij(i,3) - xi0v(i)*zij(i,3)
205 sz3 = xi0v(i)*yij(i,3) - yi0v(i)*xij(i,3)
207 xn(i,2) = y0n(i,2)*z0n(i,3) - z0n(i,2)*y0n(i,3)
208 yn(i,2) = z0n(i,2)*x0n(i,3) - x0n(i,2)*z0n(i,3)
209 zn(i,2) = x0n(i,2)*y0n(i,3) - y0n(i,2)*x0n(i,3)
211 .
max(em30,sqrt(xn(i,2)*xn(i,2)+ yn(i,2)*yn(i,2)+ zn(i,2)*zn(i,2)))
215 lb(i,2) = -(xn(i,2)*sx3 + yn(i,2)*sy3 + zn
216 lc(i,2) = (xn(i,2)*sx2 + yn(i,2)*sy2 + zn(i,2)*sz2)*nn
218 sx4 = yi0v(i)*zij(i,4) - zi0v(i)*yij(i,4)
219 sy4 = zi0v(i)*xij(i,4) - xi0v(i)*zij(i,4)
220 sz4 = xi0v(i)*yij(i,4) - yi0v(i)*xij(i,4)
222 xn(i,3) = y0n(i,3)*z0n(i,4) - z0n(i,3)*y0n(i,4)
223 yn(i,3) = z0n(i,3)*x0n(i,4) - x0n(i,3)*z0n(i,4)
224 zn(i,3) = x0n(i,3)*y0n(i,4) - y0n(i,3)*x0n(i,4)
226 .
max(em30,sqrt(xn(i,3)*xn(i,3)+ yn(i,3)*yn(i,3)+ zn(i,3)*zn(i,3)))
230 lb(i,3) = -(xn(i,3)*sx4 + yn(i,3)*sy4 + zn(i,3)*sz4)*nn
231 lc(i,3) = (xn(i,3)*sx3 + yn(i,3)*sy3 + zn(i,3)*sz3)*nn
233 xn(i,4) = y0n(i,4)*z0n(i,1) - z0n(i,4)*y0n(i,1)
234 yn(i,4) = z0n(i,4)*x0n(i,1) - x0n(i,4)*z0n(i,1)
235 zn(i,4) = x0n(i,4)*y0n(i,1) - y0n(i,4)*x0n(i,1)
237 .
max(em30,sqrt(xn(i,4)*xn(i,4)+ yn(i,4)*yn(i,4)+ zn(i,4)*zn(i,4)))
241 lb(i,4) = -(xn(i,4)*sx1 + yn(i,4)*sy1 + zn(i,4)*sz1)*nn
242 lc(i,4) = (xn(i,4)*sx4 + yn(i,4)*sy4 + zn(i,4)*sz4)*nn
246#include "vectorize.inc"
249 IF(stif(i) <= zero)cycle
251 aaa = one/
max(em30,x0n(i,1)*x0n(i,1)+y0n(i,1)*y0n(i,1)+z0n(i,1)*z0n(i,1))
252 hlc(i,1) = lc(i,1)*abs(lc(i,1))*aaa
253 hlb(i,4) = lb(i,4)*abs(lb(i,4))*aaa
254 al(i,1) = -(xi0v(i)*x0n(i,1)+yi0v(i)*y0n(i,1)+zi0v(i)*z0n(i,1))*aaa
255 al(i,1) =
max(zero,
min(one,al(i,1)))
256 aaa = one/
max(em30,x0n(i,2)*x0n(i,2)+y0n(i,2)*y0n(i,2)+z0n(i,2)*z0n(i,2))
257 hlc(i,2) = lc(i,2)*abs(lc(i,2))*aaa
258 hlb(i,1) = lb(i,1)*abs(lb(i,1))*aaa
259 al(i,2) = -(xi0v(i)*x0n(i,2)+yi0v(i)*y0n(i,2)+zi0v(i)*z0n(i,2))*aaa
260 al(i,2) =
max(zero,
min(one,al(i,2)))
261 aaa = one/
max(em30,x0n(i,3)*x0n(i,3)+y0n(i,3)*y0n(i,3)+z0n(i,3)*z0n(i,3))
262 hlc(i,3) = lc(i,3)*abs(lc(i,3))*aaa
263 hlb(i,2) = lb(i,2)*abs(lb(i,2))*aaa
264 al(i,3) = -(xi0v(i)*x0n(i,3)+yi0v(i)*y0n(i,3)+zi0v(i)*z0n(i,3))*aaa
265 al(i,3) =
max(zero,
min(one,al(i,3)))
266 aaa = one/
max(em30,x0n(i,4)*x0n(i,4)+y0n(i,4)*y0n(i,4)+z0n(i,4)*z0n(i,4))
267 hlc(i,4) = lc(i,4)*abs(lc(i,4))*aaa
268 hlb(i,3) = lb(i,3)*abs(lb(i,3))*aaa
269 al(i,4) = -(xi0v(i)*x0n(i,4)+yi0v(i)*y0n(i,4)+zi0v(i)*z0n(i,4))*aaa
270 al(i,4) =
max(zero,
min(one,al(i,4)))
277#include "vectorize.inc"
280 IF(stif(i) <= zero)cycle
282 x12(i) = xx(i,i2) - xx(i,i1)
283 y12(i) = yy(i,i2) - yy(i,i1)
284 z12(i) = zz(i,i2) - zz(i,i1)
288 lap = one-lbp(i,it)-lcp(i,it)
290 aaa = one /
max(em20,x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i))
291 hla= lap*abs(lap) * aaa
293 + hla<=hlb(i,it).AND.hla<=hlc(i,it))
THEN
294 lbp(i,it) = (xij(i,i2)*x12(i)+yij(i,i2)*y12(i)+zij(i,i2)*z12(i)) * aaa
295 lbp(i,it) =
max(zero,
min(one,lbp(i,it)))
296 lcp(i,it) = one - lbp(i,it)
297 ELSEIF(lbp(i,it)<zero.AND.
298 + hlb(i,it)<=hlc(i,it).AND.hlb(i,it)<=hla)
THEN
301 ELSEIF(lcp(i,it)<zero.AND.
302 + hlc(i,it)<=hla.AND.hlc(i,it)<=hlb(i,it))
THEN
307 lap = one-lbp(i,it)-lcp(i,it)
308 xp =lap*xx(i,5)+lbp(i,it)*xx(i,i1) + lcp(i,it)*xx(i,i2)
309 yp =lap*yy(i,5)+lbp(i,it)*yy(i,i1) + lcp(i,it)*yy(i,i2)
310 zp =lap*zz(i,5)+lbp(i,it)*zz(i,i1) + lcp(i,it)*zz(i,i2)
314 dd(i,it)=dx*dx+dy*dy+dz*dz
320#include "vectorize.inc"
323 IF (stif(i) <= zero)cycle
325 lap = one-lbp(i,it)-lcp(i,it)
327 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),
328 .
max(drad,gapmxl(i)+dgapload))
331 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))
333 IF(bb(i,it) > zero)
THEN
336 pent(i,it)=
max(zero,gap+bb(i,it))
338 pent(i,it)=
max(zero,gap-sqrt(dd(i,it)))
343#include "vectorize.inc"
346 IF (stif(i) <= zero)cycle
349 IF(lb(i,it) < -em03 .OR. lc(i,it) < -em03 .OR.
350 . lb(i,it)+lc(i,it) > one+em03) far(i,it)=1
352 xh=xi(i)+bb(i,it)*xn(i,it)
353 yh=yi(i)+bb(i,it)*yn(i,it)
354 zh=zi(i)+bb(i,it)*zn(i,it)
356 IF(ix3(i) /= ix4(i))
THEN
360 IF(mvoisn(i,it)==0)
THEN
362 IF( (xh-xx(i,i1))* nnx(i,it)+
363 . (yh-yy(i,i1))* nny(i,it)+
364 . (zh-zz(i,i1))* nnz(i,it) >= gaps(i)) far(i,it) =2
366 ELSEIF((ib1 /= 0 .AND. ib2 == 0).OR.
367 . (ib2 /= 0 .AND. ib1 == 0))
THEN
377 IF(vtx_bisector(1,1,ibx)/=zero.OR.
378 . vtx_bisector(2,1,ibx)/=zero.OR.
379 . vtx_bisector(3,1,ibx)/=zero.OR.
380 . vtx_bisector(1,2,ibx)/=zero.OR.
381 . vtx_bisector(2,2,ibx)/=zero.OR.
382 . vtx_bisector(3,2,ibx)/=zero)
THEN
383 p1 = (xh-xx(i,ix))* vtx_bisector(1,1,ibx)+
384 . (yh-yy(i,ix))* vtx_bisector(2,1,ibx)+
385 . (zh-zz(i,ix))* vtx_bisector(3,1,ibx)
386 p2 = (xh-xx(i,ix))* vtx_bisector(1,2,ibx)+
387 . (yh-yy(i,ix))* vtx_bisector(2,2,ibx)+
388 . (zh-zz(i,ix))* vtx_bisector(3,2,ibx)
389 IF(p1 >= gaps(i) .AND. p2 >= gaps(i)) far(i,it) =2
394 nn = one/
max(em20,sqrt(vx*vx+vy*vy+vz*vz))
395 pn = ((xh-xx(i,ix))*vx+(yh-yy(i,ix))*vy+(zh-zz(i,ix))*vz)*nn
396 IF(pn >= gaps(i)) far(i,it) =2
402 IF(far(i,it)==1 .OR. bb(i,it) <= zero)
THEN
404 x12(i)=xx(i,i2)-xx(i,i1)
405 y12(i)=yy(i,i2)-yy(i,i1)
406 z12(i)=zz(i,i2)-zz(i,i1)
409 px = z12(i)*nny(i,it)-y12(i)*nnz(i,it)
410 py = x12(i)*nnz(i,it)-z12(i)*nnx(i,it)
411 pz = y12(i)*nnx(i,it)-x12(i)*nny(i,it)
412 pp = px*px+py*py+pz*pz
414 ll = xij(i,i1)*xij(i,i1)+yij(i,i1)*yij(i,i1)+zij(i,i1)*zij(i,i1)
416 side(i,it)=-(xij(i,i1)*px+yij(i,i1)*py+zij(i,i1)*pz)*sqrt(one/
max(em30,ll*pp))
419 IF(side(i,it) < -zep01) far(i,it) =2
422 IF(far(i,it)==2) pent(i,it)=zero
430 d1=(xh-xx(i,1))* nnx(i,1)+
431 . (yh-yy(i,1))* nny(i,1)+
432 . (zh-zz(i,1))* nnz(i,1)
433 d2=(xh-xx(i,2))* nnx(i,2)+
434 . (yh-yy(i,2))* nny(i,2)+
435 . (zh-zz(i,2))* nnz(i,2)
436 d3=(xh-xx(i,3))* nnx(i,4)+
437 . (yh-yy(i,3))* nny(i,4)+
438 . (zh-zz(i,3))* nnz(i,4)
440 IF ( mvoisn(i,1)==0 .AND. d1 >= gaps(i))
THEN
442 ELSEIF( mvoisn(i,2)==0 .AND. d2 >= gaps(i))
THEN
444 ELSEIF( mvoisn(i,4)==0 .AND. d3 >= gaps(i))
THEN
446 ELSEIF((ib1/=0 .AND. ib2==0 .AND. ib3==0).OR.
447 . (ib2/=0 .AND. ib3==0 .AND. ib1==0).OR.
448 . (ib3/=0 .AND. ib1==0 .AND. ib2==0))
THEN
465 IF(vtx_bisector(1,1,ibx)/=zero.OR.
466 . vtx_bisector(2,1,ibx)/=zero.OR.
467 . vtx_bisector(3,1,ibx)/=zero.OR.
468 . vtx_bisector(1,2,ibx)/=zero.OR.
469 . vtx_bisector(2,2,ibx)/=zero.OR.
470 . vtx_bisector(3,2,ibx)/=zero)
THEN
471 p1 = (xh-xx(i,ix))* vtx_bisector(1,1,ibx)+
472 . (yh-yy(i,ix))* vtx_bisector(2,1,ibx)+
473 . (zh-zz(i,ix))* vtx_bisector(3,1,ibx)
474 p2 = (xh-xx(i,ix))* vtx_bisector(1,2,ibx)+
475 . (yh-yy(i,ix))* vtx_bisector(2,2,ibx)+
476 . (zh-zz(i,ix))* vtx_bisector(3,2,ibx)
477 IF(p1 >= gaps(i) .AND. p2 >= gaps(i)) far(i,it) =2
479 vx = two*xx(i,ix)-(xx(i,iy)+xx(i,iz))
480 vy = two*yy(i,ix)-(yy(i,iy)+yy(i,iz))
481 vz = two*zz(i,ix)-(zz(i,iy)+zz(i,iz))
482 nn = one/
max(em20,sqrt(vx*vx+vy*vy+vz*vz))
483 pn = ((xh-xx(i,ix))*vx+(yh-yy(i,ix))*vy+(zh-zz(i,ix))*vz)*nn
484 IF(pn >= gaps(i)) far(i,it) =2
489 IF(far(i,it)==1 .OR. bb(i,it) <= zero)
THEN
491 IF(mvoisn(i,1)/=0)
THEN
493 x12(i)=xx(i,2)-xx(i,1)
494 y12(i)=yy(i,2)-yy(i,1)
495 z12(i)=zz(i,2)-zz(i,1)
498 px = z12(i)*nny(i,1)-y12(i)*nnz(i,1)
499 py = x12(i)*nnz(i,1)-z12(i)*nnx(i,1)
500 pz = y12(i)*nnx(i,1)-x12(i)*nny(i,1)
501 pp = px*px+py*py+pz*pz
503 xi1(i) = xx(i,1) - xi(i)
504 yi1(i) = yy(i,1) - yi(i)
505 zi1(i) = zz(i,1) - zi(i)
506 ll = xi1(i)*xi1(i)+yi1(i)*yi1(i)+zi1(i)*zi1(i)
508 side(i,1)=-(xi1(i)*px+yi1(i)*py+zi1(i)*pz)*sqrt(one/
max(em30
513 IF(mvoisn(i,2)/=0)
THEN
515 x12(i)=xx(i,3)-xx(i,2)
516 y12(i)=yy(i,3)-yy(i,2)
517 z12(i)=zz(i,3)-zz(i,2)
520 px = z12(i)*nny(i,2)-y12(i)*nnz(i,2)
521 py = x12(i)*nnz(i,2)-z12(i)*nnx(i,2)
522 pz = y12(i)*nnx(i,2)-x12(i)*nny(i,2)
523 pp = px*px+py*py+pz*pz
525 xi1(i) = xx(i,2) - xi(i)
526 yi1(i) = yy(i,2) - yi(i)
527 zi1(i) = zz(i,2) - zi(i)
528 ll = xi1(i)*xi1(i)+yi1(i)*yi1(i)+zi1(i)*zi1(i)
530 side(i,2)=-(xi1(i)*px+yi1(i)*py+zi1(i)*pz)*sqrt(one/
max(em30,ll*pp))
531 IF(side(i,2) < -zep01) far(i,2)=2
535 IF(mvoisn(i,4)/=0)
THEN
537 x12(i)=xx(i,1)-xx(i,3)
538 y12(i)=yy(i,1)-yy(i,3)
539 z12(i)=zz(i,1)-zz(i,3)
542 px = z12(i)*nny(i,4)-y12(i)*nnz(i,4)
543 py = x12(i)*nnz(i,4)-z12(i)*nnx(i,4)
544 pz = y12(i)*nnx(i,4)-x12(i)*nny(i,4)
545 pp = px*px+py*py+pz*pz
547 xi1(i) = xx(i,3) - xi(i)
548 yi1(i) = yy(i,3) - yi(i)
549 zi1(i) = zz(i,3) - zi(i)
550 ll = xi1(i)*xi1(i)+yi1(i)*yi1(i)+zi1(i)*zi1(i)
552 side(i,4)=-(xi1(i)*px+yi1(i)*py+zi1(i)*pz)*sqrt(one/
max(em30,ll*pp))
553 IF(side(i,4) < -zep01) far(i,3)=2
559 IF(far(i,1)==2.OR.far(i,2)==2.OR.far(i,3)==2) pent(i,it)=zero
567 IF(stif(i) <= zero)cycle
569 IF(ix3(i)/=ix4(i))
THEN
580#include "vectorize.inc"
584 x12(i) = xx(i,i2) - xx(i,i1)
585 y12(i) = yy(i,i2) - yy(i,i1)
586 z12(i) = zz(i,i2) - zz(i,i1)
590 lap = one-lbp(i,it)-lcp(i,it)
592 aaa = one /
max(em20,x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i))
593 hla= lap*abs(lap) * aaa
595 + hla<=hlb(i,it).AND.hla<=hlc(i,it))
THEN
596 lbp(i,it) = (xij(i,i2)*x12(i)+yij(i,i2)*y12(i)+zij(i,i2)*z12(i)) * aaa
597 lbp(i,it) =
max(zero,
min(one,lbp(i,it)))
598 lcp(i,it) = one - lbp(i,it)
599 ELSEIF(lbp(i,it)<zero.AND.
600 + hlb(i,it)<=hlc(i,it).AND.hlb(i,it)<=hla)
THEN
603 ELSEIF(lcp(i,it)<zero.AND.
604 + hlc(i,it)<=hla.AND.hlc(i,it)<=hlb(i,it))
THEN
609 lap = one-lbp(i,it)-lcp(i,it)
610 xp =lap*xx(i,5)+lbp(i,it)*xx(i,i1) + lcp(i,it)*xx(i,i2)
611 yp =lap*yy(i,5)+lbp(i,it)*yy(i,i1) + lcp(i,it)*yy(i,i2)
612 zp =lap*zz(i,5)+lbp(i,it)*zz(i,i1) + lcp(i,it)*zz(i,i2)
616 dd(i,it)=dx*dx+dy*dy+dz*dz
622#include "vectorize.inc"
626 lap = one-lbp(i,it)-lcp(i,it)
628 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),
629 .
max(drad,gapmxl(i)+dgapload))
632 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))
634 IF(bb(i,it) > zero)
THEN
637 pent(i,it)=
max(zero,gap+bb(i,it))
639 pent(i,it)=
max(zero,gap-sqrt(dd(i,it)))
644#include "vectorize.inc"
649 IF(lb(i,it) < -em03 .OR. lc(i,it) < -em03 .OR.
650 . lb(i,it)+lc(i,it) > one+em03) far(i,it)=1
652 xh=xi(i)+bb(i,it)*xn(i,it)
653 yh=yi(i)+bb(i,it)*yn(i,it)
654 zh=zi(i)+bb(i,it)*zn(i,it)
658 IF(mvoisn(i,it)==0)
THEN
660 IF( (xh-xx(i,i1))* nnx(i,it)+
661 . (yh-yy(i,i1))* nny(i,it)+
662 . (zh-zz(i,i1))* nnz(i,it) >= gaps(i)) far(i,it) =2
664 ELSEIF((ib1 /= 0 .AND. ib2 == 0).OR.
665 . (ib2 /= 0 .AND. ib1 == 0))
THEN
674 IF(vtx_bisector(1,1,ibx)/=zero.OR.
675 . vtx_bisector(2,1,ibx)/=zero.OR.
676 . vtx_bisector(3,1,ibx)/=zero.OR.
677 . vtx_bisector(1,2,ibx)/=zero.OR.
678 . vtx_bisector(2,2,ibx)/=zero.OR.
679 . vtx_bisector(3,2,ibx)/=zero)
THEN
680 p1 = (xh-xx(i,ix))* vtx_bisector(1,1,ibx)+
681 . (yh-yy(i,ix))* vtx_bisector(2,1,ibx)+
682 . (zh-zz(i,ix))* vtx_bisector(3,1,ibx)
683 p2 = (xh-xx(i,ix))* vtx_bisector(1,2,ibx)+
684 . (yh-yy(i,ix))* vtx_bisector(2,2,ibx)+
685 . (zh-zz(i,ix))* vtx_bisector(3,2,ibx)
686 IF(p1 >= gaps(i) .AND. p2 >= gaps(i)) far
691 nn = one/
max(em20,sqrt(vx*vx+vy*vy+vz*vz))
692 pn = ((xh-xx(i,ix))*vx+(yh-yy(i,ix))*vy+(zh-zz(i,ix))*vz)*nn
693 IF(pn >= gaps(i)) far(i,it) =2
699 IF(far(i,it)==1 .OR. bb(i,it) <= zero)
THEN
703 x12(i)=xx(i,i2)-xx(i,i1)
704 y12(i)=yy(i,i2)-yy(i,i1)
705 z12(i)=zz(i,i2)-zz(i,i1)
708 px = z12(i)*nny(i,it)-y12(i)*nnz(i,it)
709 py = x12(i)*nnz(i,it)-z12(i
710 pz = y12(i)*nnx(i,it)-x12(i)*nny(i,it)
711 pp = px*px+py*py+pz*pz
713 ll = xij(i,i1)*xij(i,i1)+yij(i,i1)*yij(i,i1)+zij(i,i1)*zij(i,i1)
715 side(i,it)=-(xij(i,i1)*px+yij(i,i1)*py+zij(i,i1)*pz)*sqrt(one/
max(em30,ll
716 IF(side(i,it) < -zep01) far(i,it) =2
721 IF(far(i,it)==2) pent(i,it)=zero
727#include "vectorize.inc"
730 IF (stif(i) <= zero)cycle
732 subtria_n(i)= subtria(i)
734 IF(ix3(i)/=ix4(i))
THEN
735 dmin=
min(dd(i,1),dd(i,2),dd(i,3),dd(i,4))
737 IF(dd(i,jj) <= onep03*dmin)
THEN
740 lax = one-lb(i,jj)-lc(i,jj)
743 IF(lbx >= zero .AND. lcx >= zero)
THEN
746 lx =
max(zero,dd(i,jj)-bb(i,jj)*bb(i,jj))
760 IF(dd(i,1) <= dist(i))
THEN
768#include "vectorize.inc"
771 IF(stif(i) <= zero)cycle
776 IF(it/=0.AND.pent(i,it)==zero) it=0
779 IF(j /= it) pent(i,j)=zero
786 irtlm(2,n) = it+5*subtria(i)
789 time_s(1,n) = dist(i)
792 IF(ix3(i)/=ix4(i))
THEN
793 IF(pent(i,itold)==zero.OR.far(i,itold)>=2)
THEN
794 irtlm(2,n) = -irtlm(2,n)
797 IF(pent(i,itold)==zero.OR.far(i,1)>=2.OR.far(i,2)>=2.OR.far(i,3)>=2)
THEN
798 irtlm(2,n) = -irtlm(2,n)
813 irtlm_fi(nin)%P(2,n-nsn) = it+5*subtria(i)
816 time_sfi(nin)%P(2*(n-nsn-1)+1) = dist(i)
819 IF(ix3(i)/=ix4(i))
THEN
820 IF(pent(i,itold)==zero.OR.far(i,itold)>=2)
THEN
824 IF(pent(i,itold)==zero.OR.far(i,1)>=2.OR.far(i,2)>=2.OR.far(i,3)>=2)
THEN