OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i25dst3_22.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| i25dst3_22 ../engine/source/interfaces/int25/i25dst3_22.F
25!||--- called by ------------------------------------------------------
26!|| i25comp_2 ../engine/source/interfaces/int25/i25comp_2.F
27!||--- calls -----------------------------------------------------
28!|| interseca_25 ../engine/source/interfaces/int25/i25dst3_22.F
29!|| intersecb_25 ../engine/source/interfaces/int25/i25dst3_22.F
30!|| intersecv0_25 ../engine/source/interfaces/int25/i25dst3_22.F
31!||--- uses -----------------------------------------------------
32!|| tri7box ../engine/share/modules/tri7box.F
33!||====================================================================
34 SUBROUTINE i25dst3_22(
35 1 JLT ,CAND_N ,CAND_E ,ISHEL ,
36 2 XX ,YY ,ZZ ,
37 3 XI ,YI ,ZI ,
38 4 VX1 ,VX2 ,VX3 ,VX4 ,VXI ,
39 5 VY1 ,VY2 ,VY3 ,VY4 ,VYI ,
40 6 VZ1 ,VZ2 ,VZ3 ,VZ4 ,VZI ,
41 7 NIN ,NSN ,IX1 ,
42 9 IX2 ,IX3 ,IX4 ,NSVG ,STIF ,
43 A INACTI ,MSEGLO ,GAPS ,GAPM ,
44 B IRECT ,IRTLM ,TIME_S ,GAP_NM ,
45 C PENE_OLD,STIF_OLD,ITAB ,
46 D PENMIN,EPS0 ,ICONT_I,MARGE ,
47 E NAX ,NAY ,NAZ ,
48 E NBX ,NBY ,NBZ ,
49 J FAR ,PENT ,
50 L SUBTRIA ,LBS ,LCS ,LBP ,LCP ,
51 P MVOISA ,MVOISB,GAPMXL ,IBOUNDA,IBOUNDB,
52 Q VTX_BISECTOR,DRAD,DGAPLOAD)
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
1388 END
1389!||====================================================================
1390!|| interseca_25 ../engine/source/interfaces/int25/i25dst3_22.F
1391!||--- called by ------------------------------------------------------
1392!|| i25dst3_22 ../engine/source/interfaces/int25/i25dst3_22.F
1393!||====================================================================
1394 SUBROUTINE interseca_25(
1395 1 IXX3 ,IXX4 ,INTERSECA,
1396 1 V12 ,V23 ,V34 ,V41 ,
1397 2 VO12 ,VO23 ,VO34 ,VO41 ,XX1 ,
1398 3 YY1 ,ZZ1 ,XX2 ,YY2 ,ZZ2 ,
1399 4 XX3 ,YY3 ,ZZ3 ,XX4 ,YY4 ,
1400 5 ZZ4 ,XX5 ,YY5 ,ZZ5 ,
1401 6 VXI ,VYI ,VZI ,XO1 ,XO2 ,
1402 7 XO3 ,XO4 ,XO5 ,YO1 ,YO2 ,
1403 8 YO3 ,YO4 ,YO5 ,ZO1 ,ZO2 ,
1404 9 ZO3 ,ZO4 ,ZO5 ,XI ,YI ,
1405 A ZI ,XOI ,YOI ,ZOI ,ZEROM ,
1406 B UNP ,ZEROMT,UNPT )
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
1754 END
1755!||====================================================================
1756!|| intersecb_25 ../engine/source/interfaces/int25/i25dst3_22.F
1757!||--- called by ------------------------------------------------------
1758!|| i25dst3_22 ../engine/source/interfaces/int25/i25dst3_22.F
1759!||====================================================================
1760 SUBROUTINE intersecb_25(
1761 1 IXX3 ,IXX4 ,INTERSECA,INTERSECB,
1762 1 V12 ,V23 ,V34 ,V41 ,
1763 2 VO12 ,VO23 ,VO34 ,VO41 ,XX1 ,
1764 3 YY1 ,ZZ1 ,XX2 ,YY2 ,ZZ2 ,
1765 4 XX3 ,YY3 ,ZZ3 ,XX4 ,YY4 ,
1766 5 ZZ4 ,XX5 ,YY5 ,ZZ5 ,
1767 6 VXI ,VYI ,VZI ,XO1 ,XO2 ,
1768 7 XO3 ,XO4 ,XO5 ,YO1 ,YO2 ,
1769 8 YO3 ,YO4 ,YO5 ,ZO1 ,ZO2 ,
1770 9 ZO3 ,ZO4 ,ZO5 ,XI ,YI ,
1771 A ZI ,XOI ,YOI ,ZOI ,ZEROM ,
1772 B UNP ,ZEROMT,UNPT )
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
2151 END
2152!||====================================================================
2153!|| intersecv_25 ../engine/source/interfaces/int25/i25dst3_22.F
2154!||====================================================================
2155 SUBROUTINE intersecv_25(
2156 1 IXX3 ,IXX4 ,SUBTRIA,
2157 1 INTERSECT,V12 ,V23 ,V34 ,V41 ,
2158 2 VO12 ,VO23 ,VO34 ,VO41 ,XX1 ,
2159 3 YY1 ,ZZ1 ,XX2 ,YY2 ,ZZ2 ,
2160 4 XX3 ,YY3 ,ZZ3 ,XX4 ,YY4 ,
2161 5 ZZ4 ,XX5 ,YY5 ,ZZ5 ,
2162 6 VXI ,VYI ,VZI ,XO1 ,XO2 ,
2163 7 XO3 ,XO4 ,XO5 ,YO1 ,YO2 ,
2164 8 YO3 ,YO4 ,YO5 ,ZO1 ,ZO2 ,
2165 9 ZO3 ,ZO4 ,ZO5 ,XI ,YI ,
2166 A ZI ,XOI ,YOI ,ZOI ,ZEROM ,
2167 B UNP ,ZEROMT,UNPT )
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
2540 END
2541!||====================================================================
2542!|| intersecv0_25 ../engine/source/interfaces/int25/i25dst3_22.F
2543!||--- called by ------------------------------------------------------
2544!|| i25dst3_22 ../engine/source/interfaces/int25/i25dst3_22.F
2545!||====================================================================
2546 SUBROUTINE intersecv0_25(
2547 1 IXX3 ,IXX4 ,
2548 1 INTERSECT,V12 ,V23 ,V34 ,V41 ,
2549 2 VO12 ,VO23 ,VO34 ,VO41 ,XX1 ,
2550 3 YY1 ,ZZ1 ,XX2 ,YY2 ,ZZ2 ,
2551 4 XX3 ,YY3 ,ZZ3 ,XX4 ,YY4 ,
2552 5 ZZ4 ,XX5 ,YY5 ,ZZ5 ,
2553 6 XOI ,YOI ,ZOI ,XO1 ,XO2 ,
2554 7 XO3 ,XO4 ,XO5 ,YO1 ,YO2 ,
2555 8 YO3 ,YO4 ,YO5 ,ZO1 ,ZO2 ,
2556 9 ZO3 ,ZO4 ,ZO5 ,XI ,YI ,
2557 A ZI ,ZEROMT , UNPT )
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
2879 END
2880!||====================================================================
2881!|| i25glob_22 ../engine/source/interfaces/int25/i25dst3_22.F
2882!||--- called by ------------------------------------------------------
2883!|| i25comp_2 ../engine/source/interfaces/int25/i25comp_2.F
2884!||--- uses -----------------------------------------------------
2885!|| tri7box ../engine/share/modules/tri7box.F
2886!||====================================================================
2887 SUBROUTINE i25glob_22(
2888 1 JLT ,CAND_N_N ,CAND_E_N ,CAND_E ,
2889 2 NIN ,NSN ,IX1 ,IX2 ,IX3 ,
2890 3 IX4 ,NSVG ,STIF ,INACTI ,MSEGLO ,
2891 4 IRTLM ,TIME_S ,ITAB ,SUBTRIA,
2892 5 FAR ,PENT ,LB ,LC ,
2893 9 INDEX ,FARM ,PENM ,LBM ,
2894 A LCM )
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
2948 END
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
i25glob_22
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)
Definition i25dst3_22.F:2895
i25dst3_22
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)
Definition i25dst3_22.F:53
intersecv_25
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)
Definition i25dst3_22.F:2168
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::icont_i_fi
type(int_pointer), dimension(:), allocatable icont_i_fi
Definition tri7box.F:532
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