OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i25dst3_22.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "task_c.inc"
#include "comlock.inc"
#include "com08_c.inc"
#include "vectorize.inc"
#include "lockon.inc"
#include "lockoff.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine i25dst3_22 (jlt, cand_n, cand_e, ishel, xx, yy, zz, xi, yi, zi, vx1, vx2, vx3, vx4, vxi, vy1, vy2, vy3, vy4, vyi, vz1, vz2, vz3, vz4, vzi, nin, nsn, ix1, ix2, ix3, ix4, nsvg, stif, inacti, mseglo, gaps, gapm, irect, irtlm, time_s, gap_nm, pene_old, stif_old, itab, penmin, eps0, icont_i, marge, nax, nay, naz, nbx, nby, nbz, far, pent, subtria, lbs, lcs, lbp, lcp, mvoisa, mvoisb, gapmxl, ibounda, iboundb, vtx_bisector, drad, dgapload)
subroutine interseca_25 (ixx3, ixx4, interseca, v12, v23, v34, v41, vo12, vo23, vo34, vo41, xx1, yy1, zz1, xx2, yy2, zz2, xx3, yy3, zz3, xx4, yy4, zz4, xx5, yy5, zz5, vxi, vyi, vzi, xo1, xo2, xo3, xo4, xo5, yo1, yo2, yo3, yo4, yo5, zo1, zo2, zo3, zo4, zo5, xi, yi, zi, xoi, yoi, zoi, zerom, unp, zeromt, unpt)
subroutine intersecb_25 (ixx3, ixx4, interseca, intersecb, v12, v23, v34, v41, vo12, vo23, vo34, vo41, xx1, yy1, zz1, xx2, yy2, zz2, xx3, yy3, zz3, xx4, yy4, zz4, xx5, yy5, zz5, vxi, vyi, vzi, xo1, xo2, xo3, xo4, xo5, yo1, yo2, yo3, yo4, yo5, zo1, zo2, zo3, zo4, zo5, xi, yi, zi, xoi, yoi, zoi, zerom, unp, zeromt, unpt)
subroutine intersecv_25 (ixx3, ixx4, subtria, intersect, v12, v23, v34, v41, vo12, vo23, vo34, vo41, xx1, yy1, zz1, xx2, yy2, zz2, xx3, yy3, zz3, xx4, yy4, zz4, xx5, yy5, zz5, vxi, vyi, vzi, xo1, xo2, xo3, xo4, xo5, yo1, yo2, yo3, yo4, yo5, zo1, zo2, zo3, zo4, zo5, xi, yi, zi, xoi, yoi, zoi, zerom, unp, zeromt, unpt)
subroutine intersecv0_25 (ixx3, ixx4, intersect, v12, v23, v34, v41, vo12, vo23, vo34, vo41, xx1, yy1, zz1, xx2, yy2, zz2, xx3, yy3, zz3, xx4, yy4, zz4, xx5, yy5, zz5, xoi, yoi, zoi, xo1, xo2, xo3, xo4, xo5, yo1, yo2, yo3, yo4, yo5, zo1, zo2, zo3, zo4, zo5, xi, yi, zi, zeromt, unpt)
subroutine i25glob_22 (jlt, cand_n_n, cand_e_n, cand_e, nin, nsn, ix1, ix2, ix3, ix4, nsvg, stif, inacti, mseglo, irtlm, time_s, itab, subtria, far, pent, lb, lc, index, farm, penm, lbm, lcm)

Function/Subroutine Documentation

◆ i25dst3_22()

subroutine i25dst3_22 ( integer jlt,
integer, dimension(*) cand_n,
integer, dimension(*) cand_e,
integer, dimension(mvsiz) ishel,
xx,
yy,
zz,
xi,
yi,
zi,
vx1,
vx2,
vx3,
vx4,
vxi,
vy1,
vy2,
vy3,
vy4,
vyi,
vz1,
vz2,
vz3,
vz4,
vzi,
integer nin,
integer nsn,
integer, dimension(mvsiz) ix1,
integer, dimension(mvsiz) ix2,
integer, dimension(mvsiz) ix3,
integer, dimension(mvsiz) ix4,
integer, dimension(mvsiz) nsvg,
stif,
integer inacti,
integer, dimension(*) mseglo,
gaps,
gapm,
integer, dimension(4,*) irect,
integer, dimension(4,nsn) irtlm,
time_s,
gap_nm,
pene_old,
stif_old,
integer, dimension(*) itab,
penmin,
eps0,
integer, dimension(*) icont_i,
marge,
nax,
nay,
naz,
nbx,
nby,
nbz,
integer, dimension(mvsiz) far,
pent,
integer, dimension(mvsiz) subtria,
lbs,
lcs,
lbp,
lcp,
integer, dimension(mvsiz,4) mvoisa,
integer, dimension(mvsiz,4) mvoisb,
gapmxl,
integer, dimension(4,mvsiz) ibounda,
integer, dimension(4,mvsiz) iboundb,
real*4, dimension(3,2,*) vtx_bisector,
intent(in) drad,
intent(in) dgapload )

Definition at line 34 of file i25dst3_22.F.

