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

2895C-----------------------------------------------
2896C M o d u l e s
2897C-----------------------------------------------
2898 USE tri7box
2899C-----------------------------------------------
2900C I m p l i c i t T y p e s
2901C-----------------------------------------------
2902#include "implicit_f.inc"
2903C-----------------------------------------------
2904C G l o b a l P a r a m e t e r s
2905C-----------------------------------------------
2906#include "mvsiz_p.inc"
2907#include "comlock.inc"
2908C-----------------------------------------------
2909C D u m m y A r g u m e n t s
2910C-----------------------------------------------
2911 INTEGER JLT, NIN, NSN, INACTI,
2912 . CAND_N_N(*), CAND_E_N(*), CAND_E(*), ITAB(*), SUBTRIA(MVSIZ)
2913 INTEGER NSVG(MVSIZ), IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
2914 . MSEGLO(*), IRTLM(4,NSN) ,INDEX(*), FARM(4,*),
2915 . FAR(MVSIZ)
2916 my_real
2917 . stif(*), time_s(*),
2918 . pent(mvsiz),
2919 . lb(mvsiz), lc(mvsiz),
2920 . penm(4,*), lbm(4,*), lcm(4,*)
2921C-----------------------------------------------
2922C L o c a l V a r i a b l e s
2923C-----------------------------------------------
2924 INTEGER IT, I, J, L, N, MGLOB
2925C--------------------------------------------------------
2926C Back to global arrays
2927C--------------------------------------------------------
2928 DO i=1,jlt
2929 j=index(i)
2930 it=subtria(i)
2931
2932 farm(1:4,j)=0
2933 penm(1:4,j)=zero
2934 lbm(1:4,j) =zero
2935 lcm(1:4,j) =zero
2936
2937 IF(it/=0)THEN
2938 cand_e(j)=cand_e_n(i)
2939
2940 farm(it,j)=far(i)
2941 penm(it,j)=pent(i)
2942 lbm(it,j) =lb(i)
2943 lcm(it,j) =lc(i)
2944 END IF
2945
2946 END DO
2947 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 1394 of file i25dst3_22.F.

