OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i23dst3.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!|| i23dst3 ../engine/source/interfaces/int23/i23dst3.F
25!||--- called by ------------------------------------------------------
26!|| i23mainf ../engine/source/interfaces/int23/i23mainf.F
27!||====================================================================
28 SUBROUTINE i23dst3(
29 1 JLT ,CAND_N_N ,CAND_E_N ,CN_LOC ,CE_LOC ,
30 2 X1 ,X2 ,X3 ,X4 ,Y1 ,
31 3 Y2 ,Y3 ,Y4 ,Z1 ,Z2 ,
32 4 Z3 ,Z4 ,XI ,YI ,ZI ,
33 6 IX1 ,IX2 ,IX3 ,IX4 ,NSVG ,
34 7 GAPV ,INACTI ,INDEX ,
35 8 VXM ,VYM ,VZM ,H1 ,H2 ,
36 9 H3 ,H4 ,IRECT ,CAND_P ,
37 A IFPEN ,NX ,NY ,NZ ,FTXSAV ,
38 B FTYSAV ,FTZSAV ,FXT ,FYT ,FZT ,
39 C PENE ,V ,VXI ,VYI ,VZI ,
40 D MSI ,STIF ,JLT_NEW,NSMS ,KINI )
41C-----------------------------------------------
42C I m p l i c i t T y p e s
43C-----------------------------------------------
44#include "implicit_f.inc"
45C-----------------------------------------------
46C G l o b a l P a r a m e t e r s
47C-----------------------------------------------
48#include "mvsiz_p.inc"
49#include "com08_c.inc"
50#include "sms_c.inc"
51C-----------------------------------------------
52C D u m m y A r g u m e n t s
53C-----------------------------------------------
54 INTEGER JLT, INACTI, CAND_N_N(*), CAND_E_N(*), JLT_NEW
55 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
56 . NSVG(MVSIZ), CN_LOC(MVSIZ) ,CE_LOC(MVSIZ) ,INDEX(MVSIZ),
57 . IRECT(4,*), IFPEN(*), NSMS(MVSIZ), KINI(MVSIZ)
58 my_real
59 . X1(MVSIZ), X2(MVSIZ), X3(MVSIZ), X4(MVSIZ),
60 . Y1(MVSIZ), Y2(MVSIZ), Y3(MVSIZ), Y4(MVSIZ),
61 . Z1(MVSIZ), Z2(MVSIZ), Z3(MVSIZ), Z4(MVSIZ),
62 . XI(MVSIZ), YI(MVSIZ), ZI(MVSIZ),
63 . GAPV(MVSIZ), PENE(MVSIZ), CAND_P(*),
64 . NX(MVSIZ), NY(MVSIZ), NZ(MVSIZ),
65 . XG(MVSIZ), YG(MVSIZ), ZG(MVSIZ),
66 . vxm(mvsiz), vym(mvsiz), vzm(mvsiz),
67 . ftxsav(*), ftysav(*), ftzsav(*),
68 . fxt(mvsiz), fyt(mvsiz), fzt(mvsiz),
69 . vxi(mvsiz), vyi(mvsiz), vzi(mvsiz), msi(mvsiz), stif(mvsiz),
70 . h1(mvsiz), h2(mvsiz), h3(mvsiz), h4(mvsiz), v(3,*)
71C-----------------------------------------------
72C L o c a l V a r i a b l e s
73C-----------------------------------------------
74 INTEGER I, IG, J, II, INPROJ(4,MVSIZ)
75 my_real
76 . LB1(MVSIZ), LC1(MVSIZ), P1(MVSIZ),
77 . LB2(MVSIZ), LC2(MVSIZ), P2(MVSIZ),
78 . LB3(MVSIZ), LC3(MVSIZ), P3(MVSIZ),
79 . LB4(MVSIZ), LC4(MVSIZ), P4(MVSIZ),
80 . LA1, LA2, LA3, LA4,
81 . NX1(MVSIZ), NY1(MVSIZ), NZ1(MVSIZ),
82 . NX2(MVSIZ), NY2(MVSIZ), NZ2(MVSIZ),
83 . NX3(MVSIZ), NY3(MVSIZ), NZ3(MVSIZ),
84 . NX4(MVSIZ), NY4(MVSIZ), NZ4(MVSIZ),
85 . AL1(MVSIZ), AL2(MVSIZ), AL3(MVSIZ), AL4(MVSIZ),
86 . x01(mvsiz), x02(mvsiz), x03(mvsiz), x04(mvsiz),
87 . y01(mvsiz), y02(mvsiz), y03(mvsiz), y04(mvsiz),
88 . z01(mvsiz), z02(mvsiz), z03(mvsiz), z04(mvsiz),
89 . xi1(mvsiz), xi2(mvsiz), xi3(mvsiz), xi4(mvsiz),
90 . yi1(mvsiz), yi2(mvsiz), yi3(mvsiz), yi4(mvsiz),
91 . zi1(mvsiz), zi2(mvsiz), zi3(mvsiz), zi4(mvsiz),
92 . x0(mvsiz), y0(mvsiz), z0(mvsiz),
93 . xi0v(mvsiz), yi0v(mvsiz), zi0v(mvsiz),
94 . la, hla, h0, h00,
95 . hlb1(mvsiz), hlc1(mvsiz),
96 . hlb2(mvsiz), hlc2(mvsiz),
97 . hlb3(mvsiz), hlc3(mvsiz),
98 . hlb4(mvsiz), hlc4(mvsiz),
99 . xp(mvsiz), yp(mvsiz), zp(mvsiz),
100 . xp1(mvsiz), yp1(mvsiz), zp1(mvsiz),
101 . xp2(mvsiz), yp2(mvsiz), zp2(mvsiz),
102 . xp3(mvsiz), yp3(mvsiz), zp3(mvsiz),
103 . xp4(mvsiz), yp4(mvsiz), zp4(mvsiz),
104 . n11(mvsiz), n21(mvsiz), n31(mvsiz),
105 . n12(mvsiz), n22(mvsiz), n32(mvsiz),
106 . n13(mvsiz), n23(mvsiz), n33(mvsiz),
107 . n14(mvsiz), n24(mvsiz), n34(mvsiz),
108 . nvers(4,mvsiz), firstimp(4,mvsiz),
109 . xn1(mvsiz),yn1(mvsiz),zn1(mvsiz),
110 . xn2(mvsiz),yn2(mvsiz),zn2(mvsiz),
111 . xn3(mvsiz),yn3(mvsiz),zn3(mvsiz),
112 . xn4(mvsiz),yn4(mvsiz),zn4(mvsiz),
113 . ssc(mvsiz),ttc(mvsiz),
114 . a12(mvsiz),a23(mvsiz),a34(mvsiz),a41(mvsiz),
115 . b12(mvsiz),b23(mvsiz),b34(mvsiz),b41(mvsiz),
116 . ab1(mvsiz),ab2(mvsiz),
117 . n1(mvsiz),n2(mvsiz),n3(mvsiz),
118 . xx1(mvsiz),xx2(mvsiz),xx3(mvsiz),xx4(mvsiz),
119 . yy1(mvsiz),yy2(mvsiz),yy3(mvsiz),yy4(mvsiz),
120 . zz1(mvsiz),zz2(mvsiz),zz3(mvsiz),zz4(mvsiz),
121 . area(mvsiz),alp(mvsiz),var
122 my_real
123 . s2,nn,ps,
124 . x12,x23,x34,x41,sx1,sx2,sx3,sx4,sx0,
125 . y12,y23,y34,y41,sy1,sy2,sy3,sy4,sy0,
126 . z12,z23,z34,z41,sz1,sz2,sz3,sz4,sz0,
127 . gap2(mvsiz), aaa, side, dist,
128 . ll, h
129C--------------------------------------------------------
130C CAS DES PAQUETS MIXTES
131C--------------------------------------------------------
132 DO i=1,jlt
133 IF(ix3(i)/=ix4(i))THEN
134 x0(i) = fourth*(x1(i)+x2(i)+x3(i)+x4(i))
135 y0(i) = fourth*(y1(i)+y2(i)+y3(i)+y4(i))
136 z0(i) = fourth*(z1(i)+z2(i)+z3(i)+z4(i))
137 ELSE
138 x0(i) = x3(i)
139 y0(i) = y3(i)
140 z0(i) = z3(i)
141 ENDIF
142 ENDDO
143C
144 DO i=1,jlt
145C
146 x01(i) = x1(i) - x0(i)
147 y01(i) = y1(i) - y0(i)
148 z01(i) = z1(i) - z0(i)
149C
150 x02(i) = x2(i) - x0(i)
151 y02(i) = y2(i) - y0(i)
152 z02(i) = z2(i) - z0(i)
153C
154 x03(i) = x3(i) - x0(i)
155 y03(i) = y3(i) - y0(i)
156 z03(i) = z3(i) - z0(i)
157C
158 x04(i) = x4(i) - x0(i)
159 y04(i) = y4(i) - y0(i)
160 z04(i) = z4(i) - z0(i)
161C
162 xi0v(i) = x0(i) - xi(i)
163 yi0v(i) = y0(i) - yi(i)
164 zi0v(i) = z0(i) - zi(i)
165C
166 xi1(i) = x1(i) - xi(i)
167 yi1(i) = y1(i) - yi(i)
168 zi1(i) = z1(i) - zi(i)
169C
170 xi2(i) = x2(i) - xi(i)
171 yi2(i) = y2(i) - yi(i)
172 zi2(i) = z2(i) - zi(i)
173C
174 xi3(i) = x3(i) - xi(i)
175 yi3(i) = y3(i) - yi(i)
176 zi3(i) = z3(i) - zi(i)
177C
178 xi4(i) = x4(i) - xi(i)
179 yi4(i) = y4(i) - yi(i)
180 zi4(i) = z4(i) - zi(i)
181C
182 sx1 = yi0v(i)*zi1(i) - zi0v(i)*yi1(i)
183 sy1 = zi0v(i)*xi1(i) - xi0v(i)*zi1(i)
184 sz1 = xi0v(i)*yi1(i) - yi0v(i)*xi1(i)
185C
186 sx2 = yi0v(i)*zi2(i) - zi0v(i)*yi2(i)
187 sy2 = zi0v(i)*xi2(i) - xi0v(i)*zi2(i)
188 sz2 = xi0v(i)*yi2(i) - yi0v(i)*xi2(i)
189C
190 sx0 = y01(i)*z02(i) - z01(i)*y02(i)
191 sy0 = z01(i)*x02(i) - x01(i)*z02(i)
192 sz0 = x01(i)*y02(i) - y01(i)*x02(i)
193 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
194C
195 nn=sqrt(s2)
196 n11(i)=sx0*nn
197 n21(i)=sy0*nn
198 n31(i)=sz0*nn
199C
200 area(i)=nn
201C
202 lb1(i) = -(sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
203 lc1(i) = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
204C
205 sx3 = yi0v(i)*zi3(i) - zi0v(i)*yi3(i)
206 sy3 = zi0v(i)*xi3(i) - xi0v(i)*zi3(i)
207 sz3 = xi0v(i)*yi3(i) - yi0v(i)*xi3(i)
208C
209 sx0 = y02(i)*z03(i) - z02(i)*y03(i)
210 sy0 = z02(i)*x03(i) - x02(i)*z03(i)
211 sz0 = x02(i)*y03(i) - y02(i)*x03(i)
212 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
213C
214 nn=sqrt(s2)
215 n12(i)=sx0*nn
216 n22(i)=sy0*nn
217 n32(i)=sz0*nn
218C
219 area(i)=area(i)+nn
220C
221 lb2(i) = -(sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
222 lc2(i) = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
223C
224 sx4 = yi0v(i)*zi4(i) - zi0v(i)*yi4(i)
225 sy4 = zi0v(i)*xi4(i) - xi0v(i)*zi4(i)
226 sz4 = xi0v(i)*yi4(i) - yi0v(i)*xi4(i)
227C
228 sx0 = y03(i)*z04(i) - z03(i)*y04(i)
229 sy0 = z03(i)*x04(i) - x03(i)*z04(i)
230 sz0 = x03(i)*y04(i) - y03(i)*x04(i)
231 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
232C
233 nn=sqrt(s2)
234 n13(i)=sx0*nn
235 n23(i)=sy0*nn
236 n33(i)=sz0*nn
237C
238 area(i)=area(i)+nn
239C
240 lb3(i) = -(sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
241 lc3(i) = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
242C
243 sx0 = y04(i)*z01(i) - z04(i)*y01(i)
244 sy0 = z04(i)*x01(i) - x04(i)*z01(i)
245 sz0 = x04(i)*y01(i) - y04(i)*x01(i)
246 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
247C
248 nn=sqrt(s2)
249 n14(i)=sx0*nn
250 n24(i)=sy0*nn
251 n34(i)=sz0*nn
252C
253 area(i)=area(i)+nn
254 area(i)=half*area(i)
255C
256 lb4(i) = -(sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
257 lc4(i) = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
258 ENDDO
259C---------------------------------------------------------
260C
261 DO i=1,jlt
262 aaa = one/max(em30,x01(i)*x01(i)+y01(i)*y01(i)+z01(i)*z01(i))
263 hlc1(i)= lc1(i)*abs(lc1(i))*aaa
264 hlb4(i)= lb4(i)*abs(lb4(i))*aaa
265 al1(i) = -(xi0v(i)*x01(i)+yi0v(i)*y01(i)+zi0v(i)*z01(i))*aaa
266 al1(i) = max(zero,min(one,al1(i)))
267 aaa = one/max(em30,x02(i)*x02(i)+y02(i)*y02(i)+z02(i)*z02(i))
268 hlc2(i)= lc2(i)*abs(lc2(i))*aaa
269 hlb1(i)= lb1(i)*abs(lb1(i))*aaa
270 al2(i) = -(xi0v(i)*x02(i)+yi0v(i)*y02(i)+zi0v(i)*z02(i))*aaa
271 al2(i) = max(zero,min(one,al2(i)))
272 aaa = one/max(em30,x03(i)*x03(i)+y03(i)*y03(i)+z03(i)*z03(i))
273 hlc3(i)= lc3(i)*abs(lc3(i))*aaa
274 hlb2(i)= lb2(i)*abs(lb2(i))*aaa
275 al3(i) = -(xi0v(i)*x03(i)+yi0v(i)*y03(i)+zi0v(i)*z03(i))*aaa
276 al3(i) = max(zero,min(one,al3(i)))
277 aaa = one/max(em30,x04(i)*x04(i)+y04(i)*y04(i)+z04(i)*z04(i))
278 hlc4(i)= lc4(i)*abs(lc4(i))*aaa
279 hlb3(i)= lb3(i)*abs(lb3(i))*aaa
280 al4(i) = -(xi0v(i)*x04(i)+yi0v(i)*y04(i)+zi0v(i)*z04(i))*aaa
281 al4(i) = max(zero,min(one,al4(i)))
282 ENDDO
283C
284 DO i=1,jlt
285 x12 = x2(i) - x1(i)
286 y12 = y2(i) - y1(i)
287 z12 = z2(i) - z1(i)
288 la = one - lb1(i) - lc1(i)
289C HLA, HLB1, HLC1 necessaires pour triangle angle obtu
290 aaa = one / max(em20,x12*x12+y12*y12+z12*z12)
291 hla= la*abs(la) * aaa
292 inproj(1,i)=0
293 IF(la<zero.AND.
294 + hla<=hlb1(i).AND.hla<=hlc1(i))THEN
295 lb1(i) = (xi2(i)*x12+yi2(i)*y12+zi2(i)*z12) * aaa
296 lb1(i) = max(zero,min(one,lb1(i)))
297 lc1(i) = one - lb1(i)
298 inproj(1,i)=1
299 ELSEIF(lb1(i)<zero.AND.
300 + hlb1(i)<=hlc1(i).AND.hlb1(i)<=hla)THEN
301 lb1(i) = zero
302 lc1(i) = al2(i)
303 inproj(1,i)=1
304 ELSEIF(lc1(i)<zero.AND.
305 + hlc1(i)<=hla.AND.hlc1(i)<=hlb1(i))THEN
306 lc1(i) = zero
307 lb1(i) = al1(i)
308 inproj(1,i)=1
309 ENDIF
310 ENDDO
311C
312 DO i=1,jlt
313 IF(ix3(i)==ix4(i))cycle
314 x23 = x3(i) - x2(i)
315 y23 = y3(i) - y2(i)
316 z23 = z3(i) - z2(i)
317 la = one - lb2(i) - lc2(i)
318C HLA, HLB1, HLC1 necessaires pour triangle angle obtu
319 aaa = one / max(em20,x23*x23+y23*y23+z23*z23)
320 hla= la*abs(la) * aaa
321 inproj(2,i)=0
322 IF(la<zero.AND.
323 + hla<=hlb2(i).AND.hla<=hlc2(i))THEN
324 lb2(i) = (xi3(i)*x23+yi3(i)*y23+zi3(i)*z23)*aaa
325 lb2(i) = max(zero,min(one,lb2(i)))
326 lc2(i) = one - lb2(i)
327 inproj(2,i)=1
328 ELSEIF(lb2(i)<zero.AND.
329 + hlb2(i)<=hlc2(i).AND.hlb2(i)<=hla)THEN
330 lb2(i) = zero
331 lc2(i) = al3(i)
332 inproj(2,i)=1
333 ELSEIF(lc2(i)<zero.AND.
334 + hlc2(i)<=hla.AND.hlc2(i)<=hlb2(i))THEN
335 lc2(i) = zero
336 lb2(i) = al2(i)
337 inproj(2,i)=1
338 ENDIF
339 ENDDO
340C
341 DO i=1,jlt
342 IF(ix3(i)==ix4(i))cycle
343 x34 = x4(i) - x3(i)
344 y34 = y4(i) - y3(i)
345 z34 = z4(i) - z3(i)
346 la = one - lb3(i) - lc3(i)
347C HLA, HLB1, HLC1 necessaires pour triangle angle obtu
348 aaa = one / max(em20,x34*x34+y34*y34+z34*z34)
349 hla= la*abs(la) * aaa
350 inproj(3,i)=0
351 IF(la<zero.AND.
352 + hla<=hlb3(i).AND.hla<=hlc3(i))THEN
353 lb3(i) = (xi4(i)*x34+yi4(i)*y34+zi4(i)*z34)*aaa
354 lb3(i) = max(zero,min(one,lb3(i)))
355 lc3(i) = one - lb3(i)
356 inproj(3,i)=1
357 ELSEIF(lb3(i)<zero.AND.
358 + hlb3(i)<=hlc3(i).AND.hlb3(i)<=hla)THEN
359 lb3(i) = zero
360 lc3(i) = al4(i)
361 inproj(3,i)=1
362 ELSEIF(lc3(i)<zero.AND.
363 + hlc3(i)<=hla.AND.hlc3(i)<=hlb3(i))THEN
364 lc3(i) = zero
365 lb3(i) = al3(i)
366 inproj(3,i)=1
367 ENDIF
368 ENDDO
369C
370 DO i=1,jlt
371 IF(ix3(i)==ix4(i))cycle
372 x41 = x1(i) - x4(i)
373 y41 = y1(i) - y4(i)
374 z41 = z1(i) - z4(i)
375 la = one - lb4(i) - lc4(i)
376C HLA, HLB1, HLC1 necessaires pour triangle angle obtu
377 aaa = one / max(em20,x41*x41+y41*y41+z41*z41)
378 hla= la*abs(la) * aaa
379 inproj(4,i)=0
380 IF(la<zero.AND.
381 + hla<=hlb4(i).AND.hla<=hlc4(i))THEN
382 lb4(i) = (xi1(i)*x41+yi1(i)*y41+zi1(i)*z41)*aaa
383 lb4(i) = max(zero,min(one,lb4(i)))
384 lc4(i) = one - lb4(i)
385 inproj(4,i)=1
386 ELSEIF(lb4(i)<zero.AND.
387 + hlb4(i)<=hlc4(i).AND.hlb4(i)<=hla)THEN
388 lb4(i) = zero
389 lc4(i) = al1(i)
390 inproj(4,i)=1
391 ELSEIF(lc4(i)<zero.AND.
392 + hlc4(i)<=hla.AND.hlc4(i)<=hlb4(i))THEN
393 lc4(i) = zero
394 lb4(i) = al4(i)
395 inproj(4,i)=1
396 ENDIF
397 ENDDO
398C---------------------------------------------------------
399 DO i=1,jlt
400C
401 gap2(i)=gapv(i)*gapv(i)
402C
403 la1 = one - lb1(i) - lc1(i)
404 xp1(i) = lb1(i)*x1(i) + lc1(i)*x2(i) + la1*x0(i)
405 yp1(i) = lb1(i)*y1(i) + lc1(i)*y2(i) + la1*y0(i)
406 zp1(i) = lb1(i)*z1(i) + lc1(i)*z2(i) + la1*z0(i)
407C
408 nx1(i) = xi(i)-xp1(i)
409 ny1(i) = yi(i)-yp1(i)
410 nz1(i) = zi(i)-zp1(i)
411 p1(i) = nx1(i)*nx1(i) + ny1(i)*ny1(i) +nz1(i)*nz1(i)
412C
413 la2 = one - lb2(i) - lc2(i)
414 xp2(i) = lb2(i)*x2(i) + lc2(i)*x3(i) + la2*x0(i)
415 yp2(i) = lb2(i)*y2(i) + lc2(i)*y3(i) + la2*y0(i)
416 zp2(i) = lb2(i)*z2(i) + lc2(i)*z3(i) + la2*z0(i)
417C
418 nx2(i) = xi(i)-xp2(i)
419 ny2(i) = yi(i)-yp2(i)
420 nz2(i) = zi(i)-zp2(i)
421 p2(i) = nx2(i)*nx2(i) + ny2(i)*ny2(i) +nz2(i)*nz2(i)
422C
423 la3 = one - lb3(i) - lc3(i)
424 xp3(i) = lb3(i)*x3(i) + lc3(i)*x4(i) + la3*x0(i)
425 yp3(i) = lb3(i)*y3(i) + lc3(i)*y4(i) + la3*y0(i)
426 zp3(i) = lb3(i)*z3(i) + lc3(i)*z4(i) + la3*z0(i)
427C
428 nx3(i) = xi(i)-xp3(i)
429 ny3(i) = yi(i)-yp3(i)
430 nz3(i) = zi(i)-zp3(i)
431 p3(i) = nx3(i)*nx3(i) + ny3(i)*ny3(i) +nz3(i)*nz3(i)
432C
433 la4 = one - lb4(i) - lc4(i)
434 xp4(i) = lb4(i)*x4(i) + lc4(i)*x1(i) + la4*x0(i)
435 yp4(i) = lb4(i)*y4(i) + lc4(i)*y1(i) + la4*y0(i)
436 zp4(i) = lb4(i)*z4(i) + lc4(i)*z1(i) + la4*z0(i)
437C
438 nx4(i) = xi(i)-xp4(i)
439 ny4(i) = yi(i)-yp4(i)
440 nz4(i) = zi(i)-zp4(i)
441 p4(i) = nx4(i)*nx4(i) + ny4(i)*ny4(i) +nz4(i)*nz4(i)
442C
443 ENDDO
444C
445 DO i=1,jlt
446 END DO
447C
448 DO i=1,jlt
449C
450 IF(ix3(i)/=ix4(i))THEN
451C
452 h =nx1(i)*n11(i) + ny1(i)*n21(i) +nz1(i)*n31(i)
453 ll=p1(i)-h*h
454 IF(inproj(1,i)/=0.AND.ll >= gap2(i))THEN
455 p1(i)=zero
456 ELSE
457 p1(i)=max(zero,gapv(i)-abs(h))
458 END IF
459C
460 h =nx2(i)*n12(i) + ny2(i)*n22(i) +nz2(i)*n32(i)
461 ll=p2(i)-h*h
462 IF(inproj(2,i)/=0.AND.ll >= gap2(i))THEN
463 p2(i)=zero
464 ELSE
465 p2(i)=max(zero,gapv(i)-abs(h))
466 END IF
467C
468 h =nx3(i)*n13(i) + ny3(i)*n23(i) +nz3(i)*n33(i)
469 ll=p3(i)-h*h
470 IF(inproj(3,i)/=0.AND.ll >= gap2(i))THEN
471 p3(i)=zero
472 ELSE
473 p3(i)=max(zero,gapv(i)-abs(h))
474 END IF
475C
476 h =nx4(i)*n14(i) + ny4(i)*n24(i) +nz4(i)*n34(i)
477 ll=p4(i)-h*h
478 IF(inproj(4,i)/=0.AND.ll >= gap2(i))THEN
479 p4(i)=zero
480 ELSE
481 p4(i)=max(zero,gapv(i)-abs(h))
482 END IF
483C
484 pene(i) = max(p1(i),p2(i),p3(i),p4(i))
485C
486 n1(i) = p1(i)*n11(i)+p2(i)*n12(i)+p3(i)*n13(i)+p4(i)*n14(i)
487 n2(i) = p1(i)*n21(i)+p2(i)*n22(i)+p3(i)*n23(i)+p4(i)*n24(i)
488 n3(i) = p1(i)*n31(i)+p2(i)*n32(i)+p3(i)*n33(i)+p4(i)*n34(i)
489 nn=sqrt(n1(i)*n1(i)+n2(i)*n2(i)+n3(i)*n3(i))
490 nn=one/max(em20,nn)
491 n1(i)=n1(i)*nn
492 n2(i)=n2(i)*nn
493 n3(i)=n3(i)*nn
494C
495 la1 = one - lb1(i) - lc1(i)
496 la2 = one - lb2(i) - lc2(i)
497 la3 = one - lb3(i) - lc3(i)
498 la4 = one - lb4(i) - lc4(i)
499C
500 h0 = fourth *
501 . (p1(i)*la1 + p2(i)*la2 + p3(i)*la3 + p4(i)*la4)
502 h1(i) = h0 + p1(i) * lb1(i) + p4(i) * lc4(i)
503 h2(i) = h0 + p2(i) * lb2(i) + p1(i) * lc1(i)
504 h3(i) = h0 + p3(i) * lb3(i) + p2(i) * lc2(i)
505 h4(i) = h0 + p4(i) * lb4(i) + p3(i) * lc3(i)
506 h00 = one/max(em20,h1(i) + h2(i) + h3(i) + h4(i))
507 h1(i) = h1(i) * h00
508 h2(i) = h2(i) * h00
509 h3(i) = h3(i) * h00
510 h4(i) = h4(i) * h00
511C
512 ELSE ! IF(IX3(I)/=IX4(I))THEN
513C
514 h1(i) = lb1(i)
515 h2(i) = lc1(i)
516 h3(i) = one - lb1(i) - lc1(i)
517 h4(i) = zero
518C
519 n1(i) = n11(i)
520 n2(i) = n21(i)
521 n3(i) = n31(i)
522C
523 h =nx1(i)*n11(i) + ny1(i)*n21(i) +nz1(i)*n31(i)
524 ll=p1(i)-h*h
525 IF(inproj(1,i)/=0.AND.ll >= gap2(i))THEN
526 pene(i)=zero
527 ELSE
528 pene(i)=max(zero,gapv(i)-abs(h))
529 END IF
530 END IF
531 ENDDO
532C
533 DO 410 i=1,jlt
534 xp(i)=h1(i)*x1(i)+h2(i)*x2(i)+h3(i)*x3(i)+h4(i)*x4(i)
535 yp(i)=h1(i)*y1(i)+h2(i)*y2(i)+h3(i)*y3(i)+h4(i)*y4(i)
536 zp(i)=h1(i)*z1(i)+h2(i)*z2(i)+h3(i)*z3(i)+h4(i)*z4(i)
537 410 CONTINUE
538C---------------------------------------------------------
539 DO i=1,jlt
540 nx1(i) = xi(i)-xp(i)
541 ny1(i) = yi(i)-yp(i)
542 nz1(i) = zi(i)-zp(i)
543 ENDDO
544C---------------------
545 DO i=1,jlt
546C
547 nvers(1,i) = zero
548 firstimp(1,i) = zero
549C
550 nx(i)=zero
551 ny(i)=zero
552 nz(i)=zero
553 IF(pene(i)==zero)cycle
554C--------------------------------------------------------
555 side=n1(i)*nx1(i)+n2(i)*ny1(i)+n3(i)*nz1(i)
556C
557 IF(ifpen(index(i))==0.OR.tt==zero)THEN
558 firstimp(1,i) = sign(one,side)
559 IF(firstimp(1,i) < zero)THEN
560 n1(i) = -n1(i)
561 n2(i) = -n2(i)
562 n3(i) = -n3(i)
563 END IF
564 nvers(1,i) = one
565C 1st impact below gap (sorting security)
566 ELSE ! IF(IFPEN(INDEX(I))==0.OR.TT==ZERO)THEN
567 IF(ifpen(index(i)) < 0)THEN
568 n1(i) = -n1(i)
569 n2(i) = -n2(i)
570 n3(i) = -n3(i)
571C SIDE = -SIDE
572 END IF
573 nvers(1,i)= sign(one,side*ifpen(index(i)))
574 END IF
575C---------------------
576C
577C attention a la traversee de la coque => 1E20...,
578C prendre normale = normale au triangle
579 nn=one/
580 . max(em20,sqrt(n1(i)*n1(i)+n2(i)*n2(i)+n3(i)*n3(i)))
581 n1(i)=n1(i)*nn
582 n2(i)=n2(i)*nn
583 n3(i)=n3(i)*nn
584 IF(nvers(1,i)<zero)
585 . pene(i)=(gapv(i)-pene(i))+gapv(i)
586C
587 nx(i)=n1(i)
588 ny(i)=n2(i)
589 nz(i)=n3(i)
590C
591 END DO
592C---------------------
593C PENE INITIALE
594C---------------------
595 IF(tt/=zero)THEN
596 DO i=1,jlt
597 j=index(i)
598 IF(pene(i)==zero)THEN
599C
600 ftxsav(j)=zero
601 ftysav(j)=zero
602 ftzsav(j)=zero
603 cand_p(j)=zero
604 ifpen(j) =0
605 ELSE
606C
607 IF(ifpen(j)==0)THEN
608 ftxsav(j)=zero
609 ftysav(j)=zero
610 ftzsav(j)=zero
611 END IF
612C
613 IF(ifpen(j)==0.AND.pene(i)/=zero)THEN
614 cand_p(j) =pene(i)
615 END IF
616C
617 IF(ifpen(j)==0)THEN
618 IF(firstimp(1,i) > zero)THEN
619 ifpen(j) = 1
620 ELSEIF(firstimp(1,i) < zero)THEN
621 ifpen(j) =-1
622 END IF
623 ELSEIF(pene(i)/=zero)THEN
624 ifpen(j)=sign(abs(ifpen(j))+1,ifpen(j))
625 ELSE
626 ifpen(j)=0
627 END IF
628 END IF
629C
630 END DO
631 END IF
632C---------------------------------
633 DO i=1,jlt
634 vxm(i)=zero
635 vym(i)=zero
636 vzm(i)=zero
637C utile que si IFPEN/=0
638 IF(pene(i)/=0)THEN
639 vxm(i)=h1(i)*v(1,ix1(i))+h2(i)*v(1,ix2(i))+
640 . h3(i)*v(1,ix3(i))+h4(i)*v(1,ix4(i))
641 vym(i)=h1(i)*v(2,ix1(i))+h2(i)*v(2,ix2(i))+
642 . h3(i)*v(2,ix3(i))+h4(i)*v(2,ix4(i))
643 vzm(i)=h1(i)*v(3,ix1(i))+h2(i)*v(3,ix2(i))+
644 . h3(i)*v(3,ix3(i))+h4(i)*v(3,ix4(i))
645 END IF
646 END DO
647C
648 DO i=1,jlt
649 j=index(i)
650 fxt(i)=ftxsav(j)
651 fyt(i)=ftysav(j)
652 fzt(i)=ftzsav(j)
653 END DO
654C---------------------
655C JLT_NEW=JLT
656C RETURN
657C
658 DO i=1,jlt
659 IF(pene(i)/=zero.AND.stif(i)/=zero)THEN
660 jlt_new = jlt_new + 1
661 cn_loc(jlt_new) = cand_n_n(i)
662 ce_loc(jlt_new) = cand_e_n(i)
663 ix1(jlt_new) = ix1(i)
664 ix2(jlt_new) = ix2(i)
665 ix3(jlt_new) = ix3(i)
666 ix4(jlt_new) = ix4(i)
667 nsvg(jlt_new) = nsvg(i)
668 nx(jlt_new) = nx(i)
669 ny(jlt_new) = ny(i)
670 nz(jlt_new) = nz(i)
671 pene(jlt_new) = pene(i)
672 h1(jlt_new) = h1(i)
673 h2(jlt_new) = h2(i)
674 h3(jlt_new) = h3(i)
675 h4(jlt_new) = h4(i)
676 stif(jlt_new) = stif(i)
677 gapv(jlt_new) = gapv(i)
678 index(jlt_new)= index(i)
679C
680C KINI(JLT_NEW) = KINI(I)
681 vxi(jlt_new) = vxi(i)
682 vyi(jlt_new) = vyi(i)
683 vzi(jlt_new) = vzi(i)
684 msi(jlt_new) = msi(i)
685C
686 vxm(jlt_new) = vxm(i)
687 vym(jlt_new) = vym(i)
688 vzm(jlt_new) = vzm(i)
689C
690 fxt(jlt_new) = fxt(i)
691 fyt(jlt_new) = fyt(i)
692 fzt(jlt_new) = fzt(i)
693C
694 xi(jlt_new) = xi(i)
695 yi(jlt_new) = yi(i)
696 zi(jlt_new) = zi(i)
697C
698 kini(jlt_new) = kini(i)
699C
700 IF(idtmins==2.OR.idtmins_int/=0)nsms(jlt_new)=nsms(i)
701C
702 ENDIF
703 ENDDO
704C---------------------
705 RETURN
706 END
707C===============================================================================
i23dst3
subroutine i23dst3(jlt, cand_n_n, cand_e_n, cn_loc, ce_loc, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, xi, yi, zi, ix1, ix2, ix3, ix4, nsvg, gapv, inacti, index, vxm, vym, vzm, h1, h2, h3, h4, irect, cand_p, ifpen, nx, ny, nz, ftxsav, ftysav, ftzsav, fxt, fyt, fzt, pene, v, vxi, vyi, vzi, msi, stif, jlt_new, nsms, kini)
Definition i23dst3.F:41
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21