53C-----------------------------------------------
54C M o d u l e s
55C-----------------------------------------------
56 USE tri7box
57C-----------------------------------------------
58C I m p l i c i t T y p e s
59C-----------------------------------------------
60#include "implicit_f.inc"
61C-----------------------------------------------
62C G l o b a l P a r a m e t e r s
63C-----------------------------------------------
64#include "mvsiz_p.inc"
65#include "task_c.inc"
66#include "comlock.inc"
67C-----------------------------------------------
68C D u m m y A r g u m e n t s
69C-----------------------------------------------
70 INTEGER JLT, NIN, NSN, INACTI,
71 . CAND_N(*),CAND_E(*),NSVG(MVSIZ), ISHEL(MVSIZ)
72 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
73 . MSEGLO(*),IRTLM(4,NSN), SUBTRIA(MVSIZ),
74 . MVOISA(MVSIZ,4), MVOISB(MVSIZ,4),IBOUNDA(4,MVSIZ),IBOUNDB(4,MVSIZ)
75 my_real , INTENT(IN) :: dgapload ,drad
76 my_real
77 . time_s(2,nsn), gaps(*), gapm(*),
78 . pene_old(5,*),stif_old(2,nsn), penmin, eps0, marge
79 my_real
80 . xx(mvsiz,5),yy(mvsiz,5),zz(mvsiz,5),
81 . xi(mvsiz), yi(mvsiz), zi(mvsiz), stif(mvsiz),
82 . vx1(mvsiz),vx2(mvsiz),vx3(mvsiz),vx4(mvsiz),vy1(mvsiz),
83 . vy2(mvsiz),vy3(mvsiz),vy4(mvsiz),vz1(mvsiz),vz2(mvsiz),
84 . vz3(mvsiz),vz4(mvsiz),vxi(mvsiz),vyi(mvsiz),vzi(mvsiz),
85 . nax(mvsiz,5), nay(mvsiz,5), naz(mvsiz,5),
86 . nbx(mvsiz,5), nby(mvsiz,5), nbz(mvsiz,5),
87 . pent(mvsiz),
88 . lbs(mvsiz), lcs(mvsiz), lbp(mvsiz,4), lcp(mvsiz,4),
89 . gap_nm(4,mvsiz), gapmxl(mvsiz)
90 INTEGER IRECT(4,*),ITAB(*),ICONT_I(*), FAR(MVSIZ)
91 real*4 vtx_bisector(3,2,*)
92C-----------------------------------------------
93C C o m m o n B l o c k s
94C-----------------------------------------------
95#include "com08_c.inc"
96C-----------------------------------------------
97C L o c a l V a r i a b l e s
98C-----------------------------------------------
99 INTEGER I, L, N, IT, I1, I2, ITRIA(2,4),
100 . MGLOB, ITPERM(4), IA, IB, IB1, IB2, IB3, IBX, IX, IY, IZ
101 INTEGER INTERSECTA(MVSIZ), INTERSECTB(MVSIZ),
102 . INTERSECA, INTERSECB, IKEEP,
103 . FARA(MVSIZ,4), FARB(MVSIZ,4), INGAPA(MVSIZ,4), INGAPB(MVSIZ,4),
104 . SUBTRIB(MVSIZ), SUBTRIX(MVSIZ), ICONT_R(MVSIZ)
105 my_real
106 . xij(mvsiz,4), xi0v(mvsiz), xi5, xoi5,
107 . yij(mvsiz,4), yi0v(mvsiz), yi5, yoi5,
108 . zij(mvsiz,4), zi0v(mvsiz), zi5, zoi5,
109 . x51, x52, x53, x54,
110 . y51, y52, y53, y54,
111 . z51, z52, z53, z54,
112 . xo1, xo2, xo3, xo4, xo5, xoi,
113 . yo1, yo2, yo3, yo4, yo5, yoi,
114 . zo1, zo2, zo3, zo4, zo5, zoi,
115 . vo12, vo23, vo34, vo41, pene, penax, penbx
116 my_real
117 . gap_mm(mvsiz), gap, gap21,gap22,gap23,gap24,
118 . xp, yp, zp,
119 . dx, dy, dz, lx, ld, lax, lbx, lcx, dmin,
120 . dd(mvsiz,4),
121 . lap, hla, hlb(mvsiz,4), hlc(mvsiz,4),
122 . lap1,lap2,lap3,lap4,
123 . al(mvsiz,4)
124 my_real
125 . unp,zerom,eps,unpt,zeromt,aaa,
126 . pena(mvsiz,4), penb(mvsiz,4), bb(mvsiz,4),
127 . sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4,sz1,sz2,sz3,sz4,
128 . x0n(mvsiz,4), y0n(mvsiz,4), z0n(mvsiz,4),
129 . xn(mvsiz,4),yn(mvsiz,4),zn(mvsiz,4),
130 . lb(mvsiz,4), lc(mvsiz,4),
131 . sida(mvsiz,4), sidb(mvsiz,4),
132 . x12, y12, z12,
133 . px, py, pz, pp, p1, p2, xh, yh, zh, d1, d2, d3, vx, vy, vz,
134 . ll,nn,pn,
135 . v12,v23,v34,v41,
136 . nx1, nx2, nx3, nx4,
137 . ny1, ny2, ny3, ny4,
138 . nz1, nz2, nz3, nz4,
139 . sox125,sox235,sox345,sox415,
140 . soy125,soy235,soy345,soy415,
141 . soz125,soz235,soz345,soz415,
142 .
143 . eps02,tole,
144 . xx1,yy1 ,zz1,xx2,yy2,zz2,
145 . xx3 ,yy3,zz3,xx4,yy4,
146 . zz4 ,xx5,yy5,zz5
147 DATA itria/1,2,2,3,3,4,4,1/,
148 . itperm/1,4,3,2/
149 LOGICAL :: ANY_TRIANGLE
150C-----------------------------------------------------------------------
151C
152 unp = one + em4
153 zerom = zero - em4
154 eps = (two+half)/hundred
155 unpt = one + eps
156 zeromt = zero - eps
157 eps02=em3*em3
158C
159C initialisation (cf triangles)
160 fara(1:jlt,1:4) = 0
161 farb(1:jlt,1:4) = 0
162 pena(1:jlt,1:4) = zero
163 penb(1:jlt,1:4) = zero
164 dd(1:jlt,1:4) = zero
165C
166 any_triangle = .false.
167 DO i=1,jlt
168 IF(ix3(i)==ix4(i)) THEN
169 any_triangle = .true.
170 ENDIF
171 ENDDO
172
173 DO i=1,jlt
174C
175C For computing LBP, LCP, FAR=2, etc
176C
177C IF(STIF(I) == ZERO)CYCLE
178C
179 x0n(i,1) = xx(i,1) - xx(i,5)
180 y0n(i,1) = yy(i,1) - yy(i,5)
181 z0n(i,1) = zz(i,1) - zz(i,5)
182C
183 x0n(i,2) = xx(i,2) - xx(i,5)
184 y0n(i,2) = yy(i,2) - yy(i,5)
185 z0n(i,2) = zz(i,2) - zz(i,5)
186C
187 x0n(i,3) = xx(i,3) - xx(i,5)
188 y0n(i,3) = yy(i,3) - yy(i,5)
189 z0n(i,3) = zz(i,3) - zz(i,5)
190C
191 x0n(i,4) = xx(i,4) - xx(i,5)
192 y0n(i,4) = yy(i,4) - yy(i,5)
193 z0n(i,4) = zz(i,4) - zz(i,5)
194C
195 gap_mm(i)=fourth*(gap_nm(1,i)+gap_nm(2,i)+gap_nm(3,i)+gap_nm(4,i))
196 ENDDO
197
198 IF (any_triangle) THEN
199 DO i=1,jlt
200 IF(ix3(i)==ix4(i)) THEN
201 gap_mm(i)=gap_nm(3,i)
202 END IF
203 ENDDO
204 ENDIF
205
206C--------------------------------------------------------
207C
208C--------------------------------------------------------
209C
210C normals to triangles (recomputed here to avoid storing them)
211 DO i=1,jlt
212C
213C IF(STIF(I) == ZERO)CYCLE
214C
215 xn(i,1) = y0n(i,1)*z0n(i,2) - z0n(i,1)*y0n(i,2)
216 yn(i,1) = z0n(i,1)*x0n(i,2) - x0n(i,1)*z0n(i,2)
217 zn(i,1) = x0n(i,1)*y0n(i,2) - y0n(i,1)*x0n(i,2)
218
219 xn(i,2) = y0n(i,2)*z0n(i,3) - z0n(i,2)*y0n(i,3)
220 yn(i,2) = z0n(i,2)*x0n(i,3) - x0n(i,2)*z0n(i,3)
221 zn(i,2) = x0n(i,2)*y0n(i,3) - y0n(i,2)*x0n(i,3)
222
223 xn(i,3) = y0n(i,3)*z0n(i,4) - z0n(i,3)*y0n(i,4)
224 yn(i,3) = z0n(i,3)*x0n(i,4) - x0n(i,3)*z0n(i,4)
225 zn(i,3) = x0n(i,3)*y0n(i,4) - y0n(i,3)*x0n(i,4)
226
227 xn(i,4) = y0n(i,4)*z0n(i,1) - z0n(i,4)*y0n(i,1)
228 yn(i,4) = z0n(i,4)*x0n(i,1) - x0n(i,4)*z0n(i,1)
229 zn(i,4) = x0n(i,4)*y0n(i,1) - y0n(i,4)*x0n(i,1)
230
231 ENDDO
232C
233C--------------------------------------------------------
234C Search intersection
235C--------------------------------------------------------
236C
237 intersecta(1:jlt) = 0
238 intersectb(1:jlt) = 0
239
240 ingapa(1:jlt,1:4) = 0
241 ingapb(1:jlt,1:4) = 0
242
243 icont_r(1:jlt) = 0
244C
245 DO i=1,jlt
246C attention triangles
247
248 IF(ishel(i)/=0.AND.inacti==0)THEN
249 nx1 = xn(i,1)
250 ny1 = yn(i,1)
251 nz1 = zn(i,1)
252 aaa = gaps(i)+gap_nm(1,i)
253 aaa = aaa/sqrt(nx1*nx1+ny1*ny1+nz1*nz1)
254 nx1 = nx1 * aaa
255 ny1 = ny1 * aaa
256 nz1 = nz1 * aaa
257
258 nx2 = xn(i,2)
259 ny2 = yn(i,2)
260 nz2 = zn(i,2)
261 aaa = gaps(i)+gap_nm(2,i)
262 aaa = aaa/max(em30,sqrt(nx2*nx2+ny2*ny2+nz2*nz2))
263 nx2 = nx2 * aaa
264 ny2 = ny2 * aaa
265 nz2 = nz2 * aaa
266
267 IF(ix3(i)/=ix4(i))THEN
268 nx3 = xn(i,3)
269 ny3 = yn(i,3)
270 nz3 = zn(i,3)
271 aaa = gaps(i)+gap_nm(3,i)
272 aaa = aaa/max(em30,sqrt(nx3*nx3+ny3*ny3+nz3*nz3))
273 nx3 = nx3 * aaa
274 ny3 = ny3 * aaa
275 nz3 = nz3 * aaa
276
277 nx4 = xn(i,4)
278 ny4 = yn(i,4)
279 nz4 = zn(i,4)
280 aaa = gaps(i)+gap_nm(4,i)
281 aaa = aaa/max(em30,sqrt(nx4*nx4+ny4*ny4+nz4*nz4))
282 nx4 = nx4 * aaa
283 ny4 = ny4 * aaa
284 nz4 = nz4 * aaa
285
286 xx1 = xx(i,1) + nx1
287 yy1 = yy(i,1) + ny1
288 zz1 = zz(i,1) + nz1
289
290 xx2 = xx(i,2) + nx2
291 yy2 = yy(i,2) + ny2
292 zz2 = zz(i,2) + nz2
293
294 xx3 = xx(i,3) + nx3
295 yy3 = yy(i,3) + ny3
296 zz3 = zz(i,3) + nz3
297
298 xx4 = xx(i,4) + nx4
299 yy4 = yy(i,4) + ny4
300 zz4 = zz(i,4) + nz4
301
302 xx5 = fourth*(xx1+xx2+xx3+xx4)
303 yy5 = fourth*(yy1+yy2+yy3+yy4)
304 zz5 = fourth*(zz1+zz2+zz3+zz4)
305
306 ELSE
307 nx3 = xn(i,3)
308 ny3 = yn(i,3)
309 nz3 = zn(i,3)
310 aaa = gaps(i)+gap_nm(3,i)
311 aaa = aaa/max(em30,sqrt(nx3*nx3+ny3*ny3+nz3*nz3))
312 nx3 = nx3 * aaa
313 ny3 = ny3 * aaa
314 nz3 = nz3 * aaa
315
316 nx4 = xn(i,4)
317 ny4 = yn(i,4)
318 nz4 = zn(i,4)
319
320 xx1 = xx(i,1) + nx1
321 yy1 = yy(i,1) + ny1
322 zz1 = zz(i,1) + nz1
323
324 xx2 = xx(i,2) + nx2
325 yy2 = yy(i,2) + ny2
326 zz2 = zz(i,2) + nz2
327
328 xx3 = xx(i,3) + nx3
329 yy3 = yy(i,3) + ny3
330 zz3 = zz(i,3) + nz3
331
332
333
334 xx4 = xx3
335 yy4 = yy3
336 zz4 = zz3
337
338 xx5 = xx3
339 yy5 = yy3
340 zz5 = zz3
341 ENDIF
342 ELSE
343 xx1 = xx(i,1)
344 yy1 = yy(i,1)
345 zz1 = zz(i,1)
346
347 xx2 = xx(i,2)
348 yy2 = yy(i,2)
349 zz2 = zz(i,2)
350
351 xx3 = xx(i,3)
352 yy3 = yy(i,3)
353 zz3 = zz(i,3)
354
355 xx4 = xx(i,4)
356 yy4 = yy(i,4)
357 zz4 = zz(i,4)
358
359 xx5 = xx(i,5)
360 yy5 = yy(i,5)
361 zz5 = zz(i,5)
362 ENDIF
363
364 xi5 = xx5 - xi(i)
365 yi5 = yy5 - yi(i)
366 zi5 = zz5 - zi(i)
367
368 v12 = xn(i,1)*xi5 + yn(i,1)*yi5 + zn(i,1)*zi5
369 v23 = xn(i,2)*xi5 + yn(i,2)*yi5 + zn(i,2)*zi5
370 v34 = xn(i,3)*xi5 + yn(i,3)*yi5 + zn(i,3)*zi5
371 v41 = xn(i,4)*xi5 + yn(i,4)*yi5 + zn(i,4)*zi5
372
373 IF(ishel(i)==0)THEN
374 IF(v12 < zero .and. v23 < zero .and.
375 . v34 < zero .and. v41 < zero ) THEN
376C INTERSECTION IMPOSSIBLE
377 cycle
378 ENDIF
379 END IF
380C POSSIBLE INTERSECTION
381
382C--------------------------------------------------------
383C OLD COORDINATES
384C--------------------------------------------------------
385 xo1 = xx1 - dt1*vx1(i)
386 yo1 = yy1 - dt1*vy1(i)
387 zo1 = zz1 - dt1*vz1(i)
388C
389 xo2 = xx2 - dt1*vx2(i)
390 yo2 = yy2 - dt1*vy2(i)
391 zo2 = zz2 - dt1*vz2(i)
392c
393 xo3 = xx3 - dt1*vx3(i)
394 yo3 = yy3 - dt1*vy3(i)
395 zo3 = zz3 - dt1*vz3(i)
396C
397 xo4 = xx4 - dt1*vx4(i)
398 yo4 = yy4 - dt1*vy4(i)
399 zo4 = zz4 - dt1*vz4(i)
400
401 xoi = xi(i) - dt1*vxi(i)
402 yoi = yi(i) - dt1*vyi(i)
403 zoi = zi(i) - dt1*vzi(i)
404
405 IF(ix3(i) /= ix4(i))THEN
406 xo5 = fourth*(xo1+xo2+xo3+xo4)
407 yo5 = fourth*(yo1+yo2+yo3+yo4)
408 zo5 = fourth*(zo1+zo2+zo3+zo4)
409 ELSE
410 xo5 = xo3
411 yo5 = yo3
412 zo5 = zo3
413 ENDIF
414
415c compute volumes at previous time step
416
417 x51 = xo1 - xo5
418 y51 = yo1 - yo5
419 z51 = zo1 - zo5
420
421 x52 = xo2 - xo5
422 y52 = yo2 - yo5
423 z52 = zo2 - zo5
424
425 x53 = xo3 - xo5
426 y53 = yo3 - yo5
427 z53 = zo3 - zo5
428
429 x54 = xo4 - xo5
430 y54 = yo4 - yo5
431 z54 = zo4 - zo5
432
433 xoi5 = xo5 - xoi
434 yoi5 = yo5 - yoi
435 zoi5 = zo5 - zoi
436
437 sox125 = y51*z52 - z51*y52
438 soy125 = z51*x52 - x51*z52
439 soz125 = x51*y52 - y51*x52
440 vo12 = sox125*xoi5 + soy125*yoi5 + soz125*zoi5
441
442 sox235 = y52*z53 - z52*y53
443 soy235 = z52*x53 - x52*z53
444 soz235 = x52*y53 - y52*x53
445 vo23 = sox235*xoi5 + soy235*yoi5 + soz235*zoi5
446
447 sox345 = y53*z54 - z53*y54
448 soy345 = z53*x54 - x53*z54
449 soz345 = x53*y54 - y53*x54
450 vo34 = sox345*xoi5 + soy345*yoi5 + soz345*zoi5
451
452 sox415 = y54*z51 - z54*y51
453 soy415 = z54*x51 - x54*z51
454 soz415 = x54*y51 - y54*x51
455 vo41 = sox415*xoi5 + soy415*yoi5 + soz415*zoi5
456
457c compute intersection time (linear approximation)
458c compute coordinates at intersection time
459c (intersection time can be different for each sub-triangle)
460 IF(ishel(i)==0)THEN
461 CALL interseca_25(
462 1 ix3(i) ,ix4(i),interseca,
463 1 v12 ,v23 ,v34 ,v41 ,
464 2 vo12 ,vo23 ,vo34 ,vo41 ,xx(i,1),
465 3 yy(i,1) ,zz(i,1),xx(i,2),yy(i,2),zz(i,2),
466 4 xx(i,3) ,yy(i,3),zz(i,3),xx(i,4),yy(i,4),
467 5 zz(i,4) ,xx(i,5),yy(i,5),zz(i,5),
468 6 vxi(i) ,vyi(i) ,vzi(i) ,xo1 ,xo2 ,
469 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
470 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
471 9 zo3 ,zo4 ,zo5 ,xi(i) ,yi(i) ,
472 a zi(i) ,xoi ,yoi ,zoi ,zerom ,
473 b unp ,zeromt,unpt )
474
475 ikeep = 1
476 n = cand_n(i)
477 IF(n <= nsn) THEN
478 IF(icont_i(n) < 0) ikeep = 0
479 ELSE
480 IF(icont_i_fi(nin)%P(n-nsn) < 0) ikeep = 0
481 ENDIF
482 IF(ikeep == 1.AND.interseca == 0)THEN ! IF contact in the same side : intersection is missed and have to be checked
483 CALL intersecv0_25(
484 1 ix3(i) ,ix4(i) , interseca,
485 1 v12 ,v23 ,v34 ,v41 ,
486 2 vo12 ,vo23 ,vo34 ,vo41 ,xx(i,1),
487 3 yy(i,1) ,zz(i,1),xx(i,2),yy(i,2),zz(i,2),
488 4 xx(i,3) ,yy(i,3),zz(i,3),xx(i,4),yy(i,4),
489 5 zz(i,4) ,xx(i,5),yy(i,5),zz(i,5),
490 6 xoi ,yoi ,zoi ,xo1 ,xo2 ,
491 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
492 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
493 9 zo3 ,zo4 ,zo5 ,xi(i) ,yi(i) ,
494 a zi(i) ,zerom , unp )
495 icont_r(i) = interseca
496 ENDIF
497
498 intersecta(i)=interseca
499 ELSE
500 CALL intersecb_25(
501 1 ix3(i) ,ix4(i),interseca,intersecb,
502 1 v12 ,v23 ,v34 ,v41 ,
503 2 vo12 ,vo23 ,vo34 ,vo41 ,xx1,
504 3 yy1 ,zz1,xx2,yy2,zz2,
505 4 xx3 ,yy3,zz3,xx4,yy4,
506 5 zz4 ,xx5,yy5,zz5,
507 6 vxi(i) ,vyi(i) ,vzi(i) ,xo1 ,xo2 ,
508 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
509 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
510 9 zo3 ,zo4 ,zo5 ,xi(i) ,yi(i) ,
511 a zi(i) ,xoi ,yoi ,zoi ,zerom ,
512 b unp ,zeromt,unpt )
513 intersecta(i)=interseca
514 intersectb(i)=intersecb
515 END IF
516C
517 ENDDO
518C--------------------------------------------------------
519C#include "vectorize.inc"
520 DO i=1,jlt
521C
522C IF(STIF(I) == ZERO)CYCLE
523C
524 xi0v(i) = xx(i,5) - xi(i)
525 yi0v(i) = yy(i,5) - yi(i)
526 zi0v(i) = zz(i,5) - zi(i)
527C
528 xij(i,1) = xx(i,1) - xi(i)
529 yij(i,1) = yy(i,1) - yi(i)
530 zij(i,1) = zz(i,1) - zi(i)
531C
532 xij(i,2) = xx(i,2) - xi(i)
533 yij(i,2) = yy(i,2) - yi(i)
534 zij(i,2) = zz(i,2) - zi(i)
535C
536 xij(i,3) = xx(i,3) - xi(i)
537 yij(i,3) = yy(i,3) - yi(i)
538 zij(i,3) = zz(i,3) - zi(i)
539C
540 xij(i,4) = xx(i,4) - xi(i)
541 yij(i,4) = yy(i,4) - yi(i)
542 zij(i,4) = zz(i,4) - zi(i)
543C
544 sx1 = yi0v(i)*zij(i,1) - zi0v(i)*yij(i,1)
545 sy1 = zi0v(i)*xij(i,1) - xi0v(i)*zij(i,1)
546 sz1 = xi0v(i)*yij(i,1) - yi0v(i)*xij(i,1)
547C
548 sx2 = yi0v(i)*zij(i,2) - zi0v(i)*yij(i,2)
549 sy2 = zi0v(i)*xij(i,2) - xi0v(i)*zij(i,2)
550 sz2 = xi0v(i)*yij(i,2) - yi0v(i)*xij(i,2)
551C
552 sx3 = yi0v(i)*zij(i,3) - zi0v(i)*yij(i,3)
553 sy3 = zi0v(i)*xij(i,3) - xi0v(i)*zij(i,3)
554 sz3 = xi0v(i)*yij(i,3) - yi0v(i)*xij(i,3)
555C
556 sx4 = yi0v(i)*zij(i,4) - zi0v(i)*yij(i,4)
557 sy4 = zi0v(i)*xij(i,4) - xi0v(i)*zij(i,4)
558 sz4 = xi0v(i)*yij(i,4) - yi0v(i)*xij(i,4)
559C
560 nn = one/
561 . max(em30,sqrt(xn(i,1)*xn(i,1)+ yn(i,1)*yn(i,1)+ zn(i,1)*zn(i,1)))
562 xn(i,1)=xn(i,1)*nn
563 yn(i,1)=yn(i,1)*nn
564 zn(i,1)=zn(i,1)*nn
565 lb(i,1) = -(xn(i,1)*sx2 + yn(i,1)*sy2 + zn(i,1)*sz2)*nn
566 lc(i,1) = (xn(i,1)*sx1 + yn(i,1)*sy1 + zn(i,1)*sz1)*nn
567C
568 nn = one/
569 . max(em30,sqrt(xn(i,2)*xn(i,2)+ yn(i,2)*yn(i,2)+ zn(i,2)*zn(i,2)))
570 xn(i,2)=xn(i,2)*nn
571 yn(i,2)=yn(i,2)*nn
572 zn(i,2)=zn(i,2)*nn
573 lb(i,2) = -(xn(i,2)*sx3 + yn(i,2)*sy3 + zn(i,2)*sz3)*nn
574 lc(i,2) = (xn(i,2)*sx2 + yn(i,2)*sy2 + zn(i,2)*sz2)*nn
575C
576 nn = one/
577 . max(em30,sqrt(xn(i,3)*xn(i,3)+ yn(i,3)*yn(i,3)+ zn(i,3)*zn(i,3)))
578 xn(i,3)=xn(i,3)*nn
579 yn(i,3)=yn(i,3)*nn
580 zn(i,3)=zn(i,3)*nn
581 lb(i,3) = -(xn(i,3)*sx4 + yn(i,3)*sy4 + zn(i,3)*sz4)*nn
582 lc(i,3) = (xn(i,3)*sx3 + yn(i,3)*sy3 + zn(i,3)*sz3)*nn
583C
584 nn = one/
585 . max(em30,sqrt(xn(i,4)*xn(i,4)+ yn(i,4)*yn(i,4)+ zn(i,4)*zn(i,4)))
586 xn(i,4)=xn(i,4)*nn
587 yn(i,4)=yn(i,4)*nn
588 zn(i,4)=zn(i,4)*nn
589 lb(i,4) = -(xn(i,4)*sx1 + yn(i,4)*sy1 + zn(i,4)*sz1)*nn
590 lc(i,4) = (xn(i,4)*sx4 + yn(i,4)*sy4 + zn(i,4)*sz4)*nn
591C
592 END DO
593C--------------------------------------------------------
594C
595C#include "vectorize.inc"
596 DO i=1,jlt
597C
598C IF(STIF(I) == ZERO)CYCLE
599C
600 aaa = one/max(em30,x0n(i,1)*x0n(i,1)+y0n(i,1)*y0n(i,1)+z0n(i,1)*z0n(i,1))
601 hlc(i,1) = lc(i,1)*abs(lc(i,1))*aaa
602 hlb(i,4) = lb(i,4)*abs(lb(i,4))*aaa
603 al(i,1) = -(xi0v(i)*x0n(i,1)+yi0v(i)*y0n(i,1)+zi0v(i)*z0n(i,1))*aaa
604 al(i,1) = max(zero,min(one,al(i,1)))
605 aaa = one/max(em30,x0n(i,2)*x0n(i,2)+y0n(i,2)*y0n(i,2)+z0n(i,2)*z0n(i,2))
606 hlc(i,2) = lc(i,2)*abs(lc(i,2))*aaa
607 hlb(i,1) = lb(i,1)*abs(lb(i,1))*aaa
608 al(i,2) = -(xi0v(i)*x0n(i,2)+yi0v(i)*y0n(i,2)+zi0v(i)*z0n(i,2))*aaa
609 al(i,2) = max(zero,min(one,al(i,2)))
610 aaa = one/max(em30,x0n(i,3)*x0n(i,3)+y0n(i,3)*y0n(i,3)+z0n(i,3)*z0n(i,3))
611 hlc(i,3) = lc(i,3)*abs(lc(i,3))*aaa
612 hlb(i,2) = lb(i,2)*abs(lb(i,2))*aaa
613 al(i,3) = -(xi0v(i)*x0n(i,3)+yi0v(i)*y0n(i,3)+zi0v(i)*z0n(i,3))*aaa
614 al(i,3) = max(zero,min(one,al(i,3)))
615 aaa = one/max(em30,x0n(i,4)*x0n(i,4)+y0n(i,4)*y0n(i,4)+z0n(i,4)*z0n(i,4))
616 hlc(i,4) = lc(i,4)*abs(lc(i,4))*aaa
617 hlb(i,3) = lb(i,3)*abs(lb(i,3))*aaa
618 al(i,4) = -(xi0v(i)*x0n(i,4)+yi0v(i)*y0n(i,4)+zi0v(i)*z0n(i,4))*aaa
619 al(i,4) = max(zero,min(one,al(i,4)))
620C
621 END DO
622C--------------------------------------------------------
623 DO i=1,jlt
624C Projection outside the triangle
625 IF(lb(i,1) < -em03 .OR. lc(i,1) < -em03 .OR.
626 . lb(i,1)+lc(i,1) > one+em03) THEN
627 fara(i,1)=1
628 farb(i,1)=1
629 END IF
630 x12 = xx(i,2) - xx(i,1)
631 y12 = yy(i,2) - yy(i,1)
632 z12 = zz(i,2) - zz(i,1)
633 lbp(i,1) = lb(i,1)
634 lcp(i,1) = lc(i,1)
635 lap = one-lbp(i,1)-lcp(i,1)
636 aaa = one / max(em20,x12*x12+y12*y12+z12*z12)
637 hla= lap*abs(lap) * aaa
638 IF(lap<zero.AND.
639 + hla<=hlb(i,1).AND.hla<=hlc(i,1))THEN
640 lbp(i,1) = (xij(i,2)*x12+yij(i,2)*y12+zij(i,2)*z12) * aaa
641 lbp(i,1) = max(zero,min(one,lbp(i,1)))
642 lcp(i,1) = one - lbp(i,1)
643 ELSEIF(lbp(i,1)<zero.AND.
644 + hlb(i,1)<=hlc(i,1).AND.hlb(i,1)<=hla)THEN
645 lbp(i,1) = zero
646 lcp(i,1) = al(i,2)
647 ELSEIF(lcp(i,1)<zero.AND.
648 + hlc(i,1)<=hla.AND.hlc(i,1)<=hlb(i,1))THEN
649 lcp(i,1) = zero
650 lbp(i,1) = al(i,1)
651 ENDIF
652C Projection outside the triangle
653 IF(lb(i,2) < -em03 .OR. lc(i,2) < -em03 .OR.
654 . lb(i,2)+lc(i,2) > one+em03) THEN
655 fara(i,2)=1
656 farb(i,2)=1
657 ENDIF
658 x12 = xx(i,3) - xx(i,2)
659 y12 = yy(i,3) - yy(i,2)
660 z12 = zz(i,3) - zz(i,2)
661 lbp(i,2) = lb(i,2)
662 lcp(i,2) = lc(i,2)
663 lap = one-lbp(i,2)-lcp(i,2)
664 aaa = one / max(em20,x12*x12+y12*y12+z12*z12)
665 hla= lap*abs(lap) * aaa
666 IF(lap<zero.AND.
667 + hla<=hlb(i,2).AND.hla<=hlc(i,2))THEN
668 lbp(i,2) = (xij(i,3)*x12+yij(i,3)*y12+zij(i,3)*z12) * aaa
669 lbp(i,2) = max(zero,min(one,lbp(i,2)))
670 lcp(i,2) = one - lbp(i,2)
671 ELSEIF(lbp(i,2)<zero.AND.
672 + hlb(i,2)<=hlc(i,2).AND.hlb(i,2)<=hla)THEN
673 lbp(i,2) = zero
674 lcp(i,2) = al(i,3)
675 ELSEIF(lcp(i,2)<zero.AND.
676 + hlc(i,2)<=hla.AND.hlc(i,2)<=hlb(i,2))THEN
677 lcp(i,2) = zero
678 lbp(i,2) = al(i,2)
679 END IF
680C Projection outside the triangle
681 IF(lb(i,3) < -em03 .OR. lc(i,3) < -em03 .OR.
682 . lb(i,3)+lc(i,3) > one+em03) THEN
683 fara(i,3)=1
684 farb(i,3)=1
685 END IF
686 x12 = xx(i,4) - xx(i,3)
687 y12 = yy(i,4) - yy(i,3)
688 z12 = zz(i,4) - zz(i,3)
689 lbp(i,3) = lb(i,3)
690 lcp(i,3) = lc(i,3)
691 lap = one-lbp(i,3)-lcp(i,3)
692 aaa = one / max(em20,x12*x12+y12*y12+z12*z12)
693 hla= lap*abs(lap) * aaa
694 IF(lap<zero.AND.
695 + hla<=hlb(i,3).AND.hla<=hlc(i,3))THEN
696 lbp(i,3) = (xij(i,4)*x12+yij(i,4)*y12+zij(i,4)*z12) * aaa
697 lbp(i,3) = max(zero,min(one,lbp(i,3)))
698 lcp(i,3) = one - lbp(i,3)
699 ELSEIF(lbp(i,3)<zero.AND.
700 + hlb(i,3)<=hlc(i,3).AND.hlb(i,3)<=hla)THEN
701 lbp(i,3) = zero
702 lcp(i,3) = al(i,4)
703 ELSEIF(lcp(i,3)<zero.AND.
704 + hlc(i,3)<=hla.AND.hlc(i,3)<=hlb(i,3))THEN
705 lcp(i,3) = zero
706 lbp(i,3) = al(i,3)
707 ENDIF
708C Projection outside the triangle
709 IF(lb(i,4) < -em03 .OR. lc(i,4) < -em03 .OR.
710 . lb(i,4)+lc(i,4) > one+em03) THEN
711 fara(i,4)=1
712 farb(i,4)=1
713 END IF
714 x12 = xx(i,1) - xx(i,4)
715 y12 = yy(i,1) - yy(i,4)
716 z12 = zz(i,1) - zz(i,4)
717 lbp(i,4) = lb(i,4)
718 lcp(i,4) = lc(i,4)
719 lap = one-lbp(i,4)-lcp(i,4)
720 aaa = one / max(em20,x12*x12+y12*y12+z12*z12)
721 hla= lap*abs(lap) * aaa
722 IF(lap<zero.AND.
723 + hla<=hlb(i,4).AND.hla<=hlc(i,4))THEN
724 lbp(i,4) = (xij(i,1)*x12+yij(i,1)*y12+zij(i,1)*z12) * aaa
725 lbp(i,4) = max(zero,min(one,lbp(i,4)))
726 lcp(i,4) = one - lbp(i,4)
727 ELSEIF(lbp(i,4)<zero.AND.
728 + hlb(i,4)<=hlc(i,4).AND.hlb(i,4)<=hla)THEN
729 lbp(i,4) = zero
730 lcp(i,4) = al(i,1)
731 ELSEIF(lcp(i,4)<zero.AND.
732 + hlc(i,4)<=hla.AND.hlc(i,4)<=hlb(i,4))THEN
733 lcp(i,4) = zero
734 lbp(i,4) = al(i,4)
735 ENDIF
736 ENDDO
737
738
739 DO i = 1,jlt
740
741 ikeep = 1
742 IF(ishel(i) > 0)THEN
743 n = cand_n(i)
744 IF((intersecta(i) == 0.AND.intersectb(i)==0) .AND.n <= nsn)THEN
745 IF (icont_i(n) < 0 ) ikeep = 0
746 ELSEIF(intersecta(i) == 0.AND.intersectb(i)==0) THEN
747 IF (icont_i_fi(nin)%P(n-nsn) < 0 ) ikeep = 0
748 ENDIF
749 ENDIF
750
751 lap1 = one-lbp(i,1)-lcp(i,1)
752 xp =lap1*xx(i,5)+lbp(i,1)*xx(i,1) + lcp(i,1)*xx(i,2)
753 yp =lap1*yy(i,5)+lbp(i,1)*yy(i,1) + lcp(i,1)*yy(i,2)
754 zp =lap1*zz(i,5)+lbp(i,1)*zz(i,1) + lcp(i,1)*zz(i,2)
755 dx =xi(i)-xp
756 dy =yi(i)-yp
757 dz =zi(i)-zp
758 dd(i,1) =dx*dx+dy*dy+dz*dz
759C
760 lap2 = one-lbp(i,2)-lcp(i,2)
761 xp =lap2*xx(i,5)+lbp(i,2)*xx(i,2) + lcp(i,2)*xx(i,3)
762 yp =lap2*yy(i,5)+lbp(i,2)*yy(i,2) + lcp(i,2)*yy(i,3)
763 zp =lap2*zz(i,5)+lbp(i,2)*zz(i,2) + lcp(i,2)*zz(i,3)
764 dx =xi(i)-xp
765 dy =yi(i)-yp
766 dz =zi(i)-zp
767 dd(i,2) =dx*dx+dy*dy+dz*dz
768C
769 lap3 = one-lbp(i,3)-lcp(i,3)
770 xp =lap3*xx(i,5)+lbp(i,3)*xx(i,3) + lcp(i,3)*xx(i,4)
771 yp =lap3*yy(i,5)+lbp(i,3)*yy(i,3) + lcp(i,3)*yy(i,4)
772 zp =lap3*zz(i,5)+lbp(i,3)*zz(i,3) + lcp(i,3)*zz(i,4)
773 dx =xi(i)-xp
774 dy =yi(i)-yp
775 dz =zi(i)-zp
776 dd(i,3) =dx*dx+dy*dy+dz*dz
777C
778 lap4 = one-lbp(i,4)-lcp(i,4)
779 xp =lap4*xx(i,5)+lbp(i,4)*xx(i,4) + lcp(i,4)*xx(i,1)
780 yp =lap4*yy(i,5)+lbp(i,4)*yy(i,4) + lcp(i,4)*yy(i,1)
781 zp =lap4*zz(i,5)+lbp(i,4)*zz(i,4) + lcp(i,4)*zz(i,1)
782 dx =xi(i)-xp
783 dy =yi(i)-yp
784 dz =zi(i)-zp
785 dd(i,4) =dx*dx+dy*dy+dz*dz
786
787C
788 gap = min(max(drad,gaps(i)+lap1*gap_mm(i)+lbp(i,1)*gap_nm(1,i)+lcp(i,1)*gap_nm(2,i)+dgapload),
789 . max(drad,gapmxl(i)+dgapload))
790 gap21 = gap**2
791 bb(i,1) =((xx(i,5)-xi(i))*xn(i,1)+(yy(i,5)-yi(i))*yn(i,1)+(zz(i,5)-zi(i))*zn(i,1))
792C
793 gap = min(max(drad,gaps(i)+lap2*gap_mm(i)+lbp(i,2)*gap_nm(2,i)+lcp(i,2)*gap_nm(3,i)+dgapload),
794 . max(drad,gapmxl(i)+dgapload))
795 gap22 =gap**2
796 bb(i,2) =((xx(i,5)-xi(i))*xn(i,2)+(yy(i,5)-yi(i))*yn(i,2)+(zz(i,5)-zi(i))*zn(i,2))
797C
798 gap = min(max(drad,gaps(i)+lap3*gap_mm(i)+lbp(i,3)*gap_nm(3,i)+lcp(i,3)*gap_nm(4,i)+dgapload),
799 . max(drad,gapmxl(i)+dgapload))
800 gap23 =gap**2
801 bb(i,3) =((xx(i,5)-xi(i))*xn(i,3)+(yy(i,5)-yi(i))*yn(i,3)+(zz(i,5)-zi(i))*zn(i,3))
802C
803 gap = min(max(drad,gaps(i)+lap4*gap_mm(i)+lbp(i,4)*gap_nm(4,i)+lcp(i,4)*gap_nm(1,i)+dgapload),
804 . max(drad,gapmxl(i)+dgapload))
805 gap24 =gap**2
806 bb(i,4) =((xx(i,5)-xi(i))*xn(i,4)+(yy(i,5)-yi(i))*yn(i,4)+(zz(i,5)-zi(i))*zn(i,4))
807C
808 IF(dd(i,1) <= gap21) THEN
809 IF(bb(i,1) <= zero .AND. intersecta(i) /= 1.AND.ikeep==1)THEN
810 ingapa(i,1)=1
811 END IF
812 IF(bb(i,1) >= zero .AND. intersectb(i) /= 1.AND.ikeep==1)THEN
813 ingapb(i,1)=1
814 END IF
815 ENDIF
816 IF(dd(i,2) <= gap22) THEN
817 IF(bb(i,2) <= zero .AND. intersecta(i) /= 1.AND.ikeep==1)THEN
818 ingapa(i,2)=1
819 END IF
820 IF(bb(i,2) >= zero .AND. intersectb(i) /= 1.AND.ikeep==1)THEN
821 ingapb(i,2)=1
822 END IF
823 ENDIF
824 IF(dd(i,3) <= gap23) THEN
825 IF(bb(i,3) <= zero .AND. intersecta(i) /= 1.AND.ikeep==1)THEN
826 ingapa(i,3)=1
827 END IF
828 IF(bb(i,3) >= zero .AND. intersectb(i) /= 1.AND.ikeep==1)THEN
829 ingapb(i,3)=1
830 END IF
831 ENDIF
832 IF(dd(i,4) <= gap24) THEN
833 IF(bb(i,4) <= zero .AND. intersecta(i) /= 1.AND.ikeep==1)THEN
834 ingapa(i,4)=1
835 END IF
836 IF(bb(i,4) >= zero .AND. intersectb(i) /= 1.AND.ikeep==1)THEN
837 ingapb(i,4)=1
838 END IF
839 END IF
840 END DO
841C--------------------------------------------------------
842C Determination of the TQ sub-triangle its cone contains the ND does not require Far!
843C--------------------------------------------------------
844 subtria(1:jlt)=0
845 subtrib(1:jlt)=0
846 subtrix(1:jlt)=0
847 DO i=1,jlt
848C
849 IF(stif(i) <= zero)cycle
850C
851 IF(ix3(i)/=ix4(i))THEN
852
853 dmin=min(dd(i,1),dd(i,2),dd(i,3),dd(i,4))
854 ld=ep20
855 DO it=1,4
856 IF(dd(i,it) <= onep03*dmin)THEN
857 lbx = lb(i,it)
858 lcx = lc(i,it)
859 lax = one-lb(i,it)-lc(i,it)
860C
861C Privilegier sector in which we are.
862 IF(lbx >= zero .AND. lcx >= zero)THEN
863 lx=zero
864 ELSE
865 ! closest point in the angular sector == center
866 lx = max(zero,dd(i,it)-bb(i,it)*bb(i,it))
867 END IF
868
869 IF(lx < ld) THEN
870 subtrix(i)= it
871 ld = lx
872 END IF
873 END IF
874 END DO
875C
876C !
877 it=subtrix(i)
878 IF(it > 0) THEN
879 IF(intersecta(i)/=0.OR.ingapa(i,it)/=0)subtria(i)=it
880 IF (ishel(i) /= 0)THEN
881 IF(intersectb(i)/=0.OR.ingapb(i,it)/=0)subtrib(i)=it
882 END IF
883 ENDIF
884 ELSE
885
886 IF(intersecta(i)/=0.OR.ingapa(i,1)/=0)THEN
887 subtria(i)= 1
888 END IF
889C
890 IF (ishel(i) /= 0)THEN
891 IF(intersectb(i)/=0.OR.ingapb(i,1)/=0)THEN
892 subtrib(i)= 1
893 END IF
894 END IF
895 END IF
896C
897 END DO
898C--------------------------------------------------------
899C calculation of penetration
900C--------------------------------------------------------
901#include "vectorize.inc"
902 DO i =1,jlt
903C
904C IF (STIF(I) == ZERO)CYCLE
905C
906 it=subtria(i)
907 IF(it==0)cycle
908C
909 i1=itria(1,it)
910 i2=itria(2,it)
911C
912 lap = one-lbp(i,it)-lcp(i,it)
913 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,
914 . max(drad,gapmxl(i)+dgapload))
915C
916 IF(bb(i,it) > zero)THEN
917C
918C penetration approximee (cf zone limite interpolation des normales)
919 pena(i,it)=max(zero,gap+bb(i,it))
920 ELSE ! IF(BB(I,IT) > ZERO)THEN
921 pena(i,it)=max(zero,gap-sqrt(dd(i,it)))
922 END IF
923C
924C--------exclude high pene in case of ICONT_R
925 IF(icont_r(i) >0)THEN
926 aaa = max(em30,x0n(i,it)*x0n(i,1)+y0n(i,it)*y0n(i,it)+z0n(i,it)*z0n(i,it))
927 tole =eps02*aaa
928 IF (gap >zero) tole =min(tole,gap*gap)
929 IF(pena(i,it)*pena(i,it) > tole) THEN
930 pena(i,it) = zero
931 END IF
932 END if!(ICONT_R(I) >0)THEN
933C
934C
935 ENDDO
936C--------------------------------------------------------
937C Check vs bissectors ...
938C--------------------------------------------------------
939#include "vectorize.inc"
940 DO i =1,jlt
941C
942C IF (STIF(I) == ZERO)CYCLE
943C
944 it=subtria(i)
945 IF(it==0)cycle
946C
947 IF(pena(i,it)==zero) cycle
948C
949 i1=itria(1,it)
950 i2=itria(2,it)
951C
952 xh=xi(i)+bb(i,it)*xn(i,it)
953 yh=yi(i)+bb(i,it)*yn(i,it)
954 zh=zi(i)+bb(i,it)*zn(i,it)
955C
956 IF(ix3(i) /= ix4(i))THEN
957C
958 ib1=ibounda(i1,i)
959 ib2=ibounda(i2,i)
960 IF(mvoisa(i,it)==0)THEN
961C
962 IF( (xh-xx(i,i1))* nax(i,it)+
963 . (yh-yy(i,i1))* nay(i,it)+
964 . (zh-zz(i,i1))* naz(i,it) >= gaps(i)) fara(i,it) =3
965C
966 ELSEIF((ib1 /= 0 .AND. ib2 == 0).OR.
967 . (ib2 /= 0 .AND. ib1 == 0))THEN
968C
969 ibx=max(ib1,ib2)
970C
971 ibx=max(ib1,ib2)
972 IF(ib1/=0)THEN
973 ix =i1
974 ELSEIF(ib2/=0)THEN
975 ix =i2
976 END IF
977C
978 xh=xi(i)+bb(i,it)*xn(i,it)
979 yh=yi(i)+bb(i,it)*yn(i,it)
980 zh=zi(i)+bb(i,it)*zn(i,it)
981C
982 IF(vtx_bisector(1,1,ibx)/=zero.OR.
983 . vtx_bisector(2,1,ibx)/=zero.OR.
984 . vtx_bisector(3,1,ibx)/=zero.OR.
985 . vtx_bisector(1,2,ibx)/=zero.OR.
986 . vtx_bisector(2,2,ibx)/=zero.OR.
987 . vtx_bisector(3,2,ibx)/=zero)THEN
988 p1 = (xh-xx(i,ix))* vtx_bisector(1,1,ibx)+
989 . (yh-yy(i,ix))* vtx_bisector(2,1,ibx)+
990 . (zh-zz(i,ix))* vtx_bisector(3,1,ibx)
991 p2 = (xh-xx(i,ix))* vtx_bisector(1,2,ibx)+
992 . (yh-yy(i,ix))* vtx_bisector(2,2,ibx)+
993 . (zh-zz(i,ix))* vtx_bisector(3,2,ibx)
994 IF(p1 >= gaps(i) .AND. p2 >= gaps(i)) fara(i,it) =3
995 ELSE
996 vx = x0n(i,ix) ! fake bisector of angle at vertex IX
997 vy = y0n(i,ix)
998 vz = z0n(i,ix)
999 nn = one/max(em20,sqrt(vx*vx+vy*vy+vz*vz))
1000 pn = ((xh-xx(i,ix))*vx+(yh-yy(i,ix))*vy+(zh-zz(i,ix))*vz)*nn
1001 IF(pn >= gaps(i)) fara(i,it) =3
1002 END IF
1003
1004 END IF
1005C
1006C Cone
1007 IF(ingapa(i,it) == 1 .AND. (fara(i,it)==1 .OR. bb(i,it) <= zero))THEN
1008C ingap = 1 => always check if inside the cone (allows to determine which cone we are in)
1009C
1010 x12=xx(i,i2)-xx(i,i1)
1011 y12=yy(i,i2)-yy(i,i1)
1012 z12=zz(i,i2)-zz(i,i1)
1013C
1014C normal to the bisecting plane (pointing toward the inside)
1015 px = z12*nay(i,it)-y12*naz(i,it)
1016 py = x12*naz(i,it)-z12*nax(i,it)
1017 pz = y12*nax(i,it)-x12*nay(i,it)
1018 pp = px*px+py*py+pz*pz
1019C
1020 ll = xij(i,i1)*xij(i,i1)+yij(i,i1)*yij(i,i1)+zij(i,i1)*zij(i,i1)
1021C
1022 sida(i,it)=-(xij(i,i1)*px+yij(i,i1)*py+zij(i,i1)*pz)*sqrt(one/max(em30,ll*pp))
1023 IF(sida(i,it) < -zep01) fara(i,it) =2
1024C
1025 END IF
1026 ELSE
1027 ib1=ibounda(1,i)
1028 ib2=ibounda(2,i)
1029 ib3=ibounda(3,i)
1030C
1031 d1=(xh-xx(i,1))* nax(i,1)+
1032 . (yh-yy(i,1))* nay(i,1)+
1033 . (zh-zz(i,1))* naz(i,1)
1034 d2=(xh-xx(i,2))* nax(i,2)+
1035 . (yh-yy(i,2))* nay(i,2)+
1036 . (zh-zz(i,2))* naz(i,2)
1037 d3=(xh-xx(i,3))* nax(i,4)+
1038 . (yh-yy(i,3))* nay(i,4)+
1039 . (zh-zz(i,3))* naz(i,4)
1040C
1041 IF( (mvoisa(i,1) == 0 .AND. d1 >= gaps(i)).OR.
1042 . (mvoisa(i,2) == 0 .AND. d2 >= gaps(i)).OR.
1043 . (mvoisa(i,4) == 0 .AND. d3 >= gaps(i)) )THEN
1044C
1045 fara(i,it) =3
1046C
1047 ELSEIF((ib1/=0 .AND. ib2==0 .AND. ib3==0).OR.
1048 . (ib2/=0 .AND. ib3==0 .AND. ib1==0).OR.
1049 . (ib3/=0 .AND. ib1==0 .AND. ib2==0))THEN
1050C
1051 ibx=max(ib1,ib2,ib3)
1052 IF(ib1/=0)THEN
1053 ix =1
1054 iy =2
1055 iz =3
1056 ELSEIF(ib2/=0)THEN
1057 ix =2
1058 iy =3
1059 iz =1
1060 ELSEIF(ib3/=0)THEN
1061 ix =3
1062 iy =1
1063 iz =2
1064 END IF
1065C
1066 IF(vtx_bisector(1,1,ibx)/=zero.OR.
1067 . vtx_bisector(2,1,ibx)/=zero.OR.
1068 . vtx_bisector(3,1,ibx)/=zero.OR.
1069 . vtx_bisector(1,2,ibx)/=zero.OR.
1070 . vtx_bisector(2,2,ibx)/=zero.OR.
1071 . vtx_bisector(3,2,ibx)/=zero)THEN
1072 p1 = (xh-xx(i,ix))* vtx_bisector(1,1,ibx)+
1073 . (yh-yy(i,ix))* vtx_bisector(2,1,ibx)+
1074 . (zh-zz(i,ix))* vtx_bisector(3,1,ibx)
1075 p2 = (xh-xx(i,ix))* vtx_bisector(1,2,ibx)+
1076 . (yh-yy(i,ix))* vtx_bisector(2,2,ibx)+
1077 . (zh-zz(i,ix))* vtx_bisector(3,2,ibx)
1078 IF(p1 >= gaps(i) .AND. p2 >= gaps(i)) fara(i,it) =3
1079 ELSE
1080 vx = two*xx(i,ix)-(xx(i,iy)+xx(i,iz)) ! fake bisector of angle 2,1,3
1081 vy = two*yy(i,ix)-(yy(i,iy)+yy(i,iz))
1082 vz = two*zz(i,ix)-(zz(i,iy)+zz(i,iz))
1083 nn = one/max(em20,sqrt(vx*vx+vy*vy+vz*vz))
1084 pn = ((xh-xx(i,ix))*vx+(yh-yy(i,ix))*vy+(zh-zz(i,ix))*vz)*nn
1085 IF(pn >= gaps(i)) fara(i,it) =3
1086 END IF
1087C
1088 END IF
1089C
1090 IF(ingapa(i,it) == 1 .AND. (fara(i,it)==1 .OR. bb(i,it) <= zero))THEN
1091C
1092C ingap = 1 => always check if inside the cone (allows to determine which cone we are in)
1093C
1094 IF( mvoisa(i,1) /= 0 )THEN
1095C
1096 x12=xx(i,2)-xx(i,1)
1097 y12=yy(i,2)-yy(i,1)
1098 z12=zz(i,2)-zz(i,1)
1099C
1100C normal to the bisecting plane (pointing toward the inside)
1101 px = z12*nay(i,1)-y12*naz(i,1)
1102 py = x12*naz(i,1)-z12*nax(i,1)
1103 pz = y12*nax(i,1)-x12*nay(i,1)
1104 pp = px*px+py*py+pz*pz
1105C
1106 ll = xij(i,1)*xij(i,1)+yij(i,1)*yij(i,1)+zij(i,1)*zij(i,1)
1107C
1108 sida(i,1)=-(xij(i,1)*px+yij(i,1)*py+zij(i,1)*pz)*sqrt(one/max(em30,ll*pp))
1109 IF(sida(i,1) < -zep01) fara(i,it) =2
1110C
1111 END IF
1112
1113 IF( mvoisa(i,2) /= 0 )THEN
1114
1115C ingap = 1 => always check if inside the cone (allows to determine which cone we are in)
1116C
1117 x12=xx(i,3)-xx(i,2)
1118 y12=yy(i,3)-yy(i,2)
1119 z12=zz(i,3)-zz(i,2)
1120C
1121C normal to the bisecting plane (pointing toward the inside)
1122 px = z12*nay(i,2)-y12*naz(i,2)
1123 py = x12*naz(i,2)-z12*nax(i,2)
1124 pz = y12*nax(i,2)-x12*nay(i,2)
1125 pp = px*px+py*py+pz*pz
1126C
1127 ll = xij(i,2)*xij(i,2)+yij(i,2)*yij(i,2)+zij(i,2)*zij(i,2)
1128C
1129 sida(i,2)=-(xij(i,2)*px+yij(i,2)*py+zij(i,2)*pz)*sqrt(one/max(em30,ll*pp))
1130 IF(sida(i,2) < -zep01) fara(i,it) =2
1131C
1132 END IF
1133
1134 IF( mvoisa(i,4) /= 0 )THEN
1135C ingap = 1 => always check if inside the cone (allows to determine which cone we are in)
1136C
1137 x12=xx(i,1)-xx(i,3)
1138 y12=yy(i,1)-yy(i,3)
1139 z12=zz(i,1)-zz(i,3)
1140C
1141C normal to the bisecting plane (pointing toward the inside)
1142 px = z12*nay(i,4)-y12*naz(i,4)
1143 py = x12*naz(i,4)-z12*nax(i,4)
1144 pz = y12*nax(i,4)-x12*nay(i,4)
1145 pp = px*px+py*py+pz*pz
1146C
1147 ll = xij(i,3)*xij(i,3)+yij(i,3)*yij(i,3)+zij(i,3)*zij(i,3)
1148C
1149 sida(i,4)=-(xij(i,3)*px+yij(i,3)*py+zij(i,3)*pz)*sqrt(one/max(em30,ll*pp))
1150 IF(sida(i,4) < -zep01) fara(i,it) =2
1151C
1152 END IF
1153 END IF
1154 END IF
1155C
1156 IF(fara(i,it)==2 .AND. intersecta(i)==0) pena(i,it)=zero
1157 IF(fara(i,it)==3) pena(i,it)=zero
1158C
1159 ENDDO
1160C--------------------------------------------------------
1161C calculation of penetration
1162C--------------------------------------------------------
1163#include "vectorize.inc"
1164 DO i =1,jlt
1165C
1166C IF (STIF(I) == ZERO)CYCLE
1167C
1168 it=subtrib(i)
1169 IF(it==0)cycle
1170C
1171 i1=itria(1,it)
1172 i2=itria(2,it)
1173C
1174 lap = one-lbp(i,it)-lcp(i,it)
1175 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),
1176 . max(drad,gapmxl(i)+dgapload))
1177C
1178 IF(bb(i,it) < zero)THEN
1179C
1180C penetration approximee (cf zone limite interpolation des normales)
1181 penb(i,it)=max(zero,gap-bb(i,it))
1182 ELSE ! IF(BB(I,IT) < ZERO)THEN
1183 penb(i,it)=max(zero,gap-sqrt(dd(i,it)))
1184 END IF
1185C
1186 ENDDO
1187C--------------------------------------------------------
1188C Check vs bissectors ...
1189C--------------------------------------------------------
1190#include "vectorize.inc"
1191 DO i =1,jlt
1192C
1193C IF (STIF(I) == ZERO)CYCLE
1194C
1195 it=subtrib(i)
1196 IF(it==0)cycle
1197C
1198 IF(penb(i,it)==zero) cycle
1199C
1200 i1=itria(1,it)
1201 i2=itria(2,it)
1202C
1203 xh=xi(i)+bb(i,it)*xn(i,it)
1204 yh=yi(i)+bb(i,it)*yn(i,it)
1205 zh=zi(i)+bb(i,it)*zn(i,it)
1206C
1207 IF(ix3(i) /= ix4(i))THEN
1208C
1209C Cone
1210 ib1=iboundb(i1,i)
1211 ib2=iboundb(i2,i)
1212 IF(mvoisb(i,it)==0)THEN
1213C
1214 IF( (xh-xx(i,i1))* nbx(i,it)+
1215 . (yh-yy(i,i1))* nby(i,it)+
1216 . (zh-zz(i,i1))* nbz(i,it) >= gaps(i)) farb(i,it) =3
1217C
1218 ELSEIF((ib1 /= 0 .AND. ib2 == 0).OR.
1219 . (ib2 /= 0 .AND. ib1 == 0))THEN
1220C
1221 ibx=max(ib1,ib2)
1222 IF(ib1/=0)THEN
1223 ix =i1
1224 ELSEIF(ib2/=0)THEN
1225 ix =i2
1226 END IF
1227C
1228 xh=xi(i)+bb(i,it)*xn(i,it)
1229 yh=yi(i)+bb(i,it)*yn(i,it)
1230 zh=zi(i)+bb(i,it)*zn(i,it)
1231C
1232 IF(vtx_bisector(1,1,ibx)/=zero.OR.
1233 . vtx_bisector(2,1,ibx)/=zero.OR.
1234 . vtx_bisector(3,1,ibx)/=zero.OR.
1235 . vtx_bisector(1,2,ibx)/=zero.OR.
1236 . vtx_bisector(2,2,ibx)/=zero.OR.
1237 . vtx_bisector(3,2,ibx)/=zero)THEN
1238 p1 = (xh-xx(i,ix))* vtx_bisector(1,1,ibx)+
1239 . (yh-yy(i,ix))* vtx_bisector(2,1,ibx)+
1240 . (zh-zz(i,ix))* vtx_bisector(3,1,ibx)
1241 p2 = (xh-xx(i,ix))* vtx_bisector(1,2,ibx)+
1242 . (yh-yy(i,ix))* vtx_bisector(2,2,ibx)+
1243 . (zh-zz(i,ix))* vtx_bisector(3,2,ibx)
1244 IF(p1 >= gaps(i) .AND. p2 >= gaps(i)) farb(i,it) =3
1245 ELSE
1246 vx = x0n(i,ix) ! fake bisector of angle at vertex IX
1247 vy = y0n(i,ix)
1248 vz = z0n(i,ix)
1249 nn = one/max(em20,sqrt(vx*vx+vy*vy+vz*vz))
1250 pn = ((xh-xx(i,ix))*vx+(yh-yy(i,ix))*vy+(zh-zz(i,ix))*vz)*nn
1251 IF(pn >= gaps(i)) farb(i,it) =3
1252 END IF
1253
1254 END IF
1255C
1256C
1257 IF(ingapb(i,it) == 1 .AND. (farb(i,it)==1 .OR. bb(i,it) >= zero))THEN
1258C
1259C ingap = 1 => always check if inside the cone (allows to determine which cone we are in)
1260C
1261 x12=xx(i,i2)-xx(i,i1)
1262 y12=yy(i,i2)-yy(i,i1)
1263 z12=zz(i,i2)-zz(i,i1)
1264C
1265C normal to the bisecting plane (pointing toward the inside)
1266 px = z12*nby(i,it)-y12*nbz(i,it)
1267 py = x12*nbz(i,it)-z12*nbx(i,it)
1268 pz = y12*nbx(i,it)-x12*nby(i,it)
1269 pp = px*px+py*py+pz*pz
1270C
1271 ll = xij(i,i1)*xij(i,i1)+yij(i,i1)*yij(i,i1)+zij(i,i1)*zij(i,i1)
1272C
1273 sidb(i,it)= (xij(i,i1)*px+yij(i,i1)*py+zij(i,i1)*pz)*sqrt(one/max(em30,ll*pp))
1274 IF(sidb(i,it) < -zep01) farb(i,it) =2
1275C
1276 END IF
1277 ELSE
1278 ib1=iboundb(1,i)
1279 ib2=iboundb(2,i)
1280 ib3=iboundb(3,i)
1281C
1282 d1=(xh-xx(i,1))* nbx(i,1)+
1283 . (yh-yy(i,1))* nby(i,1)+
1284 . (zh-zz(i,1))* nbz(i,1)
1285 d2=(xh-xx(i,2))* nbx(i,2)+
1286 . (yh-yy(i,2))* nby(i,2)+
1287 . (zh-zz(i,2))* nbz(i,2)
1288 d3=(xh-xx(i,3))* nbx(i,4)+
1289 . (yh-yy(i,3))* nby(i,4)+
1290 . (zh-zz(i,3))* nbz(i,4)
1291C
1292 IF( (mvoisb(i,1) == 0 .AND. d1 >= gaps(i)).OR.
1293 . (mvoisb(i,2) == 0 .AND. d2 >= gaps(i)).OR.
1294 . (mvoisb(i,4) == 0 .AND. d3 >= gaps(i)) )THEN
1295C
1296 farb(i,it) =3
1297C
1298 ELSEIF((ib1/=0 .AND. ib2==0 .AND. ib3==0).OR.
1299 . (ib2/=0 .AND. ib3==0 .AND. ib1==0).OR.
1300 . (ib3/=0 .AND. ib1==0 .AND. ib2==0))THEN
1301C
1302 ibx=max(ib1,ib2,ib3)
1303 IF(ib1/=0)THEN
1304 ix =1
1305 iy =2
1306 iz =3
1307 ELSEIF(ib2/=0)THEN
1308 ix =2
1309 iy =3
1310 iz =1
1311 ELSEIF(ib3/=0)THEN
1312 ix =3
1313 iy =1
1314 iz =2
1315 END IF
1316C
1317 IF(vtx_bisector(1,1,ibx)/=zero.OR.
1318 . vtx_bisector(2,1,ibx)/=zero.OR.
1319 . vtx_bisector(3,1,ibx)/=zero.OR.
1320 . vtx_bisector(1,2,ibx)/=zero.OR.
1321 . vtx_bisector(2,2,ibx)/=zero.OR.
1322 . vtx_bisector(3,2,ibx)/=zero)THEN
1323 p1 = (xh-xx(i,ix))* vtx_bisector(1,1,ibx)+
1324 . (yh-yy(i,ix))* vtx_bisector(2,1,ibx)+
1325 . (zh-zz(i,ix))* vtx_bisector(3,1,ibx)
1326 p2 = (xh-xx(i,ix))* vtx_bisector(1,2,ibx)+
1327 . (yh-yy(i,ix))* vtx_bisector(2,2,ibx)+
1328 . (zh-zz(i,ix))* vtx_bisector(3,2,ibx)
1329 IF(p1 >= gaps(i) .AND. p2 >= gaps(i)) farb(i,it) =3
1330 ELSE
1331 vx = two*xx(i,ix)-(xx(i,iy)+xx(i,iz)) ! fake bisector of angle 2,1,3
1332 vy = two*yy(i,ix)-(yy(i,iy)+yy(i,iz))
1333 vz = two*zz(i,ix)-(zz(i,iy)+zz(i,iz))
1334 nn = one/max(em20,sqrt(vx*vx+vy*vy+vz*vz))
1335 pn = ((xh-xx(i,ix))*vx+(yh-yy(i,ix))*vy+(zh-zz(i,ix))*vz)*nn
1336 IF(pn >= gaps(i)) farb(i,it) =3
1337 END IF
1338C
1339 END IF
1340C
1341 IF(ingapb(i,it) == 1 .AND. (farb(i,it)==1 .OR. bb(i,it) >= zero))THEN
1342C
1343C ingap = 1 => always check if inside the cone (allows to determine which cone we are in)
1344C
1345 IF( mvoisb(i,1) /= 0 )THEN
1346C
1347 x12=xx(i,2)-xx(i,1)
1348 y12=yy(i,2)-yy(i,1)
1349 z12=zz(i,2)-zz(i,1)
1350C
1351C normal to the bisecting plane (pointing toward the inside)
1352 px = z12*nby(i,1)-y12*nbz(i,1)
1353 py = x12*nbz(i,1)-z12*nbx(i,1)
1354 pz = y12*nbx(i,1)-x12*nby(i,1)
1355 pp = px*px+py*py+pz*pz
1356C
1357 ll = xij(i,1)*xij(i,1)+yij(i,1)*yij(i,1)+zij(i,1)*zij(i,1)
1358C
1359 sidb(i,1)= (xij(i,1)*px+yij(i,1)*py+zij(i,1)*pz)*sqrt(one/max(em30,ll*pp))
1360 IF(sidb(i,1) < -zep01) farb(i,it) =2
1361C
1362 END IF
1363
1364 IF( mvoisb(i,2) /= 0 )THEN
1365C
1366 x12=xx(i,3)-xx(i,2)
1367 y12=yy(i,3)-yy(i,2)
1368 z12=zz(i,3)-zz(i,2)
1369C
1370C normal to the bisecting plane (pointing toward the inside)
1371 px = z12*nby(i,2)-y12*nbz(i,2)
1372 py = x12*nbz(i,2)-z12*nbx(i,2)
1373 pz = y12*nbx(i,2)-x12*nby(i,2)
1374 pp = px*px+py*py+pz*pz
1375C
1376 ll = xij(i,2)*xij(i,2)+yij(i,2)*yij(i,2)+zij(i,2)*zij(i,2)
1377C
1378 sidb(i,2)= (xij(i,2)*px+yij(i,2)*py+zij(i,2)*pz)*sqrt(one/max(em30,ll*pp))
1379 IF(sidb(i,2) < -zep01) farb(i,it) =2
1380C
1381 END IF
1382
1383 IF( mvoisb(i,4) /= 0 )THEN
1384C
1385 x12=xx(i,1)-xx(i,3)
1386 y12=yy(i,1)-yy(i,3)
1387 z12=zz(i,1)-zz(i,3)
1388C
1389C normal to the bisecting plane (pointing toward the inside)
1390 px = z12*nby(i,4)-y12*nbz(i,4)
1391 py = x12*nbz(i,4)-z12*nbx(i,4)
1392 pz = y12*nbx(i,4)-x12*nby(i,4)
1393 pp = px*px+py*py+pz*pz
1394C
1395 ll = xij(i,3)*xij(i,3)+yij(i,3)*yij(i,3)+zij(i,3)*zij(i,3)
1396C
1397 sidb(i,4)= (xij(i,3)*px+yij(i,3)*py+zij(i,3)*pz)*sqrt(one/max(em30,ll*pp))
1398 IF(sidb(i,4) < -zep01) farb(i,it) =2
1399C
1400 END IF
1401 END IF
1402 END IF
1403C
1404 IF(farb(i,it)==2 .AND. intersectb(i)==0) penb(i,it)=zero
1405 IF(farb(i,it)==3) penb(i,it)=zero
1406C
1407 ENDDO
1408C--------------------------------------------------------
1409 DO i=1,jlt
1410C
1411 IF(stif(i) <= zero)cycle
1412C
1413 ia = subtria(i)
1414 penax = zero
1415 IF(ia/=0)penax=pena(i,ia)
1416C
1417 ib = subtrib(i)
1418 penbx = zero
1419 IF(ib/=0)penbx=penb(i,ib)
1420C
1421 IF(penax > penbx .AND. ia > 0)THEN
1422 l = cand_e(i)
1423 it = ia
1424 pent(i)= penax
1425
1426 lbs(i)=lbp(i,it)
1427 lcs(i)=lcp(i,it)
1428 far(i)=fara(i,it)
1429
1430 ELSEIF(penbx > penax .AND. ib > 0)THEN
1431 l = ishel(i)
1432 cand_e(i) = l
1433
1434 it = ib
1435 pent(i) = penbx
1436
1437 lbs(i)=lcp(i,it)
1438 lcs(i)=lbp(i,it)
1439 far(i)=farb(i,it)
1440
1441 subtria(i)=itperm(it) ! IT/=0
1442 it = subtria(i)
1443
1444 ELSE
1445 it=0
1446 subtria(i) = 0
1447 pent(i) = zero
1448 END IF
1449
1450c if(itab(nsvg(i))==2421446.or.
1451c . itab(ix1(i))==2421446.or.itab(ix2(i))==2421446.or.itab(ix3(i))==2421446.or.itab(ix4(i))==2421446)
1452c . print *,'dst2 nat',ispmd+1,
1453c . itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),
1454c . it,ia,ib,penax,penbx,
1455c . ingapa(i,1:4),intersecta(i),ingapb(i,1:4),intersectb(i),
1456c . bb(i,1:4),dd(i,1:4)
1457
1458 IF(it == 0)cycle
1459 pene = pent(i)
1460
1461 n = cand_n(i)
1462
1463 mglob=mseglo(l)
1464
1465 IF(n<=nsn)THEN
1466#include "lockon.inc"
1467
1468 IF( time_s(1,n) < pene .OR.
1469 * (time_s(1,n) == pene .AND. -irtlm(1,n) < mglob ))THEN
1470 irtlm(1,n) = -mglob
1471 irtlm(2,n) = it
1472 irtlm(3,n) = cand_e(i)
1473 irtlm(4,n) = ispmd+1
1474 time_s(1,n)= pene
1475C
1476C PENE_OLD(5,N)=PENE
1477 END IF
1478c if(itab(nsvg(i))==27363)then
1479cc if(itab(nsvg(i))==27952.or.
1480cc . itab(ix1(i))==27952.or.itab(ix2(i))==27952.or.itab(ix3(i))==27952.or.itab(ix4(i))==27952)then
1481c print *,'dst2 nat',ispmd+1,
1482c . itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),
1483c . irtlm(1,n),cand_e(i),TIME_S(1,N),
1484c . it,pent(i),penax,penbx,intersecta(i),intersectb(i),lbs(i),lcs(i),lbs(i)+lcs(i)
1485c if(ia/=0)
1486c . print *,'ia',ia,bb(i,ia),ingapa(i,ia),fara(i,ia),ibounda(1:4,i)
1487c if(ib/=0)
1488c . print *,'ib',ib,bb(i,ib),ingapb(i,ib),farb(i,ib),iboundb(1:4,i)
1489c end if
1490#include "lockoff.inc"
1491 ELSE
1492#include "lockon.inc"
1493 IF( time_sfi(nin)%P(2*(n-nsn-1)+1) < pene .OR.
1494 * (time_sfi(nin)%P(2*(n-nsn-1)+1) == pene .AND. -irtlm_fi(nin)%P(1,n-nsn) < mglob ))THEN
1495 irtlm_fi(nin)%P(1,n-nsn) = -mglob
1496 irtlm_fi(nin)%P(2,n-nsn) = it
1497 irtlm_fi(nin)%P(3,n-nsn) = cand_e(i)
1498 irtlm_fi(nin)%P(4,n-nsn) = ispmd+1
1499 time_sfi(nin)%P(2*(n-nsn-1)+1) = pene
1500C
1501C PENE_OLDFI(NIN)%P(5,N-NSN)=PENE
1502 END IF
1503c if(itafi(nin)%p(n-nsn)==3817238)
1504c . print *,'dst2 rem',ispmd+1,
1505c . itafi(nin)%p(n-nsn),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),
1506c . it,IRTLM_FI(NIN)%P(1,n-nsn),mglob,cand_e(i),pent(i,it),far(i,it),
1507c . TIME_SFI(NIN)%P(N-NSN),dd(i,it),ingap(i,it),intersect(i)
1508#include "lockoff.inc"
1509 END IF
1510 END DO
1511
1512 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
if
if(complex_arithmetic) id
intersecb_25
subroutine intersecb_25(ixx3, ixx4, interseca, intersecb, v12, v23, v34, v41, vo12, vo23, vo34, vo41, xx1, yy1, zz1, xx2, yy2, zz2, xx3, yy3, zz3, xx4, yy4, zz4, xx5, yy5, zz5, vxi, vyi, vzi, xo1, xo2, xo3, xo4, xo5, yo1, yo2, yo3, yo4, yo5, zo1, zo2, zo3, zo4, zo5, xi, yi, zi, xoi, yoi, zoi, zerom, unp, zeromt, unpt)
Definition i25dst3_22.F:1898
intersecv0_25
subroutine intersecv0_25(ixx3, ixx4, intersect, v12, v23, v34, v41, vo12, vo23, vo34, vo41, xx1, yy1, zz1, xx2, yy2, zz2, xx3, yy3, zz3, xx4, yy4, zz4, xx5, yy5, zz5, xoi, yoi, zoi, xo1, xo2, xo3, xo4, xo5, yo1, yo2, yo3, yo4, yo5, zo1, zo2, zo3, zo4, zo5, xi, yi, zi, zeromt, unpt)
Definition i25dst3_22.F:2684
interseca_25
subroutine interseca_25(ixx3, ixx4, interseca, v12, v23, v34, v41, vo12, vo23, vo34, vo41, xx1, yy1, zz1, xx2, yy2, zz2, xx3, yy3, zz3, xx4, yy4, zz4, xx5, yy5, zz5, vxi, vyi, vzi, xo1, xo2, xo3, xo4, xo5, yo1, yo2, yo3, yo4, yo5, zo1, zo2, zo3, zo4, zo5, xi, yi, zi, xoi, yoi, zoi, zerom, unp, zeromt, unpt)
Definition i25dst3_22.F:1532
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
tri7box::time_sfi
type(real_pointer), dimension(:), allocatable time_sfi
Definition tri7box.F:542
tri7box::irtlm_fi
type(int_pointer2), dimension(:), allocatable irtlm_fi
Definition tri7box.F:533
tri7box::icont_i_fi
type(int_pointer), dimension(:), allocatable icont_i_fi
Definition tri7box.F:532