1407C-----------------------------------------------
1408C I m p l i c i t T y p e s
1409C-----------------------------------------------
1410#include "implicit_f.inc"
1411C-----------------------------------------------
1412C D u m m y A r g u m e n t s
1413C-----------------------------------------------
1414 INTEGER IXX3 ,IXX4 , INTERSECA
1415C REAL
1416 my_real
1417 1 v12 ,v23 ,v34 ,v41 ,
1418 2 vo12 ,vo23 ,vo34 ,vo41 ,xx1 ,
1419 3 yy1 ,zz1 ,xx2 ,yy2 ,zz2 ,
1420 4 xx3 ,yy3 ,zz3 ,xx4 ,yy4 ,
1421 5 zz4 ,xx5 ,yy5 ,zz5 ,
1422 6 vxi ,vyi ,vzi ,xo1 ,xo2 ,
1423 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
1424 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
1425 9 zo3 ,zo4 ,zo5 ,zerom ,unp ,
1426 a zeromt ,unpt ,xi ,yi ,zi ,
1427 b xoi ,yoi ,zoi
1428C-----------------------------------------------
1429C L o c a l V a r i a b l e s
1430C-----------------------------------------------
1431 my_real
1432 . x51, x52, x53, x54,
1433 . y51, y52, y53, y54,
1434 . z51, z52, z53, z54,
1435 . xi1, xi2, xi3, xi4, xi5,
1436 . yi1, yi2, yi3, yi4, yi5,
1437 . zi1, zi2, zi3, zi4, zi5,
1438 . xia, xib, xic,
1439 . yia, yib, yic,
1440 . zia, zib, zic,
1441 . xs, ys, zs,
1442 . xm1, xm2, xm3, xm4, xm5,
1443 . ym1, ym2, ym3, ym4, ym5,
1444 . zm1, zm2, zm3, zm4, zm5,la0,lb0,lc0
1445 my_real
1446 . sx,sx1,sx2,sx3,sx4,sx5,sx0,
1447 . sy,sy1,sy2,sy3,sy4,sy5,sy0,
1448 . sz,sz1,sz2,sz3,sz4,sz5,sz0,aaa,s2
1449C------------------------------------------------
1450 interseca = 0
1451
1452 IF(vo12 <= zero .AND. v12 >= zero)THEN
1453 IF(abs(vo12) < em20)THEN
1454 aaa = one
1455 ELSEIF(abs(v12) < em20)THEN
1456 aaa = zero
1457 ELSE
1458 aaa = v12 / (v12-vo12)
1459 ENDIF
1460
1461 xm1 = xx1 + aaa*(xo1-xx1)
1462 ym1 = yy1 + aaa*(yo1-yy1)
1463 zm1 = zz1 + aaa*(zo1-zz1)
1464
1465 xm2 = xx2 + aaa*(xo2-xx2)
1466 ym2 = yy2 + aaa*(yo2-yy2)
1467 zm2 = zz2 + aaa*(zo2-zz2)
1468
1469 xs = xi + aaa*(xoi-xi)
1470 ys = yi + aaa*(yoi-yi)
1471 zs = zi + aaa*(zoi-zi)
1472
1473 xm5 = xx5 + aaa*(xo5-xx5)
1474 ym5 = yy5 + aaa*(yo5-yy5)
1475 zm5 = zz5 + aaa*(zo5-zz5)
1476
1477 x51 = xm1 - xm5
1478 y51 = ym1 - ym5
1479 z51 = zm1 - zm5
1480
1481 x52 = xm2 - xm5
1482 y52 = ym2 - ym5
1483 z52 = zm2 - zm5
1484
1485 xi1 = xm1 - xs
1486 yi1 = ym1 - ys
1487 zi1 = zm1 - zs
1488
1489 xi2 = xm2 - xs
1490 yi2 = ym2 - ys
1491 zi2 = zm2 - zs
1492
1493 xi5 = xm5 - xs
1494 yi5 = ym5 - ys
1495 zi5 = zm5 - zs
1496
1497 sx0 = y51*z52 - z51*y52
1498 sy0 = z51*x52 - x51*z52
1499 sz0 = x51*y52 - y51*x52
1500
1501 sx1 = yi5*zi1 - zi5*yi1
1502 sy1 = zi5*xi1 - xi5*zi1
1503 sz1 = xi5*yi1 - yi5*xi1
1504
1505 sx5 = yi1*zi2 - zi1*yi2
1506 sy5 = zi1*xi2 - xi1*zi2
1507 sz5 = xi1*yi2 - yi1*xi2
1508
1509 sx2 = yi2*zi5 - zi2*yi5
1510 sy2 = zi2*xi5 - xi2*zi5
1511 sz2 = xi2*yi5 - yi2*xi5
1512
1513 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
1514 lb0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
1515 lc0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
1516 la0 = one-lb0-lc0
1517
1518 IF(ixx3 /= ixx4)THEN
1519 IF(lb0 >= zerom .and. lb0 <= unp .and.
1520 . lc0 >= zerom .and. lc0 <= unp .and.
1521 . la0 >= zerom .and. la0 <= unp)THEN
1522 interseca = 1
1523 ENDIF
1524C----------loose tolerance for tria
1525 ELSE
1526 IF(lb0 >= zeromt .AND. lb0 <= unpt .AND.
1527 . lc0 >= zeromt .AND. lc0 <= unpt .AND.
1528 . la0 >= zeromt .AND. la0 <= unpt)THEN
1529 interseca = 1
1530 ENDIF
1531 END if! (IXX(I,3) /= IXX(I,4))THEN
1532 ENDIF
1533 IF(ixx3 /= ixx4)THEN
1534 IF(vo23 <= zero .AND. v23 >= zero)THEN
1535 IF(abs(vo23) < em20)THEN
1536 aaa = one
1537 ELSEIF(abs(v23) < em20)THEN
1538 aaa = zero
1539 ELSE
1540 aaa = v23 / (v23-vo23)
1541 ENDIF
1542
1543 xm2 = xx2 + aaa*(xo2-xx2)
1544 ym2 = yy2 + aaa*(yo2-yy2)
1545 zm2 = zz2 + aaa*(zo2-zz2)
1546
1547 xm3 = xx3 + aaa*(xo3-xx3)
1548 ym3 = yy3 + aaa*(yo3-yy3)
1549 zm3 = zz3 + aaa*(zo3-zz3)
1550
1551 xs = xi + aaa*(xoi-xi)
1552 ys = yi + aaa*(yoi-yi)
1553 zs = zi + aaa*(zoi-zi)
1554
1555 xm5 = xx5 + aaa*(xo5-xx5)
1556 ym5 = yy5 + aaa*(yo5-yy5)
1557 zm5 = zz5 + aaa*(zo5-zz5)
1558
1559 x52 = xm2 - xm5
1560 y52 = ym2 - ym5
1561 z52 = zm2 - zm5
1562
1563 x53 = xm3 - xm5
1564 y53 = ym3 - ym5
1565 z53 = zm3 - zm5
1566
1567 xi2 = xm2 - xs
1568 yi2 = ym2 - ys
1569 zi2 = zm2 - zs
1570
1571 xi3 = xm3 - xs
1572 yi3 = ym3 - ys
1573 zi3 = zm3 - zs
1574
1575 xi5 = xm5 - xs
1576 yi5 = ym5 - ys
1577 zi5 = zm5 - zs
1578
1579 sx0 = y52*z53 - z52*y53
1580 sy0 = z52*x53 - x52*z53
1581 sz0 = x52*y53 - y52*x53
1582
1583 sx2 = yi5*zi2 - zi5*yi2
1584 sy2 = zi5*xi2 - xi5*zi2
1585 sz2 = xi5*yi2 - yi5*xi2
1586
1587 sx5 = yi2*zi3 - zi2*yi3
1588 sy5 = zi2*xi3 - xi2*zi3
1589 sz5 = xi2*yi3 - yi2*xi3
1590
1591 sx3 = yi3*zi5 - zi3*yi5
1592 sy3 = zi3*xi5 - xi3*zi5
1593 sz3 = xi3*yi5 - yi3*xi5
1594
1595 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
1596 lb0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
1597 lc0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
1598 la0 = one-lb0-lc0
1599
1600 IF(lb0 >= zerom .and. lb0 <= unp .and.
1601 . lc0 >= zerom .and. lc0 <= unp .and.
1602 . la0 >= zerom .and. la0 <= unp)THEN
1603 interseca = 1
1604 ENDIF
1605 ENDIF
1606 IF(vo34 <= zero .AND. v34 >= zero)THEN
1607 IF(abs(vo34) < em20)THEN
1608 aaa = one
1609 ELSEIF(abs(v34) < em20)THEN
1610 aaa = zero
1611 ELSE
1612 aaa = v34 / (v34-vo34)
1613 ENDIF
1614
1615 xm3 = xx3 + aaa*(xo3-xx3)
1616 ym3 = yy3 + aaa*(yo3-yy3)
1617 zm3 = zz3 + aaa*(zo3-zz3)
1618
1619 xm4 = xx4 + aaa*(xo4-xx4)
1620 ym4 = yy4 + aaa*(yo4-yy4)
1621 zm4 = zz4 + aaa*(zo4-zz4)
1622
1623 xs = xi + aaa*(xoi-xi)
1624 ys = yi + aaa*(yoi-yi)
1625 zs = zi + aaa*(zoi-zi)
1626
1627 xm5 = xx5 + aaa*(xo5-xx5)
1628 ym5 = yy5 + aaa*(yo5-yy5)
1629 zm5 = zz5 + aaa*(zo5-zz5)
1630
1631 x53 = xm3 - xm5
1632 y53 = ym3 - ym5
1633 z53 = zm3 - zm5
1634
1635 x54 = xm4 - xm5
1636 y54 = ym4 - ym5
1637 z54 = zm4 - zm5
1638
1639 xi3 = xm3 - xs
1640 yi3 = ym3 - ys
1641 zi3 = zm3 - zs
1642
1643 xi4 = xm4 - xs
1644 yi4 = ym4 - ys
1645 zi4 = zm4 - zs
1646
1647 xi5 = xm5 - xs
1648 yi5 = ym5 - ys
1649 zi5 = zm5 - zs
1650
1651 sx0 = y53*z54 - z53*y54
1652 sy0 = z53*x54 - x53*z54
1653 sz0 = x53*y54 - y53*x54
1654
1655 sx3 = yi5*zi3 - zi5*yi3
1656 sy3 = zi5*xi3 - xi5*zi3
1657 sz3 = xi5*yi3 - yi5*xi3
1658
1659 sx5 = yi3*zi4 - zi3*yi4
1660 sy5 = zi3*xi4 - xi3*zi4
1661 sz5 = xi3*yi4 - yi3*xi4
1662
1663 sx4 = yi4*zi5 - zi4*yi5
1664 sy4 = zi4*xi5 - xi4*zi5
1665 sz4 = xi4*yi5 - yi4*xi5
1666
1667 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
1668 lb0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
1669 lc0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
1670 la0 = one-lb0-lc0
1671
1672 IF(lb0 >= zerom .and. lb0 <= unp .and.
1673 . lc0 >= zerom .and. lc0 <= unp .and.
1674 . la0 >= zerom .and. la0 <= unp)THEN
1675 interseca = 1
1676 ENDIF
1677 ENDIF
1678
1679 IF(vo41 <= zero .AND. v41 >= zero)THEN
1680 IF(abs(vo41) < em20)THEN
1681 aaa = one
1682 ELSEIF(abs(v41) < em20)THEN
1683 aaa = zero
1684 ELSE
1685 aaa = v41 / (v41-vo41)
1686 ENDIF
1687
1688 xm4 = xx4 + aaa*(xo4-xx4)
1689 ym4 = yy4 + aaa*(yo4-yy4)
1690 zm4 = zz4 + aaa*(zo4-zz4)
1691
1692 xm1 = xx1 + aaa*(xo1-xx1)
1693 ym1 = yy1 + aaa*(yo1-yy1)
1694 zm1 = zz1 + aaa*(zo1-zz1)
1695
1696 xs = xi + aaa*(xoi-xi)
1697 ys = yi + aaa*(yoi-yi)
1698 zs = zi + aaa*(zoi-zi)
1699
1700 xm5 = xx5 + aaa*(xo5-xx5)
1701 ym5 = yy5 + aaa*(yo5-yy5)
1702 zm5 = zz5 + aaa*(zo5-zz5)
1703
1704 x54 = xm4 - xm5
1705 y54 = ym4 - ym5
1706 z54 = zm4 - zm5
1707
1708 x51 = xm1 - xm5
1709 y51 = ym1 - ym5
1710 z51 = zm1 - zm5
1711
1712 xi4 = xm4 - xs
1713 yi4 = ym4 - ys
1714 zi4 = zm4 - zs
1715
1716 xi1 = xm1 - xs
1717 yi1 = ym1 - ys
1718 zi1 = zm1 - zs
1719
1720 xi5 = xm5 - xs
1721 yi5 = ym5 - ys
1722 zi5 = zm5 - zs
1723
1724 sx0 = y54*z51 - z54*y51
1725 sy0 = z54*x51 - x54*z51
1726 sz0 = x54*y51 - y54*x51
1727
1728 sx4 = yi5*zi4 - zi5*yi4
1729 sy4 = zi5*xi4 - xi5*zi4
1730 sz4 = xi5*yi4 - yi5*xi4
1731
1732 sx5 = yi4*zi1 - zi4*yi1
1733 sy5 = zi4*xi1 - xi4*zi1
1734 sz5 = xi4*yi1 - yi4*xi1
1735
1736 sx1 = yi1*zi5 - zi1*yi5
1737 sy1 = zi1*xi5 - xi1*zi5
1738 sz1 = xi1*yi5 - yi1*xi5
1739
1740 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
1741 lb0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
1742 lc0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
1743 la0 = one-lb0-lc0
1744
1745 IF(lb0 >= zerom .and. lb0 <= unp .and.
1746 . lc0 >= zerom .and. lc0 <= unp .and.
1747 . la0 >= zerom .and. la0 <= unp)THEN
1748 interseca = 1
1749 ENDIF
1750 ENDIF
1751 ENDIF
1752C
1753 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 1760 of file i25dst3_22.F.

