47
48
49
50#include "implicit_f.inc"
51
52
53
54#include "mvsiz_p.inc"
55
56
57
58 INTEGER IGAP,CAND_N(*),CAND_E(*),IRECT(4,*),NLN,NLG(NLN),NSV(*),
59 . NBINFLG(*),TAG(*)
61 . gap,gap_sh(*),gapv(mvsiz),gap_s(*),gap_m(*),gapmax,gapmin,
62 . solidn_normal(3,*)
63 INTEGER, DIMENSION(MVSIZ), INTENT(IN) :: IX3,IX4
64 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: x1,x2,x3,x4
65 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: y1,y2,y3,y4
66 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: z1,z2,z3,z4
67 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: xi,yi,zi
68 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: x0,y0,z0
69 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: nx1,ny1,nz1
70 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: nx2,ny2,nz2
71 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: nx3,ny3,nz3
72 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: nx4,ny4,nz4
73 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: p1,p2,p3,p4
74 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: lb1,lb2,lb3,lb4
75 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: lc1,lc2,lc3,lc4
76
77
78
79#include "vect07_c.inc"
80
81
82
83 INTEGER I,I1,I2,I3,I4,IG,IM,IS,IB1,IB2,IB3,IB4
84
86 . al1(mvsiz), al2(mvsiz), al3(mvsiz), al4(mvsiz),
87 . x01(mvsiz), x02(mvsiz), x03(mvsiz), x04(mvsiz),
88 . y01(mvsiz), y02(mvsiz), y03(mvsiz), y04(mvsiz),
89 . z01(mvsiz), z02(mvsiz), z03(mvsiz), z04(mvsiz),
90 . xi1(mvsiz), xi2(mvsiz), xi3(mvsiz), xi4(mvsiz),
91 . yi1(mvsiz), yi2(mvsiz), yi3(mvsiz), yi4(mvsiz),
92 . zi1(mvsiz), zi2(mvsiz), zi3(mvsiz), zi4(mvsiz),
93 . hlb1(mvsiz), hlc1(mvsiz), hlb2(mvsiz),hlc2(mvsiz),
94 . hlb3(mvsiz), hlc3(mvsiz), hlb4(mvsiz),hlc4(mvsiz)
96 . s2,a1,a2,a3,a4,d1,d2,d3,d4,
97 . x12,x23,x34,x41,xi0,sx1,sx2,sx3,sx4,sx0,sxi,
98 . y12,y23,y34,y41,yi0,sy1,sy2,sy3,sy4,sy0,syi,
99 . z12,z23,z34,z41,zi0,sz1,sz2,sz3,sz4,sz0,szi,
100 . x10,y10,z10,x20,y20,z20,x30,y30,z30,x40,y40,z40,
101 . la, hla, aaa, unssqr3
102
103 INTEGER BITUNSET,BITGET,BITSET
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149 unssqr3 = one/sqr3
150 IF(igap /= 0)THEN
151 DO i=lft,llt
152 im = cand_e(i)
153 is = cand_n(i)
154
155 ig = nsv(is)
156 gapv(i) = gap_s(is) + gap_m(im)
157 gapv(i) =
min(gapv(i),gapmax)
158 gapv(i) = gapv(i) + gap_sh(im)*(one-em5)
159 gapv(i) =
max(gapmin,gapv(i))
160
161
162
163
164 i1 = irect(1,im)
165 i2 = irect(2,im)
166 i3 = irect(3,im)
167 i4 = irect(4,im)
168
169 ib1 =
bitget(nbinflg(tag(i1)),7)
170 ib2 =
bitget(nbinflg(tag(i2)),7)
171 ib3 =
bitget(nbinflg(tag(i3)),7)
172 ib4
173
174 IF(ib1+ib2+ib3+ib4 == 0)THEN
175 IF(gap_sh(im)/= zero)THEN
176
177 sx0=(y3(i)-y1(i))*(z4(i)-z2(i))-(z3(i)-z1(i))*(y4(i)-y2(i))
178 sy0=(z3(i)-z1(i))*(x4(i)-x2(i))-(x3(i)-x1(i))*(z4(i)-z2(i))
179 sz0=(x3(i)-x1(i))*(y4(i)-y2(i))-(y3(i)-y1(i))*(x4(i)-x2(i))
180 aaa = one / sqrt(
max(em20,sx0*sx0+sy0*sy0+sz0*sz0))
181 sx0=sx0*aaa
182 sy0=sy0*aaa
183 sz0=sz0*aaa
184 sxi = solidn_normal(1,ig)
185 syi = solidn_normal(2,ig)
186 szi = solidn_normal(3,ig)
187 aaa = sxi*sx0+syi*sy0+szi*sz0
188 IF(aaa > zero)THEN
189 sxi = zero
190 syi = zero
191 szi = zero
192 ENDIF
193
194 sx1 = solidn_normal(1,i1)-sxi
195 sy1 = solidn_normal(2,i1)-syi
196 sz1 = solidn_normal(3,i1)-szi
197 aaa = one / sqrt(
max(em20,sx1*sx1+sy1*sy1+sz1*sz1))
198 sx1 = sx1*aaa
199 sy1 = sy1*aaa
200 sz1 = sz1*aaa
201 aaa = sx0*sx1 + sy0*sy1 + sz0*sz1
202 IF(aaa > unssqr3)THEN
203 aaa = gap_sh(im)/aaa
204 ELSE
205 aaa = unssqr3 - aaa
206 sx1 = sx1 + aaa*sx0
207 sy1 = sy1 + aaa*sy0
208 sz1 = sz1 + aaa*sz0
209 aaa = gap_sh(im)*sqr3
210 ENDIF
211 sx1 = sx1*aaa
212 sy1 = sy1*aaa
213 sz1 = sz1*aaa
214 sx2 = solidn_normal(1,i2)-sxi
215 sy2 = solidn_normal(2,i2)-syi
216 sz2 = solidn_normal(3,i2)-szi
217 aaa = one / sqrt(
max(em20,sx2*sx2+sy2*sy2+sz2*sz2))
218 sx2 = sx2*aaa
219 sy2 = sy2*aaa
220 sz2 = sz2*aaa
221 aaa = sx0*sx2 + sy0*sy2 + sz0*sz2
222 IF(aaa > unssqr3)THEN
223 aaa = gap_sh(im)/aaa
224 ELSE
225 aaa = unssqr3 - aaa
226 sx2 = sx2 + aaa*sx0
227 sy2 = sy2 + aaa*sy0
228 sz2 = sz2 + aaa*sz0
229 aaa = gap_sh(im)*sqr3
230 ENDIF
231 sx2 = sx2*aaa
232 sy2 = sy2*aaa
233 sz2 = sz2*aaa
234 sx3 = solidn_normal(1,i3)-sxi
235 sy3 = solidn_normal(2,i3)-syi
236 sz3 = solidn_normal(3,i3)-szi
237 aaa = one / sqrt(
max(em20,sx3*sx3+sy3*sy3+sz3*sz3))
238 sx3 = sx3*aaa
239 sy3 = sy3*aaa
240 sz3 = sz3*aaa
241 aaa = sx0*sx3 + sy0*sy3 + sz0*sz3
242 IF(aaa > unssqr3)THEN
243 aaa = gap_sh(im)/aaa
244 ELSE
245 aaa = unssqr3 - aaa
246 sx3 = sx3 + aaa*sx0
247 sy3 = sy3 + aaa*sy0
248 sz3 = sz3 + aaa*sz0
249 aaa = gap_sh(im)*sqr3
250 ENDIF
251 sx3 = sx3*aaa
252 sy3 = sy3*aaa
253 sz3 = sz3*aaa
254 sx4 = solidn_normal(1,i4)-sxi
255 sy4 = solidn_normal(2,i4)-syi
256 sz4 = solidn_normal(3,i4)-szi
257 aaa = one / sqrt(
max(em20,sx4*sx4+sy4*sy4+sz4*sz4))
258 sx4 = sx4*aaa
259 sy4 = sy4*aaa
260 sz4 = sz4*aaa
261 aaa = sx0*sx4 + sy0*sy4 + sz0*sz4
262 IF(aaa > unssqr3)THEN
263 aaa = gap_sh(im)/aaa
264 ELSE
265 aaa = unssqr3 - aaa
266 sx4 = sx4 + aaa*sx0
267 sy4 = sy4 + aaa*sy0
268 sz4 = sz4 + aaa*sz0
269 aaa = gap_sh(im)*sqr3
270 ENDIF
271 sx4 = sx4*aaa
272 sy4 = sy4*aaa
273 sz4 = sz4*aaa
274 x1(i) = x1(i) - sx1
275 y1(i) = y1(i) - sy1
276 z1(i) = z1(i) - sz1
277 x2(i) = x2(i) - sx2
278 y2(i) = y2(i) - sy2
279 z2(i) = z2(i) - sz2
280 x3(i) = x3(i) - sx3
281 y3(i) = y3(i) - sy3
282 z3(i) = z3(i) - sz3
283 x4(i) = x4(i) - sx4
284 y4(i) = y4(i) - sy4
285 z4(i) = z4(i) - sz4
286 ENDIF
287 ELSEIF(i3 == i4)THEN
288
289 x10 = x1(i)
290 y10 = y1(i)
291 z10 = z1(i)
292 x20 = x2(i)
293 y20 = y2(i)
294 z20 = z2(i)
295 x30 = x3(i)
296 y30 = y3(i)
297 z30 = z3(i)
298 IF(ib1 == 1 .and. ib2 == 1)THEN
299 CALL i20cgap0(x10,y10,z10,x20,y20,z20,
300 . x30,y30,z30,x30,y30,z30,
301 . x1(i),y1(i),z1(i),x2(i),y2(i),z2(i),
302 . x3(i),y3(i),z3(i),x3(i),y3(i),z3(i),
303 . gap_m(im))
304 ENDIF
305 IF(ib2 == 1 .and. ib3 == 1)THEN
306 CALL i20cgap0(x20,y20,z20,x30,y30,z30,
307 . x10,y10,z10,x10,y10,z10,
308 . x2(i),y2(i),z2(i),x3(i),y3(i),z3(i),
309 . x1(i),y1(i),z1(i),x1(i),y1(i),z1(i),
310 . gap_m(im))
311 ENDIF
312 IF(ib3 == 1 .and. ib1 == 1)THEN
313 CALL i20cgap0(x30,y30,z30,x10,y10,z10,
314 . x20,y20,z20,x20,y20,z20,
315 . x3(i),y3(i),z3(i),x1(i),y1(i),z1(i),
316 . x2(i),y2(i),z2(i),x2(i),y2(i),z2(i),
317 . gap_m(im))
318 ENDIF
319 IF(ib1 == 1 .and. ib2+ib3 == 0)THEN
320 CALL i20cgap1(x10,y10,z10,x20,y20,z20,
321 . x30,y30,z30,
322 . x1(i),y1(i),z1(i),x2(i),y2(i),z2(i),
323 . x3(i),y3(i),z3(i),gap_m(im))
324 ENDIF
325 IF(ib2 == 1 .and. ib3+ib1 == 0)THEN
326 CALL i20cgap1(x20,y20,z20,x30,y30,z30,
327 . x10,y10,z10,
328 . x2(i),y2(i),z2(i),x3(i),y3(i),z3(i),
329 . x1(i),y1(i),z1(i),gap_m(im))
330 ENDIF
331 IF(ib3 == 1 .and. ib1+ib2 == 0)THEN
332 CALL i20cgap1(x30,y30,z30,x10,y10,z10,
333 . x20,y20,z20,
334 . x3(i),y3(i),z3(i),x1(i),y1(i),z1(i),
335 . x2(i),y2(i),z2(i),gap_m(im))
336 ENDIF
337 x4(i)=x3(i)
338 y4(i)=y3(i)
339 z4(i)=z3(i)
340 ELSE
341
342 x10 = x1(i)
343 y10 = y1(i)
344 z10 = z1(i)
345 x20 = x2(i)
346 y20 = y2(i)
347 z20 = z2(i)
348 x30 = x3(i)
349 y30 = y3(i)
350 z30 = z3(i)
351 x40 = x4(i)
352 y40 = y4(i)
353 z40 = z4(i)
354 IF(ib1 == 1 .and. ib2 == 1)THEN
355 CALL i20cgap0(x10,y10,z10,x20,y20,z20,
356 . x30,y30,z30,x40,y40,z40,
357 . x1(i),y1(i),z1(i),x2(i),y2(i),z2(i),
358 . x3(i),y3(i),z3(i),x4(i),y4(i),z4(i),
359 . gap_m(im))
360 ENDIF
361 IF(ib2 == 1 .and. ib3 == 1)THEN
362 CALL i20cgap0(x20,y20,z20,x30,y30,z30,
363 . x40,y40,z40,x10,y10,z10,
364 . x2(i),y2(i),z2(i),x3(i),y3(i),z3(i),
365 . x4(i),y4(i),z4(i),x1(i),y1(i),z1(i),
366 . gap_m(im))
367 ENDIF
368 IF(ib3 == 1 .and. ib4 == 1)THEN
369 CALL i20cgap0(x30,y30,z30,x40,y40,z40,
370 . x10,y10,z10,x20,y20,z20,
371 . x3(i),y3(i),z3(i),x4(i),y4(i),z4(i),
372 . x1(i),y1(i),z1(i),x2(i),y2(i),z2(i),
373 . gap_m(im))
374 ENDIF
375 IF(ib4 == 1 .and. ib1 == 1)THEN
376 CALL i20cgap0(x40,y40,z40,x10,y10,z10,
377 . x20,y20,z20,x30,y30,z30,
378 . x4(i),y4(i),z4(i),x1(i),y1(i),z1(i),
379 . x2(i),y2(i),z2(i),x3(i),y3(i),z3(i),
380 . gap_m(im))
381 ENDIF
382 IF(ib1 == 1 .and. ib2+ib4 == 0)THEN
383 CALL i20cgap1(x10,y10,z10,x20,y20,z20,
384 . x40,y40,z40,
385 . x1(i),y1(i),z1(i),x2(i),y2(i),z2(i),
386 . x4(i),y4(i),z4(i),gap_m(im))
387 ENDIF
388 IF(ib2 == 1 .and. ib3+ib1 == 0)THEN
389 CALL i20cgap1(x20,y20,z20,x30,y30,z30,
390 . x10,y10,z10,
391 . x2(i),y2(i),z2(i),x3(i),y3(i),z3(i),
392 . x1(i),y1(i),z1(i),gap_m(im))
393 ENDIF
394 IF(ib3 == 1 .and. ib4+ib2 == 0)THEN
395 CALL i20cgap1(x30,y30,z30,x40,y40,z40,
396 . x10,y10,z10,
397 . x3(i),y3(i),z3(i),x4(i),y4(i),z4(i),
398 . x2(i),y2(i),z2(i),gap_m(im))
399 ENDIF
400 IF(ib4 == 1 .and. ib1+ib3 == 0)THEN
401 CALL i20cgap1(x40,y40,z40,x10,y10,z10,
402 . x30,y30,z30,
403 . x4(i),y4(i),z4(i),x1(i),y1(i),z1(i),
404 . x3(i),y3(i),z3(i),gap_m(im))
405 ENDIF
406 ENDIF
407 ENDDO
408 ENDIF
409
410 DO i=lft,llt
411 x0(i) = fourth*(x1(i)+x2(i)+x3(i)+x4(i))
412 y0(i) = fourth*(y1(i)+y2(i)+y3(i)+y4(i))
413 z0(i) = fourth*(z1(i)+z2(i)+z3(i)+z4(i))
414 ENDDO
415
416 DO i=lft,llt
417 IF(ix3(i)==ix4(i))THEN
418 x0(i) = x3(i)
419 y0(i) = y3(i)
420 z0(i) = z3(i)
421 ENDIF
422 ENDDO
423
424 DO i=lft,llt
425
426 x01(i) = x1(i) - x0(i)
427 y01(i) = y1(i) - y0(i)
428 z01(i) = z1(i) - z0(i)
429
430 x02(i) = x2(i) - x0(i)
431 y02(i) = y2(i) - y0(i)
432 z02(i) = z2(i) - z0(i)
433
434 x03(i) = x3(i) - x0(i)
435 y03(i) = y3(i) - y0(i)
436 z03(i) = z3(i) - z0(i)
437
438 x04(i) = x4(i) - x0(i)
439 y04(i) = y4(i) - y0(i)
440 z04(i) = z4(i) - z0(i)
441
442 xi0 = x0(i) - xi(i)
443 yi0 = y0(i) - yi(i)
444 zi0 = z0(i) - zi(i)
445
446 xi1(i) = x1(i) - xi(i)
447 yi1(i) = y1(i) - yi(i)
448 zi1(i) = z1(i) - zi(i)
449
450 xi2(i) = x2(i) - xi(i)
451 yi2(i) = y2(i) - yi(i)
452 zi2(i) = z2(i) - zi(i)
453
454 xi3(i) = x3(i) - xi(i)
455 yi3(i) = y3(i) - yi(i)
456 zi3(i) = z3(i) - zi(i)
457
458 xi4(i) = x4(i) - xi(i)
459 yi4(i) = y4(i) - yi(i)
460 zi4(i) = z4(i) - zi(i)
461
462 sx1 = yi0*zi1(i) - zi0*yi1(i)
463 sy1 = zi0*xi1(i) - xi0*zi1(i)
464 sz1 = xi0*yi1(i) - yi0*xi1(i)
465
466 sx2 = yi0*zi2(i) - zi0*yi2(i)
467 sy2 = zi0*xi2(i) - xi0*zi2(i)
468 sz2 = xi0*yi2(i) - yi0*xi2(i)
469
470 sx0 = y01(i)*z02(i) - z01(i)*y02(i)
471 sy0 = z01(i)*x02(i) - x01(i)*z02(i)
472 sz0 = x01(i)*y02(i) - y01(i)*x02(i)
473 s2 = one/
max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
474
475 lb1(i) = -(sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
476 lc1(i) = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
477
478 sx3 = yi0*zi3(i) - zi0*yi3(i)
479 sy3 = zi0*xi3(i) - xi0*zi3(i)
480 sz3 = xi0*yi3(i) - yi0*xi3(i)
481
482 sx0 = y02(i)*z03(i) - z02(i)*y03(i)
483 sy0 = z02(i)*x03(i) - x02(i)*z03(i)
484 sz0 = x02(i)*y03(i) - y02(i)*x03(i)
485 s2 = one/
max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
486
487 lb2(i) = -(sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
488 lc2(i) = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
489
490 sx4 = yi0*zi4(i) - zi0*yi4(i)
491 sy4 = zi0*xi4(i) - xi0*zi4(i)
492 sz4 = xi0*yi4(i) - yi0*xi4(i)
493
494 sx0 = y03(i)*z04(i) - z03(i)*y04(i)
495 sy0 = z03(i)*x04(i) - x03(i)*z04(i)
496 sz0 = x03(i)*y04(i) - y03(i)*x04(i)
497 s2 = one/
max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
498
499 lb3(i) = -(sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
500 lc3(i) = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
501
502 sx0 = y04(i)*z01(i) - z04(i)*y01(i)
503 sy0 = z04(i)*x01(i) - x04(i)*z01(i)
504 sz0 = x04(i)*y01(i) - y04(i)*x01(i)
505 s2 = one/
max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
506
507 lb4(i) = -(sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
508 lc4(i) = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
509
510 aaa = one/
max(em30,x01(i)*x01(i)+y01(i)*y01(i)+z01(i)*z01(i))
511 hlc1(i)= lc1(i)*abs(lc1(i))*aaa
512 hlb4(i)= lb4(i)*abs(lb4(i))*aaa
513 al1(i) = -(xi0*x01(i)+yi0*y01(i)+zi0*z01(i))*aaa
514 al1(i) =
max(zero,
min(one,al1(i)))
515 aaa = one/
max(em30,x02(i)*x02(i)+y02(i)*y02(i)+z02(i)*z02(i))
516 hlc2(i)= lc2(i)*abs(lc2(i))*aaa
517 hlb1(i)= lb1(i)*abs(lb1(i))*aaa
518 al2(i) = -(xi0*x02(i)+yi0*y02(i)+zi0*z02(i))*aaa
519 al2(i) =
max(zero,
min(one,al2(i)))
520 aaa = one/
max(em30,x03(i)*x03(i)+y03(i)*y03(i)+z03(i)*z03(i))
521 hlc3(i)= lc3(i)*abs(lc3(i))*aaa
522 hlb2(i)= lb2(i)*abs(lb2(i))*aaa
523 al3(i) = -(xi0*x03(i)+yi0*y03(i)+zi0*z03(i))*aaa
524 al3(i) =
max(zero,
min(one,al3(i)))
525 aaa = one/
max(em30,x04(i)*x04(i)+y04(i)*y04(i)+z04(i)*z04(i))
526 hlc4(i)= lc4(i)*abs(lc4(i))*aaa
527 hlb3(i)= lb3(i)*abs(lb3(i))*aaa
528 al4(i) = -(xi0*x04(i)+yi0*y04(i)+zi0*z04(i))*aaa
529 al4(i) =
max(zero,
min(one,al4(i)))
530
531 ENDDO
532
533 DO i=lft,llt
534 x12 = x2(i) - x1(i)
535 y12 = y2(i) - y1(i)
536 z12 = z2(i) - z1(i)
537 la = one - lb1(i) - lc1(i)
538
539 aaa = one /
max(em20,x12*x12+y12*y12+z12*z12)
540 hla= la*abs(la) * aaa
541 IF(la<zero.AND.
542 + hla<=hlb1(i).AND.hla<=hlc1(i))THEN
543 lb1(i) = (xi2(i)*x12+yi2(i)*y12+zi2(i)*z12) * aaa
544 lb1(i) =
max(zero,
min(one,lb1(i)))
545 lc1(i) = one - lb1(i)
546 ELSEIF(lb1(i)<zero.AND.
547 + hlb1(i)<=hlc1(i).AND.hlb1(i)<=hla)THEN
548 lb1(i) = zero
549 lc1(i) = al2(i)
550 ELSEIF(lc1(i)<zero.AND.
551 + hlc1(i)<=hla.AND.hlc1(i)<=hlb1(i))THEN
552 lc1(i) = zero
553 lb1(i) = al1(i)
554 ENDIF
555 ENDDO
556
557 DO i=lft,llt
558 x23 = x3(i) - x2(i)
559 y23 = y3(i) - y2(i)
560 z23 = z3(i) - z2(i)
561 la = one - lb2(i) - lc2(i)
562
563 aaa = one /
max(em20,x23*x23+y23*y23+z23*z23)
564 hla= la*abs(la) * aaa
565 IF(la<zero.AND.
566 + hla<=hlb2(i).AND.hla<=hlc2(i))THEN
567 lb2(i) = (xi3(i)*x23+yi3(i)*y23+zi3(i)*z23)*aaa
568 lb2(i) =
max(zero,
min(one,lb2(i)))
569 lc2(i) = one - lb2(i)
570 ELSEIF(lb2(i)<zero.AND.
571 + hlb2(i)<=hlc2(i).AND.hlb2(i)<=hla)THEN
572 lb2(i) = zero
573 lc2(i) = al3(i)
574 ELSEIF(lc2(i)<zero.AND.
575 + hlc2(i)<=hla.AND.hlc2(i)<=hlb2(i))THEN
576 lc2(i) = zero
577 lb2(i) = al2(i)
578 ENDIF
579 ENDDO
580
581 DO i=lft,llt
582 x34 = x4(i) - x3(i)
583 y34 = y4(i) - y3(i)
584 z34 = z4(i) - z3(i)
585 la = one - lb3(i) - lc3(i)
586
587 aaa = one /
max(em20,x34*x34+y34*y34+z34*z34)
588 hla= la*abs(la) * aaa
589 IF(la<zero.AND.
590 + hla<=hlb3(i).AND.hla<=hlc3(i))THEN
591 lb3(i) = (xi4(i)*x34+yi4(i)*y34+zi4(i)*z34)*aaa
592 lb3(i) =
max(zero,
min(one,lb3(i)))
593 lc3(i) = one - lb3(i)
594 ELSEIF(lb3(i)<zero.AND.
595 + hlb3(i)<=hlc3(i).AND.hlb3(i)<=hla)THEN
596 lb3(i) = zero
597 lc3(i) = al4(i)
598 ELSEIF(lc3(i)<zero.AND.
599 + hlc3(i)<=hla.AND.hlc3(i)<=hlb3(i))THEN
600 lc3(i) = zero
601 lb3(i) = al3(i)
602 ENDIF
603 ENDDO
604
605 DO i=lft,llt
606 x41 = x1(i) - x4(i)
607 y41 = y1(i) - y4(i)
608 z41 = z1(i) - z4(i)
609 la = one - lb4(i) - lc4(i)
610
611 aaa = one /
max(em20,x41*x41+y41*y41+z41*z41)
612 hla= la*abs(la) * aaa
613 IF(la<zero.AND.
614 + hla<=hlb4(i).AND.hla<=hlc4(i))THEN
615 lb4(i) = (xi1(i)*x41+yi1(i)*y41+zi1(i)*z41)*aaa
616 lb4(i) =
max(zero,
min(one,lb4(i)))
617 lc4(i) = one - lb4(i)
618 ELSEIF(lb4(i)<zero.AND.
619 + hlb4(i)<=hlc4(i).AND.hlb4(i)<=hla)THEN
620 lb4(i) = zero
621 lc4(i) = al1(i)
622 ELSEIF(lc4(i)<zero.AND.
623 + hlc4(i)<=hla.AND.hlc4(i)<=hlb4(i))THEN
624 lc4(i) = zero
625 lb4(i) = al4(i)
626 ENDIF
627 ENDDO
628
629 DO i=lft,llt
630
631 nx1(i) = xi(i)-(x0(i) + lb1(i)*x01(i) + lc1(i)*x02(i))
632 ny1(i) = yi(i)-(y0(i) + lb1(i)*y01(i) + lc1(i)*y02(i))
633 nz1(i) = zi(i)-(z0(i) + lb1(i)*z01(i) + lc1(i)*z02(i))
634 p1(i) = nx1(i)*nx1(i) + ny1(i)*ny1(i) +nz1(i)*nz1(i)
635
636 nx2(i) = xi(i)-(x0(i) + lb2(i)*x02(i) + lc2(i)*x03(i))
637 ny2(i) = yi(i)-(y0(i) + lb2(i)*y02(i) + lc2(i)*y03(i))
638 nz2(i) = zi(i)-(z0(i) + lb2(i)*z02(i) + lc2(i)*z03(i))
639 p2(i) = nx2(i)*nx2(i) + ny2(i)*ny2(i) +nz2(i)*nz2(i)
640
641 nx3(i) = xi(i)-(x0(i) + lb3(i)*x03(i) + lc3(i)*x04(i))
642 ny3(i) = yi(i)-(y0(i) + lb3(i)*y03(i) + lc3(i)*y04(i))
643 nz3(i) = zi(i)-(z0(i) + lb3(i)*z03(i) + lc3(i)*z04(i))
644 p3(i) = nx3(i)*nx3(i) + ny3(i)*ny3(i) +nz3(i)*nz3(i)
645
646 nx4(i) = xi(i)-(x0(i) + lb4(i)*x04(i) + lc4(i)*x01(i))
647 ny4(i) = yi(i)-(y0(i) + lb4(i)*y04(i) + lc4(i)*y01(i))
648 nz4(i) = zi(i)-(z0(i) + lb4(i)*z04(i) + lc4(i)*z01(i))
649 p4(i) = nx4(i)*nx4(i) + ny4(i)*ny4(i) +nz4(i)*nz4(i)
650
651 ENDDO
652
653 RETURN
integer function bitget(i, n)
integer function bitset(i, n)
integer function bitunset(i, n)
subroutine i20cgap1(x10, y10, z10, x20, y20, z20, x40, y40, z40, x1, y1, z1, x2, y2, z2, x4, y4, z4, gap_m)
subroutine i20cgap0(x10, y10, z10, x20, y20, z20, x30, y30, z30, x40, y40, z40, x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, gap_m)