◆ i25glob_22()

subroutine i25glob_22 ( integer jlt,
integer, dimension(*) cand_n_n,
integer, dimension(*) cand_e_n,
integer, dimension(*) cand_e,
integer nin,
integer nsn,
integer, dimension(mvsiz) ix1,
integer, dimension(mvsiz) ix2,
integer, dimension(mvsiz) ix3,
integer, dimension(mvsiz) ix4,
integer, dimension(mvsiz) nsvg,
stif,
integer inacti,
integer, dimension(*) mseglo,
integer, dimension(4,nsn) irtlm,
time_s,
integer, dimension(*) itab,
integer, dimension(mvsiz) subtria,
integer, dimension(mvsiz) far,
pent,
lb,
lc,
integer, dimension(*) index,
integer, dimension(4,*) farm,
penm,
lbm,
lcm )

Definition at line 3013 of file i25dst3_22.F.

3021C-----------------------------------------------
3022C M o d u l e s
3023C-----------------------------------------------
3024 USE tri7box
3025C-----------------------------------------------
3026C I m p l i c i t T y p e s
3027C-----------------------------------------------
3028#include "implicit_f.inc"
3029C-----------------------------------------------
3030C G l o b a l P a r a m e t e r s
3031C-----------------------------------------------
3032#include "mvsiz_p.inc"
3033#include "comlock.inc"
3034C-----------------------------------------------
3035C D u m m y A r g u m e n t s
3036C-----------------------------------------------
3037 INTEGER JLT, NIN, NSN, INACTI,
3038 . CAND_N_N(*), CAND_E_N(*), CAND_E(*), ITAB(*), SUBTRIA(MVSIZ)
3039 INTEGER NSVG(MVSIZ), IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
3040 . MSEGLO(*), IRTLM(4,NSN) ,INDEX(*), FARM(4,*),
3041 . FAR(MVSIZ)
3042 my_real
3043 . stif(*), time_s(*),
3044 . pent(mvsiz),
3045 . lb(mvsiz), lc(mvsiz),
3046 . penm(4,*), lbm(4,*), lcm(4,*)
3047C-----------------------------------------------
3048C L o c a l V a r i a b l e s
3049C-----------------------------------------------
3050 INTEGER IT, I, J
3051C--------------------------------------------------------
3052C Back to global arrays
3053C--------------------------------------------------------
3054 DO i=1,jlt
3055 j=index(i)
3056 it=subtria(i)
3057
3058 farm(1:4,j)=0
3059 penm(1:4,j)=zero
3060 lbm(1:4,j) =zero
3061 lcm(1:4,j) =zero
3062
3063 IF(it/=0)THEN
3064 cand_e(j)=cand_e_n(i)
3065
3066 farm(it,j)=far(i)
3067 penm(it,j)=pent(i)
3068 lbm(it,j) =lb(i)
3069 lcm(it,j) =lc(i)
3070 END IF
3071
3072 END DO
3073 RETURN

