OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i21dst3.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!|| i21dst3 ../engine/source/interfaces/int21/i21dst3.f
25!||--- called by ------------------------------------------------------
26!|| i21optcd ../engine/source/interfaces/int21/i21optcd.F
27!||====================================================================
28 SUBROUTINE i21dst3(
29 1 JLT ,LIST ,CAND_N ,CAND_E ,IX1 ,
30 2 IX2 ,IX3 ,IX4 ,GAPV ,XI ,
31 3 YI ,ZI ,IRTLM ,CSTS ,DEPTH2 ,
32 4 NNX1 ,NNY1 ,NNZ1 ,NNX2 ,NNY2 ,
33 5 NNZ2 ,NNX3 ,NNY3 ,NNZ3 ,NNX4 ,
34 6 NNY4 ,NNZ4 ,X1 ,Y1 ,Z1 ,
35 7 X2 ,Y2 ,Z2 ,X3 ,Y3 ,
36 8 Z3 ,X4 ,Y4 ,Z4 ,DRAD2 ,
37 9 DGAPLOAD)
38C-----------------------------------------------
39C I m p l i c i t T y p e s
40C-----------------------------------------------
41#include "implicit_f.inc"
42#include "comlock.inc"
43C-----------------------------------------------
44C G l o b a l P a r a m e t e r s
45C-----------------------------------------------
46#include "mvsiz_p.inc"
47C-----------------------------------------------
48C D u m m y A r g u m e n t s
49C-----------------------------------------------
50 INTEGER JLT, LIST(*), CAND_N(*), CAND_E(*),IRTLM(2,*),
51 . IX1(*), IX2(*), IX3(*),IX4(*)
52 my_real , INTENT(IN) :: DGAPLOAD
53 my_real
54 . GAPV(*), CSTS(2,*), DEPTH2, DRAD2,
55 . XI(*) , YI(*), ZI(*),
56 . X1(*), X2(*), X3(*), X4(*),
57 . Y1(*), Y2(*), Y3(*), Y4(*),
58 . Z1(*), Z2(*), Z3(*), Z4(*),
59 . nnx1(*), nnx2(*), nnx3(*), nnx4(*),
60 . nny1(*), nny2(*), nny3(*), nny4(*),
61 . nnz1(*), nnz2(*), nnz3(*), nnz4(*)
62
63C-----------------------------------------------
64C L o c a l V a r i a b l e s
65C-----------------------------------------------
66 INTEGER I, IG, J, IS, L, N, N4N
67 INTEGER NSVG(MVSIZ), FAR(MVSIZ), I4N(MVSIZ),
68 . IRTLM_L(2,MVSIZ)
69 my_real
70 . NX1(MVSIZ), NX2(MVSIZ), NX3(MVSIZ), NX4(MVSIZ),
71 . NY1(MVSIZ), NY2(MVSIZ), NY3(MVSIZ), NY4(MVSIZ),
72 . NZ1(MVSIZ), NZ2(MVSIZ), NZ3(MVSIZ), NZ4(MVSIZ),
73 . LB1(MVSIZ), LB2(MVSIZ), LB3(MVSIZ), LB4(MVSIZ),
74 . LC1(MVSIZ), LC2(MVSIZ), LC3(MVSIZ), LC4(MVSIZ),
75
76 . x0(mvsiz), y0(mvsiz), z0(mvsiz),
77 . p1(mvsiz), p2(mvsiz), p3(mvsiz), p4(mvsiz)
78 my_real
79 . nnx0(mvsiz), nny0(mvsiz), nnz0(mvsiz)
80 my_real
81 . x01(mvsiz), x02(mvsiz), x03(mvsiz), x04(mvsiz),
82 . y01(mvsiz), y02(mvsiz), y03(mvsiz), y04(mvsiz),
83 . z01(mvsiz), z02(mvsiz), z03(mvsiz), z04(mvsiz),
84 . xi1(mvsiz), xi2(mvsiz),
85 . yi1(mvsiz), yi2(mvsiz),
86 . zi1(mvsiz), zi2(mvsiz)
87 my_real
88 . s2,d1,d2,d3,d4,
89 . x12,x23,x34,x41,xi0,
90 . y12,y23,y34,y41,yi0,
91 . z12,z23,z34,z41,zi0,
92 . xn1(mvsiz),yn1(mvsiz),zn1(mvsiz),
93 . xn2(mvsiz),yn2(mvsiz),zn2(mvsiz),
94 . xn3(mvsiz),yn3(mvsiz),zn3(mvsiz),
95 . xn4(mvsiz),yn4(mvsiz),zn4(mvsiz),
96 . sx1(mvsiz),sx2(mvsiz),
97 . sy1(mvsiz),sy2(mvsiz),
98 . sz1(mvsiz),sz2(mvsiz),
99 . xi0v(mvsiz), yi0v(mvsiz), zi0v(mvsiz),
100 . xh(mvsiz), yh(mvsiz), zh(mvsiz),
101 . nx(mvsiz), ny(mvsiz), nz(mvsiz),
102 . gap2(mvsiz), nn,
103 . x0h(mvsiz), y0h(mvsiz), z0h(mvsiz),
104 . x1h(mvsiz), y1h(mvsiz), z1h(mvsiz),
105 . x2h(mvsiz), y2h(mvsiz), z2h(mvsiz), hh(mvsiz), ll,
106 . hxi0, hyi0, hzi0, hxi1, hyi1, hzi1, hxi2, hyi2, hzi2,
107 . hx01, hy01, hz01, hx02, hy02, hz02,
108 . hxn1, hyn1, hzn1, hsx1, hsy1, hsz1, hsx2, hsy2, hsz2,
109 . side, crit, crito, lbo, lco, lb, lc, la, sl,
110 . lambda, xp, yp, zp,gapp2(mvsiz)
111 my_real
112 . cmaj(mvsiz), csts_l(2,mvsiz)
113C--------------------------------------------------------
114 DO i=1,jlt
115 is=list(i)
116 nsvg(i)= cand_n(is)
117c L = CAND_E(IS)
118C
119 END DO
120C
121C
122 DO i=1,jlt
123 irtlm_l(1,i)=0
124 irtlm_l(2,i)=0
125 csts_l(1,i) =-one
126 csts_l(2,i) =-one
127 END DO
128C--------------------------------------------------------
129C CAS DES PAQUETS MIXTES
130C--------------------------------------------------------
131 DO i=1,jlt
132 IF(ix3(i)/=ix4(i))THEN
133 x0(i) = fourth*(x1(i)+x2(i)+x3(i)+x4(i))
134 y0(i) = fourth*(y1(i)+y2(i)+y3(i)+y4(i))
135 z0(i) = fourth*(z1(i)+z2(i)+z3(i)+z4(i))
136 ELSE
137 x0(i) = x3(i)
138 y0(i) = y3(i)
139 z0(i) = z3(i)
140 ENDIF
141 ENDDO
142
143 DO i=1,jlt
144
145 IF(ix3(i)/=ix4(i))THEN
146 nnx0(i)=fourth*(nnx1(i)+nnx2(i)+nnx3(i)+nnx4(i))
147 nny0(i)=fourth*(nny1(i)+nny2(i)+nny3(i)+nny4(i))
148 nnz0(i)=fourth*(nnz1(i)+nnz2(i)+nnz3(i)+nnz4(i))
149 ELSE
150 nnx0(i)=nnx3(i)
151 nny0(i)=nny3(i)
152 nnz0(i)=nnz3(i)
153 ENDIF
154
155 ENDDO
156C
157c---------------------------------------------------------
158c courbure cubique
159c---------------------------------------------------------
160c IF(ICURV == 3)THEN
161c CALL I7CMAJ(JLT ,CMAJ ,IRECT ,NOD_NORMAL,CAND_E,
162c 2 X1 ,X2 ,X3 ,X4 ,
163c 3 Y1 ,Y2 ,Y3 ,Y4 ,
164c 4 Z1 ,Z2 ,Z3 ,Z4 ,
165c 5 NNX1 ,NNX2 ,NNX3 ,NNX4 ,
166c 6 NNY1 ,NNY2 ,NNY3 ,NNY4 ,
167c 7 NNZ1 ,NNZ2 ,NNZ3 ,NNZ4 )
168c ENDIF
169c---------------------------------------------------------
170 DO i=1,jlt
171c CMAJ(I) = ZERO
172C GAP2=(GAPV(I)+CMAJ(I))*(GAPV(I)+CMAJ(I))
173 is=list(i)
174 gap2(i)=gapv(is)*gapv(is)
175 gapp2(i)= (gapv(is)+dgapload)*(gapv(is)+dgapload)
176C
177 x01(i) = x1(i) - x0(i)
178 y01(i) = y1(i) - y0(i)
179 z01(i) = z1(i) - z0(i)
180C
181 x02(i) = x2(i) - x0(i)
182 y02(i) = y2(i) - y0(i)
183 z02(i) = z2(i) - z0(i)
184C
185 x03(i) = x3(i) - x0(i)
186 y03(i) = y3(i) - y0(i)
187 z03(i) = z3(i) - z0(i)
188C
189 x04(i) = x4(i) - x0(i)
190 y04(i) = y4(i) - y0(i)
191 z04(i) = z4(i) - z0(i)
192C
193 ENDDO
194C
195C normales aux triangles (recalculees ici pour pas les stocker).
196 DO i=1,jlt
197
198 xn1(i) = y01(i)*z02(i) - z01(i)*y02(i)
199 yn1(i) = z01(i)*x02(i) - x01(i)*z02(i)
200 zn1(i) = x01(i)*y02(i) - y01(i)*x02(i)
201
202 xn2(i) = y02(i)*z03(i) - z02(i)*y03(i)
203 yn2(i) = z02(i)*x03(i) - x02(i)*z03(i)
204 zn2(i) = x02(i)*y03(i) - y02(i)*x03(i)
205
206 xn3(i) = y03(i)*z04(i) - z03(i)*y04(i)
207 yn3(i) = z03(i)*x04(i) - x03(i)*z04(i)
208 zn3(i) = x03(i)*y04(i) - y03(i)*x04(i)
209
210 xn4(i) = y04(i)*z01(i) - z04(i)*y01(i)
211 yn4(i) = z04(i)*x01(i) - x04(i)*z01(i)
212 zn4(i) = x04(i)*y01(i) - y04(i)*x01(i)
213
214 ENDDO
215C
216C------
217C 1er triangle ...
218 DO i=1,jlt
219C
220 xi0v(i) = x0(i) - xi(i)
221 yi0v(i) = y0(i) - yi(i)
222 zi0v(i) = z0(i) - zi(i)
223C
224 xi1(i) = x1(i) - xi(i)
225 yi1(i) = y1(i) - yi(i)
226 zi1(i) = z1(i) - zi(i)
227C
228 xi2(i) = x2(i) - xi(i)
229 yi2(i) = y2(i) - yi(i)
230 zi2(i) = z2(i) - zi(i)
231C
232 sx1(i) = yi0v(i)*zi1(i) - zi0v(i)*yi1(i)
233 sy1(i) = zi0v(i)*xi1(i) - xi0v(i)*zi1(i)
234 sz1(i) = xi0v(i)*yi1(i) - yi0v(i)*xi1(i)
235C
236 sx2(i) = yi0v(i)*zi2(i) - zi0v(i)*yi2(i)
237 sy2(i) = zi0v(i)*xi2(i) - xi0v(i)*zi2(i)
238 sz2(i) = xi0v(i)*yi2(i) - yi0v(i)*xi2(i)
239C
240 s2 = one/
241 . max(em30,(xn1(i)*xn1(i)+ yn1(i)*yn1(i)+ zn1(i)*zn1(i)))
242 lb1(i) = -(xn1(i)*sx2(i) + yn1(i)*sy2(i) + zn1(i)*sz2(i))*s2
243 lc1(i) = (xn1(i)*sx1(i) + yn1(i)*sy1(i) + zn1(i)*sz1(i))*s2
244C
245 ENDDO
246C
247 DO i=1,jlt
248C
249 nx(i) = xi(i)-(x0(i) + lb1(i)*x01(i) + lc1(i)*x02(i))
250 ny(i) = yi(i)-(y0(i) + lb1(i)*y01(i) + lc1(i)*y02(i))
251 nz(i) = zi(i)-(z0(i) + lb1(i)*z01(i) + lc1(i)*z02(i))
252 hh(i) = nx(i)*xn1(i) + ny(i)*yn1(i) +nz(i)*zn1(i)
253C
254 ENDDO
255C
256C elevation : triangle //
257 DO i=1,jlt
258C
259 ll = hh(i)/
260 . max(em30,nnx1(i)*xn1(i)+nny1(i)*yn1(i)+nnz1(i)*zn1(i))
261 x1h(i)=x1(i)+ll*nnx1(i)
262 y1h(i)=y1(i)+ll*nny1(i)
263 z1h(i)=z1(i)+ll*nnz1(i)
264C
265 ll = hh(i)/
266 . max(em30,nnx2(i)*xn1(i)+nny2(i)*yn1(i)+nnz2(i)*zn1(i))
267 x2h(i)=x2(i)+ll*nnx2(i)
268 y2h(i)=y2(i)+ll*nny2(i)
269 z2h(i)=z2(i)+ll*nnz2(i)
270C
271 ll = hh(i)/
272 . max(em30,nnx0(i)*xn1(i)+nny0(i)*yn1(i)+nnz0(i)*zn1(i))
273 x0h(i)=x0(i)+ll*nnx0(i)
274 y0h(i)=y0(i)+ll*nny0(i)
275 z0h(i)=z0(i)+ll*nnz0(i)
276C
277 ENDDO
278C
279C coordonnees locale dans le triangle //
280 DO i=1,jlt
281C
282 far(i)=0
283C
284 hxi0 = x0h(i) - xi(i)
285 hyi0 = y0h(i) - yi(i)
286 hzi0 = z0h(i) - zi(i)
287C
288 hxi1 = x1h(i) - xi(i)
289 hyi1 = y1h(i) - yi(i)
290 hzi1 = z1h(i) - zi(i)
291C
292 hxi2 = x2h(i) - xi(i)
293 hyi2 = y2h(i) - yi(i)
294 hzi2 = z2h(i) - zi(i)
295C
296 hx01 = x1h(i) - x0h(i)
297 hy01 = y1h(i) - y0h(i)
298 hz01 = z1h(i) - z0h(i)
299C
300 hx02 = x2h(i) - x0h(i)
301 hy02 = y2h(i) - y0h(i)
302 hz02 = z2h(i) - z0h(i)
303C
304 hxn1 = hy01*hz02 - hz01*hy02
305 hyn1 = hz01*hx02 - hx01*hz02
306 hzn1 = hx01*hy02 - hy01*hx02
307C
308 IF(hxn1*xn1(i)+hyn1*yn1(i)+hzn1*zn1(i) <= zero) THEN
309C a optimiser
310 far(i)=1
311 cycle
312 END IF
313C
314 hsx1 = hyi0*hzi1 - hzi0*hyi1
315 hsy1 = hzi0*hxi1 - hxi0*hzi1
316 hsz1 = hxi0*hyi1 - hyi0*hxi1
317C
318 hsx2 = hyi0*hzi2 - hzi0*hyi2
319 hsy2 = hzi0*hxi2 - hxi0*hzi2
320 hsz2 = hxi0*hyi2 - hyi0*hxi2
321C
322 s2 = one/
323 . max(em30,(hxn1*hxn1+ hyn1*hyn1+ hzn1*hzn1))
324 lb1(i) = -(hxn1*hsx2 + hyn1*hsy2 + hzn1*hsz2)*s2
325 lc1(i) = (hxn1*hsx1 + hyn1*hsy1 + hzn1*hsz1)*s2
326C
327 IF(lb1(i) < -zep01 .OR. lc1(i) < -zep01 .OR.
328 . lb1(i)+lc1(i) > onep01) far(i)=1
329 ENDDO
330C
331C necessaire au calcul de distance // normale
332 DO i=1,jlt
333C
334 IF(far(i)==1)cycle
335C
336C normale interpolee dans le triangle //
337 nx(i)=(one-lb1(i)-lc1(i))*nnx0(i)+lb1(i)*nnx1(i)+lc1(i)*nnx2(i)
338 ny(i)=(one-lb1(i)-lc1(i))*nny0(i)+lb1(i)*nny1(i)+lc1(i)*nny2(i)
339 nz(i)=(one-lb1(i)-lc1(i))*nnz0(i)+lb1(i)*nnz1(i)+lc1(i)*nnz2(i)
340C
341C projection sur triangle 1
342 nn = nx(i)*xn1(i)+ ny(i)*yn1(i)+ nz(i)*zn1(i)
343 IF(nn <= zero) THEN
344C cas limite
345 far(i)=1
346 cycle
347 END IF
348 nn = one/max(em30,nn)
349C
350 lambda=(xn1(i)*xi0v(i)+yn1(i)*yi0v(i)+zn1(i)*zi0v(i))*nn
351 xp=xi(i)+lambda*nx(i)
352 yp=yi(i)+lambda*ny(i)
353 zp=zi(i)+lambda*nz(i)
354C
355 xi0v(i) = x0(i) - xp
356 yi0v(i) = y0(i) - yp
357 zi0v(i) = z0(i) - zp
358C
359 xi1(i) = x1(i) - xp
360 yi1(i) = y1(i) - yp
361 zi1(i) = z1(i) - zp
362C
363 xi2(i) = x2(i) - xp
364 yi2(i) = y2(i) - yp
365 zi2(i) = z2(i) - zp
366C
367 sx1(i) = yi0v(i)*zi1(i) - zi0v(i)*yi1(i)
368 sy1(i) = zi0v(i)*xi1(i) - xi0v(i)*zi1(i)
369 sz1(i) = xi0v(i)*yi1(i) - yi0v(i)*xi1(i)
370C
371 sx2(i) = yi0v(i)*zi2(i) - zi0v(i)*yi2(i)
372 sy2(i) = zi0v(i)*xi2(i) - xi0v(i)*zi2(i)
373 sz2(i) = xi0v(i)*yi2(i) - yi0v(i)*xi2(i)
374C
375 nn = one/
376 . max(em30,(xn1(i)*xn1(i)+ yn1(i)*yn1(i)+ zn1(i)*zn1(i)))
377 lb1(i) = -(xn1(i)*sx2(i) + yn1(i)*sy2(i) + zn1(i)*sz2(i))*nn
378 lc1(i) = (xn1(i)*sx1(i) + yn1(i)*sy1(i) + zn1(i)*sz1(i))*nn
379C
380 IF(lb1(i) < -zep01 .OR. lc1(i) < -zep01 .OR.
381 . lb1(i)+lc1(i) > onep01) far(i)=1
382C
383 ENDDO
384C
385 DO i=1,jlt
386C
387 IF(far(i)==1)cycle
388C
389C=======================================================================
390C1 on a garde -zep01 < la,b,c < 0 pour retenir un contact dans le cas
391C ou la surface main est concave
392C
393C / \
394C / \
395C / o \
396C / \
397C
398C=======================================================================
399C2 on rabat la,b,c pour revenir sur l'element lorsqu'on calcule
400C le critere:
401C vu dans le cas de surface main convexe avec grand element,
402C la projection retenue est tres eloignee et peut-etre
403C hors de la boite lors du retri => parith on, KO
404C
405C /
406C /
407C /
408C ------
409C ...
410C /
411C / o
412C
413C=======================================================================
414C
415 lb = min(one,max(lb1(i),zero))
416 lc = min(one,max(lc1(i),zero))
417 la = min(one,max(one-lb1(i)-lc1(i),zero))
418
419 sl=one/max(em20,la+lb+lc)
420 lb = lb*sl
421 lc = lc*sl
422 la = la*sl
423
424 nx1(i) = xi(i)-(la*x0(i) + lb*x1(i) + lc*x2(i))
425 ny1(i) = yi(i)-(la*y0(i) + lb*y1(i) + lc*y2(i))
426 nz1(i) = zi(i)-(la*z0(i) + lb*z1(i) + lc*z2(i))
427 p1(i) = nx1(i)*nx1(i) + ny1(i)*ny1(i) +nz1(i)*nz1(i)
428C
429C NX(I)=(1-LB1(I)-LC1(I))*NNX0(I)+LB1(I)*NNX1(I)+LC1(I)*NNX2(I)
430C NY(I)=(1-LB1(I)-LC1(I))*NNY0(I)+LB1(I)*NNY1(I)+LC1(I)*NNY2(I)
431C NZ(I)=(1-LB1(I)-LC1(I))*NNZ0(I)+LB1(I)*NNZ1(I)+LC1(I)*NNZ2(I)
432C SIDE=NX1(I)*NX(I)+NY1(I)*NY(I)+NZ1(I)*NZ(I)
433
434 side=sign(one,nx1(i)*xn1(i)+ny1(i)*yn1(i)+nz1(i)*zn1(i))
435 IF((side >= zero .AND. p1(i) < gap2(i)) .OR.
436 . (side < zero .AND. p1(i) < depth2))THEN
437 crit = abs(third-lb1(i))
438 . + abs(third-lc1(i))
439 . + abs(two_third-lb1(i)-lc1(i))
440 lbo = csts_l(1,i)
441 lco = csts_l(2,i)
442 crito = abs(third-lbo)
443 . + abs(third-lco)
444 . + abs(two_third-lbo-lco)
445 IF(crit < crito)THEN
446 is=list(i)
447 irtlm_l(1,i) = cand_e(is)
448 irtlm_l(2,i) = nint(side)
449 csts_l(1,i) = lb1(i)
450 csts_l(2,i) = lc1(i)
451 END IF
452 ELSEIF(side >= zero .AND. p1(i) < max(drad2,gapp2(i)))THEN
453 crit = abs(third-lb1(i))
454 . + abs(third-lc1(i))
455 . + abs(two_third-lb1(i)-lc1(i))
456 lbo = csts_l(1,i)
457 lco = csts_l(2,i)
458 crito = abs(third-lbo)
459 . + abs(third-lco)
460 . + abs(two_third-lbo-lco)
461 IF(crit < crito)THEN
462 is=list(i)
463 irtlm_l(1,i) = -cand_e(is)
464 irtlm_l(2,i) = nint(side)
465 csts_l(1,i) = lb1(i)
466 csts_l(2,i) = lc1(i)
467 END IF
468 END IF
469C
470 ENDDO
471C------
472 n4n=0
473 DO i=1,jlt
474 IF(ix3(i)/=ix4(i))THEN
475 n4n = n4n+1
476 i4n(n4n)=i
477 ELSE
478 ENDIF
479 ENDDO
480C------
481C 2eme triangle ...
482 DO n=1,n4n
483C
484 i=i4n(n)
485C
486 xi0v(i) = x0(i) - xi(i)
487 yi0v(i) = y0(i) - yi(i)
488 zi0v(i) = z0(i) - zi(i)
489C
490 xi1(i) = x2(i) - xi(i)
491 yi1(i) = y2(i) - yi(i)
492 zi1(i) = z2(i) - zi(i)
493C
494 xi2(i) = x3(i) - xi(i)
495 yi2(i) = y3(i) - yi(i)
496 zi2(i) = z3(i) - zi(i)
497C
498 sx1(i) = yi0v(i)*zi1(i) - zi0v(i)*yi1(i)
499 sy1(i) = zi0v(i)*xi1(i) - xi0v(i)*zi1(i)
500 sz1(i) = xi0v(i)*yi1(i) - yi0v(i)*xi1(i)
501C
502 sx2(i) = yi0v(i)*zi2(i) - zi0v(i)*yi2(i)
503 sy2(i) = zi0v(i)*xi2(i) - xi0v(i)*zi2(i)
504 sz2(i) = xi0v(i)*yi2(i) - yi0v(i)*xi2(i)
505C
506 s2 = one/
507 . max(em30,(xn2(i)*xn2(i)+ yn2(i)*yn2(i)+ zn2(i)*zn2(i)))
508 lb2(i) = -(xn2(i)*sx2(i) + yn2(i)*sy2(i) + zn2(i)*sz2(i))*s2
509 lc2(i) = (xn2(i)*sx1(i) + yn2(i)*sy1(i) + zn2(i)*sz1(i))*s2
510C
511 ENDDO
512C
513 DO n=1,n4n
514C
515 i=i4n(n)
516C
517 nx(i) = xi(i)-(x0(i) + lb2(i)*x02(i) + lc2(i)*x03(i))
518 ny(i) = yi(i)-(y0(i) + lb2(i)*y02(i) + lc2(i)*y03(i))
519 nz(i) = zi(i)-(z0(i) + lb2(i)*z02(i) + lc2(i)*z03(i))
520 hh(i) = nx(i)*xn2(i) + ny(i)*yn2(i) +nz(i)*zn2(i)
521C
522 ENDDO
523C
524C elevation : triangle //
525 DO n=1,n4n
526C
527 i=i4n(n)
528C
529 ll = hh(i)/
530 . max(em30,nnx2(i)*xn2(i)+nny2(i)*yn2(i)+nnz2(i)*zn2(i))
531 x1h(i)=x2(i)+ll*nnx2(i)
532 y1h(i)=y2(i)+ll*nny2(i)
533 z1h(i)=z2(i)+ll*nnz2(i)
534C
535 ll = hh(i)/
536 . max(em30,nnx3(i)*xn2(i)+nny3(i)*yn2(i)+nnz3(i)*zn2(i))
537 x2h(i)=x3(i)+ll*nnx3(i)
538 y2h(i)=y3(i)+ll*nny3(i)
539 z2h(i)=z3(i)+ll*nnz3(i)
540C
541 ll = hh(i)/
542 . max(em30,nnx0(i)*xn2(i)+nny0(i)*yn2(i)+nnz0(i)*zn2(i))
543 x0h(i)=x0(i)+ll*nnx0(i)
544 y0h(i)=y0(i)+ll*nny0(i)
545 z0h(i)=z0(i)+ll*nnz0(i)
546C
547 ENDDO
548C
549C coordonnees locale dans le triangle //
550 DO n=1,n4n
551C
552 i=i4n(n)
553C
554 far(i)=0
555C
556 hxi0 = x0h(i) - xi(i)
557 hyi0 = y0h(i) - yi(i)
558 hzi0 = z0h(i) - zi(i)
559C
560 hxi1 = x1h(i) - xi(i)
561 hyi1 = y1h(i) - yi(i)
562 hzi1 = z1h(i) - zi(i)
563C
564 hxi2 = x2h(i) - xi(i)
565 hyi2 = y2h(i) - yi(i)
566 hzi2 = z2h(i) - zi(i)
567C
568 hx01 = x1h(i) - x0h(i)
569 hy01 = y1h(i) - y0h(i)
570 hz01 = z1h(i) - z0h(i)
571C
572 hx02 = x2h(i) - x0h(i)
573 hy02 = y2h(i) - y0h(i)
574 hz02 = z2h(i) - z0h(i)
575C
576 hxn1 = hy01*hz02 - hz01*hy02
577 hyn1 = hz01*hx02 - hx01*hz02
578 hzn1 = hx01*hy02 - hy01*hx02
579C
580 IF(hxn1*xn2(i)+hyn1*yn2(i)+hzn1*zn2(i) <= zero) THEN
581C a optimiser
582 far(i)=1
583 cycle
584 END IF
585C
586 hsx1 = hyi0*hzi1 - hzi0*hyi1
587 hsy1 = hzi0*hxi1 - hxi0*hzi1
588 hsz1 = hxi0*hyi1 - hyi0*hxi1
589C
590 hsx2 = hyi0*hzi2 - hzi0*hyi2
591 hsy2 = hzi0*hxi2 - hxi0*hzi2
592 hsz2 = hxi0*hyi2 - hyi0*hxi2
593C
594 s2 = one/
595 . max(em30,(hxn1*hxn1+ hyn1*hyn1+ hzn1*hzn1))
596 lb2(i) = -(hxn1*hsx2 + hyn1*hsy2 + hzn1*hsz2)*s2
597 lc2(i) = (hxn1*hsx1 + hyn1*hsy1 + hzn1*hsz1)*s2
598C
599 IF(lb2(i) < -zep01 .OR. lc2(i) < -zep01 .OR.
600 . lb2(i)+lc2(i) > onep01) far(i)=1
601 ENDDO
602C
603C necessaire au calcul de distance // normale
604 DO n=1,n4n
605C
606 i=i4n(n)
607C
608 IF(far(i)==1)cycle
609C
610C normale interpolee dans le triangle //
611 nx(i)=(one-lb2(i)-lc2(i))*nnx0(i)+lb2(i)*nnx2(i)+lc2(i)*nnx3(i)
612 ny(i)=(one-lb2(i)-lc2(i))*nny0(i)+lb2(i)*nny2(i)+lc2(i)*nny3(i)
613 nz(i)=(one-lb2(i)-lc2(i))*nnz0(i)+lb2(i)*nnz2(i)+lc2(i)*nnz3(i)
614C
615C projection sur triangle 2
616 nn = nx(i)*xn2(i)+ ny(i)*yn2(i)+ nz(i)*zn2(i)
617 IF(nn <= zero) THEN
618C cas limite
619 far(i)=1
620 cycle
621 END IF
622 nn = one/max(em30,nn)
623C
624 lambda=(xn2(i)*xi0v(i)+yn2(i)*yi0v(i)+zn2(i)*zi0v(i))*nn
625 xp=xi(i)+lambda*nx(i)
626 yp=yi(i)+lambda*ny(i)
627 zp=zi(i)+lambda*nz(i)
628C
629 xi0v(i) = x0(i) - xp
630 yi0v(i) = y0(i) - yp
631 zi0v(i) = z0(i) - zp
632C
633 xi1(i) = x2(i) - xp
634 yi1(i) = y2(i) - yp
635 zi1(i) = z2(i) - zp
636C
637 xi2(i) = x3(i) - xp
638 yi2(i) = y3(i) - yp
639 zi2(i) = z3(i) - zp
640C
641 sx1(i) = yi0v(i)*zi1(i) - zi0v(i)*yi1(i)
642 sy1(i) = zi0v(i)*xi1(i) - xi0v(i)*zi1(i)
643 sz1(i) = xi0v(i)*yi1(i) - yi0v(i)*xi1(i)
644C
645 sx2(i) = yi0v(i)*zi2(i) - zi0v(i)*yi2(i)
646 sy2(i) = zi0v(i)*xi2(i) - xi0v(i)*zi2(i)
647 sz2(i) = xi0v(i)*yi2(i) - yi0v(i)*xi2(i)
648C
649 nn = one/
650 . max(em30,(xn2(i)*xn2(i)+ yn2(i)*yn2(i)+ zn2(i)*zn2(i)))
651 lb2(i) = -(xn2(i)*sx2(i) + yn2(i)*sy2(i) + zn2(i)*sz2(i))*nn
652 lc2(i) = (xn2(i)*sx1(i) + yn2(i)*sy1(i) + zn2(i)*sz1(i))*nn
653C
654 IF(lb2(i) < -zep01 .OR. lc2(i) < -zep01 .OR.
655 . lb2(i)+lc2(i) > onep01) far(i)=1
656C
657 ENDDO
658C
659 DO n=1,n4n
660C
661 i=i4n(n)
662C
663 IF(far(i)==1)cycle
664C
665 lb = min(one,max(lb2(i),zero))
666 lc = min(one,max(lc2(i),zero))
667 la = min(one,max(one-lb2(i)-lc2(i),zero))
668
669 sl=one/max(em20,la+lb+lc)
670 lb = lb*sl
671 lc = lc*sl
672 la = la*sl
673
674 nx2(i) = xi(i)-(la*x0(i) + lb*x2(i) + lc*x3(i))
675 ny2(i) = yi(i)-(la*y0(i) + lb*y2(i) + lc*y3(i))
676 nz2(i) = zi(i)-(la*z0(i) + lb*z2(i) + lc*z3(i))
677 p2(i) = nx2(i)*nx2(i) + ny2(i)*ny2(i) +nz2(i)*nz2(i)
678C
679C NX(I)=(1-LB2(I)-LC2(I))*NNX0(I)+LB2(I)*NNX2(I)+LC2(I)*NNX3(I)
680C NY(I)=(1-LB2(I)-LC2(I))*NNY0(I)+LB2(I)*NNY2(I)+LC2(I)*NNY3(I)
681C NZ(I)=(1-LB2(I)-LC2(I))*NNZ0(I)+LB2(I)*NNZ2(I)+LC2(I)*NNZ3(I)
682C SIDE=NX2(I)*NX(I)+NY2(I)*NY(I)+NZ2(I)*NZ(I)
683
684 side=sign(one,nx2(i)*xn2(i)+ny2(i)*yn2(i)+nz2(i)*zn2(i))
685 IF((side >= zero .AND. p2(i) < gap2(i)) .OR.
686 . (side < zero .AND. p2(i) < depth2))THEN
687 crit = abs(third-lb2(i))
688 . + abs(third-lc2(i))
689 . + abs(two_third-lb2(i)-lc2(i))
690 lbo = csts_l(1,i)
691 lco = csts_l(2,i)
692 crito = abs(third-lbo)
693 . + abs(third-lco)
694 . + abs(two_third-lbo-lco)
695 IF(crit < crito)THEN
696 is=list(i)
697 irtlm_l(1,i) = cand_e(is)
698 irtlm_l(2,i) = nint(two*side)
699 csts_l(1,i) = lb2(i)
700 csts_l(2,i) = lc2(i)
701 END IF
702 ELSEIF(side >= zero .AND. p2(i) < max(drad2,gapp2(i)))THEN
703 crit = abs(third-lb2(i))
704 . + abs(third-lc2(i))
705 . + abs(two_third-lb2(i)-lc2(i))
706 lbo = csts_l(1,i)
707 lco = csts_l(2,i)
708 crito = abs(third-lbo)
709 . + abs(third-lco)
710 . + abs(two_third-lbo-lco)
711 IF(crit < crito)THEN
712 is=list(i)
713 irtlm_l(1,i) = -cand_e(is)
714 irtlm_l(2,i) = nint(two*side)
715 csts_l(1,i) = lb2(i)
716 csts_l(2,i) = lc2(i)
717 END IF
718 END IF
719C
720 ENDDO
721C------
722C 3eme triangle ...
723 DO n=1,n4n
724C
725 i=i4n(n)
726C
727 xi0v(i) = x0(i) - xi(i)
728 yi0v(i) = y0(i) - yi(i)
729 zi0v(i) = z0(i) - zi(i)
730C
731 xi1(i) = x3(i) - xi(i)
732 yi1(i) = y3(i) - yi(i)
733 zi1(i) = z3(i) - zi(i)
734C
735 xi2(i) = x4(i) - xi(i)
736 yi2(i) = y4(i) - yi(i)
737 zi2(i) = z4(i) - zi(i)
738C
739 sx1(i) = yi0v(i)*zi1(i) - zi0v(i)*yi1(i)
740 sy1(i) = zi0v(i)*xi1(i) - xi0v(i)*zi1(i)
741 sz1(i) = xi0v(i)*yi1(i) - yi0v(i)*xi1(i)
742C
743 sx2(i) = yi0v(i)*zi2(i) - zi0v(i)*yi2(i)
744 sy2(i) = zi0v(i)*xi2(i) - xi0v(i)*zi2(i)
745 sz2(i) = xi0v(i)*yi2(i) - yi0v(i)*xi2(i)
746C
747 s2 = one/
748 . max(em30,(xn3(i)*xn3(i)+ yn3(i)*yn3(i)+ zn3(i)*zn3(i)))
749 lb3(i) = -(xn3(i)*sx2(i) + yn3(i)*sy2(i) + zn3(i)*sz2(i))*s2
750 lc3(i) = (xn3(i)*sx1(i) + yn3(i)*sy1(i) + zn3(i)*sz1(i))*s2
751C
752 ENDDO
753C
754 DO n=1,n4n
755C
756 i=i4n(n)
757C
758 nx(i) = xi(i)-(x0(i) + lb3(i)*x03(i) + lc3(i)*x04(i))
759 ny(i) = yi(i)-(y0(i) + lb3(i)*y03(i) + lc3(i)*y04(i))
760 nz(i) = zi(i)-(z0(i) + lb3(i)*z03(i) + lc3(i)*z04(i))
761 hh(i) = nx(i)*xn3(i) + ny(i)*yn3(i) +nz(i)*zn3(i)
762C
763 ENDDO
764C
765C elevation : triangle //
766 DO n=1,n4n
767C
768 i=i4n(n)
769C
770 ll = hh(i)/
771 . max(em30,nnx3(i)*xn3(i)+nny3(i)*yn3(i)+nnz3(i)*zn3(i))
772 x1h(i)=x3(i)+ll*nnx3(i)
773 y1h(i)=y3(i)+ll*nny3(i)
774 z1h(i)=z3(i)+ll*nnz3(i)
775C
776 ll = hh(i)/
777 . max(em30,nnx4(i)*xn3(i)+nny4(i)*yn3(i)+nnz4(i)*zn3(i))
778 x2h(i)=x4(i)+ll*nnx4(i)
779 y2h(i)=y4(i)+ll*nny4(i)
780 z2h(i)=z4(i)+ll*nnz4(i)
781C
782 ll = hh(i)/
783 . max(em30,nnx0(i)*xn3(i)+nny0(i)*yn3(i)+nnz0(i)*zn3(i))
784 x0h(i)=x0(i)+ll*nnx0(i)
785 y0h(i)=y0(i)+ll*nny0(i)
786 z0h(i)=z0(i)+ll*nnz0(i)
787C
788 ENDDO
789C
790C coordonnees locale dans le triangle //
791 DO n=1,n4n
792C
793 i=i4n(n)
794C
795 far(i)=0
796C
797 hxi0 = x0h(i) - xi(i)
798 hyi0 = y0h(i) - yi(i)
799 hzi0 = z0h(i) - zi(i)
800C
801 hxi1 = x1h(i) - xi(i)
802 hyi1 = y1h(i) - yi(i)
803 hzi1 = z1h(i) - zi(i)
804C
805 hxi2 = x2h(i) - xi(i)
806 hyi2 = y2h(i) - yi(i)
807 hzi2 = z2h(i) - zi(i)
808C
809 hx01 = x1h(i) - x0h(i)
810 hy01 = y1h(i) - y0h(i)
811 hz01 = z1h(i) - z0h(i)
812C
813 hx02 = x2h(i) - x0h(i)
814 hy02 = y2h(i) - y0h(i)
815 hz02 = z2h(i) - z0h(i)
816C
817 hxn1 = hy01*hz02 - hz01*hy02
818 hyn1 = hz01*hx02 - hx01*hz02
819 hzn1 = hx01*hy02 - hy01*hx02
820C
821 IF(hxn1*xn3(i)+hyn1*yn3(i)+hzn1*zn3(i) <= zero) THEN
822C a optimiser
823 far(i)=1
824 cycle
825 END IF
826C
827 hsx1 = hyi0*hzi1 - hzi0*hyi1
828 hsy1 = hzi0*hxi1 - hxi0*hzi1
829 hsz1 = hxi0*hyi1 - hyi0*hxi1
830C
831 hsx2 = hyi0*hzi2 - hzi0*hyi2
832 hsy2 = hzi0*hxi2 - hxi0*hzi2
833 hsz2 = hxi0*hyi2 - hyi0*hxi2
834C
835 s2 = one/
836 . max(em30,(hxn1*hxn1+ hyn1*hyn1+ hzn1*hzn1))
837 lb3(i) = -(hxn1*hsx2 + hyn1*hsy2 + hzn1*hsz2)*s2
838 lc3(i) = (hxn1*hsx1 + hyn1*hsy1 + hzn1*hsz1)*s2
839C
840 IF(lb3(i) < -zep01 .OR. lc3(i) < -zep01 .OR.
841 . lb3(i)+lc3(i) > onep01) far(i)=1
842 ENDDO
843C
844C necessaire au calcul de distance // normale
845 DO n=1,n4n
846C
847 i=i4n(n)
848C
849 IF(far(i)==1)cycle
850C
851C normale interpolee dans le triangle //
852 nx(i)=(one-lb3(i)-lc3(i))*nnx0(i)+lb3(i)*nnx3(i)+lc3(i)*nnx4(i)
853 ny(i)=(one-lb3(i)-lc3(i))*nny0(i)+lb3(i)*nny3(i)+lc3(i)*nny4(i)
854 nz(i)=(one-lb3(i)-lc3(i))*nnz0(i)+lb3(i)*nnz3(i)+lc3(i)*nnz4(i)
855C
856C projection sur triangle 3
857 nn = nx(i)*xn3(i)+ ny(i)*yn3(i)+ nz(i)*zn3(i)
858 IF(nn <= zero) THEN
859C cas limite
860 far(i)=1
861 cycle
862 END IF
863 nn = one/max(em30,nn)
864C
865 lambda=(xn3(i)*xi0v(i)+yn3(i)*yi0v(i)+zn3(i)*zi0v(i))*nn
866 xp=xi(i)+lambda*nx(i)
867 yp=yi(i)+lambda*ny(i)
868 zp=zi(i)+lambda*nz(i)
869C
870 xi0v(i) = x0(i) - xp
871 yi0v(i) = y0(i) - yp
872 zi0v(i) = z0(i) - zp
873C
874 xi1(i) = x3(i) - xp
875 yi1(i) = y3(i) - yp
876 zi1(i) = z3(i) - zp
877C
878 xi2(i) = x4(i) - xp
879 yi2(i) = y4(i) - yp
880 zi2(i) = z4(i) - zp
881C
882 sx1(i) = yi0v(i)*zi1(i) - zi0v(i)*yi1(i)
883 sy1(i) = zi0v(i)*xi1(i) - xi0v(i)*zi1(i)
884 sz1(i) = xi0v(i)*yi1(i) - yi0v(i)*xi1(i)
885C
886 sx2(i) = yi0v(i)*zi2(i) - zi0v(i)*yi2(i)
887 sy2(i) = zi0v(i)*xi2(i) - xi0v(i)*zi2(i)
888 sz2(i) = xi0v(i)*yi2(i) - yi0v(i)*xi2(i)
889C
890 nn = one/
891 . max(em30,(xn3(i)*xn3(i)+ yn3(i)*yn3(i)+ zn3(i)*zn3(i)))
892 lb3(i) = -(xn3(i)*sx2(i) + yn3(i)*sy2(i) + zn3(i)*sz2(i))*nn
893 lc3(i) = (xn3(i)*sx1(i) + yn3(i)*sy1(i) + zn3(i)*sz1(i))*nn
894C
895 IF(lb3(i) < -zep01 .OR. lc3(i) < -zep01 .OR.
896 . lb3(i)+lc3(i) > onep01) far(i)=1
897C
898 ENDDO
899C
900 DO n=1,n4n
901C
902 i=i4n(n)
903C
904 IF(far(i)==1)cycle
905
906 lb = min(one,max(lb3(i),zero))
907 lc = min(one,max(lc3(i),zero))
908 la = min(one,max(one-lb3(i)-lc3(i),zero))
909
910 sl=one/max(em20,la+lb+lc)
911 lb = lb*sl
912 lc = lc*sl
913 la = la*sl
914
915 nx3(i) = xi(i)-(la*x0(i) + lb*x3(i) + lc*x4(i))
916 ny3(i) = yi(i)-(la*y0(i) + lb*y3(i) + lc*y4(i))
917 nz3(i) = zi(i)-(la*z0(i) + lb*z3(i) + lc*z4(i))
918 p3(i) = nx3(i)*nx3(i) + ny3(i)*ny3(i) +nz3(i)*nz3(i)
919C
920C NX(I)=(1-LB3(I)-LC3(I))*NNX0(I)+LB3(I)*NNX3(I)+LC3(I)*NNX4(I)
921C NY(I)=(1-LB3(I)-LC3(I))*NNY0(I)+LB3(I)*NNY3(I)+LC3(I)*NNY4(I)
922C NZ(I)=(1-LB3(I)-LC3(I))*NNZ0(I)+LB3(I)*NNZ3(I)+LC3(I)*NNZ4(I)
923C SIDE=NX3(I)*NX(I)+NY3(I)*NY(I)+NZ3(I)*NZ(I)
924
925 side=sign(one,nx3(i)*xn3(i)+ny3(i)*yn3(i)+nz3(i)*zn3(i))
926 IF((side >= zero .AND. p3(i) < gap2(i)) .OR.
927 . (side < zero .AND. p3(i) < depth2))THEN
928 crit = abs(third-lb3(i))
929 . + abs(third-lc3(i))
930 . + abs(two_third-lb3(i)-lc3(i))
931 lbo = csts_l(1,i)
932 lco = csts_l(2,i)
933 crito = abs(third-lbo)
934 . + abs(third-lco)
935 . + abs(two_third-lbo-lco)
936 IF(crit < crito)THEN
937 is=list(i)
938 irtlm_l(1,i) = cand_e(is)
939 irtlm_l(2,i) = nint(three*side)
940 csts_l(1,i) = lb3(i)
941 csts_l(2,i) = lc3(i)
942 END IF
943 ELSEIF(side >= zero .AND. p3(i) < max(drad2,gapp2(i)))THEN
944 crit = abs(third-lb3(i))
945 . + abs(third-lc3(i))
946 . + abs(two_third-lb3(i)-lc3(i))
947 lbo = csts_l(1,i)
948 lco = csts_l(2,i)
949 crito = abs(third-lbo)
950 . + abs(third-lco)
951 . + abs(two_third-lbo-lco)
952 IF(crit < crito)THEN
953 is=list(i)
954 irtlm_l(1,i) = -cand_e(is)
955 irtlm_l(2,i) = nint(three*side)
956 csts_l(1,i) = lb3(i)
957 csts_l(2,i) = lc3(i)
958 END IF
959 END IF
960C
961 END DO
962C------
963C 4eme triangle ...
964 DO n=1,n4n
965C
966 i=i4n(n)
967C
968 xi0v(i) = x0(i) - xi(i)
969 yi0v(i) = y0(i) - yi(i)
970 zi0v(i) = z0(i) - zi(i)
971C
972 xi1(i) = x4(i) - xi(i)
973 yi1(i) = y4(i) - yi(i)
974 zi1(i) = z4(i) - zi(i)
975C
976 xi2(i) = x1(i) - xi(i)
977 yi2(i) = y1(i) - yi(i)
978 zi2(i) = z1(i) - zi(i)
979C
980 sx1(i) = yi0v(i)*zi1(i) - zi0v(i)*yi1(i)
981 sy1(i) = zi0v(i)*xi1(i) - xi0v(i)*zi1(i)
982 sz1(i) = xi0v(i)*yi1(i) - yi0v(i)*xi1(i)
983C
984 sx2(i) = yi0v(i)*zi2(i) - zi0v(i)*yi2(i)
985 sy2(i) = zi0v(i)*xi2(i) - xi0v(i)*zi2(i)
986 sz2(i) = xi0v(i)*yi2(i) - yi0v(i)*xi2(i)
987C
988 s2 = one/
989 . max(em30,(xn4(i)*xn4(i)+ yn4(i)*yn4(i)+ zn4(i)*zn4(i)))
990 lb4(i) = -(xn4(i)*sx2(i) + yn4(i)*sy2(i) + zn4(i)*sz2(i))*s2
991 lc4(i) = (xn4(i)*sx1(i) + yn4(i)*sy1(i) + zn4(i)*sz1(i))*s2
992C
993 ENDDO
994C
995 DO n=1,n4n
996C
997 i=i4n(n)
998C
999 nx(i) = xi(i)-(x0(i) + lb4(i)*x04(i) + lc4(i)*x01(i))
1000 ny(i) = yi(i)-(y0(i) + lb4(i)*y04(i) + lc4(i)*y01(i))
1001 nz(i) = zi(i)-(z0(i) + lb4(i)*z04(i) + lc4(i)*z01(i))
1002 hh(i) = nx(i)*xn4(i) + ny(i)*yn4(i) +nz(i)*zn4(i)
1003C
1004 ENDDO
1005C
1006C elevation : triangle //
1007 DO n=1,n4n
1008C
1009 i=i4n(n)
1010C
1011 ll = hh(i)/
1012 . max(em30,nnx4(i)*xn4(i)+nny4(i)*yn4(i)+nnz4(i)*zn4(i))
1013 x1h(i)=x4(i)+ll*nnx4(i)
1014 y1h(i)=y4(i)+ll*nny4(i)
1015 z1h(i)=z4(i)+ll*nnz4(i)
1016C
1017 ll = hh(i)/
1018 . max(em30,nnx1(i)*xn4(i)+nny1(i)*yn4(i)+nnz1(i)*zn4(i))
1019 x2h(i)=x1(i)+ll*nnx1(i)
1020 y2h(i)=y1(i)+ll*nny1(i)
1021 z2h(i)=z1(i)+ll*nnz1(i)
1022C
1023 ll = hh(i)/
1024 . max(em30,nnx0(i)*xn4(i)+nny0(i)*yn4(i)+nnz0(i)*zn4(i))
1025 x0h(i)=x0(i)+ll*nnx0(i)
1026 y0h(i)=y0(i)+ll*nny0(i)
1027 z0h(i)=z0(i)+ll*nnz0(i)
1028C
1029 ENDDO
1030C
1031C coordonnees locale dans le triangle //
1032 DO n=1,n4n
1033C
1034 i=i4n(n)
1035C
1036 far(i)=0
1037C
1038 hxi0 = x0h(i) - xi(i)
1039 hyi0 = y0h(i) - yi(i)
1040 hzi0 = z0h(i) - zi(i)
1041C
1042 hxi1 = x1h(i) - xi(i)
1043 hyi1 = y1h(i) - yi(i)
1044 hzi1 = z1h(i) - zi(i)
1045C
1046 hxi2 = x2h(i) - xi(i)
1047 hyi2 = y2h(i) - yi(i)
1048 hzi2 = z2h(i) - zi(i)
1049C
1050 hx01 = x1h(i) - x0h(i)
1051 hy01 = y1h(i) - y0h(i)
1052 hz01 = z1h(i) - z0h(i)
1053C
1054 hx02 = x2h(i) - x0h(i)
1055 hy02 = y2h(i) - y0h(i)
1056 hz02 = z2h(i) - z0h(i)
1057C
1058 hxn1 = hy01*hz02 - hz01*hy02
1059 hyn1 = hz01*hx02 - hx01*hz02
1060 hzn1 = hx01*hy02 - hy01*hx02
1061C
1062 IF(hxn1*xn4(i)+hyn1*yn4(i)+hzn1*zn4(i) <= zero) THEN
1063C a optimiser
1064 far(i)=1
1065 cycle
1066 END IF
1067C
1068 hsx1 = hyi0*hzi1 - hzi0*hyi1
1069 hsy1 = hzi0*hxi1 - hxi0*hzi1
1070 hsz1 = hxi0*hyi1 - hyi0*hxi1
1071C
1072 hsx2 = hyi0*hzi2 - hzi0*hyi2
1073 hsy2 = hzi0*hxi2 - hxi0*hzi2
1074 hsz2 = hxi0*hyi2 - hyi0*hxi2
1075C
1076 s2 = one/
1077 . max(em30,(hxn1*hxn1+ hyn1*hyn1+ hzn1*hzn1))
1078 lb4(i) = -(hxn1*hsx2 + hyn1*hsy2 + hzn1*hsz2)*s2
1079 lc4(i) = (hxn1*hsx1 + hyn1*hsy1 + hzn1*hsz1)*s2
1080C
1081 IF(lb4(i) < -zep01 .OR. lc4(i) < -zep01 .OR.
1082 . lb4(i)+lc4(i) > onep01) far(i)=1
1083 ENDDO
1084C
1085C necessaire au calcul de distance // normale
1086 DO n=1,n4n
1087C
1088 i=i4n(n)
1089C
1090 IF(far(i)==1)cycle
1091C
1092C normale interpolee dans le triangle //
1093 nx(i)=(one-lb4(i)-lc4(i))*nnx0(i)+lb4(i)*nnx4(i)+lc4(i)*nnx1(i)
1094 ny(i)=(one-lb4(i)-lc4(i))*nny0(i)+lb4(i)*nny4(i)+lc4(i)*nny1(i)
1095 nz(i)=(one-lb4(i)-lc4(i))*nnz0(i)+lb4(i)*nnz4(i)+lc4(i)*nnz1(i)
1096C
1097C projection sur triangle 4
1098 nn = nx(i)*xn4(i)+ ny(i)*yn4(i)+ nz(i)*zn4(i)
1099 IF(nn <= zero) THEN
1100C cas limite
1101 far(i)=1
1102 cycle
1103 END IF
1104 nn = one/max(em30,nn)
1105C
1106 lambda=(xn4(i)*xi0v(i)+yn4(i)*yi0v(i)+zn4(i)*zi0v(i))*nn
1107 xp=xi(i)+lambda*nx(i)
1108 yp=yi(i)+lambda*ny(i)
1109 zp=zi(i)+lambda*nz(i)
1110C
1111 xi0v(i) = x0(i) - xp
1112 yi0v(i) = y0(i) - yp
1113 zi0v(i) = z0(i) - zp
1114C
1115 xi1(i) = x4(i) - xp
1116 yi1(i) = y4(i) - yp
1117 zi1(i) = z4(i) - zp
1118C
1119 xi2(i) = x1(i) - xp
1120 yi2(i) = y1(i) - yp
1121 zi2(i) = z1(i) - zp
1122C
1123 sx1(i) = yi0v(i)*zi1(i) - zi0v(i)*yi1(i)
1124 sy1(i) = zi0v(i)*xi1(i) - xi0v(i)*zi1(i)
1125 sz1(i) = xi0v(i)*yi1(i) - yi0v(i)*xi1(i)
1126C
1127 sx2(i) = yi0v(i)*zi2(i) - zi0v(i)*yi2(i)
1128 sy2(i) = zi0v(i)*xi2(i) - xi0v(i)*zi2(i)
1129 sz2(i) = xi0v(i)*yi2(i) - yi0v(i)*xi2(i)
1130C
1131 nn = one/
1132 . max(em30,(xn4(i)*xn4(i)+ yn4(i)*yn4(i)+ zn4(i)*zn4(i)))
1133 lb4(i) = -(xn4(i)*sx2(i) + yn4(i)*sy2(i) + zn4(i)*sz2(i))*nn
1134 lc4(i) = (xn4(i)*sx1(i) + yn4(i)*sy1(i) + zn4(i)*sz1(i))*nn
1135C
1136 IF(lb4(i) < -zep01 .OR. lc4(i) < -zep01 .OR.
1137 . lb4(i)+lc4(i) > onep01) far(i)=1
1138C
1139 ENDDO
1140
1141
1142 DO n=1,n4n
1143C
1144 i=i4n(n)
1145C
1146 IF(far(i)==1)cycle
1147C
1148 lb = min(one,max(lb4(i),zero))
1149 lc = min(one,max(lc4(i),zero))
1150 la = min(one,max(one-lb4(i)-lc4(i),zero))
1151
1152 sl=one/max(em20,la+lb+lc)
1153 lb = lb*sl
1154 lc = lc*sl
1155 la = la*sl
1156
1157 nx4(i) = xi(i)-(la*x0(i) + lb*x4(i) + lc*x1(i))
1158 ny4(i) = yi(i)-(la*y0(i) + lb*y4(i) + lc*y1(i))
1159 nz4(i) = zi(i)-(la*z0(i) + lb*z4(i) + lc*z1(i))
1160 p4(i) = nx4(i)*nx4(i) + ny4(i)*ny4(i) +nz4(i)*nz4(i)
1161C
1162C NX(I)=(1-LB4(I)-LC4(I))*NNX0(I)+LB4(I)*NNX4(I)+LC4(I)*NNX1(I)
1163C NY(I)=(1-LB4(I)-LC4(I))*NNY0(I)+LB4(I)*NNY4(I)+LC4(I)*NNY1(I)
1164C NZ(I)=(1-LB4(I)-LC4(I))*NNZ0(I)+LB4(I)*NNZ4(I)+LC4(I)*NNZ1(I)
1165C SIDE=NX4(I)*NX(I)+NY4(I)*NY(I)+NZ4(I)*NZ(I)
1166
1167 side=sign(one,nx4(i)*xn4(i)+ny4(i)*yn4(i)+nz4(i)*zn4(i))
1168 IF((side >= zero .AND. p4(i) < gap2(i)) .OR.
1169 . (side < zero .AND. p4(i) < depth2))THEN
1170 crit = abs(third-lb4(i))
1171 . + abs(third-lc4(i))
1172 . + abs(two_third-lb4(i)-lc4(i))
1173 lbo = csts_l(1,i)
1174 lco = csts_l(2,i)
1175 crito = abs(third-lbo)
1176 . + abs(third-lco)
1177 . + abs(two_third-lbo-lco)
1178 IF(crit < crito)THEN
1179 is=list(i)
1180 irtlm_l(1,i) = cand_e(is)
1181 irtlm_l(2,i) = nint(four*side)
1182 csts_l(1,i) = lb4(i)
1183 csts_l(2,i) = lc4(i)
1184 END IF
1185 ELSEIF(side >= zero .AND. p4(i) < max(drad2,gapp2(i)))THEN
1186 crit = abs(third-lb4(i))
1187 . + abs(third-lc4(i))
1188 . + abs(two_third-lb4(i)-lc4(i))
1189 lbo = csts_l(1,i)
1190 lco = csts_l(2,i)
1191 crito = abs(third-lbo)
1192 . + abs(third-lco)
1193 . + abs(two_third-lbo-lco)
1194 IF(crit < crito)THEN
1195 is=list(i)
1196 irtlm_l(1,i) = -cand_e(is)
1197 irtlm_l(2,i) = nint(four*side)
1198 csts_l(1,i) = lb4(i)
1199 csts_l(2,i) = lc4(i)
1200 END IF
1201 END IF
1202C
1203 ENDDO
1204C---------------------
1205C
1206C---------------------
1207 DO i=1,jlt
1208 IF(irtlm_l(1,i)/=0)THEN
1209 lb=csts_l(1,i)
1210 lc=csts_l(2,i)
1211 crit = abs(third-lb)
1212 . + abs(third-lc)
1213 . + abs(two_third-lb-lc)
1214 is=list(i)
1215#include "lockon.inc"
1216 lbo = csts(1,cand_n(is))
1217 lco = csts(2,cand_n(is))
1218 crito = abs(third-lbo)
1219 . + abs(third-lco)
1220 . + abs(two_third-lbo-lco)
1221 IF(crit < crito .OR.
1222 . crit==crito.AND.
1223 . abs(irtlm_l(1,i)) < abs(irtlm(1,cand_n(is))))THEN
1224 irtlm(1,cand_n(is)) = irtlm_l(1,i)
1225 irtlm(2,cand_n(is)) = irtlm_l(2,i)
1226 csts(1,cand_n(is)) = lb
1227 csts(2,cand_n(is)) = lc
1228 END IF
1229
1230#include "lockoff.inc"
1231 END IF
1232C
1233 ENDDO
1234C---------------------
1235 RETURN
1236 END
i21dst3
subroutine i21dst3(jlt, list, cand_n, cand_e, ix1, ix2, ix3, ix4, gapv, xi, yi, zi, irtlm, csts, depth2, nnx1, nny1, nnz1, nnx2, nny2, nnz2, nnx3, nny3, nnz3, nnx4, nny4, nnz4, x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, drad2, dgapload)
Definition i21dst3.F:38
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21