1773C-----------------------------------------------
1774C I m p l i c i t T y p e s
1775C-----------------------------------------------
1776#include "implicit_f.inc"
1777C-----------------------------------------------
1778C D u m m y A r g u m e n t s
1779C-----------------------------------------------
1780 INTEGER IXX3 ,IXX4 , INTERSECA, INTERSECB
1781C REAL
1782 my_real
1783 1 v12 ,v23 ,v34 ,v41 ,
1784 2 vo12 ,vo23 ,vo34 ,vo41 ,xx1 ,
1785 3 yy1 ,zz1 ,xx2 ,yy2 ,zz2 ,
1786 4 xx3 ,yy3 ,zz3 ,xx4 ,yy4 ,
1787 5 zz4 ,xx5 ,yy5 ,zz5 ,
1788 6 vxi ,vyi ,vzi ,xo1 ,xo2 ,
1789 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
1790 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
1791 9 zo3 ,zo4 ,zo5 ,zerom ,unp ,
1792 a zeromt ,unpt ,xi ,yi ,zi ,
1793 b xoi ,yoi ,zoi
1794C-----------------------------------------------
1795C L o c a l V a r i a b l e s
1796C-----------------------------------------------
1797 my_real
1798 . x51, x52, x53, x54,
1799 . y51, y52, y53, y54,
1800 . z51, z52, z53, z54,
1801 . xi1, xi2, xi3, xi4, xi5,
1802 . yi1, yi2, yi3, yi4, yi5,
1803 . zi1, zi2, zi3, zi4, zi5,
1804 . xia, xib, xic,
1805 . yia, yib, yic,
1806 . zia, zib, zic,
1807 . xs, ys, zs,
1808 . xm1, xm2, xm3, xm4, xm5,
1809 . ym1, ym2, ym3, ym4, ym5,
1810 . zm1, zm2, zm3, zm4, zm5,la0,lb0,lc0
1811 my_real
1812 . sx,sx1,sx2,sx3,sx4,sx5,sx0,
1813 . sy,sy1,sy2,sy3,sy4,sy5,sy0,
1814 . sz,sz1,sz2,sz3,sz4,sz5,sz0,aaa,s2
1815C------------------------------------------------
1816 interseca = 0
1817 intersecb = 0
1818
1819 IF(vo12 * v12 <= zero)THEN
1820 IF(abs(vo12) < em20)THEN
1821 aaa = one
1822 ELSEIF(abs(v12) < em20)THEN
1823 aaa = zero
1824 ELSE
1825 aaa = v12 / (v12-vo12)
1826 ENDIF
1827
1828 xm1 = xx1 + aaa*(xo1-xx1)
1829 ym1 = yy1 + aaa*(yo1-yy1)
1830 zm1 = zz1 + aaa*(zo1-zz1)
1831
1832 xm2 = xx2 + aaa*(xo2-xx2)
1833 ym2 = yy2 + aaa*(yo2-yy2)
1834 zm2 = zz2 + aaa*(zo2-zz2)
1835
1836 xs = xi + aaa*(xoi-xi)
1837 ys = yi + aaa*(yoi-yi)
1838 zs = zi + aaa*(zoi-zi)
1839
1840 xm5 = xx5 + aaa*(xo5-xx5)
1841 ym5 = yy5 + aaa*(yo5-yy5)
1842 zm5 = zz5 + aaa*(zo5-zz5)
1843
1844 x51 = xm1 - xm5
1845 y51 = ym1 - ym5
1846 z51 = zm1 - zm5
1847
1848 x52 = xm2 - xm5
1849 y52 = ym2 - ym5
1850 z52 = zm2 - zm5
1851
1852 xi1 = xm1 - xs
1853 yi1 = ym1 - ys
1854 zi1 = zm1 - zs
1855
1856 xi2 = xm2 - xs
1857 yi2 = ym2 - ys
1858 zi2 = zm2 - zs
1859
1860 xi5 = xm5 - xs
1861 yi5 = ym5 - ys
1862 zi5 = zm5 - zs
1863
1864 sx0 = y51*z52 - z51*y52
1865 sy0 = z51*x52 - x51*z52
1866 sz0 = x51*y52 - y51*x52
1867
1868 sx1 = yi5*zi1 - zi5*yi1
1869 sy1 = zi5*xi1 - xi5*zi1
1870 sz1 = xi5*yi1 - yi5*xi1
1871
1872 sx5 = yi1*zi2 - zi1*yi2
1873 sy5 = zi1*xi2 - xi1*zi2
1874 sz5 = xi1*yi2 - yi1*xi2
1875
1876 sx2 = yi2*zi5 - zi2*yi5
1877 sy2 = zi2*xi5 - xi2*zi5
1878 sz2 = xi2*yi5 - yi2*xi5
1879
1880 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
1881 lb0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
1882 lc0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
1883 la0 = one-lb0-lc0
1884
1885 IF(ixx3 /= ixx4)THEN
1886 IF(lb0 >= zerom .and. lb0 <= unp .and.
1887 . lc0 >= zerom .and. lc0 <= unp .and.
1888 . la0 >= zerom .and. la0 <= unp)THEN
1889 IF(vo12 <= zero .AND. v12 >= zero)THEN
1890 interseca = 1
1891 ELSE ! (VO12 > ZERO .AND. V12 < ZERO) .OR.
1892 ! (VO12 > ZERO .AND. V12 == ZERO) .OR.
1893 ! (VO12 == ZERO .AND. V12 < ZERO)
1894 intersecb = 1
1895 END IF
1896 ENDIF
1897C----------loose tolerance for tria
1898 ELSE
1899 IF(lb0 >= zeromt .AND. lb0 <= unpt .AND.
1900 . lc0 >= zeromt .AND. lc0 <= unpt .AND.
1901 . la0 >= zeromt .AND. la0 <= unpt)THEN
1902 IF(vo12 <= zero .AND. v12 >= zero)THEN
1903 interseca = 1
1904 ELSE ! (VO12 > ZERO .AND. V12 < ZERO) .OR.
1905 ! (VO12 > ZERO .AND. V12 == ZERO) .OR.
1906 ! (VO12 == ZERO .AND. V12 < ZERO)
1907 intersecb = 1
1908 ENDIF
1909 ENDIF
1910 END if! (IXX(I,3) /= IXX(I,4))THEN
1911 ENDIF
1912 IF(ixx3 /= ixx4)THEN
1913 IF(vo23 * v23 <= zero)THEN
1914 IF(abs(vo23) < em20)THEN
1915 aaa = one
1916 ELSEIF(abs(v23) < em20)THEN
1917 aaa = zero
1918 ELSE
1919 aaa = v23 / (v23-vo23)
1920 ENDIF
1921
1922 xm2 = xx2 + aaa*(xo2-xx2)
1923 ym2 = yy2 + aaa*(yo2-yy2)
1924 zm2 = zz2 + aaa*(zo2-zz2)
1925
1926 xm3 = xx3 + aaa*(xo3-xx3)
1927 ym3 = yy3 + aaa*(yo3-yy3)
1928 zm3 = zz3 + aaa*(zo3-zz3)
1929
1930 xs = xi + aaa*(xoi-xi)
1931 ys = yi + aaa*(yoi-yi)
1932 zs = zi + aaa*(zoi-zi)
1933
1934 xm5 = xx5 + aaa*(xo5-xx5)
1935 ym5 = yy5 + aaa*(yo5-yy5)
1936 zm5 = zz5 + aaa*(zo5-zz5)
1937
1938 x52 = xm2 - xm5
1939 y52 = ym2 - ym5
1940 z52 = zm2 - zm5
1941
1942 x53 = xm3 - xm5
1943 y53 = ym3 - ym5
1944 z53 = zm3 - zm5
1945
1946 xi2 = xm2 - xs
1947 yi2 = ym2 - ys
1948 zi2 = zm2 - zs
1949
1950 xi3 = xm3 - xs
1951 yi3 = ym3 - ys
1952 zi3 = zm3 - zs
1953
1954 xi5 = xm5 - xs
1955 yi5 = ym5 - ys
1956 zi5 = zm5 - zs
1957
1958 sx0 = y52*z53 - z52*y53
1959 sy0 = z52*x53 - x52*z53
1960 sz0 = x52*y53 - y52*x53
1961
1962 sx2 = yi5*zi2 - zi5*yi2
1963 sy2 = zi5*xi2 - xi5*zi2
1964 sz2 = xi5*yi2 - yi5*xi2
1965
1966 sx5 = yi2*zi3 - zi2*yi3
1967 sy5 = zi2*xi3 - xi2*zi3
1968 sz5 = xi2*yi3 - yi2*xi3
1969
1970 sx3 = yi3*zi5 - zi3*yi5
1971 sy3 = zi3*xi5 - xi3*zi5
1972 sz3 = xi3*yi5 - yi3*xi5
1973
1974 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
1975 lb0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
1976 lc0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
1977 la0 = one-lb0-lc0
1978
1979 IF(lb0 >= zerom .and. lb0 <= unp .and.
1980 . lc0 >= zerom .and. lc0 <= unp .and.
1981 . la0 >= zerom .and. la0 <= unp)THEN
1982 IF(vo23 <= zero .AND. v23 >= zero)THEN
1983 interseca = 1
1984 ELSE ! (VO23 > ZERO .AND. V23 < ZERO) .OR.
1985 ! (VO23 > ZERO .AND. V23 == ZERO) .OR.
1986 ! (VO23 == ZERO .AND. V23 < ZERO)
1987 intersecb = 1
1988 END IF
1989 ENDIF
1990 ENDIF
1991 IF(vo34 * v34 <= zero)THEN
1992 IF(abs(vo34) < em20)THEN
1993 aaa = one
1994 ELSEIF(abs(v34) < em20)THEN
1995 aaa = zero
1996 ELSE
1997 aaa = v34 / (v34-vo34)
1998 ENDIF
1999
2000 xm3 = xx3 + aaa*(xo3-xx3)
2001 ym3 = yy3 + aaa*(yo3-yy3)
2002 zm3 = zz3 + aaa*(zo3-zz3)
2003
2004 xm4 = xx4 + aaa*(xo4-xx4)
2005 ym4 = yy4 + aaa*(yo4-yy4)
2006 zm4 = zz4 + aaa*(zo4-zz4)
2007
2008 xs = xi + aaa*(xoi-xi)
2009 ys = yi + aaa*(yoi-yi)
2010 zs = zi + aaa*(zoi-zi)
2011
2012 xm5 = xx5 + aaa*(xo5-xx5)
2013 ym5 = yy5 + aaa*(yo5-yy5)
2014 zm5 = zz5 + aaa*(zo5-zz5)
2015
2016 x53 = xm3 - xm5
2017 y53 = ym3 - ym5
2018 z53 = zm3 - zm5
2019
2020 x54 = xm4 - xm5
2021 y54 = ym4 - ym5
2022 z54 = zm4 - zm5
2023
2024 xi3 = xm3 - xs
2025 yi3 = ym3 - ys
2026 zi3 = zm3 - zs
2027
2028 xi4 = xm4 - xs
2029 yi4 = ym4 - ys
2030 zi4 = zm4 - zs
2031
2032 xi5 = xm5 - xs
2033 yi5 = ym5 - ys
2034 zi5 = zm5 - zs
2035
2036 sx0 = y53*z54 - z53*y54
2037 sy0 = z53*x54 - x53*z54
2038 sz0 = x53*y54 - y53*x54
2039
2040 sx3 = yi5*zi3 - zi5*yi3
2041 sy3 = zi5*xi3 - xi5*zi3
2042 sz3 = xi5*yi3 - yi5*xi3
2043
2044 sx5 = yi3*zi4 - zi3*yi4
2045 sy5 = zi3*xi4 - xi3*zi4
2046 sz5 = xi3*yi4 - yi3*xi4
2047
2048 sx4 = yi4*zi5 - zi4*yi5
2049 sy4 = zi4*xi5 - xi4*zi5
2050 sz4 = xi4*yi5 - yi4*xi5
2051
2052 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2053 lb0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
2054 lc0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
2055 la0 = one-lb0-lc0
2056
2057 IF(lb0 >= zerom .and. lb0 <= unp .and.
2058 . lc0 >= zerom .and. lc0 <= unp .and.
2059 . la0 >= zerom .and. la0 <= unp)THEN
2060 IF(vo34 <= zero .AND. v34 >= zero)THEN
2061 interseca = 1
2062 ELSE ! (VO34 > ZERO .AND. V34 < ZERO) .OR.
2063 ! (VO34 > ZERO .AND. V34 == ZERO) .OR.
2064 ! (VO34 == ZERO .AND. V34 < ZERO)
2065 intersecb = 1
2066 END IF
2067 ENDIF
2068 ENDIF
2069
2070 IF(vo41 * v41 <= zero)THEN
2071 IF(abs(vo41) < em20)THEN
2072 aaa = one
2073 ELSEIF(abs(v41) < em20)THEN
2074 aaa = zero
2075 ELSE
2076 aaa = v41 / (v41-vo41)
2077 ENDIF
2078
2079 xm4 = xx4 + aaa*(xo4-xx4)
2080 ym4 = yy4 + aaa*(yo4-yy4)
2081 zm4 = zz4 + aaa*(zo4-zz4)
2082
2083 xm1 = xx1 + aaa*(xo1-xx1)
2084 ym1 = yy1 + aaa*(yo1-yy1)
2085 zm1 = zz1 + aaa*(zo1-zz1)
2086
2087 xs = xi + aaa*(xoi-xi)
2088 ys = yi + aaa*(yoi-yi)
2089 zs = zi + aaa*(zoi-zi)
2090
2091 xm5 = xx5 + aaa*(xo5-xx5)
2092 ym5 = yy5 + aaa*(yo5-yy5)
2093 zm5 = zz5 + aaa*(zo5-zz5)
2094
2095 x54 = xm4 - xm5
2096 y54 = ym4 - ym5
2097 z54 = zm4 - zm5
2098
2099 x51 = xm1 - xm5
2100 y51 = ym1 - ym5
2101 z51 = zm1 - zm5
2102
2103 xi4 = xm4 - xs
2104 yi4 = ym4 - ys
2105 zi4 = zm4 - zs
2106
2107 xi1 = xm1 - xs
2108 yi1 = ym1 - ys
2109 zi1 = zm1 - zs
2110
2111 xi5 = xm5 - xs
2112 yi5 = ym5 - ys
2113 zi5 = zm5 - zs
2114
2115 sx0 = y54*z51 - z54*y51
2116 sy0 = z54*x51 - x54*z51
2117 sz0 = x54*y51 - y54*x51
2118
2119 sx4 = yi5*zi4 - zi5*yi4
2120 sy4 = zi5*xi4 - xi5*zi4
2121 sz4 = xi5*yi4 - yi5*xi4
2122
2123 sx5 = yi4*zi1 - zi4*yi1
2124 sy5 = zi4*xi1 - xi4*zi1
2125 sz5 = xi4*yi1 - yi4*xi1
2126
2127 sx1 = yi1*zi5 - zi1*yi5
2128 sy1 = zi1*xi5 - xi1*zi5
2129 sz1 = xi1*yi5 - yi1*xi5
2130
2131 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2132 lb0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
2133 lc0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
2134 la0 = one-lb0-lc0
2135
2136 IF(lb0 >= zerom .and. lb0 <= unp .and.
2137 . lc0 >= zerom .and. lc0 <= unp .and.
2138 . la0 >= zerom .and. la0 <= unp)THEN
2139 IF(vo41 <= zero .AND. v41 >= zero)THEN
2140 interseca = 1
2141 ELSE ! (VO41 > ZERO .AND. V41 < ZERO) .OR.
2142 ! (VO41 > ZERO .AND. V41 == ZERO) .OR.
2143 ! (VO41 == ZERO .AND. V41 < ZERO)
2144 intersecb = 1
2145 END IF
2146 ENDIF
2147 ENDIF
2148 ENDIF
2149C
2150 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 2546 of file i25dst3_22.F.