◆ interseca_25()

subroutine interseca_25 ( integer ixx3,
integer ixx4,
integer interseca,
v12,
v23,
v34,
v41,
vo12,
vo23,
vo34,
vo41,
xx1,
yy1,
zz1,
xx2,
yy2,
zz2,
xx3,
yy3,
zz3,
xx4,
yy4,
zz4,
xx5,
yy5,
zz5,
vxi,
vyi,
vzi,
xo1,
xo2,
xo3,
xo4,
xo5,
yo1,
yo2,
yo3,
yo4,
yo5,
zo1,
zo2,
zo3,
zo4,
zo5,
xi,
yi,
zi,
xoi,
yoi,
zoi,
zerom,
unp,
zeromt,
unpt )

Definition at line 1519 of file i25dst3_22.F.

1532C-----------------------------------------------
1533C I m p l i c i t T y p e s
1534C-----------------------------------------------
1535#include "implicit_f.inc"
1536C-----------------------------------------------
1537C D u m m y A r g u m e n t s
1538C-----------------------------------------------
1539 INTEGER IXX3 ,IXX4 , INTERSECA
1540C REAL
1541 my_real
1542 1 v12 ,v23 ,v34 ,v41 ,
1543 2 vo12 ,vo23 ,vo34 ,vo41 ,xx1 ,
1544 3 yy1 ,zz1 ,xx2 ,yy2 ,zz2 ,
1545 4 xx3 ,yy3 ,zz3 ,xx4 ,yy4 ,
1546 5 zz4 ,xx5 ,yy5 ,zz5 ,
1547 6 vxi ,vyi ,vzi ,xo1 ,xo2 ,
1548 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
1549 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
1550 9 zo3 ,zo4 ,zo5 ,zerom ,unp ,
1551 a zeromt ,unpt ,xi ,yi ,zi ,
1552 b xoi ,yoi ,zoi
1553C-----------------------------------------------
1554C L o c a l V a r i a b l e s
1555C-----------------------------------------------
1556 my_real
1557 . x51, x52, x53, x54,
1558 . y51, y52, y53, y54,
1559 . z51, z52, z53, z54,
1560 . xi1, xi2, xi3, xi4, xi5,
1561 . yi1, yi2, yi3, yi4, yi5,
1562 . zi1, zi2, zi3, zi4, zi5,
1563 .
1564 .
1565 .
1566 . xs, ys, zs,
1567 . xm1, xm2, xm3, xm4, xm5,
1568 . ym1, ym2, ym3, ym4, ym5,
1569 . zm1, zm2, zm3, zm4, zm5,la0,lb0,lc0
1570 my_real
1571 . sx1,sx2,sx3,sx4,sx5,sx0,
1572 . sy1,sy2,sy3,sy4,sy5,sy0,
1573 . sz1,sz2,sz3,sz4,sz5,sz0,aaa,s2
1574C------------------------------------------------
1575 interseca = 0
1576
1577 IF(vo12 <= zero .AND. v12 >= zero)THEN
1578 IF(abs(vo12) < em20)THEN
1579 aaa = one
1580 ELSEIF(abs(v12) < em20)THEN
1581 aaa = zero
1582 ELSE
1583 aaa = v12 / (v12-vo12)
1584 ENDIF
1585
1586 xm1 = xx1 + aaa*(xo1-xx1)
1587 ym1 = yy1 + aaa*(yo1-yy1)
1588 zm1 = zz1 + aaa*(zo1-zz1)
1589
1590 xm2 = xx2 + aaa*(xo2-xx2)
1591 ym2 = yy2 + aaa*(yo2-yy2)
1592 zm2 = zz2 + aaa*(zo2-zz2)
1593
1594 xs = xi + aaa*(xoi-xi)
1595 ys = yi + aaa*(yoi-yi)
1596 zs = zi + aaa*(zoi-zi)
1597
1598 xm5 = xx5 + aaa*(xo5-xx5)
1599 ym5 = yy5 + aaa*(yo5-yy5)
1600 zm5 = zz5 + aaa*(zo5-zz5)
1601
1602 x51 = xm1 - xm5
1603 y51 = ym1 - ym5
1604 z51 = zm1 - zm5
1605
1606 x52 = xm2 - xm5
1607 y52 = ym2 - ym5
1608 z52 = zm2 - zm5
1609
1610 xi1 = xm1 - xs
1611 yi1 = ym1 - ys
1612 zi1 = zm1 - zs
1613
1614 xi2 = xm2 - xs
1615 yi2 = ym2 - ys
1616 zi2 = zm2 - zs
1617
1618 xi5 = xm5 - xs
1619 yi5 = ym5 - ys
1620 zi5 = zm5 - zs
1621
1622 sx0 = y51*z52 - z51*y52
1623 sy0 = z51*x52 - x51*z52
1624 sz0 = x51*y52 - y51*x52
1625
1626 sx1 = yi5*zi1 - zi5*yi1
1627 sy1 = zi5*xi1 - xi5*zi1
1628 sz1 = xi5*yi1 - yi5*xi1
1629
1630 sx5 = yi1*zi2 - zi1*yi2
1631 sy5 = zi1*xi2 - xi1*zi2
1632 sz5 = xi1*yi2 - yi1*xi2
1633
1634 sx2 = yi2*zi5 - zi2*yi5
1635 sy2 = zi2*xi5 - xi2*zi5
1636 sz2 = xi2*yi5 - yi2*xi5
1637
1638 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
1639 lb0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
1640 lc0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
1641 la0 = one-lb0-lc0
1642
1643 IF(ixx3 /= ixx4)THEN
1644 IF(lb0 >= zerom .and. lb0 <= unp .and.
1645 . lc0 >= zerom .and. lc0 <= unp .and.
1646 . la0 >= zerom .and. la0 <= unp)THEN
1647 interseca = 1
1648 ENDIF
1649C----------loose tolerance for tria
1650 ELSE
1651 IF(lb0 >= zeromt .AND. lb0 <= unpt .AND.
1652 . lc0 >= zeromt .AND. lc0 <= unpt .AND.
1653 . la0 >= zeromt .AND. la0 <= unpt)THEN
1654 interseca = 1
1655 ENDIF
1656 END if! (IXX(I,3) /= IXX(I,4))THEN
1657 ENDIF
1658 IF(ixx3 /= ixx4)THEN
1659 IF(vo23 <= zero .AND. v23 >= zero)THEN
1660 IF(abs(vo23) < em20)THEN
1661 aaa = one
1662 ELSEIF(abs(v23) < em20)THEN
1663 aaa = zero
1664 ELSE
1665 aaa = v23 / (v23-vo23)
1666 ENDIF
1667
1668 xm2 = xx2 + aaa*(xo2-xx2)
1669 ym2 = yy2 + aaa*(yo2-yy2)
1670 zm2 = zz2 + aaa*(zo2-zz2)
1671
1672 xm3 = xx3 + aaa*(xo3-xx3)
1673 ym3 = yy3 + aaa*(yo3-yy3)
1674 zm3 = zz3 + aaa*(zo3-zz3)
1675
1676 xs = xi + aaa*(xoi-xi)
1677 ys = yi + aaa*(yoi-yi)
1678 zs = zi + aaa*(zoi-zi)
1679
1680 xm5 = xx5 + aaa*(xo5-xx5)
1681 ym5 = yy5 + aaa*(yo5-yy5)
1682 zm5 = zz5 + aaa*(zo5-zz5)
1683
1684 x52 = xm2 - xm5
1685 y52 = ym2 - ym5
1686 z52 = zm2 - zm5
1687
1688 x53 = xm3 - xm5
1689 y53 = ym3 - ym5
1690 z53 = zm3 - zm5
1691
1692 xi2 = xm2 - xs
1693 yi2 = ym2 - ys
1694 zi2 = zm2 - zs
1695
1696 xi3 = xm3 - xs
1697 yi3 = ym3 - ys
1698 zi3 = zm3 - zs
1699
1700 xi5 = xm5 - xs
1701 yi5 = ym5 - ys
1702 zi5 = zm5 - zs
1703
1704 sx0 = y52*z53 - z52*y53
1705 sy0 = z52*x53 - x52*z53
1706 sz0 = x52*y53 - y52*x53
1707
1708 sx2 = yi5*zi2 - zi5*yi2
1709 sy2 = zi5*xi2 - xi5*zi2
1710 sz2 = xi5*yi2 - yi5*xi2
1711
1712 sx5 = yi2*zi3 - zi2*yi3
1713 sy5 = zi2*xi3 - xi2*zi3
1714 sz5 = xi2*yi3 - yi2*xi3
1715
1716 sx3 = yi3*zi5 - zi3*yi5
1717 sy3 = zi3*xi5 - xi3*zi5
1718 sz3 = xi3*yi5 - yi3*xi5
1719
1720 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
1721 lb0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
1722 lc0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
1723 la0 = one-lb0-lc0
1724
1725 IF(lb0 >= zerom .and. lb0 <= unp .and.
1726 . lc0 >= zerom .and. lc0 <= unp .and.
1727 . la0 >= zerom .and. la0 <= unp)THEN
1728 interseca = 1
1729 ENDIF
1730 ENDIF
1731 IF(vo34 <= zero .AND. v34 >= zero)THEN
1732 IF(abs(vo34) < em20)THEN
1733 aaa = one
1734 ELSEIF(abs(v34) < em20)THEN
1735 aaa = zero
1736 ELSE
1737 aaa = v34 / (v34-vo34)
1738 ENDIF
1739
1740 xm3 = xx3 + aaa*(xo3-xx3)
1741 ym3 = yy3 + aaa*(yo3-yy3)
1742 zm3 = zz3 + aaa*(zo3-zz3)
1743
1744 xm4 = xx4 + aaa*(xo4-xx4)
1745 ym4 = yy4 + aaa*(yo4-yy4)
1746 zm4 = zz4 + aaa*(zo4-zz4)
1747
1748 xs = xi + aaa*(xoi-xi)
1749 ys = yi + aaa*(yoi-yi)
1750 zs = zi + aaa*(zoi-zi)
1751
1752 xm5 = xx5 + aaa*(xo5-xx5)
1753 ym5 = yy5 + aaa*(yo5-yy5)
1754 zm5 = zz5 + aaa*(zo5-zz5)
1755
1756 x53 = xm3 - xm5
1757 y53 = ym3 - ym5
1758 z53 = zm3 - zm5
1759
1760 x54 = xm4 - xm5
1761 y54 = ym4 - ym5
1762 z54 = zm4 - zm5
1763
1764 xi3 = xm3 - xs
1765 yi3 = ym3 - ys
1766 zi3 = zm3 - zs
1767
1768 xi4 = xm4 - xs
1769 yi4 = ym4 - ys
1770 zi4 = zm4 - zs
1771
1772 xi5 = xm5 - xs
1773 yi5 = ym5 - ys
1774 zi5 = zm5 - zs
1775
1776 sx0 = y53*z54 - z53*y54
1777 sy0 = z53*x54 - x53*z54
1778 sz0 = x53*y54 - y53*x54
1779
1780 sx3 = yi5*zi3 - zi5*yi3
1781 sy3 = zi5*xi3 - xi5*zi3
1782 sz3 = xi5*yi3 - yi5*xi3
1783
1784 sx5 = yi3*zi4 - zi3*yi4
1785 sy5 = zi3*xi4 - xi3*zi4
1786 sz5 = xi3*yi4 - yi3*xi4
1787
1788 sx4 = yi4*zi5 - zi4*yi5
1789 sy4 = zi4*xi5 - xi4*zi5
1790 sz4 = xi4*yi5 - yi4*xi5
1791
1792 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
1793 lb0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
1794 lc0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
1795 la0 = one-lb0-lc0
1796
1797 IF(lb0 >= zerom .and. lb0 <= unp .and.
1798 . lc0 >= zerom .and. lc0 <= unp .and.
1799 . la0 >= zerom .and. la0 <= unp)THEN
1800 interseca = 1
1801 ENDIF
1802 ENDIF
1803
1804 IF(vo41 <= zero .AND. v41 >= zero)THEN
1805 IF(abs(vo41) < em20)THEN
1806 aaa = one
1807 ELSEIF(abs(v41) < em20)THEN
1808 aaa = zero
1809 ELSE
1810 aaa = v41 / (v41-vo41)
1811 ENDIF
1812
1813 xm4 = xx4 + aaa*(xo4-xx4)
1814 ym4 = yy4 + aaa*(yo4-yy4)
1815 zm4 = zz4 + aaa*(zo4-zz4)
1816
1817 xm1 = xx1 + aaa*(xo1-xx1)
1818 ym1 = yy1 + aaa*(yo1-yy1)
1819 zm1 = zz1 + aaa*(zo1-zz1)
1820
1821 xs = xi + aaa*(xoi-xi)
1822 ys = yi + aaa*(yoi-yi)
1823 zs = zi + aaa*(zoi-zi)
1824
1825 xm5 = xx5 + aaa*(xo5-xx5)
1826 ym5 = yy5 + aaa*(yo5-yy5)
1827 zm5 = zz5 + aaa*(zo5-zz5)
1828
1829 x54 = xm4 - xm5
1830 y54 = ym4 - ym5
1831 z54 = zm4 - zm5
1832
1833 x51 = xm1 - xm5
1834 y51 = ym1 - ym5
1835 z51 = zm1 - zm5
1836
1837 xi4 = xm4 - xs
1838 yi4 = ym4 - ys
1839 zi4 = zm4 - zs
1840
1841 xi1 = xm1 - xs
1842 yi1 = ym1 - ys
1843 zi1 = zm1 - zs
1844
1845 xi5 = xm5 - xs
1846 yi5 = ym5 - ys
1847 zi5 = zm5 - zs
1848
1849 sx0 = y54*z51 - z54*y51
1850 sy0 = z54*x51 - x54*z51
1851 sz0 = x54*y51 - y54*x51
1852
1853 sx4 = yi5*zi4 - zi5*yi4
1854 sy4 = zi5*xi4 - xi5*zi4
1855 sz4 = xi5*yi4 - yi5*xi4
1856
1857 sx5 = yi4*zi1 - zi4*yi1
1858 sy5 = zi4*xi1 - xi4*zi1
1859 sz5 = xi4*yi1 - yi4*xi1
1860
1861 sx1 = yi1*zi5 - zi1*yi5
1862 sy1 = zi1*xi5 - xi1*zi5
1863 sz1 = xi1*yi5 - yi1*xi5
1864
1865 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
1866 lb0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
1867 lc0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
1868 la0 = one-lb0-lc0
1869
1870 IF(lb0 >= zerom .and. lb0 <= unp .and.
1871 . lc0 >= zerom .and. lc0 <= unp .and.
1872 . la0 >= zerom .and. la0 <= unp)THEN
1873 interseca = 1
1874 ENDIF
1875 ENDIF
1876 ENDIF
1877C
1878 RETURN

