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 . MSEGLO(*),IRTLM(4,NSN), SUBTRIA(MVSIZ),
74 . mvoisa(mvsiz,4), mvoisb(mvsiz,4),ibounda(4,mvsiz),iboundb(4,mvsiz)
75 my_real , INTENT(IN) :: dgapload ,drad
76 my_real
77 . time_s(2,nsn), gaps(*), gapm(*),
78 . pene_old(5,*),stif_old(2,nsn), penmin, eps0, marge
79 my_real
80 . xx(mvsiz,5),yy(mvsiz,5),zz(mvsiz,5),
81 . xi(mvsiz), yi(mvsiz), zi(mvsiz), stif(mvsiz),
82 . vx1(mvsiz),vx2(mvsiz),vx3(mvsiz),vx4(mvsiz),vy1(mvsiz),
83 . vy2(mvsiz),vy3(mvsiz),vy4(mvsiz),vz1(mvsiz),vz2(mvsiz),
84 . vz3(mvsiz),vz4(mvsiz),vxi(mvsiz),vyi(mvsiz),vzi(mvsiz),
85 . nax(mvsiz,5), nay(mvsiz,5), naz(mvsiz,5),
86 . nbx(mvsiz,5), nby(mvsiz,5), nbz(mvsiz,5),
87 . pent(mvsiz),
88 . lbs(mvsiz), lcs(mvsiz), lbp(mvsiz,4), lcp(mvsiz,4),
89 . gap_nm(4,mvsiz), gapmxl(mvsiz)
90 INTEGER IRECT(4,*),ITAB(*),ICONT_I(*), FAR(MVSIZ)
91 REAL*4 VTX_BISECTOR(3,2,*)
92C-----------------------------------------------
93C C o m m o n B l o c k s
94C-----------------------------------------------
95#include "com08_c.inc"
96C-----------------------------------------------
97C L o c a l V a r i a b l e s
98C-----------------------------------------------
99 INTEGER I, L, N, IT, I1, I2, ITRIA(2,4),
100 . MGLOB, ITPERM(4), IA, IB, IB1, IB2, IB3, IBX, IX, IY, IZ
101 INTEGER INTERSECTA(MVSIZ), INTERSECTB(MVSIZ),
102 . INTERSECA, INTERSECB, IKEEP,
103 . FARA(MVSIZ,4), FARB(MVSIZ,4), INGAPA(MVSIZ,4), INGAPB(MVSIZ,4),
104 . SUBTRIB(MVSIZ), SUBTRIX(MVSIZ), ICONT_R(MVSIZ)
105 my_real
106 . XIJ(MVSIZ,4), XI0V(MVSIZ), XI5, XOI5,
107 . YIJ(MVSIZ,4), YI0V(MVSIZ), YI5, YOI5,
108 . zij(mvsiz,4), zi0v(mvsiz), zi5, zoi5,
109 . x51, x52, x53, x54,
110 . y51, y52, y53, y54,
111 . z51, z52, z53, z54,
112 . xo1, xo2, xo3, xo4, xo5, xoi,
113 . yo1, yo2, yo3, yo4, yo5, yoi,
114 . zo1, zo2, zo3, zo4, zo5, zoi,
115 . vo12, vo23, vo34, vo41, pene, penax, penbx
116 my_real
117 . gap_mm(mvsiz), gap, gap21,gap22,gap23,gap24,
118 . xp, yp, zp,
119 . dx, dy, dz, lx, ld, lax, lbx, lcx, dmin,
120 . dd(mvsiz,4),
121 . lap, hla, hlb(mvsiz,4), hlc(mvsiz,4),
122 . lap1,lap2,lap3,lap4,
123 . al(mvsiz,4)
124 my_real
125 . unp,zerom,eps,unpt,zeromt,aaa,
126 . pena(mvsiz,4), penb(mvsiz,4), bb(mvsiz,4),
127 . sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4,sz1,sz2,sz3,sz4,
128 . x0n(mvsiz,4), y0n(mvsiz,4), z0n(mvsiz,4),
129 . xn(mvsiz,4),yn(mvsiz,4),zn(mvsiz,4),
130 . lb(mvsiz,4), lc(mvsiz,4),
131 . sida(mvsiz,4), sidb(mvsiz,4),
132 . x12, y12, z12,
133 . px, py, pz, pp, p1, p2, xh, yh, zh, d1, d2, d3, vx, vy, vz,
134 . ll,nn,pn,
135 . v12,v23,v34,v41,
136 . nx1, nx2, nx3, nx4,
137 . ny1, ny2, ny3, ny4,
138 . nz1, nz2, nz3, nz4,
139 . sox125,sox235,sox345,sox415,
140 . soy125,soy235,soy345,soy415,
141 . soz125,soz235,soz345,soz415,
142 .
143 . eps02,tole,
144 . xx1,yy1 ,zz1,xx2,yy2,zz2,
145 . xx3 ,yy3,zz3,xx4,yy4,
146 . zz4 ,xx5,yy5,zz5
147 DATA itria/1,2,2,3,3,4,4,1/,
148 . itperm/1,4,3,2/
149 LOGICAL :: ANY_TRIANGLE
150C-----------------------------------------------------------------------
151C
152 UNP = one + em4
153 zerom = zero - em4
154 eps = (two+half)/hundred
155 unpt = one + eps
156 zeromt = zero - eps
157 eps02=em3*em3
158C
159C initialisation (cf triangles)
160 fara(1:jlt,1:4) = 0
161 farb(1:jlt,1:4) = 0
162 pena(1:jlt,1:4) = zero
163 penb(1:jlt,1:4) = zero
164 dd(1:jlt,1:4) = zero
165C
166 any_triangle = .false.
167 DO i=1,jlt
168 IF(ix3(i)==ix4(i)) THEN
169 any_triangle = .true.
170 ENDIF
171 ENDDO
172
173 DO i=1,jlt
174C
175C For computing LBP, LCP, FAR=2, etc
176C
177C IF(STIF(I) == ZERO)CYCLE
178C
179 x0n(i,1) = xx(i,1) - xx(i,5)
180 y0n(i,1) = yy(i,1) - yy(i,5)
181 z0n(i,1) = zz(i,1) - zz(i,5)
182C
183 x0n(i,2) = xx(i,2) - xx(i,5)
184 y0n(i,2) = yy(i,2) - yy(i,5)
185 z0n(i,2) = zz(i,2) - zz(i,5)
186C
187 x0n(i,3) = xx(i,3) - xx(i,5)
188 y0n(i,3) = yy(i,3) - yy(i,5)
189 z0n(i,3) = zz(i,3) - zz(i,5)
190C
191 x0n(i,4) = xx(i,4) - xx(i,5)
192 y0n(i,4) = yy(i,4) - yy(i,5)
193 z0n(i,4) = zz(i,4) - zz(i,5)
194C
195 gap_mm(i)=fourth*(gap_nm(1,i)+gap_nm(2,i)+gap_nm(3,i)+gap_nm(4,i))
196 ENDDO
197
198 IF (any_triangle) THEN
199 DO i=1,jlt
200 IF(ix3(i)==ix4(i)) THEN
201 gap_mm(i)=gap_nm(3,i)
202 END IF
203 ENDDO
204 ENDIF
205
206C--------------------------------------------------------
207C
208C--------------------------------------------------------
209C
210C normals to triangles (recomputed here to avoid storing them)
211 DO i=1,jlt
212C
213C IF(STIF(I) == ZERO)CYCLE
214C
215 xn(i,1) = y0n(i,1)*z0n(i,2) - z0n(i,1)*y0n(i,2)
216 yn(i,1) = z0n(i,1)*x0n(i,2) - x0n(i,1)*z0n(i,2)
217 zn(i,1) = x0n(i,1)*y0n(i,2) - y0n(i,1)*x0n(i,2)
218
219 xn(i,2) = y0n(i,2)*z0n(i,3) - z0n(i,2)*y0n(i,3)
220 yn(i,2) = z0n(i,2)*x0n(i,3) - x0n(i,2)*z0n(i,3)
221 zn(i,2) = x0n(i,2)*y0n(i,3) - y0n(i,2)*x0n(i,3)
222
223 xn(i,3) = y0n(i,3)*z0n(i,4) - z0n(i,3)*y0n(i,4)
224 yn(i,3) = z0n(i,3)*x0n(i,4) - x0n(i,3)*z0n(i,4)
225 zn(i,3) = x0n(i,3)*y0n(i,4) - y0n(i,3)*x0n(i,4)
226
227 xn(i,4) = y0n(i,4)*z0n(i,1) - z0n(i,4)*y0n(i,1)
228 yn(i,4) = z0n(i,4)*x0n(i,1) - x0n(i,4)*z0n(i,1)
229 zn(i,4) = x0n(i,4)*y0n(i,1) - y0n(i,4)*x0n(i,1)
230
231 ENDDO
232C
233C--------------------------------------------------------
234C Search intersection
235C--------------------------------------------------------
236C
237 intersecta(1:jlt) = 0
238 intersectb(1:jlt) = 0
239
240 ingapa(1:jlt,1:4) = 0
241 ingapb(1:jlt,1:4) = 0
242
243 icont_r(1:jlt) = 0
244C
245 DO i=1,jlt
246C attention triangles
247
248 IF(ishel(i)/=0.AND.inacti==0)THEN
249 nx1 = xn(i,1)
250 ny1 = yn(i,1)
251 nz1 = zn(i,1)
252 aaa = gaps(i)+gap_nm(1,i)
253 aaa = aaa/sqrt(nx1*nx1+ny1*ny1+nz1*nz1)
254 nx1 = nx1 * aaa
255 ny1 = ny1 * aaa
256 nz1 = nz1 * aaa
257
258 nx2 = xn(i,2)
259 ny2 = yn(i,2)
260 nz2 = zn(i,2)
261 aaa = gaps(i)+gap_nm(2,i)
262 aaa = aaa/max(em30,sqrt(nx2*nx2+ny2*ny2+nz2*nz2))
263 nx2 = nx2 * aaa
264 ny2 = ny2 * aaa
265 nz2 = nz2 * aaa
266
267 IF(ix3(i)/=ix4(i))THEN
268 nx3 = xn(i,3)
269 ny3 = yn(i,3)
270 nz3 = zn(i,3)
271 aaa = gaps(i)+gap_nm(3,i)
272 aaa = aaa/max(em30,sqrt(nx3*nx3+ny3*ny3+nz3*nz3))
273 nx3 = nx3 * aaa
274 ny3 = ny3 * aaa
275 nz3 = nz3 * aaa
276
277 nx4 = xn(i,4)
278 ny4 = yn(i,4)
279 nz4 = zn(i,4)
280 aaa = gaps(i)+gap_nm(4,i)
281 aaa = aaa/max(em30,sqrt(nx4*nx4+ny4*ny4+nz4*nz4))
282 nx4 = nx4 * aaa
283 ny4 = ny4 * aaa
284 nz4 = nz4 * aaa
285
286 xx1 = xx(i,1) + nx1
287 yy1 = yy(i,1) + ny1
288 zz1 = zz(i,1) + nz1
289
290 xx2 = xx(i,2) + nx2
291 yy2 = yy(i,2) + ny2
292 zz2 = zz(i,2) + nz2
293
294 xx3 = xx(i,3) + nx3
295 yy3 = yy(i,3) + ny3
296 zz3 = zz(i,3) + nz3
297
298 xx4 = xx(i,4) + nx4
299 yy4 = yy(i,4) + ny4
300 zz4 = zz(i,4) + nz4
301
302 xx5 = fourth*(xx1+xx2+xx3+xx4)
303 yy5 = fourth*(yy1+yy2+yy3+yy4)
304 zz5 = fourth*(zz1+zz2+zz3+zz4)
305
306 ELSE
307 nx3 = xn(i,3)
308 ny3 = yn(i,3)
309 nz3 = zn(i,3)
310 aaa = gaps(i)+gap_nm(3,i)
311 aaa = aaa/max(em30,sqrt(nx3*nx3+ny3*ny3+nz3*nz3))
312 nx3 = nx3 * aaa
313 ny3 = ny3 * aaa
314 nz3 = nz3 * aaa
315
316 nx4 = xn(i,4)
317 ny4 = yn(i,4)
318 nz4 = zn(i,4)
319
320 xx1 = xx(i,1) + nx1
321 yy1 = yy(i,1) + ny1
322 zz1 = zz(i,1) + nz1
323
324 xx2 = xx(i,2) + nx2
325 yy2 = yy(i,2) + ny2
326 zz2 = zz(i,2) + nz2
327
328 xx3 = xx(i,3) + nx3
329 yy3 = yy(i,3) + ny3
330 zz3 = zz(i,3) + nz3
331
332
333
334 xx4 = xx3
335 yy4 = yy3
336 zz4 = zz3
337
338 xx5 = xx3
339 yy5 = yy3
340 zz5 = zz3
341 ENDIF
342 ELSE
343 xx1 = xx(i,1)
344 yy1 = yy(i,1)
345 zz1 = zz(i,1)
346
347 xx2 = xx(i,2)
348 yy2 = yy(i,2)
349 zz2 = zz(i,2)
350
351 xx3 = xx(i,3)
352 yy3 = yy(i,3)
353 zz3 = zz(i,3)
354
355 xx4 = xx(i,4)
356 yy4 = yy(i,4)
357 zz4 = zz(i,4)
358
359 xx5 = xx(i,5)
360 yy5 = yy(i,5)
361 zz5 = zz(i,5)
362 ENDIF
363
364 xi5 = xx5 - xi(i)
365 yi5 = yy5 - yi(i)
366 zi5 = zz5 - zi(i)
367
368 v12 = xn(i,1)*xi5 + yn(i,1)*yi5 + zn(i,1)*zi5
369 v23 = xn(i,2)*xi5 + yn(i,2)*yi5 + zn(i,2)*zi5
370 v34 = xn(i,3)*xi5 + yn(i,3)*yi5 + zn(i,3)*zi5
371 v41 = xn(i,4)*xi5 + yn(i,4)*yi5 + zn(i,4)*zi5
372
373 IF(ishel(i)==0)THEN
374 IF(v12 < zero .and. v23 < zero .and.
375 . v34 < zero .and. v41 < zero ) THEN
376C INTERSECTION IMPOSSIBLE
377 cycle
378 ENDIF
379 END IF
380C POSSIBLE INTERSECTION
381
382C--------------------------------------------------------
383C OLD COORDINATES
384C--------------------------------------------------------
385 xo1 = xx1 - dt1*vx1(i)
386 yo1 = yy1 - dt1*vy1(i)
387 zo1 = zz1 - dt1*vz1(i)
388C
389 xo2 = xx2 - dt1*vx2(i)
390 yo2 = yy2 - dt1*vy2(i)
391 zo2 = zz2 - dt1*vz2(i)
392c
393 xo3 = xx3 - dt1*vx3(i)
394 yo3 = yy3 - dt1*vy3(i)
395 zo3 = zz3 - dt1*vz3(i)
396C
397 xo4 = xx4 - dt1*vx4(i)
398 yo4 = yy4 - dt1*vy4(i)
399 zo4 = zz4 - dt1*vz4(i)
400
401 xoi = xi(i) - dt1*vxi(i)
402 yoi = yi(i) - dt1*vyi(i)
403 zoi = zi(i) - dt1*vzi(i)
404
405 IF(ix3(i) /= ix4(i))THEN
406 xo5 = fourth*(xo1+xo2+xo3+xo4)
407 yo5 = fourth*(yo1+yo2+yo3+yo4)
408 zo5 = fourth*(zo1+zo2+zo3+zo4)
409 ELSE
410 xo5 = xo3
411 yo5 = yo3
412 zo5 = zo3
413 ENDIF
414
415c compute volumes at previous time step
416
417 x51 = xo1 - xo5
418 y51 = yo1 - yo5
419 z51 = zo1 - zo5
420
421 x52 = xo2 - xo5
422 y52 = yo2 - yo5
423 z52 = zo2 - zo5
424
425 x53 = xo3 - xo5
426 y53 = yo3 - yo5
427 z53 = zo3 - zo5
428
429 x54 = xo4 - xo5
430 y54 = yo4 - yo5
431 z54 = zo4 - zo5
432
433 xoi5 = xo5 - xoi
434 yoi5 = yo5 - yoi
435 zoi5 = zo5 - zoi
436
437 sox125 = y51*z52 - z51*y52
438 soy125 = z51*x52 - x51*z52
439 soz125 = x51*y52 - y51*x52
440 vo12 = sox125*xoi5 + soy125*yoi5 + soz125*zoi5
441
442 sox235 = y52*z53 - z52*y53
443 soy235 = z52*x53 - x52*z53
444 soz235 = x52*y53 - y52*x53
445 vo23 = sox235*xoi5 + soy235*yoi5 + soz235*zoi5
446
447 sox345 = y53*z54 - z53*y54
448 soy345 = z53*x54 - x53*z54
449 soz345 = x53*y54 - y53*x54
450 vo34 = sox345*xoi5 + soy345*yoi5 + soz345*zoi5
451
452 sox415 = y54*z51 - z54*y51
453 soy415 = z54*x51 - x54*z51
454 soz415 = x54*y51 - y54*x51
455 vo41 = sox415*xoi5 + soy415*yoi5 + soz415*zoi5
456
457c compute intersection time (linear approximation)
458c compute coordinates at intersection time
459c (intersection time can be different for each sub-triangle)
460 IF(ishel(i)==0)THEN
461 CALL interseca_25(
462 1 ix3(i) ,ix4(i),interseca,
463 1 v12 ,v23 ,v34 ,v41 ,
464 2 vo12 ,vo23 ,vo34 ,vo41 ,xx(i,1),
465 3 yy(i,1) ,zz(i,1),xx(i,2),yy(i,2),zz(i,2),
466 4 xx(i,3) ,yy(i,3),zz(i,3),xx(i,4),yy(i,4),
467 5 zz(i,4) ,xx(i,5),yy(i,5),zz(i,5),
468 6 vxi(i) ,vyi(i) ,vzi(i) ,xo1 ,xo2 ,
469 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
470 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
471 9 zo3 ,zo4 ,zo5 ,xi(i) ,yi(i) ,
472 a zi(i) ,xoi ,yoi ,zoi ,zerom ,
473 b unp ,zeromt,unpt )
474
475 ikeep = 1
476 n = cand_n(i)
477 IF(n <= nsn) THEN
478 IF(icont_i(n) < 0) ikeep = 0
479 ELSE
480 IF(icont_i_fi(nin)%P(n-nsn) < 0) ikeep = 0
481 ENDIF
482 IF(ikeep == 1.AND.interseca == 0)THEN ! IF contact in the same side : intersection is missed and have to be checked
483 CALL intersecv0_25(
484 1 ix3(i) ,ix4(i) , interseca,
485 1 v12 ,v23 ,v34 ,v41 ,
486 2 vo12 ,vo23 ,vo34 ,vo41 ,xx(i,1),
487 3 yy(i,1) ,zz(i,1),xx(i,2),yy(i,2),zz(i,2),
488 4 xx(i,3) ,yy(i,3),zz(i,3),xx(i,4),yy(i,4),
489 5 zz(i,4) ,xx(i,5),yy(i,5),zz(i,5),
490 6 xoi ,yoi ,zoi ,xo1 ,xo2 ,
491 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
492 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
493 9 zo3 ,zo4 ,zo5 ,xi(i) ,yi(i) ,
494 a zi(i) ,zerom , unp )
495 icont_r(i) = interseca
496 ENDIF
497
498 intersecta(i)=interseca
499 ELSE
500 CALL intersecb_25(
501 1 ix3(i) ,ix4(i),interseca,intersecb,
502 1 v12 ,v23 ,v34 ,v41 ,
503 2 vo12 ,vo23 ,vo34 ,vo41 ,xx1,
504 3 yy1 ,zz1,xx2,yy2,zz2,
505 4 xx3 ,yy3,zz3,xx4,yy4,
506 5 zz4 ,xx5,yy5,zz5,
507 6 vxi(i) ,vyi(i) ,vzi(i) ,xo1 ,xo2 ,
508 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
509 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
510 9 zo3 ,zo4 ,zo5 ,xi(i) ,yi(i) ,
511 a zi(i) ,xoi ,yoi ,zoi ,zerom ,
512 b unp ,zeromt,unpt )
513 intersecta(i)=interseca
514 intersectb(i)=intersecb
515 END IF
516C
517 ENDDO
518C--------------------------------------------------------
519C#include "vectorize.inc"
520 DO i=1,jlt
521C
522C IF(STIF(I) == ZERO)CYCLE
523C
524 xi0v(i) = xx(i,5) - xi(i)
525 yi0v(i) = yy(i,5) - yi(i)
526 zi0v(i) = zz(i,5) - zi(i)
527C
528 xij(i,1) = xx(i,1) - xi(i)
529 yij(i,1) = yy(i,1) - yi(i)
530 zij(i,1) = zz(i,1) - zi(i)
531C
532 xij(i,2) = xx(i,2) - xi(i)
533 yij(i,2) = yy(i,2) - yi(i)
534 zij(i,2) = zz(i,2) - zi(i)
535C
536 xij(i,3) = xx(i,3) - xi(i)
537 yij(i,3) = yy(i,3) - yi(i)
538 zij(i,3) = zz(i,3) - zi(i)
539C
540 xij(i,4) = xx(i,4) - xi(i)
541 yij(i,4) = yy(i,4) - yi(i)
542 zij(i,4) = zz(i,4) - zi(i)
543C
544 sx1 = yi0v(i)*zij(i,1) - zi0v(i)*yij(i,1)
545 sy1 = zi0v(i)*xij(i,1) - xi0v(i)*zij(i,1)
546 sz1 = xi0v(i)*yij(i,1) - yi0v(i)*xij(i,1)
547C
548 sx2 = yi0v(i)*zij(i,2) - zi0v(i)*yij(i,2)
549 sy2 = zi0v(i)*xij(i,2) - xi0v(i)*zij(i,2)
550 sz2 = xi0v(i)*yij(i,2) - yi0v(i)*xij(i,2)
551C
552 sx3 = yi0v(i)*zij(i,3) - zi0v(i)*yij(i,3)
553 sy3 = zi0v(i)*xij(i,3) - xi0v(i)*zij(i,3)
554 sz3 = xi0v(i)*yij(i,3) - yi0v(i)*xij(i,3)
555C
556 sx4 = yi0v(i)*zij(i,4) - zi0v(i)*yij(i,4)
557 sy4 = zi0v(i)*xij(i,4) - xi0v(i)*zij(i,4)
558 sz4 = xi0v(i)*yij(i,4) - yi0v(i)*xij(i,4)
559C
560 nn = one/
561 . max(em30,sqrt(xn(i,1)*xn(i,1)+ yn(i,1)*yn(i,1)+ zn(i,1)*zn(i,1)))
562 xn(i,1)=xn(i,1)*nn
563 yn(i,1)=yn(i,1)*nn
564 zn(i,1)=zn(i,1)*nn
565 lb(i,1) = -(xn(i,1)*sx2 + yn(i,1)*sy2 + zn(i,1)*sz2)*nn
566 lc(i,1) = (xn(i,1)*sx1 + yn(i,1)*sy1 + zn(i,1)*sz1)*nn
567C
568 nn = one/
569 . max(em30,sqrt(xn(i,2)*xn(i,2)+ yn(i,2)*yn(i,2)+ zn(i,2)*zn(i,2)))
570 xn(i,2)=xn(i,2)*nn
571 yn(i,2)=yn(i,2)*nn
572 zn(i,2)=zn(i,2)*nn
573 lb(i,2) = -(xn(i,2)*sx3 + yn(i,2)*sy3 + zn(i,2)*sz3)*nn
574 lc(i,2) = (xn(i,2)*sx2 + yn(i,2)*sy2 + zn(i,2)*sz2)*nn
575C
576 nn = one/
577 . max(em30,sqrt(xn(i,3)*xn(i,3)+ yn(i,3)*yn(i,3)+ zn(i,3)*zn(i,3)))
578 xn(i,3)=xn(i,3)*nn
579 yn(i,3)=yn(i,3)*nn
580 zn(i,3)=zn(i,3)*nn
581 lb(i,3) = -(xn(i,3)*sx4 + yn(i,3)*sy4 + zn(i,3)*sz4)*nn
582 lc(i,3) = (xn(i,3)*sx3 + yn(i,3)*sy3 + zn(i,3)*sz3)*nn
583C
584 nn = one/
585 . max(em30,sqrt(xn(i,4)*xn(i,4)+ yn(i,4)*yn(i,4)+ zn(i,4)*zn(i,4)))
586 xn(i,4)=xn(i,4)*nn
587 yn(i,4)=yn(i,4)*nn
588 zn(i,4)=zn(i,4)*nn
589 lb(i,4) = -(xn(i,4)*sx1 + yn(i,4)*sy1 + zn(i,4)*sz1)*nn
590 lc(i,4) = (xn(i,4)*sx4 + yn(i,4)*sy4 + zn(i,4)*sz4)*nn
591C
592 END DO
593C--------------------------------------------------------
594C
595C#include "vectorize.inc"
596 DO i=1,jlt
597C
598C IF(STIF(I) == ZERO)CYCLE
599C
600 aaa = one/max(em30,x0n(i,1)*x0n(i,1)+y0n(i,1)*y0n(i,1)+z0n(i,1)*z0n(i,1))
601 hlc(i,1) = lc(i,1)*abs(lc(i,1))*aaa
602 hlb(i,4) = lb(i,4)*abs(lb(i,4))*aaa
603 al(i,1) = -(xi0v(i)*x0n(i,1)+yi0v(i)*y0n(i,1)+zi0v(i)*z0n(i,1))*aaa
604 al(i,1) = max(zero,min(one,al(i,1)))
605 aaa = one/max(em30,x0n(i,2)*x0n(i,2)+y0n(i,2)*y0n(i,2)+z0n(i,2)*z0n(i,2))
606 hlc(i,2) = lc(i,2)*abs(lc(i,2))*aaa
607 hlb(i,1) = lb(i,1)*abs(lb(i,1))*aaa
608 al(i,2) = -(xi0v(i)*x0n(i,2)+yi0v(i)*y0n(i,2)+zi0v(i)*z0n(i,2))*aaa
609 al(i,2) = max(zero,min(one,al(i,2)))
610 aaa = one/max(em30,x0n(i,3)*x0n(i,3)+y0n(i,3)*y0n(i,3)+z0n(i,3)*z0n(i,3))
611 hlc(i,3) = lc(i,3)*abs(lc(i,3))*aaa
612 hlb(i,2) = lb(i,2)*abs(lb(i,2))*aaa
613 al(i,3) = -(xi0v(i)*x0n(i,3)+yi0v(i)*y0n(i,3)+zi0v(i)*z0n(i,3))*aaa
614 al(i,3) = max(zero,min(one,al(i,3)))
615 aaa = one/max(em30,x0n(i,4)*x0n(i,4)+y0n(i,4)*y0n(i,4)+z0n(i,4)*z0n(i,4))
616 hlc(i,4) = lc(i,4)*abs(lc(i,4))*aaa
617 hlb(i,3) = lb(i,3)*abs(lb(i,3))*aaa
618 al(i,4) = -(xi0v(i)*x0n(i,4)+yi0v(i)*y0n(i,4)+zi0v(i)*z0n(i,4))*aaa
619 al(i,4) = max(zero,min(one,al(i,4)))
620C
621 END DO
622C--------------------------------------------------------
623 DO i=1,jlt
624C Projection outside the triangle
625 IF(lb(i,1) < -em03 .OR. lc(i,1) < -em03 .OR.
626 . lb(i,1)+lc(i,1) > one+em03) THEN
627 fara(i,1)=1
628 farb(i,1)=1
629 END IF
630 x12 = xx(i,2) - xx(i,1)
631 y12 = yy(i,2) - yy(i,1)
632 z12 = zz(i,2) - zz(i,1)
633 lbp(i,1) = lb(i,1)
634 lcp(i,1) = lc(i,1)
635 lap = one-lbp(i,1)-lcp(i,1)
636 aaa = one / max(em20,x12*x12+y12*y12+z12*z12)
637 hla= lap*abs(lap) * aaa
638 IF(lap<zero.AND.
639 + hla<=hlb(i,1).AND.hla<=hlc(i,1))THEN
640 lbp(i,1) = (xij(i,2)*x12+yij(i,2)*y12+zij(i,2)*z12) * aaa
641 lbp(i,1) = max(zero,min(one,lbp(i,1)))
642 lcp(i,1) = one - lbp(i,1)
643 ELSEIF(lbp(i,1)<zero.AND.
644 + hlb(i,1)<=hlc(i,1).AND.hlb(i,1)<=hla)THEN
645 lbp(i,1) = zero
646 lcp(i,1) = al(i,2)
647 ELSEIF(lcp(i,1)<zero.AND.
648 + hlc(i,1)<=hla.AND.hlc(i,1)<=hlb(i,1))THEN
649 lcp(i,1) = zero
650 lbp(i,1) = al(i,1)
651 ENDIF
652C Projection outside the triangle
653 IF(lb(i,2) < -em03 .OR. lc(i,2) < -em03 .OR.
654 . lb(i,2)+lc(i,2) > one+em03) THEN
655 fara(i,2)=1
656 farb(i,2)=1
657 ENDIF
658 x12 = xx(i,3) - xx(i,2)
659 y12 = yy(i,3) - yy(i,2)
660 z12 = zz(i,3) - zz(i,2)
661 lbp(i,2) = lb(i,2)
662 lcp(i,2) = lc(i,2)
663 lap = one-lbp(i,2)-lcp(i,2)
664 aaa = one / max(em20,x12*x12+y12*y12+z12*z12)
665 hla= lap*abs(lap) * aaa
666 IF(lap<zero.AND.
667 + hla<=hlb(i,2).AND.hla<=hlc(i,2))THEN
668 lbp(i,2) = (xij(i,3)*x12+yij(i,3)*y12+zij(i,3)*z12) * aaa
669 lbp(i,2) = max(zero,min(one,lbp(i,2)))
670 lcp(i,2) = one - lbp(i,2)
671 ELSEIF(lbp(i,2)<zero.AND.
672 + hlb(i,2)<=hlc(i,2).AND.hlb(i,2)<=hla)THEN
673 lbp(i,2) = zero
674 lcp(i,2) = al(i,3)
675 ELSEIF(lcp(i,2)<zero.AND.
676 + hlc(i,2)<=hla.AND.hlc(i,2)<=hlb(i,2))THEN
677 lcp(i,2) = zero
678 lbp(i,2) = al(i,2)
679 END IF
680C Projection outside the triangle
681 IF(lb(i,3) < -em03 .OR. lc(i,3) < -em03 .OR.
682 . lb(i,3)+lc(i,3) > one+em03) THEN
683 fara(i,3)=1
684 farb(i,3)=1
685 END IF
686 x12 = xx(i,4) - xx(i,3)
687 y12 = yy(i,4) - yy(i,3)
688 z12 = zz(i,4) - zz(i,3)
689 lbp(i,3) = lb(i,3)
690 lcp(i,3) = lc(i,3)
691 lap = one-lbp(i,3)-lcp(i,3)
692 aaa = one / max(em20,x12*x12+y12*y12+z12*z12)
693 hla= lap*abs(lap) * aaa
694 IF(lap<zero.AND.
695 + hla<=hlb(i,3).AND.hla<=hlc(i,3))THEN
696 lbp(i,3) = (xij(i,4)*x12+yij(i,4)*y12+zij(i,4)*z12) * aaa
697 lbp(i,3) = max(zero,min(one,lbp(i,3)))
698 lcp(i,3) = one - lbp(i,3)
699 ELSEIF(lbp(i,3)<zero.AND.
700 + hlb(i,3)<=hlc(i,3).AND.hlb(i,3)<=hla)THEN
701 lbp(i,3) = zero
702 lcp(i,3) = al(i,4)
703 ELSEIF(lcp(i,3)<zero.AND.
704 + hlc(i,3)<=hla.AND.hlc(i,3)<=hlb(i,3))THEN
705 lcp(i,3) = zero
706 lbp(i,3) = al(i,3)
707 ENDIF
708C Projection outside the triangle
709 IF(lb(i,4) < -em03 .OR. lc(i,4) < -em03 .OR.
710 . lb(i,4)+lc(i,4) > one+em03) THEN
711 fara(i,4)=1
712 farb(i,4)=1
713 END IF
714 x12 = xx(i,1) - xx(i,4)
715 y12 = yy(i,1) - yy(i,4)
716 z12 = zz(i,1) - zz(i,4)
717 lbp(i,4) = lb(i,4)
718 lcp(i,4) = lc(i,4)
719 lap = one-lbp(i,4)-lcp(i,4)
720 aaa = one / max(em20,x12*x12+y12*y12+z12*z12)
721 hla= lap*abs(lap) * aaa
722 IF(lap<zero.AND.
723 + hla<=hlb(i,4).AND.hla<=hlc(i,4))THEN
724 lbp(i,4) = (xij(i,1)*x12+yij(i,1)*y12+zij(i,1)*z12) * aaa
725 lbp(i,4) = max(zero,min(one,lbp(i,4)))
726 lcp(i,4) = one - lbp(i,4)
727 ELSEIF(lbp(i,4)<zero.AND.
728 + hlb(i,4)<=hlc(i,4).AND.hlb(i,4)<=hla)THEN
729 lbp(i,4) = zero
730 lcp(i,4) = al(i,1)
731 ELSEIF(lcp(i,4)<zero.AND.
732 + hlc(i,4)<=hla.AND.hlc(i,4)<=hlb(i,4))THEN
733 lcp(i,4) = zero
734 lbp(i,4) = al(i,4)
735 ENDIF
736 ENDDO
737
738
739 DO i = 1,jlt
740
741 ikeep = 1
742 IF(ishel(i) > 0)THEN
743 n = cand_n(i)
744 IF((intersecta(i) == 0.AND.intersectb(i)==0) .AND.n <= nsn)THEN
745 IF (icont_i(n) < 0 ) ikeep = 0
746 ELSEIF(intersecta(i) == 0.AND.intersectb(i)==0) THEN
747 IF (icont_i_fi(nin)%P(n-nsn) < 0 ) ikeep = 0
748 ENDIF
749 ENDIF
750
751 lap1 = one-lbp(i,1)-lcp(i,1)
752 xp =lap1*xx(i,5)+lbp(i,1)*xx(i,1) + lcp(i,1)*xx(i,2)
753 yp =lap1*yy(i,5)+lbp(i,1)*yy(i,1) + lcp(i,1)*yy(i,2)
754 zp =lap1*zz(i,5)+lbp(i,1)*zz(i,1) + lcp(i,1)*zz(i,2)
755 dx =xi(i)-xp
756 dy =yi(i)-yp
757 dz =zi(i)-zp
758 dd(i,1) =dx*dx+dy*dy+dz*dz
759C
760 lap2 = one-lbp(i,2)-lcp(i,2)
761 xp =lap2*xx(i,5)+lbp(i,2)*xx(i,2) + lcp(i,2)*xx(i,3)
762 yp =lap2*yy(i,5)+lbp(i,2)*yy(i,2) + lcp(i,2)*yy(i,3)
763 zp =lap2*zz(i,5)+lbp(i,2)*zz(i,2) + lcp(i,2)*zz(i,3)
764 dx =xi(i)-xp
765 dy =yi(i)-yp
766 dz =zi(i)-zp
767 dd(i,2) =dx*dx+dy*dy+dz*dz
768C
769 lap3 = one-lbp(i,3)-lcp(i,3)
770 xp =lap3*xx(i,5)+lbp(i,3)*xx(i,3) + lcp(i,3)*xx(i,4)
771 yp =lap3*yy(i,5)+lbp(i,3)*yy(i,3) + lcp(i,3)*yy(i,4)
772 zp =lap3*zz(i,5)+lbp(i,3)*zz(i,3) + lcp(i,3)*zz(i,4)
773 dx =xi(i)-xp
774 dy =yi(i)-yp
775 dz =zi(i)-zp
776 dd(i,3) =dx*dx+dy*dy+dz*dz
777C
778 lap4 = one-lbp(i,4)-lcp(i,4)
779 xp =lap4*xx(i,5)+lbp(i,4)*xx(i,4) + lcp(i,4)*xx(i,1)
780 yp =lap4*yy(i,5)+lbp(i,4)*yy(i,4) + lcp(i,4)*yy(i,1)
781 zp =lap4*zz(i,5)+lbp(i,4)*zz(i,4) + lcp(i,4)*zz(i,1)
782 dx =xi(i)-xp
783 dy =yi(i)-yp
784 dz =zi(i)-zp
785 dd(i,4) =dx*dx+dy*dy+dz*dz
786
787C
788 gap = min(max(drad,gaps(i)+lap1*gap_mm(i)+lbp(i,1)*gap_nm(1,i)+lcp(i,1)*gap_nm(2,i)+dgapload),
789 . max(drad,gapmxl(i)+dgapload))
790 gap21 = gap**2
791 bb(i,1) =((xx(i,5)-xi(i))*xn(i,1)+(yy(i,5)-yi(i))*yn(i,1)+(zz(i,5)-zi(i))*zn(i,1))
792C
793 gap = min(max(drad,gaps(i)+lap2*gap_mm(i)+lbp(i,2)*gap_nm(2,i)+lcp(i,2)*gap_nm(3,i)+dgapload),
794 . max(drad,gapmxl(i)+dgapload))
795 gap22 =gap**2
796 bb(i,2) =((xx(i,5)-xi(i))*xn(i,2)+(yy(i,5)-yi(i))*yn(i,2)+(zz(i,5)-zi(i))*zn(i,2))
797C
798 gap = min(max(drad,gaps(i)+lap3*gap_mm(i)+lbp(i,3)*gap_nm(3,i)+lcp(i,3)*gap_nm(4,i)+dgapload),
799 . max(drad,gapmxl(i)+dgapload))
800 gap23 =gap**2
801 bb(i,3) =((xx(i,5)-xi(i))*xn(i,3)+(yy(i,5)-yi(i))*yn(i,3)+(zz(i,5)-zi(i))*zn(i,3))
802C
803 gap = min(max(drad,gaps(i)+lap4*gap_mm(i)+lbp(i,4)*gap_nm(4,i)+lcp(i,4)*gap_nm(1,i)+dgapload),
804 . max(drad,gapmxl(i)+dgapload))
805 gap24 =gap**2
806 bb(i,4) =((xx(i,5)-xi(i))*xn(i,4)+(yy(i,5)-yi(i))*yn(i,4)+(zz(i,5)-zi(i))*zn(i,4))
807C
808 IF(dd(i,1) <= gap21) THEN
809 IF(bb(i,1) <= zero .AND. intersecta(i) /= 1.AND.ikeep==1)THEN
810 ingapa(i,1)=1
811 END IF
812 IF(bb(i,1) >= zero .AND. intersectb(i) /= 1.AND.ikeep==1)THEN
813 ingapb(i,1)=1
814 END IF
815 ENDIF
816 IF(dd(i,2) <= gap22) THEN
817 IF(bb(i,2) <= zero .AND. intersecta(i) /= 1.AND.ikeep==1)THEN
818 ingapa(i,2)=1
819 END IF
820 IF(bb(i,2) >= zero .AND. intersectb(i) /= 1.AND.ikeep==1)THEN
821 ingapb(i,2)=1
822 END IF
823 ENDIF
824 IF(dd(i,3) <= gap23) THEN
825 IF(bb(i,3) <= zero .AND. intersecta(i) /= 1.AND.ikeep==1)THEN
826 ingapa(i,3)=1
827 END IF
828 IF(bb(i,3) >= zero .AND. intersectb(i) /= 1.AND.ikeep==1)THEN
829 ingapb(i,3)=1
830 END IF
831 ENDIF
832 IF(dd(i,4) <= gap24) THEN
833 IF(bb(i,4) <= zero .AND. intersecta(i) /= 1.AND.ikeep==1)THEN
834 ingapa(i,4)=1
835 END IF
836 IF(bb(i,4) >= zero .AND. intersectb(i) /= 1.AND.ikeep==1)THEN
837 ingapb(i,4)=1
838 END IF
839 END IF
840 END DO
841C--------------------------------------------------------
842C Determination of the TQ sub-triangle its cone contains the ND does not require Far!
843C--------------------------------------------------------
844 subtria(1:jlt)=0
845 subtrib(1:jlt)=0
846 subtrix(1:jlt)=0
847 DO i=1,jlt
848C
849 IF(stif(i) <= zero)cycle
850C
851 IF(ix3(i)/=ix4(i))THEN
852
853 dmin=min(dd(i,1),dd(i,2),dd(i,3),dd(i,4))
854 ld=ep20
855 DO it=1,4
856 IF(dd(i,it) <= onep03*dmin)THEN
857 lbx = lb(i,it)
858 lcx = lc(i,it)
859 lax = one-lb(i,it)-lc(i,it)
860C
861C Privilegier sector in which we are.
862 IF(lbx >= zero .AND. lcx >= zero)THEN
863 lx=zero
864 ELSE
865 ! closest point in the angular sector == center
866 lx = max(zero,dd(i,it)-bb(i,it)*bb(i,it))
867 END IF
868
869 IF(lx < ld) THEN
870 subtrix(i)= it
871 ld = lx
872 END IF
873 END IF
874 END DO
875C
876C !
877 it=subtrix(i)
878 IF(it > 0) THEN
879 IF(intersecta(i)/=0.OR.ingapa(i,it)/=0)subtria(i)=it
880 IF (ishel(i) /= 0)THEN
881 IF(intersectb(i)/=0.OR.ingapb(i,it)/=0)subtrib(i)=it
882 END IF
883 ENDIF
884 ELSE
885
886 IF(intersecta(i)/=0.OR.ingapa(i,1)/=0)THEN
887 subtria(i)= 1
888 END IF
889C
890 IF (ishel(i) /= 0)THEN
891 IF(intersectb(i)/=0.OR.ingapb(i,1)/=0)THEN
892 subtrib(i)= 1
893 END IF
894 END IF
895 END IF
896C
897 END DO
898C--------------------------------------------------------
899C calculation of penetration
900C--------------------------------------------------------
901#include "vectorize.inc"
902 DO i =1,jlt
903C
904C IF (STIF(I) == ZERO)CYCLE
905C
906 it=subtria(i)
907 IF(it==0)cycle
908C
909 i1=itria(1,it)
910 i2=itria(2,it)
911C
912 lap = one-lbp(i,it)-lcp(i,it)
913 gap = min(max(drad,gaps(i)+lap*gap_mm(i)+lbp(i,it)*gap_nm(i1,i)+lcp(i,it)*gap_nm(i2,i))+dgapload,
914 . max(drad,gapmxl(i)+dgapload))
915C
916 IF(bb(i,it) > zero)THEN
917C
918C penetration approximee (cf zone limite interpolation des normales)
919 pena(i,it)=max(zero,gap+bb(i,it))
920 ELSE ! IF(BB(I,IT) > ZERO)THEN
921 pena(i,it)=max(zero,gap-sqrt(dd(i,it)))
922 END IF
923C
924C--------exclude high pene in case of ICONT_R
925 IF(icont_r(i) >0)THEN
926 aaa = max(em30,x0n(i,it)*x0n(i,1)+y0n(i,it)*y0n(i,it)+z0n(i,it)*z0n(i,it))
927 tole =eps02*aaa
928 IF (gap >zero) tole =min(tole,gap*gap)
929 IF(pena(i,it)*pena(i,it) > tole) THEN
930 pena(i,it) = zero
931 END IF
932 END if!(ICONT_R(I) >0)THEN
933C
934C
935 ENDDO
936C--------------------------------------------------------
937C Check vs bissectors ...
938C--------------------------------------------------------
939#include "vectorize.inc"
940 DO i =1,jlt
941C
942C IF (STIF(I) == ZERO)CYCLE
943C
944 it=subtria(i)
945 IF(it==0)cycle
946C
947 IF(pena(i,it)==zero) cycle
948C
949 i1=itria(1,it)
950 i2=itria(2,it)
951C
952 xh=xi(i)+bb(i,it)*xn(i,it)
953 yh=yi(i)+bb(i,it)*yn(i,it)
954 zh=zi(i)+bb(i,it)*zn(i,it)
955C
956 IF(ix3(i) /= ix4(i))THEN
957C
958 ib1=ibounda(i1,i)
959 ib2=ibounda(i2,i)
960 IF(mvoisa(i,it)==0)THEN
961C
962 IF( (xh-xx(i,i1))* nax(i,it)+
963 . (yh-yy(i,i1))* nay(i,it)+
964 . (zh-zz(i,i1))* naz(i,it) >= gaps(i)) fara(i,it) =3
965C
966 ELSEIF((ib1 /= 0 .AND. ib2 == 0).OR.
967 . (ib2 /= 0 .AND. ib1 == 0))THEN
968C
969 ibx=max(ib1,ib2)
970C
971 ibx=max(ib1,ib2)
972 IF(ib1/=0)THEN
973 ix =i1
974 ELSEIF(ib2/=0)THEN
975 ix =i2
976 END IF
977C
978 xh=xi(i)+bb(i,it)*xn(i,it)
979 yh=yi(i)+bb(i,it)*yn(i,it)
980 zh=zi(i)+bb(i,it)*zn(i,it)
981C
982 IF(vtx_bisector(1,1,ibx)/=zero.OR.
983 . vtx_bisector(2,1,ibx)/=zero.OR.
984 . vtx_bisector(3,1,ibx)/=zero.OR.
985 . vtx_bisector(1,2,ibx)/=zero.OR.
986 . vtx_bisector(2,2,ibx)/=zero.OR.
987 . vtx_bisector(3,2,ibx)/=zero)THEN
988 p1 = (xh-xx(i,ix))* vtx_bisector(1,1,ibx)+
989 . (yh-yy(i,ix))* vtx_bisector(2,1,ibx)+
990 . (zh-zz(i,ix))* vtx_bisector(3,1,ibx)
991 p2 = (xh-xx(i,ix))* vtx_bisector(1,2,ibx)+
992 . (yh-yy(i,ix))* vtx_bisector(2,2,ibx)+
993 . (zh-zz(i,ix))* vtx_bisector(3,2,ibx)
994 IF(p1 >= gaps(i) .AND. p2 >= gaps(i)) fara(i,it) =3
995 ELSE
996 vx = x0n(i,ix) ! fake bisector of angle at vertex IX
997 vy = y0n(i,ix)
998 vz = z0n(i,ix)
999 nn = one/max(em20,sqrt(vx*vx+vy*vy+vz*vz))
1000 pn = ((xh-xx(i,ix))*vx+(yh-yy(i,ix))*vy+(zh-zz(i,ix))*vz)*nn
1001 IF(pn >= gaps(i)) fara(i,it) =3
1002 END IF
1003
1004 END IF
1005C
1006C Cone
1007 IF(ingapa(i,it) == 1 .AND. (fara(i,it)==1 .OR. bb(i,it) <= zero))THEN
1008C ingap = 1 => always check if inside the cone (allows to determine which cone we are in)
1009C
1010 x12=xx(i,i2)-xx(i,i1)
1011 y12=yy(i,i2)-yy(i,i1)
1012 z12=zz(i,i2)-zz(i,i1)
1013C
1014C normal to the bisecting plane (pointing toward the inside)
1015 px = z12*nay(i,it)-y12*naz(i,it)
1016 py = x12*naz(i,it)-z12*nax(i,it)
1017 pz = y12*nax(i,it)-x12*nay(i,it)
1018 pp = px*px+py*py+pz*pz
1019C
1020 ll = xij(i,i1)*xij(i,i1)+yij(i,i1)*yij(i,i1)+zij(i,i1)*zij(i,i1)
1021C
1022 sida(i,it)=-(xij(i,i1)*px+yij(i,i1)*py+zij(i,i1)*pz)*sqrt(one/max(em30,ll*pp))
1023 IF(sida(i,it) < -zep01) fara(i,it) =2
1024C
1025 END IF
1026 ELSE
1027 ib1=ibounda(1,i)
1028 ib2=ibounda(2,i)
1029 ib3=ibounda(3,i)
1030C
1031 d1=(xh-xx(i,1))* nax(i,1)+
1032 . (yh-yy(i,1))* nay(i,1)+
1033 . (zh-zz(i,1))* naz(i,1)
1034 d2=(xh-xx(i,2))* nax(i,2)+
1035 . (yh-yy(i,2))* nay(i,2)+
1036 . (zh-zz(i,2))* naz(i,2)
1037 d3=(xh-xx(i,3))* nax(i,4)+
1038 . (yh-yy(i,3))* nay(i,4)+
1039 . (zh-zz(i,3))* naz(i,4)
1040C
1041 IF( (mvoisa(i,1) == 0 .AND. d1 >= gaps(i)).OR.
1042 . (mvoisa(i,2) == 0 .AND. d2 >= gaps(i)).OR.
1043 . (mvoisa(i,4) == 0 .AND. d3 >= gaps(i)) )THEN
1044C
1045 fara(i,it) =3
1046C
1047 ELSEIF((ib1/=0 .AND. ib2==0 .AND. ib3==0).OR.
1048 . (ib2/=0 .AND. ib3==0 .AND. ib1==0).OR.
1049 . (ib3/=0 .AND. ib1==0 .AND. ib2==0))THEN
1050C
1051 ibx=max(ib1,ib2,ib3)
1052 IF(ib1/=0)THEN
1053 ix =1
1054 iy =2
1055 iz =3
1056 ELSEIF(ib2/=0)THEN
1057 ix =2
1058 iy =3
1059 iz =1
1060 ELSEIF(ib3/=0)THEN
1061 ix =3
1062 iy =1
1063 iz =2
1064 END IF
1065C
1066 IF(vtx_bisector(1,1,ibx)/=zero.OR.
1067 . vtx_bisector(2,1,ibx)/=zero.OR.
1068 . vtx_bisector(3,1,ibx)/=zero.OR.
1069 . vtx_bisector(1,2,ibx)/=zero.OR.
1070 . vtx_bisector(2,2,ibx)/=zero.OR.
1071 . vtx_bisector(3,2,ibx)/=zero)THEN
1072 p1 = (xh-xx(i,ix))* vtx_bisector(1,1,ibx)+
1073 . (yh-yy(i,ix))* vtx_bisector(2,1,ibx)+
1074 . (zh-zz(i,ix))* vtx_bisector(3,1,ibx)
1075 p2 = (xh-xx(i,ix))* vtx_bisector(1,2,ibx)+
1076 . (yh-yy(i,ix))* vtx_bisector(2,2,ibx)+
1077 . (zh-zz(i,ix))* vtx_bisector(3,2,ibx)
1078 IF(p1 >= gaps(i) .AND. p2 >= gaps(i)) fara(i,it) =3
1079 ELSE
1080 vx = two*xx(i,ix)-(xx(i,iy)+xx(i,iz)) ! fake bisector of angle 2,1,3
1081 vy = two*yy(i,ix)-(yy(i,iy)+yy(i,iz))
1082 vz = two*zz(i,ix)-(zz(i,iy)+zz(i,iz))
1083 nn = one/max(em20,sqrt(vx*vx+vy*vy+vz*vz))
1084 pn = ((xh-xx(i,ix))*vx+(yh-yy(i,ix))*vy+(zh-zz(i,ix))*vz)*nn
1085 IF(pn >= gaps(i)) fara(i,it) =3
1086 END IF
1087C
1088 END IF
1089C
1090 IF(ingapa(i,it) == 1 .AND. (fara(i,it)==1 .OR. bb(i,it) <= zero))THEN
1091C
1092C ingap = 1 => always check if inside the cone (allows to determine which cone we are in)
1093C
1094 IF( mvoisa(i,1) /= 0 )THEN
1095C
1096 x12=xx(i,2)-xx(i,1)
1097 y12=yy(i,2)-yy(i,1)
1098 z12=zz(i,2)-zz(i,1)
1099C
1100C normal to the bisecting plane (pointing toward the inside)
1101 px = z12*nay(i,1)-y12*naz(i,1)
1102 py = x12*naz(i,1)-z12*nax(i,1)
1103 pz = y12*nax(i,1)-x12*nay(i,1)
1104 pp = px*px+py*py+pz*pz
1105C
1106 ll = xij(i,1)*xij(i,1)+yij(i,1)*yij(i,1)+zij(i,1)*zij(i,1)
1107C
1108 sida(i,1)=-(xij(i,1)*px+yij(i,1)*py+zij(i,1)*pz)*sqrt(one/max(em30,ll*pp))
1109 IF(sida(i,1) < -zep01) fara(i,it) =2
1110C
1111 END IF
1112
1113 IF( mvoisa(i,2) /= 0 )THEN
1114
1115C ingap = 1 => always check if inside the cone (allows to determine which cone we are in)
1116C
1117 x12=xx(i,3)-xx(i,2)
1118 y12=yy(i,3)-yy(i,2)
1119 z12=zz(i,3)-zz(i,2)
1120C
1121C normal to the bisecting plane (pointing toward the inside)
1122 px = z12*nay(i,2)-y12*naz(i,2)
1123 py = x12*naz(i,2)-z12*nax(i,2)
1124 pz = y12*nax(i,2)-x12*nay(i,2)
1125 pp = px*px+py*py+pz*pz
1126C
1127 ll = xij(i,2)*xij(i,2)+yij(i,2)*yij(i,2)+zij(i,2)*zij(i,2)
1128C
1129 sida(i,2)=-(xij(i,2)*px+yij(i,2)*py+zij(i,2)*pz)*sqrt(one/max(em30,ll*pp))
1130 IF(sida(i,2) < -zep01) fara(i,it) =2
1131C
1132 END IF
1133
1134 IF( mvoisa(i,4) /= 0 )THEN
1135C ingap = 1 => always check if inside the cone (allows to determine which cone we are in)
1136C
1137 x12=xx(i,1)-xx(i,3)
1138 y12=yy(i,1)-yy(i,3)
1139 z12=zz(i,1)-zz(i,3)
1140C
1141C normal to the bisecting plane (pointing toward the inside)
1142 px = z12*nay(i,4)-y12*naz(i,4)
1143 py = x12*naz(i,4)-z12*nax(i,4)
1144 pz = y12*nax(i,4)-x12*nay(i,4)
1145 pp = px*px+py*py+pz*pz
1146C
1147 ll = xij(i,3)*xij(i,3)+yij(i,3)*yij(i,3)+zij(i,3)*zij(i,3)
1148C
1149 sida(i,4)=-(xij(i,3)*px+yij(i,3)*py+zij(i,3)*pz)*sqrt(one/max(em30,ll*pp))
1150 IF(sida(i,4) < -zep01) fara(i,it) =2
1151C
1152 END IF
1153 END IF
1154 END IF
1155C
1156 IF(fara(i,it)==2 .AND. intersecta(i)==0) pena(i,it)=zero
1157 IF(fara(i,it)==3) pena(i,it)=zero
1158C
1159 ENDDO
1160C--------------------------------------------------------
1161C calculation of penetration
1162C--------------------------------------------------------
1163#include "vectorize.inc"
1164 DO i =1,jlt
1165C
1166C IF (STIF(I) == ZERO)CYCLE
1167C
1168 it=subtrib(i)
1169 IF(it==0)cycle
1170C
1171 i1=itria(1,it)
1172 i2=itria(2,it)
1173C
1174 lap = one-lbp(i,it)-lcp(i,it)
1175 gap = min(max(drad,gaps(i)+lap*gap_mm(i)+lbp(i,it)*gap_nm(i1,i)+lcp(i,it)*gap_nm(i2,i)+dgapload),
1176 . max(drad,gapmxl(i)+dgapload))
1177C
1178 IF(bb(i,it) < zero)THEN
1179C
1180C penetration approximee (cf zone limite interpolation des normales)
1181 penb(i,it)=max(zero,gap-bb(i,it))
1182 ELSE ! IF(BB(I,IT) < ZERO)THEN
1183 penb(i,it)=max(zero,gap-sqrt(dd(i,it)))
1184 END IF
1185C
1186 ENDDO
1187C--------------------------------------------------------
1188C Check vs bissectors ...
1189C--------------------------------------------------------
1190#include "vectorize.inc"
1191 DO i =1,jlt
1192C
1193C IF (STIF(I) == ZERO)CYCLE
1194C
1195 it=subtrib(i)
1196 IF(it==0)cycle
1197C
1198 IF(penb(i,it)==zero) cycle
1199C
1200 i1=itria(1,it)
1201 i2=itria(2,it)
1202C
1203 xh=xi(i)+bb(i,it)*xn(i,it)
1204 yh=yi(i)+bb(i,it)*yn(i,it)
1205 zh=zi(i)+bb(i,it)*zn(i,it)
1206C
1207 IF(ix3(i) /= ix4(i))THEN
1208C
1209C Cone
1210 ib1=iboundb(i1,i)
1211 ib2=iboundb(i2,i)
1212 IF(mvoisb(i,it)==0)THEN
1213C
1214 IF( (xh-xx(i,i1))* nbx(i,it)+
1215 . (yh-yy(i,i1))* nby(i,it)+
1216 . (zh-zz(i,i1))* nbz(i,it) >= gaps(i)) farb(i,it) =3
1217C
1218 ELSEIF((ib1 /= 0 .AND. ib2 == 0).OR.
1219 . (ib2 /= 0 .AND. ib1 == 0))THEN
1220C
1221 ibx=max(ib1,ib2)
1222 IF(ib1/=0)THEN
1223 ix =i1
1224 ELSEIF(ib2/=0)THEN
1225 ix =i2
1226 END IF
1227C
1228 xh=xi(i)+bb(i,it)*xn(i,it)
1229 yh=yi(i)+bb(i,it)*yn(i,it)
1230 zh=zi(i)+bb(i,it)*zn(i,it)
1231C
1232 IF(vtx_bisector(1,1,ibx)/=zero.OR.
1233 . vtx_bisector(2,1,ibx)/=zero.OR.
1234 . vtx_bisector(3,1,ibx)/=zero.OR.
1235 . vtx_bisector(1,2,ibx)/=zero.OR.
1236 . vtx_bisector(2,2,ibx)/=zero.OR.
1237 . vtx_bisector(3,2,ibx)/=zero)THEN
1238 p1 = (xh-xx(i,ix))* vtx_bisector(1,1,ibx)+
1239 . (yh-yy(i,ix))* vtx_bisector(2,1,ibx)+
1240 . (zh-zz(i,ix))* vtx_bisector(3,1,ibx)
1241 p2 = (xh-xx(i,ix))* vtx_bisector(1,2,ibx)+
1242 . (yh-yy(i,ix))* vtx_bisector(2,2,ibx)+
1243 . (zh-zz(i,ix))* vtx_bisector(3,2,ibx)
1244 IF(p1 >= gaps(i) .AND. p2 >= gaps(i)) farb(i,it) =3
1245 ELSE
1246 vx = x0n(i,ix) ! fake bisector of angle at vertex IX
1247 vy = y0n(i,ix)
1248 vz = z0n(i,ix)
1249 nn = one/max(em20,sqrt(vx*vx+vy*vy+vz*vz))
1250 pn = ((xh-xx(i,ix))*vx+(yh-yy(i,ix))*vy+(zh-zz(i,ix))*vz)*nn
1251 IF(pn >= gaps(i)) farb(i,it) =3
1252 END IF
1253
1254 END IF
1255C
1256C
1257 IF(ingapb(i,it) == 1 .AND. (farb(i,it)==1 .OR. bb(i,it) >= zero))THEN
1258C
1259C ingap = 1 => always check if inside the cone (allows to determine which cone we are in)
1260C
1261 x12=xx(i,i2)-xx(i,i1)
1262 y12=yy(i,i2)-yy(i,i1)
1263 z12=zz(i,i2)-zz(i,i1)
1264C
1265C normal to the bisecting plane (pointing toward the inside)
1266 px = z12*nby(i,it)-y12*nbz(i,it)
1267 py = x12*nbz(i,it)-z12*nbx(i,it)
1268 pz = y12*nbx(i,it)-x12*nby(i,it)
1269 pp = px*px+py*py+pz*pz
1270C
1271 ll = xij(i,i1)*xij(i,i1)+yij(i,i1)*yij(i,i1)+zij(i,i1)*zij(i,i1)
1272C
1273 sidb(i,it)= (xij(i,i1)*px+yij(i,i1)*py+zij(i,i1)*pz)*sqrt(one/max(em30,ll*pp))
1274 IF(sidb(i,it) < -zep01) farb(i,it) =2
1275C
1276 END IF
1277 ELSE
1278 ib1=iboundb(1,i)
1279 ib2=iboundb(2,i)
1280 ib3=iboundb(3,i)
1281C
1282 d1=(xh-xx(i,1))* nbx(i,1)+
1283 . (yh-yy(i,1))* nby(i,1)+
1284 . (zh-zz(i,1))* nbz(i,1)
1285 d2=(xh-xx(i,2))* nbx(i,2)+
1286 . (yh-yy(i,2))* nby(i,2)+
1287 . (zh-zz(i,2))* nbz(i,2)
1288 d3=(xh-xx(i,3))* nbx(i,4)+
1289 . (yh-yy(i,3))* nby(i,4)+
1290 . (zh-zz(i,3))* nbz(i,4)
1291C
1292 IF( (mvoisb(i,1) == 0 .AND. d1 >= gaps(i)).OR.
1293 . (mvoisb(i,2) == 0 .AND. d2 >= gaps(i)).OR.
1294 . (mvoisb(i,4) == 0 .AND. d3 >= gaps(i)) )THEN
1295C
1296 farb(i,it) =3
1297C
1298 ELSEIF((ib1/=0 .AND. ib2==0 .AND. ib3==0).OR.
1299 . (ib2/=0 .AND. ib3==0 .AND. ib1==0).OR.
1300 . (ib3/=0 .AND. ib1==0 .AND. ib2==0))THEN
1301C
1302 ibx=max(ib1,ib2,ib3)
1303 IF(ib1/=0)THEN
1304 ix =1
1305 iy =2
1306 iz =3
1307 ELSEIF(ib2/=0)THEN
1308 ix =2
1309 iy =3
1310 iz =1
1311 ELSEIF(ib3/=0)THEN
1312 ix =3
1313 iy =1
1314 iz =2
1315 END IF
1316C
1317 IF(vtx_bisector(1,1,ibx)/=zero.OR.
1318 . vtx_bisector(2,1,ibx)/=zero.OR.
1319 . vtx_bisector(3,1,ibx)/=zero.OR.
1320 . vtx_bisector(1,2,ibx)/=zero.OR.
1321 . vtx_bisector(2,2,ibx)/=zero.OR.
1322 . vtx_bisector(3,2,ibx)/=zero)THEN
1323 p1 = (xh-xx(i,ix))* vtx_bisector(1,1,ibx)+
1324 . (yh-yy(i,ix))* vtx_bisector(2,1,ibx)+
1325 . (zh-zz(i,ix))* vtx_bisector(3,1,ibx)
1326 p2 = (xh-xx(i,ix))* vtx_bisector(1,2,ibx)+
1327 . (yh-yy(i,ix))* vtx_bisector(2,2,ibx)+
1328 . (zh-zz(i,ix))* vtx_bisector(3,2,ibx)
1329 IF(p1 >= gaps(i) .AND. p2 >= gaps(i)) farb(i,it) =3
1330 ELSE
1331 vx = two*xx(i,ix)-(xx(i,iy)+xx(i,iz)) ! fake bisector of angle 2,1,3
1332 vy = two*yy(i,ix)-(yy(i,iy)+yy(i,iz))
1333 vz = two*zz(i,ix)-(zz(i,iy)+zz(i,iz))
1334 nn = one/max(em20,sqrt(vx*vx+vy*vy+vz*vz))
1335 pn = ((xh-xx(i,ix))*vx+(yh-yy(i,ix))*vy+(zh-zz(i,ix))*vz)*nn
1336 IF(pn >= gaps(i)) farb(i,it) =3
1337 END IF
1338C
1339 END IF
1340C
1341 IF(ingapb(i,it) == 1 .AND. (farb(i,it)==1 .OR. bb(i,it) >= zero))THEN
1342C
1343C ingap = 1 => always check if inside the cone (allows to determine which cone we are in)
1344C
1345 IF( mvoisb(i,1) /= 0 )THEN
1346C
1347 x12=xx(i,2)-xx(i,1)
1348 y12=yy(i,2)-yy(i,1)
1349 z12=zz(i,2)-zz(i,1)
1350C
1351C normal to the bisecting plane (pointing toward the inside)
1352 px = z12*nby(i,1)-y12*nbz(i,1)
1353 py = x12*nbz(i,1)-z12*nbx(i,1)
1354 pz = y12*nbx(i,1)-x12*nby(i,1)
1355 pp = px*px+py*py+pz*pz
1356C
1357 ll = xij(i,1)*xij(i,1)+yij(i,1)*yij(i,1)+zij(i,1)*zij(i,1)
1358C
1359 sidb(i,1)= (xij(i,1)*px+yij(i,1)*py+zij(i,1)*pz)*sqrt(one/max(em30,ll*pp))
1360 IF(sidb(i,1) < -zep01) farb(i,it) =2
1361C
1362 END IF
1363
1364 IF( mvoisb(i,2) /= 0 )THEN
1365C
1366 x12=xx(i,3)-xx(i,2)
1367 y12=yy(i,3)-yy(i,2)
1368 z12=zz(i,3)-zz(i,2)
1369C
1370C normal to the bisecting plane (pointing toward the inside)
1371 px = z12*nby(i,2)-y12*nbz(i,2)
1372 py = x12*nbz(i,2)-z12*nbx(i,2)
1373 pz = y12*nbx(i,2)-x12*nby(i,2)
1374 pp = px*px+py*py+pz*pz
1375C
1376 ll = xij(i,2)*xij(i,2)+yij(i,2)*yij(i,2)+zij(i,2)*zij(i,2)
1377C
1378 sidb(i,2)= (xij(i,2)*px+yij(i,2)*py+zij(i,2)*pz)*sqrt(one/max(em30,ll*pp))
1379 IF(sidb(i,2) < -zep01) farb(i,it) =2
1380C
1381 END IF
1382
1383 IF( mvoisb(i,4) /= 0 )THEN
1384C
1385 x12=xx(i,1)-xx(i,3)
1386 y12=yy(i,1)-yy(i,3)
1387 z12=zz(i,1)-zz(i,3)
1388C
1389C normal to the bisecting plane (pointing toward the inside)
1390 px = z12*nby(i,4)-y12*nbz(i,4)
1391 py = x12*nbz(i,4)-z12*nbx(i,4)
1392 pz = y12*nbx(i,4)-x12*nby(i,4)
1393 pp = px*px+py*py+pz*pz
1394C
1395 ll = xij(i,3)*xij(i,3)+yij(i,3)*yij(i,3)+zij(i,3)*zij(i,3)
1396C
1397 sidb(i,4)= (xij(i,3)*px+yij(i,3)*py+zij(i,3)*pz)*sqrt(one/max(em30,ll*pp))
1398 IF(sidb(i,4) < -zep01) farb(i,it) =2
1399C
1400 END IF
1401 END IF
1402 END IF
1403C
1404 IF(farb(i,it)==2 .AND. intersectb(i)==0) penb(i,it)=zero
1405 IF(farb(i,it)==3) penb(i,it)=zero
1406C
1407 ENDDO
1408C--------------------------------------------------------
1409 DO i=1,jlt
1410C
1411 IF(stif(i) <= zero)cycle
1412C
1413 ia = subtria(i)
1414 penax = zero
1415 IF(ia/=0)penax=pena(i,ia)
1416C
1417 ib = subtrib(i)
1418 penbx = zero
1419 IF(ib/=0)penbx=penb(i,ib)
1420C
1421 IF(penax > penbx .AND. ia > 0)THEN
1422 l = cand_e(i)
1423 it = ia
1424 pent(i)= penax
1425
1426 lbs(i)=lbp(i,it)
1427 lcs(i)=lcp(i,it)
1428 far(i)=fara(i,it)
1429
1430 ELSEIF(penbx > penax .AND. ib > 0)THEN
1431 l = ishel(i)
1432 cand_e(i) = l
1433
1434 it = ib
1435 pent(i) = penbx
1436
1437 lbs(i)=lcp(i,it)
1438 lcs(i)=lbp(i,it)
1439 far(i)=farb(i,it)
1440
1441 subtria(i)=itperm(it) ! IT/=0
1442 it = subtria(i)
1443
1444 ELSE
1445 it=0
1446 subtria(i) = 0
1447 pent(i) = zero
1448 END IF
1449
1450c if(itab(nsvg(i))==2421446.or.
1451c . itab(ix1(i))==2421446.or.itab(ix2(i))==2421446.or.itab(ix3(i))==2421446.or.itab(ix4(i))==2421446)
1452c . print *,'dst2 nat',ispmd+1,
1453c . itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),
1454c . it,ia,ib,penax,penbx,
1455c . ingapa(i,1:4),intersecta(i),ingapb(i,1:4),intersectb(i),
1456c . bb(i,1:4),dd(i,1:4)
1457
1458 IF(it == 0)cycle
1459 pene = pent(i)
1460
1461 n = cand_n(i)
1462
1463 mglob=mseglo(l)
1464
1465 IF(n<=nsn)THEN
1466#include "lockon.inc"
1467
1468 IF( time_s(1,n) < pene .OR.
1469 * (time_s(1,n) == pene .AND. -irtlm(1,n) < mglob ))THEN
1470 irtlm(1,n) = -mglob
1471 irtlm(2,n) = it
1472 irtlm(3,n) = cand_e(i)
1473 irtlm(4,n) = ispmd+1
1474 time_s(1,n)= pene
1475C
1476C PENE_OLD(5,N)=PENE
1477 END IF
1478c if(itab(nsvg(i))==27363)then
1479cc if(itab(nsvg(i))==27952.or.
1480cc . itab(ix1(i))==27952.or.itab(ix2(i))==27952.or.itab(ix3(i))==27952.or.itab(ix4(i))==27952)then
1481c print *,'dst2 nat',ispmd+1,
1482c . itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),
1483c . irtlm(1,n),cand_e(i),TIME_S(1,N),
1484c . it,pent(i),penax,penbx,intersecta(i),intersectb(i),lbs(i),lcs(i),lbs(i)+lcs(i)
1485c if(ia/=0)
1486c . print *,'ia',ia,bb(i,ia),ingapa(i,ia),fara(i,ia),ibounda(1:4,i)
1487c if(ib/=0)
1488c . print *,'ib',ib,bb(i,ib),ingapb(i,ib),farb(i,ib),iboundb(1:4,i)
1489c end if
1490#include "lockoff.inc"
1491 ELSE
1492#include "lockon.inc"
1493 IF( time_sfi(nin)%P(2*(n-nsn-1)+1) < pene .OR.
1494 * (time_sfi(nin)%P(2*(n-nsn-1)+1) == pene .AND. -irtlm_fi(nin)%P(1,n-nsn) < mglob ))THEN
1495 irtlm_fi(nin)%P(1,n-nsn) = -mglob
1496 irtlm_fi(nin)%P(2,n-nsn) = it
1497 irtlm_fi(nin)%P(3,n-nsn) = cand_e(i)
1498 irtlm_fi(nin)%P(4,n-nsn) = ispmd+1
1499 time_sfi(nin)%P(2*(n-nsn-1)+1) = pene
1500C
1501C PENE_OLDFI(NIN)%P(5,N-NSN)=PENE
1502 END IF
1503c if(itafi(nin)%p(n-nsn)==3817238)
1504c . print *,'dst2 rem',ispmd+1,
1505c . itafi(nin)%p(n-nsn),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),
1506c . it,IRTLM_FI(NIN)%P(1,n-nsn),mglob,cand_e(i),pent(i,it),far(i,it),
1507c . TIME_SFI(NIN)%P(N-NSN),dd(i,it),ingap(i,it),intersect(i)
1508#include "lockoff.inc"
1509 END IF
1510 END DO
1511
1512 RETURN
1513 END
1514!||====================================================================
1515!|| interseca_25 ../engine/source/interfaces/int25/i25dst3_22.F
1516!||--- called by ------------------------------------------------------
1517!|| i25dst3_22 ../engine/source/interfaces/int25/i25dst3_22.F
1518!||====================================================================
1519 SUBROUTINE interseca_25(
1520 1 IXX3 ,IXX4 ,INTERSECA,
1521 1 V12 ,V23 ,V34 ,V41 ,
1522 2 VO12 ,VO23 ,VO34 ,VO41 ,XX1 ,
1523 3 YY1 ,ZZ1 ,XX2 ,YY2 ,ZZ2 ,
1524 4 XX3 ,YY3 ,ZZ3 ,XX4 ,YY4 ,
1525 5 ZZ4 ,XX5 ,YY5 ,ZZ5 ,
1526 6 VXI ,VYI ,VZI ,XO1 ,XO2 ,
1527 7 XO3 ,XO4 ,XO5 ,YO1 ,YO2 ,
1528 8 YO3 ,YO4 ,YO5 ,ZO1 ,ZO2 ,
1529 9 ZO3 ,ZO4 ,ZO5 ,XI ,YI ,
1530 A ZI ,XOI ,YOI ,ZOI ,ZEROM ,
1531 B UNP ,ZEROMT,UNPT )
1532C-----------------------------------------------
1533C I m p l i c i t T y p e s
1534C-----------------------------------------------
1535#include "implicit_f.inc"
1536C-----------------------------------------------
1537C D u m m y A r g u m e n t s
1538C-----------------------------------------------
1539 INTEGER IXX3 ,IXX4 , INTERSECA
1540C REAL
1541 my_real
1542 1 v12 ,v23 ,v34 ,v41 ,
1543 2 vo12 ,vo23 ,vo34 ,vo41 ,xx1 ,
1544 3 yy1 ,zz1 ,xx2 ,yy2 ,zz2 ,
1545 4 xx3 ,yy3 ,zz3 ,xx4 ,yy4 ,
1546 5 zz4 ,xx5 ,yy5 ,zz5 ,
1547 6 vxi ,vyi ,vzi ,xo1 ,xo2 ,
1548 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
1549 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
1550 9 zo3 ,zo4 ,zo5 ,zerom ,unp ,
1551 a zeromt ,unpt ,xi ,yi ,zi ,
1552 b xoi ,yoi ,zoi
1553C-----------------------------------------------
1554C L o c a l V a r i a b l e s
1555C-----------------------------------------------
1556 my_real
1557 . x51, x52, x53, x54,
1558 . y51, y52, y53, y54,
1559 . z51, z52, z53, z54,
1560 . xi1, xi2, xi3, xi4, xi5,
1561 . yi1, yi2, yi3, yi4, yi5,
1562 . zi1, zi2, zi3, zi4, zi5,
1563 .
1564 .
1565 .
1566 . xs, ys, zs,
1567 . xm1, xm2, xm3, xm4, xm5,
1568 . ym1, ym2, ym3, ym4, ym5,
1569 . zm1, zm2, zm3, zm4, zm5,la0,lb0,lc0
1570 my_real
1571 . sx1,sx2,sx3,sx4,sx5,sx0,
1572 . sy1,sy2,sy3,sy4,sy5,sy0,
1573 . sz1,sz2,sz3,sz4,sz5,sz0,aaa,s2
1574C------------------------------------------------
1575 interseca = 0
1576
1577 IF(vo12 <= zero .AND. v12 >= zero)THEN
1578 IF(abs(vo12) < em20)THEN
1579 aaa = one
1580 ELSEIF(abs(v12) < em20)THEN
1581 aaa = zero
1582 ELSE
1583 aaa = v12 / (v12-vo12)
1584 ENDIF
1585
1586 xm1 = xx1 + aaa*(xo1-xx1)
1587 ym1 = yy1 + aaa*(yo1-yy1)
1588 zm1 = zz1 + aaa*(zo1-zz1)
1589
1590 xm2 = xx2 + aaa*(xo2-xx2)
1591 ym2 = yy2 + aaa*(yo2-yy2)
1592 zm2 = zz2 + aaa*(zo2-zz2)
1593
1594 xs = xi + aaa*(xoi-xi)
1595 ys = yi + aaa*(yoi-yi)
1596 zs = zi + aaa*(zoi-zi)
1597
1598 xm5 = xx5 + aaa*(xo5-xx5)
1599 ym5 = yy5 + aaa*(yo5-yy5)
1600 zm5 = zz5 + aaa*(zo5-zz5)
1601
1602 x51 = xm1 - xm5
1603 y51 = ym1 - ym5
1604 z51 = zm1 - zm5
1605
1606 x52 = xm2 - xm5
1607 y52 = ym2 - ym5
1608 z52 = zm2 - zm5
1609
1610 xi1 = xm1 - xs
1611 yi1 = ym1 - ys
1612 zi1 = zm1 - zs
1613
1614 xi2 = xm2 - xs
1615 yi2 = ym2 - ys
1616 zi2 = zm2 - zs
1617
1618 xi5 = xm5 - xs
1619 yi5 = ym5 - ys
1620 zi5 = zm5 - zs
1621
1622 sx0 = y51*z52 - z51*y52
1623 sy0 = z51*x52 - x51*z52
1624 sz0 = x51*y52 - y51*x52
1625
1626 sx1 = yi5*zi1 - zi5*yi1
1627 sy1 = zi5*xi1 - xi5*zi1
1628 sz1 = xi5*yi1 - yi5*xi1
1629
1630 sx5 = yi1*zi2 - zi1*yi2
1631 sy5 = zi1*xi2 - xi1*zi2
1632 sz5 = xi1*yi2 - yi1*xi2
1633
1634 sx2 = yi2*zi5 - zi2*yi5
1635 sy2 = zi2*xi5 - xi2*zi5
1636 sz2 = xi2*yi5 - yi2*xi5
1637
1638 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
1639 lb0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
1640 lc0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
1641 la0 = one-lb0-lc0
1642
1643 IF(ixx3 /= ixx4)THEN
1644 IF(lb0 >= zerom .and. lb0 <= unp .and.
1645 . lc0 >= zerom .and. lc0 <= unp .and.
1646 . la0 >= zerom .and. la0 <= unp)THEN
1647 interseca = 1
1648 ENDIF
1649C----------loose tolerance for tria
1650 ELSE
1651 IF(lb0 >= zeromt .AND. lb0 <= unpt .AND.
1652 . lc0 >= zeromt .AND. lc0 <= unpt .AND.
1653 . la0 >= zeromt .AND. la0 <= unpt)THEN
1654 interseca = 1
1655 ENDIF
1656 END if! (IXX(I,3) /= IXX(I,4))THEN
1657 ENDIF
1658 IF(ixx3 /= ixx4)THEN
1659 IF(vo23 <= zero .AND. v23 >= zero)THEN
1660 IF(abs(vo23) < em20)THEN
1661 aaa = one
1662 ELSEIF(abs(v23) < em20)THEN
1663 aaa = zero
1664 ELSE
1665 aaa = v23 / (v23-vo23)
1666 ENDIF
1667
1668 xm2 = xx2 + aaa*(xo2-xx2)
1669 ym2 = yy2 + aaa*(yo2-yy2)
1670 zm2 = zz2 + aaa*(zo2-zz2)
1671
1672 xm3 = xx3 + aaa*(xo3-xx3)
1673 ym3 = yy3 + aaa*(yo3-yy3)
1674 zm3 = zz3 + aaa*(zo3-zz3)
1675
1676 xs = xi + aaa*(xoi-xi)
1677 ys = yi + aaa*(yoi-yi)
1678 zs = zi + aaa*(zoi-zi)
1679
1680 xm5 = xx5 + aaa*(xo5-xx5)
1681 ym5 = yy5 + aaa*(yo5-yy5)
1682 zm5 = zz5 + aaa*(zo5-zz5)
1683
1684 x52 = xm2 - xm5
1685 y52 = ym2 - ym5
1686 z52 = zm2 - zm5
1687
1688 x53 = xm3 - xm5
1689 y53 = ym3 - ym5
1690 z53 = zm3 - zm5
1691
1692 xi2 = xm2 - xs
1693 yi2 = ym2 - ys
1694 zi2 = zm2 - zs
1695
1696 xi3 = xm3 - xs
1697 yi3 = ym3 - ys
1698 zi3 = zm3 - zs
1699
1700 xi5 = xm5 - xs
1701 yi5 = ym5 - ys
1702 zi5 = zm5 - zs
1703
1704 sx0 = y52*z53 - z52*y53
1705 sy0 = z52*x53 - x52*z53
1706 sz0 = x52*y53 - y52*x53
1707
1708 sx2 = yi5*zi2 - zi5*yi2
1709 sy2 = zi5*xi2 - xi5*zi2
1710 sz2 = xi5*yi2 - yi5*xi2
1711
1712 sx5 = yi2*zi3 - zi2*yi3
1713 sy5 = zi2*xi3 - xi2*zi3
1714 sz5 = xi2*yi3 - yi2*xi3
1715
1716 sx3 = yi3*zi5 - zi3*yi5
1717 sy3 = zi3*xi5 - xi3*zi5
1718 sz3 = xi3*yi5 - yi3*xi5
1719
1720 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
1721 lb0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
1722 lc0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
1723 la0 = one-lb0-lc0
1724
1725 IF(lb0 >= zerom .and. lb0 <= unp .and.
1726 . lc0 >= zerom .and. lc0 <= unp .and.
1727 . la0 >= zerom .and. la0 <= unp)THEN
1728 interseca = 1
1729 ENDIF
1730 ENDIF
1731 IF(vo34 <= zero .AND. v34 >= zero)THEN
1732 IF(abs(vo34) < em20)THEN
1733 aaa = one
1734 ELSEIF(abs(v34) < em20)THEN
1735 aaa = zero
1736 ELSE
1737 aaa = v34 / (v34-vo34)
1738 ENDIF
1739
1740 xm3 = xx3 + aaa*(xo3-xx3)
1741 ym3 = yy3 + aaa*(yo3-yy3)
1742 zm3 = zz3 + aaa*(zo3-zz3)
1743
1744 xm4 = xx4 + aaa*(xo4-xx4)
1745 ym4 = yy4 + aaa*(yo4-yy4)
1746 zm4 = zz4 + aaa*(zo4-zz4)
1747
1748 xs = xi + aaa*(xoi-xi)
1749 ys = yi + aaa*(yoi-yi)
1750 zs = zi + aaa*(zoi-zi)
1751
1752 xm5 = xx5 + aaa*(xo5-xx5)
1753 ym5 = yy5 + aaa*(yo5-yy5)
1754 zm5 = zz5 + aaa*(zo5-zz5)
1755
1756 x53 = xm3 - xm5
1757 y53 = ym3 - ym5
1758 z53 = zm3 - zm5
1759
1760 x54 = xm4 - xm5
1761 y54 = ym4 - ym5
1762 z54 = zm4 - zm5
1763
1764 xi3 = xm3 - xs
1765 yi3 = ym3 - ys
1766 zi3 = zm3 - zs
1767
1768 xi4 = xm4 - xs
1769 yi4 = ym4 - ys
1770 zi4 = zm4 - zs
1771
1772 xi5 = xm5 - xs
1773 yi5 = ym5 - ys
1774 zi5 = zm5 - zs
1775
1776 sx0 = y53*z54 - z53*y54
1777 sy0 = z53*x54 - x53*z54
1778 sz0 = x53*y54 - y53*x54
1779
1780 sx3 = yi5*zi3 - zi5*yi3
1781 sy3 = zi5*xi3 - xi5*zi3
1782 sz3 = xi5*yi3 - yi5*xi3
1783
1784 sx5 = yi3*zi4 - zi3*yi4
1785 sy5 = zi3*xi4 - xi3*zi4
1786 sz5 = xi3*yi4 - yi3*xi4
1787
1788 sx4 = yi4*zi5 - zi4*yi5
1789 sy4 = zi4*xi5 - xi4*zi5
1790 sz4 = xi4*yi5 - yi4*xi5
1791
1792 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
1793 lb0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
1794 lc0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
1795 la0 = one-lb0-lc0
1796
1797 IF(lb0 >= zerom .and. lb0 <= unp .and.
1798 . lc0 >= zerom .and. lc0 <= unp .and.
1799 . la0 >= zerom .and. la0 <= unp)THEN
1800 interseca = 1
1801 ENDIF
1802 ENDIF
1803
1804 IF(vo41 <= zero .AND. v41 >= zero)THEN
1805 IF(abs(vo41) < em20)THEN
1806 aaa = one
1807 ELSEIF(abs(v41) < em20)THEN
1808 aaa = zero
1809 ELSE
1810 aaa = v41 / (v41-vo41)
1811 ENDIF
1812
1813 xm4 = xx4 + aaa*(xo4-xx4)
1814 ym4 = yy4 + aaa*(yo4-yy4)
1815 zm4 = zz4 + aaa*(zo4-zz4)
1816
1817 xm1 = xx1 + aaa*(xo1-xx1)
1818 ym1 = yy1 + aaa*(yo1-yy1)
1819 zm1 = zz1 + aaa*(zo1-zz1)
1820
1821 xs = xi + aaa*(xoi-xi)
1822 ys = yi + aaa*(yoi-yi)
1823 zs = zi + aaa*(zoi-zi)
1824
1825 xm5 = xx5 + aaa*(xo5-xx5)
1826 ym5 = yy5 + aaa*(yo5-yy5)
1827 zm5 = zz5 + aaa*(zo5-zz5)
1828
1829 x54 = xm4 - xm5
1830 y54 = ym4 - ym5
1831 z54 = zm4 - zm5
1832
1833 x51 = xm1 - xm5
1834 y51 = ym1 - ym5
1835 z51 = zm1 - zm5
1836
1837 xi4 = xm4 - xs
1838 yi4 = ym4 - ys
1839 zi4 = zm4 - zs
1840
1841 xi1 = xm1 - xs
1842 yi1 = ym1 - ys
1843 zi1 = zm1 - zs
1844
1845 xi5 = xm5 - xs
1846 yi5 = ym5 - ys
1847 zi5 = zm5 - zs
1848
1849 sx0 = y54*z51 - z54*y51
1850 sy0 = z54*x51 - x54*z51
1851 sz0 = x54*y51 - y54*x51
1852
1853 sx4 = yi5*zi4 - zi5*yi4
1854 sy4 = zi5*xi4 - xi5*zi4
1855 sz4 = xi5*yi4 - yi5*xi4
1856
1857 sx5 = yi4*zi1 - zi4*yi1
1858 sy5 = zi4*xi1 - xi4*zi1
1859 sz5 = xi4*yi1 - yi4*xi1
1860
1861 sx1 = yi1*zi5 - zi1*yi5
1862 sy1 = zi1*xi5 - xi1*zi5
1863 sz1 = xi1*yi5 - yi1*xi5
1864
1865 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
1866 lb0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
1867 lc0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
1868 la0 = one-lb0-lc0
1869
1870 IF(lb0 >= zerom .and. lb0 <= unp .and.
1871 . lc0 >= zerom .and. lc0 <= unp .and.
1872 . la0 >= zerom .and. la0 <= unp)THEN
1873 interseca = 1
1874 ENDIF
1875 ENDIF
1876 ENDIF
1877C
1878 RETURN
1879 END
1880!||====================================================================
1881!|| intersecb_25 ../engine/source/interfaces/int25/i25dst3_22.F
1882!||--- called by ------------------------------------------------------
1883!|| i25dst3_22 ../engine/source/interfaces/int25/i25dst3_22.F
1884!||====================================================================
1885 SUBROUTINE intersecb_25(
1886 1 IXX3 ,IXX4 ,INTERSECA,INTERSECB,
1887 1 V12 ,V23 ,V34 ,V41 ,
1888 2 VO12 ,VO23 ,VO34 ,VO41 ,XX1 ,
1889 3 YY1 ,ZZ1 ,XX2 ,YY2 ,ZZ2 ,
1890 4 XX3 ,YY3 ,ZZ3 ,XX4 ,YY4 ,
1891 5 ZZ4 ,XX5 ,YY5 ,ZZ5 ,
1892 6 VXI ,VYI ,VZI ,XO1 ,XO2 ,
1893 7 XO3 ,XO4 ,XO5 ,YO1 ,YO2 ,
1894 8 YO3 ,YO4 ,YO5 ,ZO1 ,ZO2 ,
1895 9 ZO3 ,ZO4 ,ZO5 ,XI ,YI ,
1896 A ZI ,XOI ,YOI ,ZOI ,ZEROM ,
1897 B UNP ,ZEROMT,UNPT )
1898C-----------------------------------------------
1899C I m p l i c i t T y p e s
1900C-----------------------------------------------
1901#include "implicit_f.inc"
1902C-----------------------------------------------
1903C D u m m y A r g u m e n t s
1904C-----------------------------------------------
1905 INTEGER IXX3 ,IXX4 , INTERSECA, INTERSECB
1906C REAL
1907 my_real
1908 1 v12 ,v23 ,v34 ,v41 ,
1909 2 vo12 ,vo23 ,vo34 ,vo41 ,xx1 ,
1910 3 yy1 ,zz1 ,xx2 ,yy2 ,zz2 ,
1911 4 xx3 ,yy3 ,zz3 ,xx4 ,yy4 ,
1912 5 zz4 ,xx5 ,yy5 ,zz5 ,
1913 6 vxi ,vyi ,vzi ,xo1 ,xo2 ,
1914 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
1915 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
1916 9 zo3 ,zo4 ,zo5 ,zerom ,unp ,
1917 a zeromt ,unpt ,xi ,yi ,zi ,
1918 b xoi ,yoi ,zoi
1919C-----------------------------------------------
1920C L o c a l V a r i a b l e s
1921C-----------------------------------------------
1922 my_real
1923 . x51, x52, x53, x54,
1924 . y51, y52, y53, y54,
1925 . z51, z52, z53, z54,
1926 . xi1, xi2, xi3, xi4, xi5,
1927 . yi1, yi2, yi3, yi4, yi5,
1928 . zi1, zi2, zi3, zi4, zi5,
1929 .
1930 .
1931 .
1932 . xs, ys, zs,
1933 . xm1, xm2, xm3, xm4, xm5,
1934 . ym1, ym2, ym3, ym4, ym5,
1935 . zm1, zm2, zm3, zm4, zm5,la0,lb0,lc0
1936 my_real
1937 . sx1,sx2,sx3,sx4,sx5,sx0,
1938 . sy1,sy2,sy3,sy4,sy5,sy0,
1939 . sz1,sz2,sz3,sz4,sz5,sz0,aaa,s2
1940C------------------------------------------------
1941 interseca = 0
1942 intersecb = 0
1943
1944
1945 IF(vo12 * v12 <= zero)THEN
1946 IF(abs(vo12) < em20)THEN
1947 aaa = one
1948 ELSEIF(abs(v12) < em20)THEN
1949 aaa = zero
1950 ELSE
1951 aaa = v12 / (v12-vo12)
1952 ENDIF
1953
1954 xm1 = xx1 + aaa*(xo1-xx1)
1955 ym1 = yy1 + aaa*(yo1-yy1)
1956 zm1 = zz1 + aaa*(zo1-zz1)
1957
1958 xm2 = xx2 + aaa*(xo2-xx2)
1959 ym2 = yy2 + aaa*(yo2-yy2)
1960 zm2 = zz2 + aaa*(zo2-zz2)
1961
1962 xs = xi + aaa*(xoi-xi)
1963 ys = yi + aaa*(yoi-yi)
1964 zs = zi + aaa*(zoi-zi)
1965
1966 xm5 = xx5 + aaa*(xo5-xx5)
1967 ym5 = yy5 + aaa*(yo5-yy5)
1968 zm5 = zz5 + aaa*(zo5-zz5)
1969
1970 x51 = xm1 - xm5
1971 y51 = ym1 - ym5
1972 z51 = zm1 - zm5
1973
1974 x52 = xm2 - xm5
1975 y52 = ym2 - ym5
1976 z52 = zm2 - zm5
1977
1978 xi1 = xm1 - xs
1979 yi1 = ym1 - ys
1980 zi1 = zm1 - zs
1981
1982 xi2 = xm2 - xs
1983 yi2 = ym2 - ys
1984 zi2 = zm2 - zs
1985
1986 xi5 = xm5 - xs
1987 yi5 = ym5 - ys
1988 zi5 = zm5 - zs
1989
1990 sx0 = y51*z52 - z51*y52
1991 sy0 = z51*x52 - x51*z52
1992 sz0 = x51*y52 - y51*x52
1993
1994 sx1 = yi5*zi1 - zi5*yi1
1995 sy1 = zi5*xi1 - xi5*zi1
1996 sz1 = xi5*yi1 - yi5*xi1
1997
1998 sx5 = yi1*zi2 - zi1*yi2
1999 sy5 = zi1*xi2 - xi1*zi2
2000 sz5 = xi1*yi2 - yi1*xi2
2001
2002 sx2 = yi2*zi5 - zi2*yi5
2003 sy2 = zi2*xi5 - xi2*zi5
2004 sz2 = xi2*yi5 - yi2*xi5
2005
2006 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2007 lb0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
2008 lc0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
2009 la0 = one-lb0-lc0
2010
2011 IF(ixx3 /= ixx4)THEN
2012 IF(lb0 >= zerom .and. lb0 <= unp .and.
2013 . lc0 >= zerom .and. lc0 <= unp .and.
2014 . la0 >= zerom .and. la0 <= unp)THEN
2015 IF(vo12 <= zero .AND. v12 >= zero)THEN
2016 interseca = 1
2017 ELSE !((VO12 > ZERO .AND. V12 < ZERO) .OR.
2018 . ! (VO12 > ZERO .AND. V12 == ZERO) .OR.
2019 . !(VO12 == ZERO .AND. V12 < ZERO)) THEN
2020 intersecb = 1
2021 END IF
2022 ENDIF
2023C----------loose tolerance for tria
2024 ELSE
2025 IF(lb0 >= zeromt .AND. lb0 <= unpt .AND.
2026 . lc0 >= zeromt .AND. lc0 <= unpt .AND.
2027 . la0 >= zeromt .AND. la0 <= unpt)THEN
2028 IF(vo12 <= zero .AND. v12 >= zero)THEN
2029 interseca = 1
2030 ELSE !((VO12 > ZERO .AND. V12 < ZERO) .OR.
2031 . !(VO12 > ZERO .AND. V12 == ZERO) .OR.
2032 . !(VO12 == ZERO .AND. V12 < ZERO)) THEN
2033 intersecb = 1
2034 ENDIF
2035 ENDIF
2036 END if! (IXX(I,3) /= IXX(I,4))THEN
2037 ENDIF
2038 IF(ixx3 /= ixx4)THEN
2039 IF(vo23 * v23 <= zero)THEN
2040 IF(abs(vo23) < em20)THEN
2041 aaa = one
2042 ELSEIF(abs(v23) < em20)THEN
2043 aaa = zero
2044 ELSE
2045 aaa = v23 / (v23-vo23)
2046 ENDIF
2047
2048 xm2 = xx2 + aaa*(xo2-xx2)
2049 ym2 = yy2 + aaa*(yo2-yy2)
2050 zm2 = zz2 + aaa*(zo2-zz2)
2051
2052 xm3 = xx3 + aaa*(xo3-xx3)
2053 ym3 = yy3 + aaa*(yo3-yy3)
2054 zm3 = zz3 + aaa*(zo3-zz3)
2055
2056 xs = xi + aaa*(xoi-xi)
2057 ys = yi + aaa*(yoi-yi)
2058 zs = zi + aaa*(zoi-zi)
2059
2060 xm5 = xx5 + aaa*(xo5-xx5)
2061 ym5 = yy5 + aaa*(yo5-yy5)
2062 zm5 = zz5 + aaa*(zo5-zz5)
2063
2064 x52 = xm2 - xm5
2065 y52 = ym2 - ym5
2066 z52 = zm2 - zm5
2067
2068 x53 = xm3 - xm5
2069 y53 = ym3 - ym5
2070 z53 = zm3 - zm5
2071
2072 xi2 = xm2 - xs
2073 yi2 = ym2 - ys
2074 zi2 = zm2 - zs
2075
2076 xi3 = xm3 - xs
2077 yi3 = ym3 - ys
2078 zi3 = zm3 - zs
2079
2080 xi5 = xm5 - xs
2081 yi5 = ym5 - ys
2082 zi5 = zm5 - zs
2083
2084 sx0 = y52*z53 - z52*y53
2085 sy0 = z52*x53 - x52*z53
2086 sz0 = x52*y53 - y52*x53
2087
2088 sx2 = yi5*zi2 - zi5*yi2
2089 sy2 = zi5*xi2 - xi5*zi2
2090 sz2 = xi5*yi2 - yi5*xi2
2091
2092 sx5 = yi2*zi3 - zi2*yi3
2093 sy5 = zi2*xi3 - xi2*zi3
2094 sz5 = xi2*yi3 - yi2*xi3
2095
2096 sx3 = yi3*zi5 - zi3*yi5
2097 sy3 = zi3*xi5 - xi3*zi5
2098 sz3 = xi3*yi5 - yi3*xi5
2099
2100 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2101 lb0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
2102 lc0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
2103 la0 = one-lb0-lc0
2104
2105 IF(lb0 >= zerom .and. lb0 <= unp .and.
2106 . lc0 >= zerom .and. lc0 <= unp .and.
2107 . la0 >= zerom .and. la0 <= unp)THEN
2108 IF(vo23 <= zero .AND. v23 >= zero)THEN
2109 interseca = 1
2110 ELSE !((VO23 > ZERO .AND. V23 < ZERO) .OR.
2111 . !(VO23 > ZERO .AND. V23 == ZERO) .OR.
2112 . !(VO23 == ZERO .AND. V23 < ZERO)) THEN
2113 intersecb = 1
2114 END IF
2115 ENDIF
2116 ENDIF
2117 IF(vo34 * v34 <= zero)THEN
2118 IF(abs(vo34) < em20)THEN
2119 aaa = one
2120 ELSEIF(abs(v34) < em20)THEN
2121 aaa = zero
2122 ELSE
2123 aaa = v34 / (v34-vo34)
2124 ENDIF
2125
2126 xm3 = xx3 + aaa*(xo3-xx3)
2127 ym3 = yy3 + aaa*(yo3-yy3)
2128 zm3 = zz3 + aaa*(zo3-zz3)
2129
2130 xm4 = xx4 + aaa*(xo4-xx4)
2131 ym4 = yy4 + aaa*(yo4-yy4)
2132 zm4 = zz4 + aaa*(zo4-zz4)
2133
2134 xs = xi + aaa*(xoi-xi)
2135 ys = yi + aaa*(yoi-yi)
2136 zs = zi + aaa*(zoi-zi)
2137
2138 xm5 = xx5 + aaa*(xo5-xx5)
2139 ym5 = yy5 + aaa*(yo5-yy5)
2140 zm5 = zz5 + aaa*(zo5-zz5)
2141
2142 x53 = xm3 - xm5
2143 y53 = ym3 - ym5
2144 z53 = zm3 - zm5
2145
2146 x54 = xm4 - xm5
2147 y54 = ym4 - ym5
2148 z54 = zm4 - zm5
2149
2150 xi3 = xm3 - xs
2151 yi3 = ym3 - ys
2152 zi3 = zm3 - zs
2153
2154 xi4 = xm4 - xs
2155 yi4 = ym4 - ys
2156 zi4 = zm4 - zs
2157
2158 xi5 = xm5 - xs
2159 yi5 = ym5 - ys
2160 zi5 = zm5 - zs
2161
2162 sx0 = y53*z54 - z53*y54
2163 sy0 = z53*x54 - x53*z54
2164 sz0 = x53*y54 - y53*x54
2165
2166 sx3 = yi5*zi3 - zi5*yi3
2167 sy3 = zi5*xi3 - xi5*zi3
2168 sz3 = xi5*yi3 - yi5*xi3
2169
2170 sx5 = yi3*zi4 - zi3*yi4
2171 sy5 = zi3*xi4 - xi3*zi4
2172 sz5 = xi3*yi4 - yi3*xi4
2173
2174 sx4 = yi4*zi5 - zi4*yi5
2175 sy4 = zi4*xi5 - xi4*zi5
2176 sz4 = xi4*yi5 - yi4*xi5
2177
2178 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2179 lb0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
2180 lc0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
2181 la0 = one-lb0-lc0
2182
2183 IF(lb0 >= zerom .and. lb0 <= unp .and.
2184 . lc0 >= zerom .and. lc0 <= unp .and.
2185 . la0 >= zerom .and. la0 <= unp)THEN
2186 IF(vo34 <= zero .AND. v34 >= zero)THEN
2187 interseca = 1
2188 ELSE !((VO34 > ZERO .AND. V34 < ZERO) .OR.
2189 . !(VO34 > ZERO .AND. V34 == ZERO) .OR.
2190 . !(VO34 == ZERO .AND. V34 < ZERO)) THEN
2191 intersecb = 1
2192 END IF
2193 ENDIF
2194 ENDIF
2195
2196 IF(vo41 * v41 <= zero)THEN
2197 IF(abs(vo41) < em20)THEN
2198 aaa = one
2199 ELSEIF(abs(v41) < em20)THEN
2200 aaa = zero
2201 ELSE
2202 aaa = v41 / (v41-vo41)
2203 ENDIF
2204
2205 xm4 = xx4 + aaa*(xo4-xx4)
2206 ym4 = yy4 + aaa*(yo4-yy4)
2207 zm4 = zz4 + aaa*(zo4-zz4)
2208
2209 xm1 = xx1 + aaa*(xo1-xx1)
2210 ym1 = yy1 + aaa*(yo1-yy1)
2211 zm1 = zz1 + aaa*(zo1-zz1)
2212
2213 xs = xi + aaa*(xoi-xi)
2214 ys = yi + aaa*(yoi-yi)
2215 zs = zi + aaa*(zoi-zi)
2216
2217 xm5 = xx5 + aaa*(xo5-xx5)
2218 ym5 = yy5 + aaa*(yo5-yy5)
2219 zm5 = zz5 + aaa*(zo5-zz5)
2220
2221 x54 = xm4 - xm5
2222 y54 = ym4 - ym5
2223 z54 = zm4 - zm5
2224
2225 x51 = xm1 - xm5
2226 y51 = ym1 - ym5
2227 z51 = zm1 - zm5
2228
2229 xi4 = xm4 - xs
2230 yi4 = ym4 - ys
2231 zi4 = zm4 - zs
2232
2233 xi1 = xm1 - xs
2234 yi1 = ym1 - ys
2235 zi1 = zm1 - zs
2236
2237 xi5 = xm5 - xs
2238 yi5 = ym5 - ys
2239 zi5 = zm5 - zs
2240
2241 sx0 = y54*z51 - z54*y51
2242 sy0 = z54*x51 - x54*z51
2243 sz0 = x54*y51 - y54*x51
2244
2245 sx4 = yi5*zi4 - zi5*yi4
2246 sy4 = zi5*xi4 - xi5*zi4
2247 sz4 = xi5*yi4 - yi5*xi4
2248
2249 sx5 = yi4*zi1 - zi4*yi1
2250 sy5 = zi4*xi1 - xi4*zi1
2251 sz5 = xi4*yi1 - yi4*xi1
2252
2253 sx1 = yi1*zi5 - zi1*yi5
2254 sy1 = zi1*xi5 - xi1*zi5
2255 sz1 = xi1*yi5 - yi1*xi5
2256
2257 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2258 lb0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
2259 lc0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
2260 la0 = one-lb0-lc0
2261
2262 IF(lb0 >= zerom .and. lb0 <= unp .and.
2263 . lc0 >= zerom .and. lc0 <= unp .and.
2264 . la0 >= zerom .and. la0 <= unp)THEN
2265 IF(vo41 <= zero .AND. v41 >= zero)THEN
2266 interseca = 1
2267 ELSE !((VO41 > ZERO .AND. V41 < ZERO) .OR.
2268 . !(VO41 > ZERO .AND. V41 == ZERO) .OR.
2269 . !(VO41 == ZERO .AND. V41 < ZERO)) THEN
2270 intersecb = 1
2271 END IF
2272 ENDIF
2273 ENDIF
2274 ENDIF
2275C
2276 RETURN
2277 END
2278!||====================================================================
2279!|| intersecv_25 ../engine/source/interfaces/int25/i25dst3_22.F
2280!||====================================================================
2281 SUBROUTINE intersecv_25(
2282 1 IXX3 ,IXX4 ,SUBTRIA,
2283 1 INTERSECT,V12 ,V23 ,V34 ,V41 ,
2284 2 VO12 ,VO23 ,VO34 ,VO41 ,XX1 ,
2285 3 YY1 ,ZZ1 ,XX2 ,YY2 ,ZZ2 ,
2286 4 XX3 ,YY3 ,ZZ3 ,XX4 ,YY4 ,
2287 5 ZZ4 ,XX5 ,YY5 ,ZZ5 ,
2288 6 VXI ,VYI ,VZI ,XO1 ,XO2 ,
2289 7 XO3 ,XO4 ,XO5 ,YO1 ,YO2 ,
2290 8 YO3 ,YO4 ,YO5 ,ZO1 ,ZO2 ,
2291 9 ZO3 ,ZO4 ,ZO5 ,XI ,YI ,
2292 A ZI ,XOI ,YOI ,ZOI ,ZEROM ,
2293 B UNP ,ZEROMT,UNPT )
2294C-----------------------------------------------
2295C I m p l i c i t T y p e s
2296C-----------------------------------------------
2297#include "implicit_f.inc"
2298C-----------------------------------------------
2299C D u m m y A r g u m e n t s
2300C-----------------------------------------------
2301 INTEGER IXX3 ,IXX4 ,INTERSECT, SUBTRIA
2302C REAL
2303 my_real
2304 1 v12 ,v23 ,v34 ,v41 ,
2305 2 vo12 ,vo23 ,vo34 ,vo41 ,xx1 ,
2306 3 yy1 ,zz1 ,xx2 ,yy2 ,zz2 ,
2307 4 xx3 ,yy3 ,zz3 ,xx4 ,yy4 ,
2308 5 zz4 ,xx5 ,yy5 ,zz5 ,
2309 6 vxi ,vyi ,vzi ,xo1 ,xo2 ,
2310 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
2311 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
2312 9 zo3 ,zo4 ,zo5 ,zerom ,unp ,
2313 a zeromt ,unpt ,xi ,yi ,zi ,
2314 b xoi ,yoi ,zoi
2315C-----------------------------------------------
2316C L o c a l V a r i a b l e s
2317C-----------------------------------------------
2318 my_real
2319 . x51, x52, x53, x54,
2320 . y51, y52, y53, y54,
2321 . z51, z52, z53, z54,
2322 . xi1, xi2, xi3, xi4, xi5,
2323 . yi1, yi2, yi3, yi4, yi5,
2324 . zi1, zi2, zi3, zi4, zi5,
2325 .
2326 .
2327 .
2328 . xs, ys, zs,
2329 . xm1, xm2, xm3, xm4, xm5,
2330 . ym1, ym2, ym3, ym4, ym5,
2331 . zm1, zm2, zm3, zm4, zm5,la0,lb0,lc0
2332 my_real
2333 . sx1,sx2,sx3,sx4,sx5,sx0,
2334 . sy1,sy2,sy3,sy4,sy5,sy0,
2335 . sz1,sz2,sz3,sz4,sz5,sz0,aaa,s2
2336C------------------------------------------------
2337 IF(vo12 <= zero .AND. v12 >= zero)THEN
2338 IF(abs(vo12) < em20)THEN
2339 aaa = one
2340 ELSEIF(abs(v12) < em20)THEN
2341 aaa = zero
2342 ELSE
2343 aaa = v12 / (v12-vo12)
2344 ENDIF
2345
2346 xm1 = xx1 + aaa*(xo1-xx1)
2347 ym1 = yy1 + aaa*(yo1-yy1)
2348 zm1 = zz1 + aaa*(zo1-zz1)
2349
2350 xm2 = xx2 + aaa*(xo2-xx2)
2351 ym2 = yy2 + aaa*(yo2-yy2)
2352 zm2 = zz2 + aaa*(zo2-zz2)
2353
2354 xs = xi + aaa*(xoi-xi)
2355 ys = yi + aaa*(yoi-yi)
2356 zs = zi + aaa*(zoi-zi)
2357
2358 xm5 = xx5 + aaa*(xo5-xx5)
2359 ym5 = yy5 + aaa*(yo5-yy5)
2360 zm5 = zz5 + aaa*(zo5-zz5)
2361
2362 x51 = xm1 - xm5
2363 y51 = ym1 - ym5
2364 z51 = zm1 - zm5
2365
2366 x52 = xm2 - xm5
2367 y52 = ym2 - ym5
2368 z52 = zm2 - zm5
2369
2370 xi1 = xm1 - xs
2371 yi1 = ym1 - ys
2372 zi1 = zm1 - zs
2373
2374 xi2 = xm2 - xs
2375 yi2 = ym2 - ys
2376 zi2 = zm2 - zs
2377
2378 xi5 = xm5 - xs
2379 yi5 = ym5 - ys
2380 zi5 = zm5 - zs
2381
2382 sx0 = y51*z52 - z51*y52
2383 sy0 = z51*x52 - x51*z52
2384 sz0 = x51*y52 - y51*x52
2385
2386 sx1 = yi5*zi1 - zi5*yi1
2387 sy1 = zi5*xi1 - xi5*zi1
2388 sz1 = xi5*yi1 - yi5*xi1
2389
2390 sx5 = yi1*zi2 - zi1*yi2
2391 sy5 = zi1*xi2 - xi1*zi2
2392 sz5 = xi1*yi2 - yi1*xi2
2393
2394 sx2 = yi2*zi5 - zi2*yi5
2395 sy2 = zi2*xi5 - xi2*zi5
2396 sz2 = xi2*yi5 - yi2*xi5
2397
2398 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2399 lb0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
2400 lc0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
2401 la0 = one-lb0-lc0
2402
2403 IF(ixx3 /= ixx4)THEN
2404 IF(lb0 >= zerom .and. lb0 <= unp .and.
2405 . lc0 >= zerom .and. lc0 <= unp .and.
2406 . la0 >= zerom .and. la0 <= unp)THEN
2407 IF(vo12 <= zero .AND. v12 >= zero)THEN
2408 intersect = 1
2409 ELSE ! (VO12 > ZERO .AND. V12 < ZERO) .OR.
2410 ! (VO12 > ZERO .AND. V12 == ZERO) .OR.
2411 ! (VO12 == ZERO .AND. V12 < ZERO)
2412 intersect = -1
2413 END IF
2414 subtria= 1 ! subtriangle
2415 ENDIF
2416C----------loose tolerance for tria
2417 ELSE
2418 IF(lb0 >= zeromt .AND. lb0 <= unpt .AND.
2419 . lc0 >= zeromt .AND. lc0 <= unpt .AND.
2420 . la0 >= zeromt .AND. la0 <= unpt)THEN
2421 intersect = 1
2422 subtria= 1 ! subtriangle
2423 ENDIF
2424 END if! (IXX(I,3) /= IXX(I,4))THEN
2425 ENDIF
2426 IF(ixx3 /= ixx4)THEN
2427 IF(vo23 <= zero .AND. v23 >= zero)THEN
2428 IF(abs(vo23) < em20)THEN
2429 aaa = one
2430 ELSEIF(abs(v23) < em20)THEN
2431 aaa = zero
2432 ELSE
2433 aaa = v23 / (v23-vo23)
2434 ENDIF
2435
2436 xm2 = xx2 + aaa*(xo2-xx2)
2437 ym2 = yy2 + aaa*(yo2-yy2)
2438 zm2 = zz2 + aaa*(zo2-zz2)
2439
2440 xm3 = xx3 + aaa*(xo3-xx3)
2441 ym3 = yy3 + aaa*(yo3-yy3)
2442 zm3 = zz3 + aaa*(zo3-zz3)
2443
2444 xs = xi + aaa*(xoi-xi)
2445 ys = yi + aaa*(yoi-yi)
2446 zs = zi + aaa*(zoi-zi)
2447
2448 xm5 = xx5 + aaa*(xo5-xx5)
2449 ym5 = yy5 + aaa*(yo5-yy5)
2450 zm5 = zz5 + aaa*(zo5-zz5)
2451
2452 x52 = xm2 - xm5
2453 y52 = ym2 - ym5
2454 z52 = zm2 - zm5
2455
2456 x53 = xm3 - xm5
2457 y53 = ym3 - ym5
2458 z53 = zm3 - zm5
2459
2460 xi2 = xm2 - xs
2461 yi2 = ym2 - ys
2462 zi2 = zm2 - zs
2463
2464 xi3 = xm3 - xs
2465 yi3 = ym3 - ys
2466 zi3 = zm3 - zs
2467
2468 xi5 = xm5 - xs
2469 yi5 = ym5 - ys
2470 zi5 = zm5 - zs
2471
2472 sx0 = y52*z53 - z52*y53
2473 sy0 = z52*x53 - x52*z53
2474 sz0 = x52*y53 - y52*x53
2475
2476 sx2 = yi5*zi2 - zi5*yi2
2477 sy2 = zi5*xi2 - xi5*zi2
2478 sz2 = xi5*yi2 - yi5*xi2
2479
2480 sx5 = yi2*zi3 - zi2*yi3
2481 sy5 = zi2*xi3 - xi2*zi3
2482 sz5 = xi2*yi3 - yi2*xi3
2483
2484 sx3 = yi3*zi5 - zi3*yi5
2485 sy3 = zi3*xi5 - xi3*zi5
2486 sz3 = xi3*yi5 - yi3*xi5
2487
2488 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2489 lb0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
2490 lc0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
2491 la0 = one-lb0-lc0
2492
2493 IF(lb0 >= zerom .and. lb0 <= unp .and.
2494 . lc0 >= zerom .and. lc0 <= unp .and.
2495 . la0 >= zerom .and. la0 <= unp)THEN
2496 IF(vo23 <= zero .AND. v23 >= zero)THEN
2497 intersect = 1
2498 ELSE ! (VO23 > ZERO .AND. V23 < ZERO) .OR.
2499 ! (VO23 > ZERO .AND. V23 == ZERO) .OR.
2500 ! (VO23 == ZERO .AND. V23 < ZERO)
2501 intersect = -1
2502 END IF
2503
2504 ENDIF
2505 ENDIF
2506 IF(vo34 <= zero .AND. v34 >= zero)THEN
2507 IF(abs(vo34) < em20)THEN
2508 aaa = one
2509 ELSEIF(abs(v34) < em20)THEN
2510 aaa = zero
2511 ELSE
2512 aaa = v34 / (v34-vo34)
2513 ENDIF
2514
2515 xm3 = xx3 + aaa*(xo3-xx3)
2516 ym3 = yy3 + aaa*(yo3-yy3)
2517 zm3 = zz3 + aaa*(zo3-zz3)
2518
2519 xm4 = xx4 + aaa*(xo4-xx4)
2520 ym4 = yy4 + aaa*(yo4-yy4)
2521 zm4 = zz4 + aaa*(zo4-zz4)
2522
2523 xs = xi + aaa*(xoi-xi)
2524 ys = yi + aaa*(yoi-yi)
2525 zs = zi + aaa*(zoi-zi)
2526
2527 xm5 = xx5 + aaa*(xo5-xx5)
2528 ym5 = yy5 + aaa*(yo5-yy5)
2529 zm5 = zz5 + aaa*(zo5-zz5)
2530
2531 x53 = xm3 - xm5
2532 y53 = ym3 - ym5
2533 z53 = zm3 - zm5
2534
2535 x54 = xm4 - xm5
2536 y54 = ym4 - ym5
2537 z54 = zm4 - zm5
2538
2539 xi3 = xm3 - xs
2540 yi3 = ym3 - ys
2541 zi3 = zm3 - zs
2542
2543 xi4 = xm4 - xs
2544 yi4 = ym4 - ys
2545 zi4 = zm4 - zs
2546
2547 xi5 = xm5 - xs
2548 yi5 = ym5 - ys
2549 zi5 = zm5 - zs
2550
2551 sx0 = y53*z54 - z53*y54
2552 sy0 = z53*x54 - x53*z54
2553 sz0 = x53*y54 - y53*x54
2554
2555 sx3 = yi5*zi3 - zi5*yi3
2556 sy3 = zi5*xi3 - xi5*zi3
2557 sz3 = xi5*yi3 - yi5*xi3
2558
2559 sx5 = yi3*zi4 - zi3*yi4
2560 sy5 = zi3*xi4 - xi3*zi4
2561 sz5 = xi3*yi4 - yi3*xi4
2562
2563 sx4 = yi4*zi5 - zi4*yi5
2564 sy4 = zi4*xi5 - xi4*zi5
2565 sz4 = xi4*yi5 - yi4*xi5
2566
2567 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2568 lb0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
2569 lc0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
2570 la0 = one-lb0-lc0
2571
2572 IF(lb0 >= zerom .and. lb0 <= unp .and.
2573 . lc0 >= zerom .and. lc0 <= unp .and.
2574 . la0 >= zerom .and. la0 <= unp)THEN
2575 IF(vo34 <= zero .AND. v34 >= zero)THEN
2576 intersect = 1
2577 ELSE ! (VO34 > ZERO .AND. V34 < ZERO) .OR.
2578 ! (VO34 > ZERO .AND. V34 == ZERO) .OR.
2579 ! (VO34 == ZERO .AND. V34 < ZERO)
2580 intersect = -1
2581 END IF
2582 ENDIF
2583 ENDIF
2584
2585 IF(vo41 <= zero .AND. v41 >= zero)THEN
2586 IF(abs(vo41) < em20)THEN
2587 aaa = one
2588 ELSEIF(abs(v41) < em20)THEN
2589 aaa = zero
2590 ELSE
2591 aaa = v41 / (v41-vo41)
2592 ENDIF
2593
2594 xm4 = xx4 + aaa*(xo4-xx4)
2595 ym4 = yy4 + aaa*(yo4-yy4)
2596 zm4 = zz4 + aaa*(zo4-zz4)
2597
2598 xm1 = xx1 + aaa*(xo1-xx1)
2599 ym1 = yy1 + aaa*(yo1-yy1)
2600 zm1 = zz1 + aaa*(zo1-zz1)
2601
2602 xs = xi + aaa*(xoi-xi)
2603 ys = yi + aaa*(yoi-yi)
2604 zs = zi + aaa*(zoi-zi)
2605
2606 xm5 = xx5 + aaa*(xo5-xx5)
2607 ym5 = yy5 + aaa*(yo5-yy5)
2608 zm5 = zz5 + aaa*(zo5-zz5)
2609
2610 x54 = xm4 - xm5
2611 y54 = ym4 - ym5
2612 z54 = zm4 - zm5
2613
2614 x51 = xm1 - xm5
2615 y51 = ym1 - ym5
2616 z51 = zm1 - zm5
2617
2618 xi4 = xm4 - xs
2619 yi4 = ym4 - ys
2620 zi4 = zm4 - zs
2621
2622 xi1 = xm1 - xs
2623 yi1 = ym1 - ys
2624 zi1 = zm1 - zs
2625
2626 xi5 = xm5 - xs
2627 yi5 = ym5 - ys
2628 zi5 = zm5 - zs
2629
2630 sx0 = y54*z51 - z54*y51
2631 sy0 = z54*x51 - x54*z51
2632 sz0 = x54*y51 - y54*x51
2633
2634 sx4 = yi5*zi4 - zi5*yi4
2635 sy4 = zi5*xi4 - xi5*zi4
2636 sz4 = xi5*yi4 - yi5*xi4
2637
2638 sx5 = yi4*zi1 - zi4*yi1
2639 sy5 = zi4*xi1 - xi4*zi1
2640 sz5 = xi4*yi1 - yi4*xi1
2641
2642 sx1 = yi1*zi5 - zi1*yi5
2643 sy1 = zi1*xi5 - xi1*zi5
2644 sz1 = xi1*yi5 - yi1*xi5
2645
2646 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2647 lb0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
2648 lc0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
2649 la0 = one-lb0-lc0
2650
2651 IF(lb0 >= zerom .and. lb0 <= unp .and.
2652 . lc0 >= zerom .and. lc0 <= unp .and.
2653 . la0 >= zerom .and. la0 <= unp)THEN
2654 IF(vo41 <= zero .AND. v41 >= zero)THEN
2655 intersect = 1
2656 ELSE ! (VO41 > ZERO .AND. V41 < ZERO) .OR.
2657 ! (VO41 > ZERO .AND. V41 == ZERO) .OR.
2658 ! (VO41 == ZERO .AND. V41 < ZERO)
2659 intersect = -1
2660 END IF
2661 ENDIF
2662 ENDIF
2663 ENDIF
2664C
2665 RETURN
2666 END
2667!||====================================================================
2668!|| intersecv0_25 ../engine/source/interfaces/int25/i25dst3_22.F
2669!||--- called by ------------------------------------------------------
2670!|| i25dst3_22 ../engine/source/interfaces/int25/i25dst3_22.F
2671!||====================================================================
2672 SUBROUTINE intersecv0_25(
2673 1 IXX3 ,IXX4 ,
2674 1 INTERSECT,V12 ,V23 ,V34 ,V41 ,
2675 2 VO12 ,VO23 ,VO34 ,VO41 ,XX1 ,
2676 3 YY1 ,ZZ1 ,XX2 ,YY2 ,ZZ2 ,
2677 4 XX3 ,YY3 ,ZZ3 ,XX4 ,YY4 ,
2678 5 ZZ4 ,XX5 ,YY5 ,ZZ5 ,
2679 6 XOI ,YOI ,ZOI ,XO1 ,XO2 ,
2680 7 XO3 ,XO4 ,XO5 ,YO1 ,YO2 ,
2681 8 YO3 ,YO4 ,YO5 ,ZO1 ,ZO2 ,
2682 9 ZO3 ,ZO4 ,ZO5 ,XI ,YI ,
2683 A ZI ,ZEROMT , UNPT )
2684C-----------------------------------------------
2685C I m p l i c i t T y p e s
2686C-----------------------------------------------
2687#include "implicit_f.inc"
2688C-----------------------------------------------
2689C D u m m y A r g u m e n t s
2690C-----------------------------------------------
2691 INTEGER IXX3 ,IXX4 ,INTERSECT
2692C REAL
2693 my_real
2694 1 v12 ,v23 ,v34 ,v41 ,
2695 2 vo12 ,vo23 ,vo34 ,vo41 ,xx1 ,
2696 3 yy1 ,zz1 ,xx2 ,yy2 ,zz2 ,
2697 4 xx3 ,yy3 ,zz3 ,xx4 ,yy4 ,
2698 5 zz4 ,xx5 ,yy5 ,zz5 ,
2699 6 xoi ,yoi ,zoi ,xo1 ,xo2 ,
2700 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
2701 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
2702 9 zo3 ,zo4 ,zo5 ,zeromt,unpt ,
2703 a xi ,yi ,zi
2704C-----------------------------------------------
2705C C o m m o n B l o c k s
2706C-----------------------------------------------
2707#include "com08_c.inc"
2708C-----------------------------------------------
2709C L o c a l V a r i a b l e s
2710C-----------------------------------------------
2711 my_real
2712 . x51, x52, x53, x54,
2713 . y51, y52, y53, y54,
2714 . z51, z52, z53, z54,
2715 . xi1, xi2, xi3, xi4, xi5,
2716 . yi1, yi2, yi3, yi4, yi5,
2717 . zi1, zi2, zi3, zi4, zi5,
2718 .
2719 .
2720 .
2721 .
2722 . xm1, xm2, xm3, xm4, xm5,
2723 . ym1, ym2, ym3, ym4, ym5,
2724 . zm1, zm2, zm3, zm4, zm5,la0,lb0,lc0
2725 my_real
2726 . sx1,sx2,sx3,sx4,sx5,sx0,
2727 . sy1,sy2,sy3,sy4,sy5,sy0,
2728 . sz1,sz2,sz3,sz4,sz5,sz0,aaa,adt,s2
2729C------------------------------------------------
2730C---------------- cas V>0 V_old>0
2731 IF(vo12 > zero .and. v12 >= zero )THEN
2732 aaa = one
2733 adt = dt1
2734
2735 xm1 = xo1
2736 ym1 = yo1
2737 zm1 = zo1
2738
2739 xm2 = xo2
2740 ym2 = yo2
2741 zm2 = zo2
2742
2743c XOI = XI - ADT*VXI
2744c YOI = YI - ADT*VYI
2745c ZOI = ZI - ADT*VZI
2746
2747 xm5 = xo5
2748 ym5 = yo5
2749 zm5 = zo5
2750
2751 x51 = xm1 - xm5
2752 y51 = ym1 - ym5
2753 z51 = zm1 - zm5
2754
2755 x52 = xm2 - xm5
2756 y52 = ym2 - ym5
2757 z52 = zm2 - zm5
2758
2759 xi1 = xm1 - xoi
2760 yi1 = ym1 - yoi
2761 zi1 = zm1 - zoi
2762
2763 xi2 = xm2 - xoi
2764 yi2 = ym2 - yoi
2765 zi2 = zm2 - zoi
2766
2767 xi5 = xm5 - xoi
2768 yi5 = ym5 - yoi
2769 zi5 = zm5 - zoi
2770
2771 sx0 = y51*z52 - z51*y52
2772 sy0 = z51*x52 - x51*z52
2773 sz0 = x51*y52 - y51*x52
2774
2775 sx1 = yi5*zi1 - zi5*yi1
2776 sy1 = zi5*xi1 - xi5*zi1
2777 sz1 = xi5*yi1 - yi5*xi1
2778
2779 sx5 = yi1*zi2 - zi1*yi2
2780 sy5 = zi1*xi2 - xi1*zi2
2781 sz5 = xi1*yi2 - yi1*xi2
2782
2783 sx2 = yi2*zi5 - zi2*yi5
2784 sy2 = zi2*xi5 - xi2*zi5
2785 sz2 = xi2*yi5 - yi2*xi5
2786
2787 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2788 lb0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
2789 lc0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
2790 la0 = one-lb0-lc0
2791
2792 IF(lb0 >= zeromt .AND. lb0 <= unpt .AND.
2793 . lc0 >= zeromt .AND. lc0 <= unpt .AND.
2794 . la0 >= zeromt .AND. la0 <= unpt)THEN
2795 intersect = 1
2796 ENDIF
2797 ENDIF
2798
2799 IF(ixx3 /= ixx4)THEN
2800 IF(vo23 > zero .AND. v23 >= zero)THEN
2801 aaa = one
2802 adt = dt1
2803
2804 xm2 = xo2
2805 ym2 = yo2
2806 zm2 = zo2
2807
2808 xm3 = xo3
2809 ym3 = yo3
2810 zm3 = zo3
2811
2812c XOI = XI - ADT*VXI
2813c YOI = YI - ADT*VYI
2814c ZOI = ZI - ADT*VZI
2815
2816 xm5 = xo5
2817 ym5 = yo5
2818 zm5 = zo5
2819
2820 x52 = xm2 - xm5
2821 y52 = ym2 - ym5
2822 z52 = zm2 - zm5
2823
2824 x53 = xm3 - xm5
2825 y53 = ym3 - ym5
2826 z53 = zm3 - zm5
2827
2828 xi2 = xm2 - xoi
2829 yi2 = ym2 - yoi
2830 zi2 = zm2 - zoi
2831
2832 xi3 = xm3 - xoi
2833 yi3 = ym3 - yoi
2834 zi3 = zm3 - zoi
2835
2836 xi5 = xm5 - xoi
2837 yi5 = ym5 - yoi
2838 zi5 = zm5 - zoi
2839
2840 sx0 = y52*z53 - z52*y53
2841 sy0 = z52*x53 - x52*z53
2842 sz0 = x52*y53 - y52*x53
2843
2844 sx2 = yi5*zi2 - zi5*yi2
2845 sy2 = zi5*xi2 - xi5*zi2
2846 sz2 = xi5*yi2 - yi5*xi2
2847
2848 sx5 = yi2*zi3 - zi2*yi3
2849 sy5 = zi2*xi3 - xi2*zi3
2850 sz5 = xi2*yi3 - yi2*xi3
2851
2852 sx3 = yi3*zi5 - zi3*yi5
2853 sy3 = zi3*xi5 - xi3*zi5
2854 sz3 = xi3*yi5 - yi3*xi5
2855
2856 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2857 lb0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
2858 lc0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
2859 la0 = one-lb0-lc0
2860
2861 IF(lb0 >= zeromt.and. lb0 <= unpt .and.
2862 . lc0 >= zeromt.and. lc0 <= unpt .and.
2863 . la0 >= zeromt.and. la0 <= unpt)THEN
2864 intersect = 1
2865 ENDIF
2866 ENDIF
2867 IF(vo34 > zero .and. v34 >= zero )THEN
2868 aaa = one
2869 adt = dt1
2870
2871 xm3 = xo3
2872 ym3 = yo3
2873 zm3 = zo3
2874
2875 xm4 = xo4
2876 ym4 = yo4
2877 zm4 = zo4
2878
2879c XOI = XI - ADT*VXI
2880c YOI = YI - ADT*VYI
2881c ZOI = ZI - ADT*VZI
2882
2883 xm5 = xo5
2884 ym5 = yo5
2885 zm5 = zo5
2886
2887 x53 = xm3 - xm5
2888 y53 = ym3 - ym5
2889 z53 = zm3 - zm5
2890
2891 x54 = xm4 - xm5
2892 y54 = ym4 - ym5
2893 z54 = zm4 - zm5
2894
2895 xi3 = xm3 - xoi
2896 yi3 = ym3 - yoi
2897 zi3 = zm3 - zoi
2898
2899 xi4 = xm4 - xoi
2900 yi4 = ym4 - yoi
2901 zi4 = zm4 - zoi
2902
2903 xi5 = xm5 - xoi
2904 yi5 = ym5 - yoi
2905 zi5 = zm5 - zoi
2906
2907 sx0 = y53*z54 - z53*y54
2908 sy0 = z53*x54 - x53*z54
2909 sz0 = x53*y54 - y53*x54
2910
2911 sx3 = yi5*zi3 - zi5*yi3
2912 sy3 = zi5*xi3 - xi5*zi3
2913 sz3 = xi5*yi3 - yi5*xi3
2914
2915 sx5 = yi3*zi4 - zi3*yi4
2916 sy5 = zi3*xi4 - xi3*zi4
2917 sz5 = xi3*yi4 - yi3*xi4
2918
2919 sx4 = yi4*zi5 - zi4*yi5
2920 sy4 = zi4*xi5 - xi4*zi5
2921 sz4 = xi4*yi5 - yi4*xi5
2922
2923 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2924 lb0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
2925 lc0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
2926 la0 = one-lb0-lc0
2927
2928 IF(lb0 >= zeromt.and. lb0 <= unpt .and.
2929 . lc0 >= zeromt.and. lc0 <= unpt .and.
2930 . la0 >= zeromt.and. la0 <= unpt)THEN
2931 intersect = 1
2932 ENDIF
2933 ENDIF
2934
2935 IF(vo41 > zero .and. v41 >= zero )THEN
2936 aaa = one
2937 adt = dt1
2938
2939 xm4 = xo4
2940 ym4 = yo4
2941 zm4 = zo4
2942
2943 xm1 = xo1
2944 ym1 = yo1
2945 zm1 = zo1
2946
2947c XOI = XI - ADT*VXI
2948c YOI = YI - ADT*VYI
2949c ZOI = ZI - ADT*VZI
2950
2951 xm5 = xo5
2952 ym5 = yo5
2953 zm5 = zo5
2954
2955 x54 = xm4 - xm5
2956 y54 = ym4 - ym5
2957 z54 = zm4 - zm5
2958
2959 x51 = xm1 - xm5
2960 y51 = ym1 - ym5
2961 z51 = zm1 - zm5
2962
2963 xi4 = xm4 - xoi
2964 yi4 = ym4 - yoi
2965 zi4 = zm4 - zoi
2966
2967 xi1 = xm1 - xoi
2968 yi1 = ym1 - yoi
2969 zi1 = zm1 - zoi
2970
2971 xi5 = xm5 - xoi
2972 yi5 = ym5 - yoi
2973 zi5 = zm5 - zoi
2974
2975 sx0 = y54*z51 - z54*y51
2976 sy0 = z54*x51 - x54*z51
2977 sz0 = x54*y51 - y54*x51
2978
2979 sx4 = yi5*zi4 - zi5*yi4
2980 sy4 = zi5*xi4 - xi5*zi4
2981 sz4 = xi5*yi4 - yi5*xi4
2982
2983 sx5 = yi4*zi1 - zi4*yi1
2984 sy5 = zi4*xi1 - xi4*zi1
2985 sz5 = xi4*yi1 - yi4*xi1
2986
2987 sx1 = yi1*zi5 - zi1*yi5
2988 sy1 = zi1*xi5 - xi1*zi5
2989 sz1 = xi1*yi5 - yi1*xi5
2990
2991 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
2992 lb0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
2993 lc0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
2994 la0 = one-lb0-lc0
2995
2996 IF(lb0 >= zeromt.and. lb0 <= unpt .and.
2997 . lc0 >= zeromt.and. lc0 <= unpt .and.
2998 . la0 >= zeromt.and. la0 <= unpt)THEN
2999 intersect = 1
3000 ENDIF
3001 ENDIF
3002 ENDIF
3003C
3004 RETURN
3005 END
3006!||====================================================================
3007!|| i25glob_22 ../engine/source/interfaces/int25/i25dst3_22.F
3008!||--- called by ------------------------------------------------------
3009!|| i25comp_2 ../engine/source/interfaces/int25/i25comp_2.F
3010!||--- uses -----------------------------------------------------
3011!|| tri7box ../engine/share/modules/tri7box.F
3012!||====================================================================
3013 SUBROUTINE i25glob_22(
3014 1 JLT ,CAND_N_N ,CAND_E_N ,CAND_E ,
3015 2 NIN ,NSN ,IX1 ,IX2 ,IX3 ,
3016 3 IX4 ,NSVG ,STIF ,INACTI ,MSEGLO ,
3017 4 IRTLM ,TIME_S ,ITAB ,SUBTRIA,
3018 5 FAR ,PENT ,LB ,LC ,
3019 9 INDEX ,FARM ,PENM ,LBM ,
3020 A LCM )
3021C-----------------------------------------------
3022C M o d u l e s
3023C-----------------------------------------------
3024 USE tri7box
3025C-----------------------------------------------
3026C I m p l i c i t T y p e s
3027C-----------------------------------------------
3028#include "implicit_f.inc"
3029C-----------------------------------------------
3030C G l o b a l P a r a m e t e r s
3031C-----------------------------------------------
3032#include "mvsiz_p.inc"
3033#include "comlock.inc"
3034C-----------------------------------------------
3035C D u m m y A r g u m e n t s
3036C-----------------------------------------------
3037 INTEGER JLT, NIN, NSN, INACTI,
3038 . CAND_N_N(*), CAND_E_N(*), CAND_E(*), ITAB(*), SUBTRIA(MVSIZ)
3039 INTEGER NSVG(MVSIZ), IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
3040 . mseglo(*), irtlm(4,nsn) ,index(*), farm(4,*),
3041 . far(mvsiz)
3042 my_real
3043 . stif(*), time_s(*),
3044 . pent(mvsiz),
3045 . lb(mvsiz), lc(mvsiz),
3046 . penm(4,*), lbm(4,*), lcm(4,*)
3047C-----------------------------------------------
3048C L o c a l V a r i a b l e s
3049C-----------------------------------------------
3050 INTEGER IT, I, J
3051C--------------------------------------------------------
3052C Back to global arrays
3053C--------------------------------------------------------
3054 DO i=1,jlt
3055 j=index(i)
3056 it=subtria(i)
3057
3058 farm(1:4,j)=0
3059 penm(1:4,j)=zero
3060 lbm(1:4,j) =zero
3061 lcm(1:4,j) =zero
3062
3063 IF(it/=0)THEN
3064 cand_e(j)=cand_e_n(i)
3065
3066 farm(it,j)=far(i)
3067 penm(it,j)=pent(i)
3068 lbm(it,j) =lb(i)
3069 lcm(it,j) =lc(i)
3070 END IF
3071
3072 END DO
3073 RETURN
3074 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:1898
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:3021
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:2294
intersecv0_25
subroutine intersecv0_25(ixx3, ixx4, intersect, v12, v23, v34, v41, vo12, vo23, vo34, vo41, xx1, yy1, zz1, xx2, yy2, zz2, xx3, yy3, zz3, xx4, yy4, zz4, xx5, yy5, zz5, xoi, yoi, zoi, xo1, xo2, xo3, xo4, xo5, yo1, yo2, yo3, yo4, yo5, zo1, zo2, zo3, zo4, zo5, xi, yi, zi, zeromt, unpt)
Definition i25dst3_22.F:2684
interseca_25
subroutine interseca_25(ixx3, ixx4, interseca, v12, v23, v34, v41, vo12, vo23, vo34, vo41, xx1, yy1, zz1, xx2, yy2, zz2, xx3, yy3, zz3, xx4, yy4, zz4, xx5, yy5, zz5, vxi, vyi, vzi, xo1, xo2, xo3, xo4, xo5, yo1, yo2, yo3, yo4, yo5, zo1, zo2, zo3, zo4, zo5, xi, yi, zi, xoi, yoi, zoi, zerom, unp, zeromt, unpt)
Definition i25dst3_22.F:1532
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
tri7box::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