2558C-----------------------------------------------
2559C I m p l i c i t T y p e s
2560C-----------------------------------------------
2561#include "implicit_f.inc"
2562C-----------------------------------------------
2563C D u m m y A r g u m e n t s
2564C-----------------------------------------------
2565 INTEGER IXX3 ,IXX4 ,INTERSECT
2566C REAL
2567 my_real
2568 1 v12 ,v23 ,v34 ,v41 ,
2569 2 vo12 ,vo23 ,vo34 ,vo41 ,xx1 ,
2570 3 yy1 ,zz1 ,xx2 ,yy2 ,zz2 ,
2571 4 xx3 ,yy3 ,zz3 ,xx4 ,yy4 ,
2572 5 zz4 ,xx5 ,yy5 ,zz5 ,
2573 6 xoi ,yoi ,zoi ,xo1 ,xo2 ,
2574 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
2575 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
2576 9 zo3 ,zo4 ,zo5 ,zeromt,unpt ,
2577 a xi ,yi ,zi
2578C-----------------------------------------------
2579C C o m m o n B l o c k s
2580C-----------------------------------------------
2581#include "com08_c.inc"
2582C-----------------------------------------------
2583C L o c a l V a r i a b l e s
2584C-----------------------------------------------
2585 my_real
2586 . x51, x52, x53, x54,
2587 . y51, y52, y53, y54,
2588 . z51, z52, z53, z54,
2589 . xi1, xi2, xi3, xi4, xi5,
2590 . yi1, yi2, yi3, yi4, yi5,
2591 . zi1, zi2, zi3, zi4, zi5,
2592 . xia, xib, xic,
2593 . yia, yib, yic,
2594 . zia, zib, zic,
2595 . xs,ys,zs,
2596 . xm1, xm2, xm3, xm4, xm5,
2597 . ym1, ym2, ym3, ym4, ym5,
2598 . zm1, zm2, zm3, zm4, zm5,la0,lb0,lc0
2599 my_real
2600 . sx,sx1,sx2,sx3,sx4,sx5,sx0,
2601 . sy,sy1,sy2,sy3,sy4,sy5,sy0,
2602 . sz,sz1,sz2,sz3,sz4,sz5,sz0,aaa,adt,s2
2603C------------------------------------------------
2604C---------------- cas V>0 V_old>0
2605 IF(vo12 > zero .and. v12 >= zero )THEN
2606 aaa = one
2607 adt = dt1
2608
2609 xm1 = xo1
2610 ym1 = yo1
2611 zm1 = zo1
2612
2613 xm2 = xo2
2614 ym2 = yo2
2615 zm2 = zo2
2616
2617c XOI = XI - ADT*VXI
2618c YOI = YI - ADT*VYI
2619c ZOI = ZI - ADT*VZI
2620
2621 xm5 = xo5
2622 ym5 = yo5
2623 zm5 = zo5
2624
2625 x51 = xm1 - xm5
2626 y51 = ym1 - ym5
2627 z51 = zm1 - zm5
2628
2629 x52 = xm2 - xm5
2630 y52 = ym2 - ym5
2631 z52 = zm2 - zm5
2632
2633 xi1 = xm1 - xoi
2634 yi1 = ym1 - yoi
2635 zi1 = zm1 - zoi
2636
2637 xi2 = xm2 - xoi
2638 yi2 = ym2 - yoi
2639 zi2 = zm2 - zoi
2640
2641 xi5 = xm5 - xoi
2642 yi5 = ym5 - yoi
2643 zi5 = zm5 - zoi
2644
2645 sx0 = y51*z52 - z51*y52
2646 sy0 = z51*x52 - x51*z52
2647 sz0 = x51*y52 - y51*x52
2648
2649 sx1 = yi5*zi1 - zi5*yi1
2650 sy1 = zi5*xi1 - xi5*zi1
2651 sz1 = xi5*yi1 - yi5*xi1
2652
2653 sx5 = yi1*zi2 - zi1*yi2
2654 sy5 = zi1*xi2 - xi1*zi2
2655 sz5 = xi1*yi2 - yi1*xi2
2656
2657 sx2 = yi2*zi5 - zi2*yi5
2658 sy2 = zi2*xi5 - xi2*zi5
2659 sz2 = xi2*yi5 - yi2*xi5
2660
2661 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2662 lb0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
2663 lc0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
2664 la0 = one-lb0-lc0
2665
2666 IF(lb0 >= zeromt .AND. lb0 <= unpt .AND.
2667 . lc0 >= zeromt .AND. lc0 <= unpt .AND.
2668 . la0 >= zeromt .AND. la0 <= unpt)THEN
2669 intersect = 1
2670 ENDIF
2671 ENDIF
2672
2673 IF(ixx3 /= ixx4)THEN
2674 IF(vo23 > zero .AND. v23 >= zero)THEN
2675 aaa = one
2676 adt = dt1
2677
2678 xm2 = xo2
2679 ym2 = yo2
2680 zm2 = zo2
2681
2682 xm3 = xo3
2683 ym3 = yo3
2684 zm3 = zo3
2685
2686c XOI = XI - ADT*VXI
2687c YOI = YI - ADT*VYI
2688c ZOI = ZI - ADT*VZI
2689
2690 xm5 = xo5
2691 ym5 = yo5
2692 zm5 = zo5
2693
2694 x52 = xm2 - xm5
2695 y52 = ym2 - ym5
2696 z52 = zm2 - zm5
2697
2698 x53 = xm3 - xm5
2699 y53 = ym3 - ym5
2700 z53 = zm3 - zm5
2701
2702 xi2 = xm2 - xoi
2703 yi2 = ym2 - yoi
2704 zi2 = zm2 - zoi
2705
2706 xi3 = xm3 - xoi
2707 yi3 = ym3 - yoi
2708 zi3 = zm3 - zoi
2709
2710 xi5 = xm5 - xoi
2711 yi5 = ym5 - yoi
2712 zi5 = zm5 - zoi
2713
2714 sx0 = y52*z53 - z52*y53
2715 sy0 = z52*x53 - x52*z53
2716 sz0 = x52*y53 - y52*x53
2717
2718 sx2 = yi5*zi2 - zi5*yi2
2719 sy2 = zi5*xi2 - xi5*zi2
2720 sz2 = xi5*yi2 - yi5*xi2
2721
2722 sx5 = yi2*zi3 - zi2*yi3
2723 sy5 = zi2*xi3 - xi2*zi3
2724 sz5 = xi2*yi3 - yi2*xi3
2725
2726 sx3 = yi3*zi5 - zi3*yi5
2727 sy3 = zi3*xi5 - xi3*zi5
2728 sz3 = xi3*yi5 - yi3*xi5
2729
2730 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2731 lb0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
2732 lc0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
2733 la0 = one-lb0-lc0
2734
2735 IF(lb0 >= zeromt.and. lb0 <= unpt .and.
2736 . lc0 >= zeromt.and. lc0 <= unpt .and.
2737 . la0 >= zeromt.and. la0 <= unpt)THEN
2738 intersect = 1
2739 ENDIF
2740 ENDIF
2741 IF(vo34 > zero .and. v34 >= zero )THEN
2742 aaa = one
2743 adt = dt1
2744
2745 xm3 = xo3
2746 ym3 = yo3
2747 zm3 = zo3
2748
2749 xm4 = xo4
2750 ym4 = yo4
2751 zm4 = zo4
2752
2753c XOI = XI - ADT*VXI
2754c YOI = YI - ADT*VYI
2755c ZOI = ZI - ADT*VZI
2756
2757 xm5 = xo5
2758 ym5 = yo5
2759 zm5 = zo5
2760
2761 x53 = xm3 - xm5
2762 y53 = ym3 - ym5
2763 z53 = zm3 - zm5
2764
2765 x54 = xm4 - xm5
2766 y54 = ym4 - ym5
2767 z54 = zm4 - zm5
2768
2769 xi3 = xm3 - xoi
2770 yi3 = ym3 - yoi
2771 zi3 = zm3 - zoi
2772
2773 xi4 = xm4 - xoi
2774 yi4 = ym4 - yoi
2775 zi4 = zm4 - zoi
2776
2777 xi5 = xm5 - xoi
2778 yi5 = ym5 - yoi
2779 zi5 = zm5 - zoi
2780
2781 sx0 = y53*z54 - z53*y54
2782 sy0 = z53*x54 - x53*z54
2783 sz0 = x53*y54 - y53*x54
2784
2785 sx3 = yi5*zi3 - zi5*yi3
2786 sy3 = zi5*xi3 - xi5*zi3
2787 sz3 = xi5*yi3 - yi5*xi3
2788
2789 sx5 = yi3*zi4 - zi3*yi4
2790 sy5 = zi3*xi4 - xi3*zi4
2791 sz5 = xi3*yi4 - yi3*xi4
2792
2793 sx4 = yi4*zi5 - zi4*yi5
2794 sy4 = zi4*xi5 - xi4*zi5
2795 sz4 = xi4*yi5 - yi4*xi5
2796
2797 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2798 lb0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
2799 lc0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
2800 la0 = one-lb0-lc0
2801
2802 IF(lb0 >= zeromt.and. lb0 <= unpt .and.
2803 . lc0 >= zeromt.and. lc0 <= unpt .and.
2804 . la0 >= zeromt.and. la0 <= unpt)THEN
2805 intersect = 1
2806 ENDIF
2807 ENDIF
2808
2809 IF(vo41 > zero .and. v41 >= zero )THEN
2810 aaa = one
2811 adt = dt1
2812
2813 xm4 = xo4
2814 ym4 = yo4
2815 zm4 = zo4
2816
2817 xm1 = xo1
2818 ym1 = yo1
2819 zm1 = zo1
2820
2821c XOI = XI - ADT*VXI
2822c YOI = YI - ADT*VYI
2823c ZOI = ZI - ADT*VZI
2824
2825 xm5 = xo5
2826 ym5 = yo5
2827 zm5 = zo5
2828
2829 x54 = xm4 - xm5
2830 y54 = ym4 - ym5
2831 z54 = zm4 - zm5
2832
2833 x51 = xm1 - xm5
2834 y51 = ym1 - ym5
2835 z51 = zm1 - zm5
2836
2837 xi4 = xm4 - xoi
2838 yi4 = ym4 - yoi
2839 zi4 = zm4 - zoi
2840
2841 xi1 = xm1 - xoi
2842 yi1 = ym1 - yoi
2843 zi1 = zm1 - zoi
2844
2845 xi5 = xm5 - xoi
2846 yi5 = ym5 - yoi
2847 zi5 = zm5 - zoi
2848
2849 sx0 = y54*z51 - z54*y51
2850 sy0 = z54*x51 - x54*z51
2851 sz0 = x54*y51 - y54*x51
2852
2853 sx4 = yi5*zi4 - zi5*yi4
2854 sy4 = zi5*xi4 - xi5*zi4
2855 sz4 = xi5*yi4 - yi5*xi4
2856
2857 sx5 = yi4*zi1 - zi4*yi1
2858 sy5 = zi4*xi1 - xi4*zi1
2859 sz5 = xi4*yi1 - yi4*xi1
2860
2861 sx1 = yi1*zi5 - zi1*yi5
2862 sy1 = zi1*xi5 - xi1*zi5
2863 sz1 = xi1*yi5 - yi1*xi5
2864
2865 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2866 lb0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
2867 lc0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
2868 la0 = one-lb0-lc0
2869
2870 IF(lb0 >= zeromt.and. lb0 <= unpt .and.
2871 . lc0 >= zeromt.and. lc0 <= unpt .and.
2872 . la0 >= zeromt.and. la0 <= unpt)THEN
2873 intersect = 1
2874 ENDIF
2875 ENDIF
2876 ENDIF
2877C
2878 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 2155 of file i25dst3_22.F.