◆ intersecb_25()

subroutine intersecb_25 ( integer ixx3,
integer ixx4,
integer interseca,
integer intersecb,
v12,
v23,
v34,
v41,
vo12,
vo23,
vo34,
vo41,
xx1,
yy1,
zz1,
xx2,
yy2,
zz2,
xx3,
yy3,
zz3,
xx4,
yy4,
zz4,
xx5,
yy5,
zz5,
vxi,
vyi,
vzi,
xo1,
xo2,
xo3,
xo4,
xo5,
yo1,
yo2,
yo3,
yo4,
yo5,
zo1,
zo2,
zo3,
zo4,
zo5,
xi,
yi,
zi,
xoi,
yoi,
zoi,
zerom,
unp,
zeromt,
unpt )

Definition at line 1885 of file i25dst3_22.F.

1898C-----------------------------------------------
1899C I m p l i c i t T y p e s
1900C-----------------------------------------------
1901#include "implicit_f.inc"
1902C-----------------------------------------------
1903C D u m m y A r g u m e n t s
1904C-----------------------------------------------
1905 INTEGER IXX3 ,IXX4 , INTERSECA, INTERSECB
1906C REAL
1907 my_real
1908 1 v12 ,v23 ,v34 ,v41 ,
1909 2 vo12 ,vo23 ,vo34 ,vo41 ,xx1 ,
1910 3 yy1 ,zz1 ,xx2 ,yy2 ,zz2 ,
1911 4 xx3 ,yy3 ,zz3 ,xx4 ,yy4 ,
1912 5 zz4 ,xx5 ,yy5 ,zz5 ,
1913 6 vxi ,vyi ,vzi ,xo1 ,xo2 ,
1914 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
1915 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
1916 9 zo3 ,zo4 ,zo5 ,zerom ,unp ,
1917 a zeromt ,unpt ,xi ,yi ,zi ,
1918 b xoi ,yoi ,zoi
1919C-----------------------------------------------
1920C L o c a l V a r i a b l e s
1921C-----------------------------------------------
1922 my_real
1923 . x51, x52, x53, x54,
1924 . y51, y52, y53, y54,
1925 . z51, z52, z53, z54,
1926 . xi1, xi2, xi3, xi4, xi5,
1927 . yi1, yi2, yi3, yi4, yi5,
1928 . zi1, zi2, zi3, zi4, zi5,
1929 .
1930 .
1931 .
1932 . xs, ys, zs,
1933 . xm1, xm2, xm3, xm4, xm5,
1934 . ym1, ym2, ym3, ym4, ym5,
1935 . zm1, zm2, zm3, zm4, zm5,la0,lb0,lc0
1936 my_real
1937 . sx1,sx2,sx3,sx4,sx5,sx0,
1938 . sy1,sy2,sy3,sy4,sy5,sy0,
1939 . sz1,sz2,sz3,sz4,sz5,sz0,aaa,s2
1940C------------------------------------------------
1941 interseca = 0
1942 intersecb = 0
1943
1944
1945 IF(vo12 * v12 <= zero)THEN
1946 IF(abs(vo12) < em20)THEN
1947 aaa = one
1948 ELSEIF(abs(v12) < em20)THEN
1949 aaa = zero
1950 ELSE
1951 aaa = v12 / (v12-vo12)
1952 ENDIF
1953
1954 xm1 = xx1 + aaa*(xo1-xx1)
1955 ym1 = yy1 + aaa*(yo1-yy1)
1956 zm1 = zz1 + aaa*(zo1-zz1)
1957
1958 xm2 = xx2 + aaa*(xo2-xx2)
1959 ym2 = yy2 + aaa*(yo2-yy2)
1960 zm2 = zz2 + aaa*(zo2-zz2)
1961
1962 xs = xi + aaa*(xoi-xi)
1963 ys = yi + aaa*(yoi-yi)
1964 zs = zi + aaa*(zoi-zi)
1965
1966 xm5 = xx5 + aaa*(xo5-xx5)
1967 ym5 = yy5 + aaa*(yo5-yy5)
1968 zm5 = zz5 + aaa*(zo5-zz5)
1969
1970 x51 = xm1 - xm5
1971 y51 = ym1 - ym5
1972 z51 = zm1 - zm5
1973
1974 x52 = xm2 - xm5
1975 y52 = ym2 - ym5
1976 z52 = zm2 - zm5
1977
1978 xi1 = xm1 - xs
1979 yi1 = ym1 - ys
1980 zi1 = zm1 - zs
1981
1982 xi2 = xm2 - xs
1983 yi2 = ym2 - ys
1984 zi2 = zm2 - zs
1985
1986 xi5 = xm5 - xs
1987 yi5 = ym5 - ys
1988 zi5 = zm5 - zs
1989
1990 sx0 = y51*z52 - z51*y52
1991 sy0 = z51*x52 - x51*z52
1992 sz0 = x51*y52 - y51*x52
1993
1994 sx1 = yi5*zi1 - zi5*yi1
1995 sy1 = zi5*xi1 - xi5*zi1
1996 sz1 = xi5*yi1 - yi5*xi1
1997
1998 sx5 = yi1*zi2 - zi1*yi2
1999 sy5 = zi1*xi2 - xi1*zi2
2000 sz5 = xi1*yi2 - yi1*xi2
2001
2002 sx2 = yi2*zi5 - zi2*yi5
2003 sy2 = zi2*xi5 - xi2*zi5
2004 sz2 = xi2*yi5 - yi2*xi5
2005
2006 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2007 lb0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
2008 lc0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
2009 la0 = one-lb0-lc0
2010
2011 IF(ixx3 /= ixx4)THEN
2012 IF(lb0 >= zerom .and. lb0 <= unp .and.
2013 . lc0 >= zerom .and. lc0 <= unp .and.
2014 . la0 >= zerom .and. la0 <= unp)THEN
2015 IF(vo12 <= zero .AND. v12 >= zero)THEN
2016 interseca = 1
2017 ELSE !((VO12 > ZERO .AND. V12 < ZERO) .OR.
2018 . ! (VO12 > ZERO .AND. V12 == ZERO) .OR.
2019 . !(VO12 == ZERO .AND. V12 < ZERO)) THEN
2020 intersecb = 1
2021 END IF
2022 ENDIF
2023C----------loose tolerance for tria
2024 ELSE
2025 IF(lb0 >= zeromt .AND. lb0 <= unpt .AND.
2026 . lc0 >= zeromt .AND. lc0 <= unpt .AND.
2027 . la0 >= zeromt .AND. la0 <= unpt)THEN
2028 IF(vo12 <= zero .AND. v12 >= zero)THEN
2029 interseca = 1
2030 ELSE !((VO12 > ZERO .AND. V12 < ZERO) .OR.
2031 . !(VO12 > ZERO .AND. V12 == ZERO) .OR.
2032 . !(VO12 == ZERO .AND. V12 < ZERO)) THEN
2033 intersecb = 1
2034 ENDIF
2035 ENDIF
2036 END if! (IXX(I,3) /= IXX(I,4))THEN
2037 ENDIF
2038 IF(ixx3 /= ixx4)THEN
2039 IF(vo23 * v23 <= zero)THEN
2040 IF(abs(vo23) < em20)THEN
2041 aaa = one
2042 ELSEIF(abs(v23) < em20)THEN
2043 aaa = zero
2044 ELSE
2045 aaa = v23 / (v23-vo23)
2046 ENDIF
2047
2048 xm2 = xx2 + aaa*(xo2-xx2)
2049 ym2 = yy2 + aaa*(yo2-yy2)
2050 zm2 = zz2 + aaa*(zo2-zz2)
2051
2052 xm3 = xx3 + aaa*(xo3-xx3)
2053 ym3 = yy3 + aaa*(yo3-yy3)
2054 zm3 = zz3 + aaa*(zo3-zz3)
2055
2056 xs = xi + aaa*(xoi-xi)
2057 ys = yi + aaa*(yoi-yi)
2058 zs = zi + aaa*(zoi-zi)
2059
2060 xm5 = xx5 + aaa*(xo5-xx5)
2061 ym5 = yy5 + aaa*(yo5-yy5)
2062 zm5 = zz5 + aaa*(zo5-zz5)
2063
2064 x52 = xm2 - xm5
2065 y52 = ym2 - ym5
2066 z52 = zm2 - zm5
2067
2068 x53 = xm3 - xm5
2069 y53 = ym3 - ym5
2070 z53 = zm3 - zm5
2071
2072 xi2 = xm2 - xs
2073 yi2 = ym2 - ys
2074 zi2 = zm2 - zs
2075
2076 xi3 = xm3 - xs
2077 yi3 = ym3 - ys
2078 zi3 = zm3 - zs
2079
2080 xi5 = xm5 - xs
2081 yi5 = ym5 - ys
2082 zi5 = zm5 - zs
2083
2084 sx0 = y52*z53 - z52*y53
2085 sy0 = z52*x53 - x52*z53
2086 sz0 = x52*y53 - y52*x53
2087
2088 sx2 = yi5*zi2 - zi5*yi2
2089 sy2 = zi5*xi2 - xi5*zi2
2090 sz2 = xi5*yi2 - yi5*xi2
2091
2092 sx5 = yi2*zi3 - zi2*yi3
2093 sy5 = zi2*xi3 - xi2*zi3
2094 sz5 = xi2*yi3 - yi2*xi3
2095
2096 sx3 = yi3*zi5 - zi3*yi5
2097 sy3 = zi3*xi5 - xi3*zi5
2098 sz3 = xi3*yi5 - yi3*xi5
2099
2100 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2101 lb0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
2102 lc0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
2103 la0 = one-lb0-lc0
2104
2105 IF(lb0 >= zerom .and. lb0 <= unp .and.
2106 . lc0 >= zerom .and. lc0 <= unp .and.
2107 . la0 >= zerom .and. la0 <= unp)THEN
2108 IF(vo23 <= zero .AND. v23 >= zero)THEN
2109 interseca = 1
2110 ELSE !((VO23 > ZERO .AND. V23 < ZERO) .OR.
2111 . !(VO23 > ZERO .AND. V23 == ZERO) .OR.
2112 . !(VO23 == ZERO .AND. V23 < ZERO)) THEN
2113 intersecb = 1
2114 END IF
2115 ENDIF
2116 ENDIF
2117 IF(vo34 * v34 <= zero)THEN
2118 IF(abs(vo34) < em20)THEN
2119 aaa = one
2120 ELSEIF(abs(v34) < em20)THEN
2121 aaa = zero
2122 ELSE
2123 aaa = v34 / (v34-vo34)
2124 ENDIF
2125
2126 xm3 = xx3 + aaa*(xo3-xx3)
2127 ym3 = yy3 + aaa*(yo3-yy3)
2128 zm3 = zz3 + aaa*(zo3-zz3)
2129
2130 xm4 = xx4 + aaa*(xo4-xx4)
2131 ym4 = yy4 + aaa*(yo4-yy4)
2132 zm4 = zz4 + aaa*(zo4-zz4)
2133
2134 xs = xi + aaa*(xoi-xi)
2135 ys = yi + aaa*(yoi-yi)
2136 zs = zi + aaa*(zoi-zi)
2137
2138 xm5 = xx5 + aaa*(xo5-xx5)
2139 ym5 = yy5 + aaa*(yo5-yy5)
2140 zm5 = zz5 + aaa*(zo5-zz5)
2141
2142 x53 = xm3 - xm5
2143 y53 = ym3 - ym5
2144 z53 = zm3 - zm5
2145
2146 x54 = xm4 - xm5
2147 y54 = ym4 - ym5
2148 z54 = zm4 - zm5
2149
2150 xi3 = xm3 - xs
2151 yi3 = ym3 - ys
2152 zi3 = zm3 - zs
2153
2154 xi4 = xm4 - xs
2155 yi4 = ym4 - ys
2156 zi4 = zm4 - zs
2157
2158 xi5 = xm5 - xs
2159 yi5 = ym5 - ys
2160 zi5 = zm5 - zs
2161
2162 sx0 = y53*z54 - z53*y54
2163 sy0 = z53*x54 - x53*z54
2164 sz0 = x53*y54 - y53*x54
2165
2166 sx3 = yi5*zi3 - zi5*yi3
2167 sy3 = zi5*xi3 - xi5*zi3
2168 sz3 = xi5*yi3 - yi5*xi3
2169
2170 sx5 = yi3*zi4 - zi3*yi4
2171 sy5 = zi3*xi4 - xi3*zi4
2172 sz5 = xi3*yi4 - yi3*xi4
2173
2174 sx4 = yi4*zi5 - zi4*yi5
2175 sy4 = zi4*xi5 - xi4*zi5
2176 sz4 = xi4*yi5 - yi4*xi5
2177
2178 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2179 lb0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
2180 lc0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
2181 la0 = one-lb0-lc0
2182
2183 IF(lb0 >= zerom .and. lb0 <= unp .and.
2184 . lc0 >= zerom .and. lc0 <= unp .and.
2185 . la0 >= zerom .and. la0 <= unp)THEN
2186 IF(vo34 <= zero .AND. v34 >= zero)THEN
2187 interseca = 1
2188 ELSE !((VO34 > ZERO .AND. V34 < ZERO) .OR.
2189 . !(VO34 > ZERO .AND. V34 == ZERO) .OR.
2190 . !(VO34 == ZERO .AND. V34 < ZERO)) THEN
2191 intersecb = 1
2192 END IF
2193 ENDIF
2194 ENDIF
2195
2196 IF(vo41 * v41 <= zero)THEN
2197 IF(abs(vo41) < em20)THEN
2198 aaa = one
2199 ELSEIF(abs(v41) < em20)THEN
2200 aaa = zero
2201 ELSE
2202 aaa = v41 / (v41-vo41)
2203 ENDIF
2204
2205 xm4 = xx4 + aaa*(xo4-xx4)
2206 ym4 = yy4 + aaa*(yo4-yy4)
2207 zm4 = zz4 + aaa*(zo4-zz4)
2208
2209 xm1 = xx1 + aaa*(xo1-xx1)
2210 ym1 = yy1 + aaa*(yo1-yy1)
2211 zm1 = zz1 + aaa*(zo1-zz1)
2212
2213 xs = xi + aaa*(xoi-xi)
2214 ys = yi + aaa*(yoi-yi)
2215 zs = zi + aaa*(zoi-zi)
2216
2217 xm5 = xx5 + aaa*(xo5-xx5)
2218 ym5 = yy5 + aaa*(yo5-yy5)
2219 zm5 = zz5 + aaa*(zo5-zz5)
2220
2221 x54 = xm4 - xm5
2222 y54 = ym4 - ym5
2223 z54 = zm4 - zm5
2224
2225 x51 = xm1 - xm5
2226 y51 = ym1 - ym5
2227 z51 = zm1 - zm5
2228
2229 xi4 = xm4 - xs
2230 yi4 = ym4 - ys
2231 zi4 = zm4 - zs
2232
2233 xi1 = xm1 - xs
2234 yi1 = ym1 - ys
2235 zi1 = zm1 - zs
2236
2237 xi5 = xm5 - xs
2238 yi5 = ym5 - ys
2239 zi5 = zm5 - zs
2240
2241 sx0 = y54*z51 - z54*y51
2242 sy0 = z54*x51 - x54*z51
2243 sz0 = x54*y51 - y54*x51
2244
2245 sx4 = yi5*zi4 - zi5*yi4
2246 sy4 = zi5*xi4 - xi5*zi4
2247 sz4 = xi5*yi4 - yi5*xi4
2248
2249 sx5 = yi4*zi1 - zi4*yi1
2250 sy5 = zi4*xi1 - xi4*zi1
2251 sz5 = xi4*yi1 - yi4*xi1
2252
2253 sx1 = yi1*zi5 - zi1*yi5
2254 sy1 = zi1*xi5 - xi1*zi5
2255 sz1 = xi1*yi5 - yi1*xi5
2256
2257 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2258 lb0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
2259 lc0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
2260 la0 = one-lb0-lc0
2261
2262 IF(lb0 >= zerom .and. lb0 <= unp .and.
2263 . lc0 >= zerom .and. lc0 <= unp .and.
2264 . la0 >= zerom .and. la0 <= unp)THEN
2265 IF(vo41 <= zero .AND. v41 >= zero)THEN
2266 interseca = 1
2267 ELSE !((VO41 > ZERO .AND. V41 < ZERO) .OR.
2268 . !(VO41 > ZERO .AND. V41 == ZERO) .OR.
2269 . !(VO41 == ZERO .AND. V41 < ZERO)) THEN
2270 intersecb = 1
2271 END IF
2272 ENDIF
2273 ENDIF
2274 ENDIF
2275C
2276 RETURN

◆ intersecv0_25()

subroutine intersecv0_25 ( integer ixx3,
integer ixx4,
integer intersect,
v12,
v23,
v34,
v41,
vo12,
vo23,
vo34,
vo41,
xx1,
yy1,
zz1,
xx2,
yy2,
zz2,
xx3,
yy3,
zz3,
xx4,
yy4,
zz4,
xx5,
yy5,
zz5,
xoi,
yoi,
zoi,
xo1,
xo2,
xo3,
xo4,
xo5,
yo1,
yo2,
yo3,
yo4,
yo5,
zo1,
zo2,
zo3,
zo4,
zo5,
xi,
yi,
zi,
zeromt,
unpt )

Definition at line 2672 of file i25dst3_22.F.

2684C-----------------------------------------------
2685C I m p l i c i t T y p e s
2686C-----------------------------------------------
2687#include "implicit_f.inc"
2688C-----------------------------------------------
2689C D u m m y A r g u m e n t s
2690C-----------------------------------------------
2691 INTEGER IXX3 ,IXX4 ,INTERSECT
2692C REAL
2693 my_real
2694 1 v12 ,v23 ,v34 ,v41 ,
2695 2 vo12 ,vo23 ,vo34 ,vo41 ,xx1 ,
2696 3 yy1 ,zz1 ,xx2 ,yy2 ,zz2 ,
2697 4 xx3 ,yy3 ,zz3 ,xx4 ,yy4 ,
2698 5 zz4 ,xx5 ,yy5 ,zz5 ,
2699 6 xoi ,yoi ,zoi ,xo1 ,xo2 ,
2700 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
2701 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
2702 9 zo3 ,zo4 ,zo5 ,zeromt,unpt ,
2703 a xi ,yi ,zi
2704C-----------------------------------------------
2705C C o m m o n B l o c k s
2706C-----------------------------------------------
2707#include "com08_c.inc"
2708C-----------------------------------------------
2709C L o c a l V a r i a b l e s
2710C-----------------------------------------------
2711 my_real
2712 . x51, x52, x53, x54,
2713 . y51, y52, y53, y54,
2714 . z51, z52, z53, z54,
2715 . xi1, xi2, xi3, xi4, xi5,
2716 . yi1, yi2, yi3, yi4, yi5,
2717 . zi1, zi2, zi3, zi4, zi5,
2718 .
2719 .
2720 .
2721 .
2722 . xm1, xm2, xm3, xm4, xm5,
2723 . ym1, ym2, ym3, ym4, ym5,
2724 . zm1, zm2, zm3, zm4, zm5,la0,lb0,lc0
2725 my_real
2726 . sx1,sx2,sx3,sx4,sx5,sx0,
2727 . sy1,sy2,sy3,sy4,sy5,sy0,
2728 . sz1,sz2,sz3,sz4,sz5,sz0,aaa,adt,s2
2729C------------------------------------------------
2730C---------------- cas V>0 V_old>0
2731 IF(vo12 > zero .and. v12 >= zero )THEN
2732 aaa = one
2733 adt = dt1
2734
2735 xm1 = xo1
2736 ym1 = yo1
2737 zm1 = zo1
2738
2739 xm2 = xo2
2740 ym2 = yo2
2741 zm2 = zo2
2742
2743c XOI = XI - ADT*VXI
2744c YOI = YI - ADT*VYI
2745c ZOI = ZI - ADT*VZI
2746
2747 xm5 = xo5
2748 ym5 = yo5
2749 zm5 = zo5
2750
2751 x51 = xm1 - xm5
2752 y51 = ym1 - ym5
2753 z51 = zm1 - zm5
2754
2755 x52 = xm2 - xm5
2756 y52 = ym2 - ym5
2757 z52 = zm2 - zm5
2758
2759 xi1 = xm1 - xoi
2760 yi1 = ym1 - yoi
2761 zi1 = zm1 - zoi
2762
2763 xi2 = xm2 - xoi
2764 yi2 = ym2 - yoi
2765 zi2 = zm2 - zoi
2766
2767 xi5 = xm5 - xoi
2768 yi5 = ym5 - yoi
2769 zi5 = zm5 - zoi
2770
2771 sx0 = y51*z52 - z51*y52
2772 sy0 = z51*x52 - x51*z52
2773 sz0 = x51*y52 - y51*x52
2774
2775 sx1 = yi5*zi1 - zi5*yi1
2776 sy1 = zi5*xi1 - xi5*zi1
2777 sz1 = xi5*yi1 - yi5*xi1
2778
2779 sx5 = yi1*zi2 - zi1*yi2
2780 sy5 = zi1*xi2 - xi1*zi2
2781 sz5 = xi1*yi2 - yi1*xi2
2782
2783 sx2 = yi2*zi5 - zi2*yi5
2784 sy2 = zi2*xi5 - xi2*zi5
2785 sz2 = xi2*yi5 - yi2*xi5
2786
2787 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2788 lb0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
2789 lc0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
2790 la0 = one-lb0-lc0
2791
2792 IF(lb0 >= zeromt .AND. lb0 <= unpt .AND.
2793 . lc0 >= zeromt .AND. lc0 <= unpt .AND.
2794 . la0 >= zeromt .AND. la0 <= unpt)THEN
2795 intersect = 1
2796 ENDIF
2797 ENDIF
2798
2799 IF(ixx3 /= ixx4)THEN
2800 IF(vo23 > zero .AND. v23 >= zero)THEN
2801 aaa = one
2802 adt = dt1
2803
2804 xm2 = xo2
2805 ym2 = yo2
2806 zm2 = zo2
2807
2808 xm3 = xo3
2809 ym3 = yo3
2810 zm3 = zo3
2811
2812c XOI = XI - ADT*VXI
2813c YOI = YI - ADT*VYI
2814c ZOI = ZI - ADT*VZI
2815
2816 xm5 = xo5
2817 ym5 = yo5
2818 zm5 = zo5
2819
2820 x52 = xm2 - xm5
2821 y52 = ym2 - ym5
2822 z52 = zm2 - zm5
2823
2824 x53 = xm3 - xm5
2825 y53 = ym3 - ym5
2826 z53 = zm3 - zm5
2827
2828 xi2 = xm2 - xoi
2829 yi2 = ym2 - yoi
2830 zi2 = zm2 - zoi
2831
2832 xi3 = xm3 - xoi
2833 yi3 = ym3 - yoi
2834 zi3 = zm3 - zoi
2835
2836 xi5 = xm5 - xoi
2837 yi5 = ym5 - yoi
2838 zi5 = zm5 - zoi
2839
2840 sx0 = y52*z53 - z52*y53
2841 sy0 = z52*x53 - x52*z53
2842 sz0 = x52*y53 - y52*x53
2843
2844 sx2 = yi5*zi2 - zi5*yi2
2845 sy2 = zi5*xi2 - xi5*zi2
2846 sz2 = xi5*yi2 - yi5*xi2
2847
2848 sx5 = yi2*zi3 - zi2*yi3
2849 sy5 = zi2*xi3 - xi2*zi3
2850 sz5 = xi2*yi3 - yi2*xi3
2851
2852 sx3 = yi3*zi5 - zi3*yi5
2853 sy3 = zi3*xi5 - xi3*zi5
2854 sz3 = xi3*yi5 - yi3*xi5
2855
2856 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2857 lb0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
2858 lc0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
2859 la0 = one-lb0-lc0
2860
2861 IF(lb0 >= zeromt.and. lb0 <= unpt .and.
2862 . lc0 >= zeromt.and. lc0 <= unpt .and.
2863 . la0 >= zeromt.and. la0 <= unpt)THEN
2864 intersect = 1
2865 ENDIF
2866 ENDIF
2867 IF(vo34 > zero .and. v34 >= zero )THEN
2868 aaa = one
2869 adt = dt1
2870
2871 xm3 = xo3
2872 ym3 = yo3
2873 zm3 = zo3
2874
2875 xm4 = xo4
2876 ym4 = yo4
2877 zm4 = zo4
2878
2879c XOI = XI - ADT*VXI
2880c YOI = YI - ADT*VYI
2881c ZOI = ZI - ADT*VZI
2882
2883 xm5 = xo5
2884 ym5 = yo5
2885 zm5 = zo5
2886
2887 x53 = xm3 - xm5
2888 y53 = ym3 - ym5
2889 z53 = zm3 - zm5
2890
2891 x54 = xm4 - xm5
2892 y54 = ym4 - ym5
2893 z54 = zm4 - zm5
2894
2895 xi3 = xm3 - xoi
2896 yi3 = ym3 - yoi
2897 zi3 = zm3 - zoi
2898
2899 xi4 = xm4 - xoi
2900 yi4 = ym4 - yoi
2901 zi4 = zm4 - zoi
2902
2903 xi5 = xm5 - xoi
2904 yi5 = ym5 - yoi
2905 zi5 = zm5 - zoi
2906
2907 sx0 = y53*z54 - z53*y54
2908 sy0 = z53*x54 - x53*z54
2909 sz0 = x53*y54 - y53*x54
2910
2911 sx3 = yi5*zi3 - zi5*yi3
2912 sy3 = zi5*xi3 - xi5*zi3
2913 sz3 = xi5*yi3 - yi5*xi3
2914
2915 sx5 = yi3*zi4 - zi3*yi4
2916 sy5 = zi3*xi4 - xi3*zi4
2917 sz5 = xi3*yi4 - yi3*xi4
2918
2919 sx4 = yi4*zi5 - zi4*yi5
2920 sy4 = zi4*xi5 - xi4*zi5
2921 sz4 = xi4*yi5 - yi4*xi5
2922
2923 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2924 lb0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
2925 lc0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
2926 la0 = one-lb0-lc0
2927
2928 IF(lb0 >= zeromt.and. lb0 <= unpt .and.
2929 . lc0 >= zeromt.and. lc0 <= unpt .and.
2930 . la0 >= zeromt.and. la0 <= unpt)THEN
2931 intersect = 1
2932 ENDIF
2933 ENDIF
2934
2935 IF(vo41 > zero .and. v41 >= zero )THEN
2936 aaa = one
2937 adt = dt1
2938
2939 xm4 = xo4
2940 ym4 = yo4
2941 zm4 = zo4
2942
2943 xm1 = xo1
2944 ym1 = yo1
2945 zm1 = zo1
2946
2947c XOI = XI - ADT*VXI
2948c YOI = YI - ADT*VYI
2949c ZOI = ZI - ADT*VZI
2950
2951 xm5 = xo5
2952 ym5 = yo5
2953 zm5 = zo5
2954
2955 x54 = xm4 - xm5
2956 y54 = ym4 - ym5
2957 z54 = zm4 - zm5
2958
2959 x51 = xm1 - xm5
2960 y51 = ym1 - ym5
2961 z51 = zm1 - zm5
2962
2963 xi4 = xm4 - xoi
2964 yi4 = ym4 - yoi
2965 zi4 = zm4 - zoi
2966
2967 xi1 = xm1 - xoi
2968 yi1 = ym1 - yoi
2969 zi1 = zm1 - zoi
2970
2971 xi5 = xm5 - xoi
2972 yi5 = ym5 - yoi
2973 zi5 = zm5 - zoi
2974
2975 sx0 = y54*z51 - z54*y51
2976 sy0 = z54*x51 - x54*z51
2977 sz0 = x54*y51 - y54*x51
2978
2979 sx4 = yi5*zi4 - zi5*yi4
2980 sy4 = zi5*xi4 - xi5*zi4
2981 sz4 = xi5*yi4 - yi5*xi4
2982
2983 sx5 = yi4*zi1 - zi4*yi1
2984 sy5 = zi4*xi1 - xi4*zi1
2985 sz5 = xi4*yi1 - yi4*xi1
2986
2987 sx1 = yi1*zi5 - zi1*yi5
2988 sy1 = zi1*xi5 - xi1*zi5
2989 sz1 = xi1*yi5 - yi1*xi5
2990
2991 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2992 lb0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
2993 lc0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
2994 la0 = one-lb0-lc0
2995
2996 IF(lb0 >= zeromt.and. lb0 <= unpt .and.
2997 . lc0 >= zeromt.and. lc0 <= unpt .and.
2998 . la0 >= zeromt.and. la0 <= unpt)THEN
2999 intersect = 1
3000 ENDIF
3001 ENDIF
3002 ENDIF
3003C
3004 RETURN
intersect
static2 void intersect(char *uplo, char *diag, Int j, Int start, Int end, Int action, Int *ptrsizebuff, complex **pptrbuff, complex *ptrblock, Int m, Int n, MDESC *ma, Int ia, Int ja, Int templateheight0, Int templatewidth0, MDESC *mb, Int ib, Int jb, Int templateheight1, Int templatewidth1)
Definition pctrmr2.c:176