2168C-----------------------------------------------
2169C I m p l i c i t T y p e s
2170C-----------------------------------------------
2171#include "implicit_f.inc"
2172C-----------------------------------------------
2173C D u m m y A r g u m e n t s
2174C-----------------------------------------------
2175 INTEGER IXX3 ,IXX4 ,INTERSECT, SUBTRIA
2176C REAL
2177 my_real
2178 1 v12 ,v23 ,v34 ,v41 ,
2179 2 vo12 ,vo23 ,vo34 ,vo41 ,xx1 ,
2180 3 yy1 ,zz1 ,xx2 ,yy2 ,zz2 ,
2181 4 xx3 ,yy3 ,zz3 ,xx4 ,yy4 ,
2182 5 zz4 ,xx5 ,yy5 ,zz5 ,
2183 6 vxi ,vyi ,vzi ,xo1 ,xo2 ,
2184 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
2185 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
2186 9 zo3 ,zo4 ,zo5 ,zerom ,unp ,
2187 a zeromt ,unpt ,xi ,yi ,zi ,
2188 b xoi ,yoi ,zoi
2189C-----------------------------------------------
2190C L o c a l V a r i a b l e s
2191C-----------------------------------------------
2192 my_real
2193 . x51, x52, x53, x54,
2194 . y51, y52, y53, y54,
2195 . z51, z52, z53, z54,
2196 . xi1, xi2, xi3, xi4, xi5,
2197 . yi1, yi2, yi3, yi4, yi5,
2198 . zi1, zi2, zi3, zi4, zi5,
2199 . xia, xib, xic,
2200 . yia, yib, yic,
2201 . zia, zib, zic,
2202 . xs, ys, zs,
2203 . xm1, xm2, xm3, xm4, xm5,
2204 . ym1, ym2, ym3, ym4, ym5,
2205 . zm1, zm2, zm3, zm4, zm5,la0,lb0,lc0
2206 my_real
2207 . sx,sx1,sx2,sx3,sx4,sx5,sx0,
2208 . sy,sy1,sy2,sy3,sy4,sy5,sy0,
2209 . sz,sz1,sz2,sz3,sz4,sz5,sz0,aaa,s2
2210C------------------------------------------------
2211 IF(vo12 <= zero .AND. v12 >= zero)THEN
2212 IF(abs(vo12) < em20)THEN
2213 aaa = one
2214 ELSEIF(abs(v12) < em20)THEN
2215 aaa = zero
2216 ELSE
2217 aaa = v12 / (v12-vo12)
2218 ENDIF
2219
2220 xm1 = xx1 + aaa*(xo1-xx1)
2221 ym1 = yy1 + aaa*(yo1-yy1)
2222 zm1 = zz1 + aaa*(zo1-zz1)
2223
2224 xm2 = xx2 + aaa*(xo2-xx2)
2225 ym2 = yy2 + aaa*(yo2-yy2)
2226 zm2 = zz2 + aaa*(zo2-zz2)
2227
2228 xs = xi + aaa*(xoi-xi)
2229 ys = yi + aaa*(yoi-yi)
2230 zs = zi + aaa*(zoi-zi)
2231
2232 xm5 = xx5 + aaa*(xo5-xx5)
2233 ym5 = yy5 + aaa*(yo5-yy5)
2234 zm5 = zz5 + aaa*(zo5-zz5)
2235
2236 x51 = xm1 - xm5
2237 y51 = ym1 - ym5
2238 z51 = zm1 - zm5
2239
2240 x52 = xm2 - xm5
2241 y52 = ym2 - ym5
2242 z52 = zm2 - zm5
2243
2244 xi1 = xm1 - xs
2245 yi1 = ym1 - ys
2246 zi1 = zm1 - zs
2247
2248 xi2 = xm2 - xs
2249 yi2 = ym2 - ys
2250 zi2 = zm2 - zs
2251
2252 xi5 = xm5 - xs
2253 yi5 = ym5 - ys
2254 zi5 = zm5 - zs
2255
2256 sx0 = y51*z52 - z51*y52
2257 sy0 = z51*x52 - x51*z52
2258 sz0 = x51*y52 - y51*x52
2259
2260 sx1 = yi5*zi1 - zi5*yi1
2261 sy1 = zi5*xi1 - xi5*zi1
2262 sz1 = xi5*yi1 - yi5*xi1
2263
2264 sx5 = yi1*zi2 - zi1*yi2
2265 sy5 = zi1*xi2 - xi1*zi2
2266 sz5 = xi1*yi2 - yi1*xi2
2267
2268 sx2 = yi2*zi5 - zi2*yi5
2269 sy2 = zi2*xi5 - xi2*zi5
2270 sz2 = xi2*yi5 - yi2*xi5
2271
2272 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2273 lb0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
2274 lc0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
2275 la0 = one-lb0-lc0
2276
2277 IF(ixx3 /= ixx4)THEN
2278 IF(lb0 >= zerom .and. lb0 <= unp .and.
2279 . lc0 >= zerom .and. lc0 <= unp .and.
2280 . la0 >= zerom .and. la0 <= unp)THEN
2281 IF(vo12 <= zero .AND. v12 >= zero)THEN
2282 intersect = 1
2283 ELSE ! (VO12 > ZERO .AND. V12 < ZERO) .OR.
2284 ! (VO12 > ZERO .AND. V12 == ZERO) .OR.
2285 ! (VO12 == ZERO .AND. V12 < ZERO)
2286 intersect = -1
2287 END IF
2288 subtria= 1 ! subtriangle
2289 ENDIF
2290C----------loose tolerance for tria
2291 ELSE
2292 IF(lb0 >= zeromt .AND. lb0 <= unpt .AND.
2293 . lc0 >= zeromt .AND. lc0 <= unpt .AND.
2294 . la0 >= zeromt .AND. la0 <= unpt)THEN
2295 intersect = 1
2296 subtria= 1 ! subtriangle
2297 ENDIF
2298 END if! (IXX(I,3) /= IXX(I,4))THEN
2299 ENDIF
2300 IF(ixx3 /= ixx4)THEN
2301 IF(vo23 <= zero .AND. v23 >= zero)THEN
2302 IF(abs(vo23) < em20)THEN
2303 aaa = one
2304 ELSEIF(abs(v23) < em20)THEN
2305 aaa = zero
2306 ELSE
2307 aaa = v23 / (v23-vo23)
2308 ENDIF
2309
2310 xm2 = xx2 + aaa*(xo2-xx2)
2311 ym2 = yy2 + aaa*(yo2-yy2)
2312 zm2 = zz2 + aaa*(zo2-zz2)
2313
2314 xm3 = xx3 + aaa*(xo3-xx3)
2315 ym3 = yy3 + aaa*(yo3-yy3)
2316 zm3 = zz3 + aaa*(zo3-zz3)
2317
2318 xs = xi + aaa*(xoi-xi)
2319 ys = yi + aaa*(yoi-yi)
2320 zs = zi + aaa*(zoi-zi)
2321
2322 xm5 = xx5 + aaa*(xo5-xx5)
2323 ym5 = yy5 + aaa*(yo5-yy5)
2324 zm5 = zz5 + aaa*(zo5-zz5)
2325
2326 x52 = xm2 - xm5
2327 y52 = ym2 - ym5
2328 z52 = zm2 - zm5
2329
2330 x53 = xm3 - xm5
2331 y53 = ym3 - ym5
2332 z53 = zm3 - zm5
2333
2334 xi2 = xm2 - xs
2335 yi2 = ym2 - ys
2336 zi2 = zm2 - zs
2337
2338 xi3 = xm3 - xs
2339 yi3 = ym3 - ys
2340 zi3 = zm3 - zs
2341
2342 xi5 = xm5 - xs
2343 yi5 = ym5 - ys
2344 zi5 = zm5 - zs
2345
2346 sx0 = y52*z53 - z52*y53
2347 sy0 = z52*x53 - x52*z53
2348 sz0 = x52*y53 - y52*x53
2349
2350 sx2 = yi5*zi2 - zi5*yi2
2351 sy2 = zi5*xi2 - xi5*zi2
2352 sz2 = xi5*yi2 - yi5*xi2
2353
2354 sx5 = yi2*zi3 - zi2*yi3
2355 sy5 = zi2*xi3 - xi2*zi3
2356 sz5 = xi2*yi3 - yi2*xi3
2357
2358 sx3 = yi3*zi5 - zi3*yi5
2359 sy3 = zi3*xi5 - xi3*zi5
2360 sz3 = xi3*yi5 - yi3*xi5
2361
2362 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2363 lb0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
2364 lc0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
2365 la0 = one-lb0-lc0
2366
2367 IF(lb0 >= zerom .and. lb0 <= unp .and.
2368 . lc0 >= zerom .and. lc0 <= unp .and.
2369 . la0 >= zerom .and. la0 <= unp)THEN
2370 IF(vo23 <= zero .AND. v23 >= zero)THEN
2371 intersect = 1
2372 ELSE ! (VO23 > ZERO .AND. V23 < ZERO) .OR.
2373 ! (VO23 > ZERO .AND. V23 == ZERO) .OR.
2374 ! (VO23 == ZERO .AND. V23 < ZERO)
2375 intersect = -1
2376 END IF
2377
2378 ENDIF
2379 ENDIF
2380 IF(vo34 <= zero .AND. v34 >= zero)THEN
2381 IF(abs(vo34) < em20)THEN
2382 aaa = one
2383 ELSEIF(abs(v34) < em20)THEN
2384 aaa = zero
2385 ELSE
2386 aaa = v34 / (v34-vo34)
2387 ENDIF
2388
2389 xm3 = xx3 + aaa*(xo3-xx3)
2390 ym3 = yy3 + aaa*(yo3-yy3)
2391 zm3 = zz3 + aaa*(zo3-zz3)
2392
2393 xm4 = xx4 + aaa*(xo4-xx4)
2394 ym4 = yy4 + aaa*(yo4-yy4)
2395 zm4 = zz4 + aaa*(zo4-zz4)
2396
2397 xs = xi + aaa*(xoi-xi)
2398 ys = yi + aaa*(yoi-yi)
2399 zs = zi + aaa*(zoi-zi)
2400
2401 xm5 = xx5 + aaa*(xo5-xx5)
2402 ym5 = yy5 + aaa*(yo5-yy5)
2403 zm5 = zz5 + aaa*(zo5-zz5)
2404
2405 x53 = xm3 - xm5
2406 y53 = ym3 - ym5
2407 z53 = zm3 - zm5
2408
2409 x54 = xm4 - xm5
2410 y54 = ym4 - ym5
2411 z54 = zm4 - zm5
2412
2413 xi3 = xm3 - xs
2414 yi3 = ym3 - ys
2415 zi3 = zm3 - zs
2416
2417 xi4 = xm4 - xs
2418 yi4 = ym4 - ys
2419 zi4 = zm4 - zs
2420
2421 xi5 = xm5 - xs
2422 yi5 = ym5 - ys
2423 zi5 = zm5 - zs
2424
2425 sx0 = y53*z54 - z53*y54
2426 sy0 = z53*x54 - x53*z54
2427 sz0 = x53*y54 - y53*x54
2428
2429 sx3 = yi5*zi3 - zi5*yi3
2430 sy3 = zi5*xi3 - xi5*zi3
2431 sz3 = xi5*yi3 - yi5*xi3
2432
2433 sx5 = yi3*zi4 - zi3*yi4
2434 sy5 = zi3*xi4 - xi3*zi4
2435 sz5 = xi3*yi4 - yi3*xi4
2436
2437 sx4 = yi4*zi5 - zi4*yi5
2438 sy4 = zi4*xi5 - xi4*zi5
2439 sz4 = xi4*yi5 - yi4*xi5
2440
2441 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2442 lb0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
2443 lc0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
2444 la0 = one-lb0-lc0
2445
2446 IF(lb0 >= zerom .and. lb0 <= unp .and.
2447 . lc0 >= zerom .and. lc0 <= unp .and.
2448 . la0 >= zerom .and. la0 <= unp)THEN
2449 IF(vo34 <= zero .AND. v34 >= zero)THEN
2450 intersect = 1
2451 ELSE ! (VO34 > ZERO .AND. V34 < ZERO) .OR.
2452 ! (VO34 > ZERO .AND. V34 == ZERO) .OR.
2453 ! (VO34 == ZERO .AND. V34 < ZERO)
2454 intersect = -1
2455 END IF
2456 ENDIF
2457 ENDIF
2458
2459 IF(vo41 <= zero .AND. v41 >= zero)THEN
2460 IF(abs(vo41) < em20)THEN
2461 aaa = one
2462 ELSEIF(abs(v41) < em20)THEN
2463 aaa = zero
2464 ELSE
2465 aaa = v41 / (v41-vo41)
2466 ENDIF
2467
2468 xm4 = xx4 + aaa*(xo4-xx4)
2469 ym4 = yy4 + aaa*(yo4-yy4)
2470 zm4 = zz4 + aaa*(zo4-zz4)
2471
2472 xm1 = xx1 + aaa*(xo1-xx1)
2473 ym1 = yy1 + aaa*(yo1-yy1)
2474 zm1 = zz1 + aaa*(zo1-zz1)
2475
2476 xs = xi + aaa*(xoi-xi)
2477 ys = yi + aaa*(yoi-yi)
2478 zs = zi + aaa*(zoi-zi)
2479
2480 xm5 = xx5 + aaa*(xo5-xx5)
2481 ym5 = yy5 + aaa*(yo5-yy5)
2482 zm5 = zz5 + aaa*(zo5-zz5)
2483
2484 x54 = xm4 - xm5
2485 y54 = ym4 - ym5
2486 z54 = zm4 - zm5
2487
2488 x51 = xm1 - xm5
2489 y51 = ym1 - ym5
2490 z51 = zm1 - zm5
2491
2492 xi4 = xm4 - xs
2493 yi4 = ym4 - ys
2494 zi4 = zm4 - zs
2495
2496 xi1 = xm1 - xs
2497 yi1 = ym1 - ys
2498 zi1 = zm1 - zs
2499
2500 xi5 = xm5 - xs
2501 yi5 = ym5 - ys
2502 zi5 = zm5 - zs
2503
2504 sx0 = y54*z51 - z54*y51
2505 sy0 = z54*x51 - x54*z51
2506 sz0 = x54*y51 - y54*x51
2507
2508 sx4 = yi5*zi4 - zi5*yi4
2509 sy4 = zi5*xi4 - xi5*zi4
2510 sz4 = xi5*yi4 - yi5*xi4
2511
2512 sx5 = yi4*zi1 - zi4*yi1
2513 sy5 = zi4*xi1 - xi4*zi1
2514 sz5 = xi4*yi1 - yi4*xi1
2515
2516 sx1 = yi1*zi5 - zi1*yi5
2517 sy1 = zi1*xi5 - xi1*zi5
2518 sz1 = xi1*yi5 - yi1*xi5
2519
2520 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2521 lb0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
2522 lc0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
2523 la0 = one-lb0-lc0
2524
2525 IF(lb0 >= zerom .and. lb0 <= unp .and.
2526 . lc0 >= zerom .and. lc0 <= unp .and.
2527 . la0 >= zerom .and. la0 <= unp)THEN
2528 IF(vo41 <= zero .AND. v41 >= zero)THEN
2529 intersect = 1
2530 ELSE ! (VO41 > ZERO .AND. V41 < ZERO) .OR.
2531 ! (VO41 > ZERO .AND. V41 == ZERO) .OR.
2532 ! (VO41 == ZERO .AND. V41 < ZERO)
2533 intersect = -1
2534 END IF
2535 ENDIF
2536 ENDIF
2537 ENDIF
2538C
2539 RETURN