◆ intersecv_25()

subroutine intersecv_25 ( integer ixx3,
integer ixx4,
integer subtria,
integer intersect,
v12,
v23,
v34,
v41,
vo12,
vo23,
vo34,
vo41,
xx1,
yy1,
zz1,
xx2,
yy2,
zz2,
xx3,
yy3,
zz3,
xx4,
yy4,
zz4,
xx5,
yy5,
zz5,
vxi,
vyi,
vzi,
xo1,
xo2,
xo3,
xo4,
xo5,
yo1,
yo2,
yo3,
yo4,
yo5,
zo1,
zo2,
zo3,
zo4,
zo5,
xi,
yi,
zi,
xoi,
yoi,
zoi,
zerom,
unp,
zeromt,
unpt )

Definition at line 2281 of file i25dst3_22.F.

2294C-----------------------------------------------
2295C I m p l i c i t T y p e s
2296C-----------------------------------------------
2297#include "implicit_f.inc"
2298C-----------------------------------------------
2299C D u m m y A r g u m e n t s
2300C-----------------------------------------------
2301 INTEGER IXX3 ,IXX4 ,INTERSECT, SUBTRIA
2302C REAL
2303 my_real
2304 1 v12 ,v23 ,v34 ,v41 ,
2305 2 vo12 ,vo23 ,vo34 ,vo41 ,xx1 ,
2306 3 yy1 ,zz1 ,xx2 ,yy2 ,zz2 ,
2307 4 xx3 ,yy3 ,zz3 ,xx4 ,yy4 ,
2308 5 zz4 ,xx5 ,yy5 ,zz5 ,
2309 6 vxi ,vyi ,vzi ,xo1 ,xo2 ,
2310 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
2311 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
2312 9 zo3 ,zo4 ,zo5 ,zerom ,unp ,
2313 a zeromt ,unpt ,xi ,yi ,zi ,
2314 b xoi ,yoi ,zoi
2315C-----------------------------------------------
2316C L o c a l V a r i a b l e s
2317C-----------------------------------------------
2318 my_real
2319 . x51, x52, x53, x54,
2320 . y51, y52, y53, y54,
2321 . z51, z52, z53, z54,
2322 . xi1, xi2, xi3, xi4, xi5,
2323 . yi1, yi2, yi3, yi4, yi5,
2324 . zi1, zi2, zi3, zi4, zi5,
2325 .
2326 .
2327 .
2328 . xs, ys, zs,
2329 . xm1, xm2, xm3, xm4, xm5,
2330 . ym1, ym2, ym3, ym4, ym5,
2331 . zm1, zm2, zm3, zm4, zm5,la0,lb0,lc0
2332 my_real
2333 . sx1,sx2,sx3,sx4,sx5,sx0,
2334 . sy1,sy2,sy3,sy4,sy5,sy0,
2335 . sz1,sz2,sz3,sz4,sz5,sz0,aaa,s2
2336C------------------------------------------------
2337 IF(vo12 <= zero .AND. v12 >= zero)THEN
2338 IF(abs(vo12) < em20)THEN
2339 aaa = one
2340 ELSEIF(abs(v12) < em20)THEN
2341 aaa = zero
2342 ELSE
2343 aaa = v12 / (v12-vo12)
2344 ENDIF
2345
2346 xm1 = xx1 + aaa*(xo1-xx1)
2347 ym1 = yy1 + aaa*(yo1-yy1)
2348 zm1 = zz1 + aaa*(zo1-zz1)
2349
2350 xm2 = xx2 + aaa*(xo2-xx2)
2351 ym2 = yy2 + aaa*(yo2-yy2)
2352 zm2 = zz2 + aaa*(zo2-zz2)
2353
2354 xs = xi + aaa*(xoi-xi)
2355 ys = yi + aaa*(yoi-yi)
2356 zs = zi + aaa*(zoi-zi)
2357
2358 xm5 = xx5 + aaa*(xo5-xx5)
2359 ym5 = yy5 + aaa*(yo5-yy5)
2360 zm5 = zz5 + aaa*(zo5-zz5)
2361
2362 x51 = xm1 - xm5
2363 y51 = ym1 - ym5
2364 z51 = zm1 - zm5
2365
2366 x52 = xm2 - xm5
2367 y52 = ym2 - ym5
2368 z52 = zm2 - zm5
2369
2370 xi1 = xm1 - xs
2371 yi1 = ym1 - ys
2372 zi1 = zm1 - zs
2373
2374 xi2 = xm2 - xs
2375 yi2 = ym2 - ys
2376 zi2 = zm2 - zs
2377
2378 xi5 = xm5 - xs
2379 yi5 = ym5 - ys
2380 zi5 = zm5 - zs
2381
2382 sx0 = y51*z52 - z51*y52
2383 sy0 = z51*x52 - x51*z52
2384 sz0 = x51*y52 - y51*x52
2385
2386 sx1 = yi5*zi1 - zi5*yi1
2387 sy1 = zi5*xi1 - xi5*zi1
2388 sz1 = xi5*yi1 - yi5*xi1
2389
2390 sx5 = yi1*zi2 - zi1*yi2
2391 sy5 = zi1*xi2 - xi1*zi2
2392 sz5 = xi1*yi2 - yi1*xi2
2393
2394 sx2 = yi2*zi5 - zi2*yi5
2395 sy2 = zi2*xi5 - xi2*zi5
2396 sz2 = xi2*yi5 - yi2*xi5
2397
2398 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2399 lb0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
2400 lc0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
2401 la0 = one-lb0-lc0
2402
2403 IF(ixx3 /= ixx4)THEN
2404 IF(lb0 >= zerom .and. lb0 <= unp .and.
2405 . lc0 >= zerom .and. lc0 <= unp .and.
2406 . la0 >= zerom .and. la0 <= unp)THEN
2407 IF(vo12 <= zero .AND. v12 >= zero)THEN
2408 intersect = 1
2409 ELSE ! (VO12 > ZERO .AND. V12 < ZERO) .OR.
2410 ! (VO12 > ZERO .AND. V12 == ZERO) .OR.
2411 ! (VO12 == ZERO .AND. V12 < ZERO)
2412 intersect = -1
2413 END IF
2414 subtria= 1 ! subtriangle
2415 ENDIF
2416C----------loose tolerance for tria
2417 ELSE
2418 IF(lb0 >= zeromt .AND. lb0 <= unpt .AND.
2419 . lc0 >= zeromt .AND. lc0 <= unpt .AND.
2420 . la0 >= zeromt .AND. la0 <= unpt)THEN
2421 intersect = 1
2422 subtria= 1 ! subtriangle
2423 ENDIF
2424 END if! (IXX(I,3) /= IXX(I,4))THEN
2425 ENDIF
2426 IF(ixx3 /= ixx4)THEN
2427 IF(vo23 <= zero .AND. v23 >= zero)THEN
2428 IF(abs(vo23) < em20)THEN
2429 aaa = one
2430 ELSEIF(abs(v23) < em20)THEN
2431 aaa = zero
2432 ELSE
2433 aaa = v23 / (v23-vo23)
2434 ENDIF
2435
2436 xm2 = xx2 + aaa*(xo2-xx2)
2437 ym2 = yy2 + aaa*(yo2-yy2)
2438 zm2 = zz2 + aaa*(zo2-zz2)
2439
2440 xm3 = xx3 + aaa*(xo3-xx3)
2441 ym3 = yy3 + aaa*(yo3-yy3)
2442 zm3 = zz3 + aaa*(zo3-zz3)
2443
2444 xs = xi + aaa*(xoi-xi)
2445 ys = yi + aaa*(yoi-yi)
2446 zs = zi + aaa*(zoi-zi)
2447
2448 xm5 = xx5 + aaa*(xo5-xx5)
2449 ym5 = yy5 + aaa*(yo5-yy5)
2450 zm5 = zz5 + aaa*(zo5-zz5)
2451
2452 x52 = xm2 - xm5
2453 y52 = ym2 - ym5
2454 z52 = zm2 - zm5
2455
2456 x53 = xm3 - xm5
2457 y53 = ym3 - ym5
2458 z53 = zm3 - zm5
2459
2460 xi2 = xm2 - xs
2461 yi2 = ym2 - ys
2462 zi2 = zm2 - zs
2463
2464 xi3 = xm3 - xs
2465 yi3 = ym3 - ys
2466 zi3 = zm3 - zs
2467
2468 xi5 = xm5 - xs
2469 yi5 = ym5 - ys
2470 zi5 = zm5 - zs
2471
2472 sx0 = y52*z53 - z52*y53
2473 sy0 = z52*x53 - x52*z53
2474 sz0 = x52*y53 - y52*x53
2475
2476 sx2 = yi5*zi2 - zi5*yi2
2477 sy2 = zi5*xi2 - xi5*zi2
2478 sz2 = xi5*yi2 - yi5*xi2
2479
2480 sx5 = yi2*zi3 - zi2*yi3
2481 sy5 = zi2*xi3 - xi2*zi3
2482 sz5 = xi2*yi3 - yi2*xi3
2483
2484 sx3 = yi3*zi5 - zi3*yi5
2485 sy3 = zi3*xi5 - xi3*zi5
2486 sz3 = xi3*yi5 - yi3*xi5
2487
2488 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2489 lb0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
2490 lc0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
2491 la0 = one-lb0-lc0
2492
2493 IF(lb0 >= zerom .and. lb0 <= unp .and.
2494 . lc0 >= zerom .and. lc0 <= unp .and.
2495 . la0 >= zerom .and. la0 <= unp)THEN
2496 IF(vo23 <= zero .AND. v23 >= zero)THEN
2497 intersect = 1
2498 ELSE ! (VO23 > ZERO .AND. V23 < ZERO) .OR.
2499 ! (VO23 > ZERO .AND. V23 == ZERO) .OR.
2500 ! (VO23 == ZERO .AND. V23 < ZERO)
2501 intersect = -1
2502 END IF
2503
2504 ENDIF
2505 ENDIF
2506 IF(vo34 <= zero .AND. v34 >= zero)THEN
2507 IF(abs(vo34) < em20)THEN
2508 aaa = one
2509 ELSEIF(abs(v34) < em20)THEN
2510 aaa = zero
2511 ELSE
2512 aaa = v34 / (v34-vo34)
2513 ENDIF
2514
2515 xm3 = xx3 + aaa*(xo3-xx3)
2516 ym3 = yy3 + aaa*(yo3-yy3)
2517 zm3 = zz3 + aaa*(zo3-zz3)
2518
2519 xm4 = xx4 + aaa*(xo4-xx4)
2520 ym4 = yy4 + aaa*(yo4-yy4)
2521 zm4 = zz4 + aaa*(zo4-zz4)
2522
2523 xs = xi + aaa*(xoi-xi)
2524 ys = yi + aaa*(yoi-yi)
2525 zs = zi + aaa*(zoi-zi)
2526
2527 xm5 = xx5 + aaa*(xo5-xx5)
2528 ym5 = yy5 + aaa*(yo5-yy5)
2529 zm5 = zz5 + aaa*(zo5-zz5)
2530
2531 x53 = xm3 - xm5
2532 y53 = ym3 - ym5
2533 z53 = zm3 - zm5
2534
2535 x54 = xm4 - xm5
2536 y54 = ym4 - ym5
2537 z54 = zm4 - zm5
2538
2539 xi3 = xm3 - xs
2540 yi3 = ym3 - ys
2541 zi3 = zm3 - zs
2542
2543 xi4 = xm4 - xs
2544 yi4 = ym4 - ys
2545 zi4 = zm4 - zs
2546
2547 xi5 = xm5 - xs
2548 yi5 = ym5 - ys
2549 zi5 = zm5 - zs
2550
2551 sx0 = y53*z54 - z53*y54
2552 sy0 = z53*x54 - x53*z54
2553 sz0 = x53*y54 - y53*x54
2554
2555 sx3 = yi5*zi3 - zi5*yi3
2556 sy3 = zi5*xi3 - xi5*zi3
2557 sz3 = xi5*yi3 - yi5*xi3
2558
2559 sx5 = yi3*zi4 - zi3*yi4
2560 sy5 = zi3*xi4 - xi3*zi4
2561 sz5 = xi3*yi4 - yi3*xi4
2562
2563 sx4 = yi4*zi5 - zi4*yi5
2564 sy4 = zi4*xi5 - xi4*zi5
2565 sz4 = xi4*yi5 - yi4*xi5
2566
2567 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2568 lb0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
2569 lc0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
2570 la0 = one-lb0-lc0
2571
2572 IF(lb0 >= zerom .and. lb0 <= unp .and.
2573 . lc0 >= zerom .and. lc0 <= unp .and.
2574 . la0 >= zerom .and. la0 <= unp)THEN
2575 IF(vo34 <= zero .AND. v34 >= zero)THEN
2576 intersect = 1
2577 ELSE ! (VO34 > ZERO .AND. V34 < ZERO) .OR.
2578 ! (VO34 > ZERO .AND. V34 == ZERO) .OR.
2579 ! (VO34 == ZERO .AND. V34 < ZERO)
2580 intersect = -1
2581 END IF
2582 ENDIF
2583 ENDIF
2584
2585 IF(vo41 <= zero .AND. v41 >= zero)THEN
2586 IF(abs(vo41) < em20)THEN
2587 aaa = one
2588 ELSEIF(abs(v41) < em20)THEN
2589 aaa = zero
2590 ELSE
2591 aaa = v41 / (v41-vo41)
2592 ENDIF
2593
2594 xm4 = xx4 + aaa*(xo4-xx4)
2595 ym4 = yy4 + aaa*(yo4-yy4)
2596 zm4 = zz4 + aaa*(zo4-zz4)
2597
2598 xm1 = xx1 + aaa*(xo1-xx1)
2599 ym1 = yy1 + aaa*(yo1-yy1)
2600 zm1 = zz1 + aaa*(zo1-zz1)
2601
2602 xs = xi + aaa*(xoi-xi)
2603 ys = yi + aaa*(yoi-yi)
2604 zs = zi + aaa*(zoi-zi)
2605
2606 xm5 = xx5 + aaa*(xo5-xx5)
2607 ym5 = yy5 + aaa*(yo5-yy5)
2608 zm5 = zz5 + aaa*(zo5-zz5)
2609
2610 x54 = xm4 - xm5
2611 y54 = ym4 - ym5
2612 z54 = zm4 - zm5
2613
2614 x51 = xm1 - xm5
2615 y51 = ym1 - ym5
2616 z51 = zm1 - zm5
2617
2618 xi4 = xm4 - xs
2619 yi4 = ym4 - ys
2620 zi4 = zm4 - zs
2621
2622 xi1 = xm1 - xs
2623 yi1 = ym1 - ys
2624 zi1 = zm1 - zs
2625
2626 xi5 = xm5 - xs
2627 yi5 = ym5 - ys
2628 zi5 = zm5 - zs
2629
2630 sx0 = y54*z51 - z54*y51
2631 sy0 = z54*x51 - x54*z51
2632 sz0 = x54*y51 - y54*x51
2633
2634 sx4 = yi5*zi4 - zi5*yi4
2635 sy4 = zi5*xi4 - xi5*zi4
2636 sz4 = xi5*yi4 - yi5*xi4
2637
2638 sx5 = yi4*zi1 - zi4*yi1
2639 sy5 = zi4*xi1 - xi4*zi1
2640 sz5 = xi4*yi1 - yi4*xi1
2641
2642 sx1 = yi1*zi5 - zi1*yi5
2643 sy1 = zi1*xi5 - xi1*zi5
2644 sz1 = xi1*yi5 - yi1*xi5
2645
2646 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2647 lb0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
2648 lc0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
2649 la0 = one-lb0-lc0
2650
2651 IF(lb0 >= zerom .and. lb0 <= unp .and.
2652 . lc0 >= zerom .and. lc0 <= unp .and.
2653 . la0 >= zerom .and. la0 <= unp)THEN
2654 IF(vo41 <= zero .AND. v41 >= zero)THEN
2655 intersect = 1
2656 ELSE ! (VO41 > ZERO .AND. V41 < ZERO) .OR.
2657 ! (VO41 > ZERO .AND. V41 == ZERO) .OR.
2658 ! (VO41 == ZERO .AND. V41 < ZERO)
2659 intersect = -1
2660 END IF
2661 ENDIF
2662 ENDIF
2663 ENDIF
2664C
2665 RETURN