OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i23dst3.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine i23dst3 (jlt, cand_n, cand_e, irect, nsv, gap_s, x, msr, pene, ifpen, igap, gap, gapmax, gapmin, gapv, gap_m)

Function/Subroutine Documentation

◆ i23dst3()

subroutine i23dst3 ( integer jlt,
integer, dimension(*) cand_n,
integer, dimension(*) cand_e,
integer, dimension(4,*) irect,
integer, dimension(*) nsv,
gap_s,
x,
integer, dimension(*) msr,
pene,
integer, dimension(*) ifpen,
integer igap,
gap,
gapmax,
gapmin,
gapv,
gap_m )

Definition at line 28 of file i23dst3.F.

33C-----------------------------------------------
34C I m p l i c i t T y p e s
35C-----------------------------------------------
36#include "implicit_f.inc"
37C-----------------------------------------------
38C G l o b a l P a r a m e t e r s
39C-----------------------------------------------
40#include "mvsiz_p.inc"
41C-----------------------------------------------
42C D u m m y A r g u m e n t s
43C-----------------------------------------------
44 INTEGER JLT, CAND_N(*), CAND_E(*), IRECT(4,*),
45 . NSV(*), IFPEN(*), IGAP, MSR(*)
46 my_real
47 . x(3,*), pene(*), gap_s(*),
48 . gap, gapmin, gapmax, gapv(mvsiz), gap_m(*)
49C-----------------------------------------------
50C L o c a l V a r i a b l e s
51C-----------------------------------------------
52 INTEGER I, IG, J, L, INPROJ(4,MVSIZ)
53 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
54 . NSVG(MVSIZ)
55 my_real
56 . x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz),
57 . y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz),
58 . z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz),
59 . xi(mvsiz), yi(mvsiz), zi(mvsiz),
60 . lb1(mvsiz), lc1(mvsiz), p1(mvsiz),
61 . lb2(mvsiz), lc2(mvsiz), p2(mvsiz),
62 . lb3(mvsiz), lc3(mvsiz), p3(mvsiz),
63 . lb4(mvsiz), lc4(mvsiz), p4(mvsiz),
64 . la1, la2, la3, la4,
65 . nx1(mvsiz), ny1(mvsiz), nz1(mvsiz),
66 . nx2(mvsiz), ny2(mvsiz), nz2(mvsiz),
67 . nx3(mvsiz), ny3(mvsiz), nz3(mvsiz),
68 . nx4(mvsiz), ny4(mvsiz), nz4(mvsiz),
69 . al1(mvsiz), al2(mvsiz), al3(mvsiz), al4(mvsiz),
70 . x01(mvsiz), x02(mvsiz), x03(mvsiz), x04(mvsiz),
71 . y01(mvsiz), y02(mvsiz), y03(mvsiz), y04(mvsiz),
72 . z01(mvsiz), z02(mvsiz), z03(mvsiz), z04(mvsiz),
73 . xi1(mvsiz), xi2(mvsiz), xi3(mvsiz), xi4(mvsiz),
74 . yi1(mvsiz), yi2(mvsiz), yi3(mvsiz), yi4(mvsiz),
75 . zi1(mvsiz), zi2(mvsiz), zi3(mvsiz), zi4(mvsiz),
76 . x0(mvsiz), y0(mvsiz), z0(mvsiz),
77 . xi0v(mvsiz), yi0v(mvsiz), zi0v(mvsiz),
78 . la, hla, h0, h00,
79 . hlb1(mvsiz), hlc1(mvsiz),
80 . hlb2(mvsiz), hlc2(mvsiz),
81 . hlb3(mvsiz), hlc3(mvsiz),
82 . hlb4(mvsiz), hlc4(mvsiz),
83 . xp(mvsiz), yp(mvsiz), zp(mvsiz),
84 . xp1(mvsiz), yp1(mvsiz), zp1(mvsiz),
85 . xp2(mvsiz), yp2(mvsiz), zp2(mvsiz),
86 . xp3(mvsiz), yp3(mvsiz), zp3(mvsiz),
87 . xp4(mvsiz), yp4(mvsiz), zp4(mvsiz),
88 . n11(mvsiz), n21(mvsiz), n31(mvsiz),
89 . n12(mvsiz), n22(mvsiz), n32(mvsiz),
90 . n13(mvsiz), n23(mvsiz), n33(mvsiz),
91 . n14(mvsiz), n24(mvsiz), n34(mvsiz),
92 . nx(mvsiz), ny(mvsiz), nz(mvsiz),
93 . xg(mvsiz), yg(mvsiz), zg(mvsiz),
94 . h1(mvsiz), h2(mvsiz), h3(mvsiz), h4(mvsiz),
95 . nvers(4,mvsiz), firstimp(4,mvsiz),
96 . xn1(mvsiz),yn1(mvsiz),zn1(mvsiz),
97 . xn2(mvsiz),yn2(mvsiz),zn2(mvsiz),
98 . xn3(mvsiz),yn3(mvsiz),zn3(mvsiz),
99 . xn4(mvsiz),yn4(mvsiz),zn4(mvsiz),
100 . ssc(mvsiz),ttc(mvsiz),
101 . a12(mvsiz),a23(mvsiz),a34(mvsiz),a41(mvsiz),
102 . b12(mvsiz),b23(mvsiz),b34(mvsiz),b41(mvsiz),
103 . ab1(mvsiz),ab2(mvsiz),
104 . n1(mvsiz),n2(mvsiz),n3(mvsiz),
105 . xx1(mvsiz),xx2(mvsiz),xx3(mvsiz),xx4(mvsiz),
106 . yy1(mvsiz),yy2(mvsiz),yy3(mvsiz),yy4(mvsiz),
107 . zz1(mvsiz),zz2(mvsiz),zz3(mvsiz),zz4(mvsiz),
108 . area(mvsiz),alp(mvsiz),var
109 my_real
110 . s2,nn,ps,
111 . x12,x23,x34,x41,sx1,sx2,sx3,sx4,sx0,
112 . y12,y23,y34,y41,sy1,sy2,sy3,sy4,sy0,
113 . z12,z23,z34,z41,sz1,sz2,sz3,sz4,sz0,
114 . gap2(mvsiz), aaa, side, dist,
115 . ll, h
116C--------------------------------------------------------
117 IF(igap==0)THEN
118 DO i=1,jlt
119 gapv(i)=gap
120 END DO
121 ELSE
122 DO i=1,jlt
123 gapv(i)=gap_s(cand_n(i))+gap_m(cand_e(i))
124 gapv(i)=min(gapv(i),gapmax)
125 gapv(i)=max(gapmin,gapv(i))
126 END DO
127 END IF
128C--------------------------------------------------------
129C
130 DO i=1,jlt
131 nsvg(i)= nsv(cand_n(i))
132 l = cand_e(i)
133C
134 ix1(i)=irect(1,l)
135 ix2(i)=irect(2,l)
136 ix3(i)=irect(3,l)
137 ix4(i)=irect(4,l)
138 END DO
139C
140 DO i=1,jlt
141 xi(i)=x(1,nsvg(i))
142 yi(i)=x(2,nsvg(i))
143 zi(i)=x(3,nsvg(i))
144 x1(i)=x(1,ix1(i))
145 y1(i)=x(2,ix1(i))
146 z1(i)=x(3,ix1(i))
147 x2(i)=x(1,ix2(i))
148 y2(i)=x(2,ix2(i))
149 z2(i)=x(3,ix2(i))
150 x3(i)=x(1,ix3(i))
151 y3(i)=x(2,ix3(i))
152 z3(i)=x(3,ix3(i))
153 x4(i)=x(1,ix4(i))
154 y4(i)=x(2,ix4(i))
155 z4(i)=x(3,ix4(i))
156C
157 pene(i)=zero
158 END DO
159C--------------------------------------------------------
160C CAS DES PAQUETS MIXTES
161C--------------------------------------------------------
162 DO i=1,jlt
163 IF(ix3(i)/=ix4(i))THEN
164 x0(i) = fourth*(x1(i)+x2(i)+x3(i)+x4(i))
165 y0(i) = fourth*(y1(i)+y2(i)+y3(i)+y4(i))
166 z0(i) = fourth*(z1(i)+z2(i)+z3(i)+z4(i))
167 ELSE
168 x0(i) = x3(i)
169 y0(i) = y3(i)
170 z0(i) = z3(i)
171 ENDIF
172 ENDDO
173C
174 DO i=1,jlt
175C
176 x01(i) = x1(i) - x0(i)
177 y01(i) = y1(i) - y0(i)
178 z01(i) = z1(i) - z0(i)
179C
180 x02(i) = x2(i) - x0(i)
181 y02(i) = y2(i) - y0(i)
182 z02(i) = z2(i) - z0(i)
183C
184 x03(i) = x3(i) - x0(i)
185 y03(i) = y3(i) - y0(i)
186 z03(i) = z3(i) - z0(i)
187C
188 x04(i) = x4(i) - x0(i)
189 y04(i) = y4(i) - y0(i)
190 z04(i) = z4(i) - z0(i)
191C
192 xi0v(i) = x0(i) - xi(i)
193 yi0v(i) = y0(i) - yi(i)
194 zi0v(i) = z0(i) - zi(i)
195C
196 xi1(i) = x1(i) - xi(i)
197 yi1(i) = y1(i) - yi(i)
198 zi1(i) = z1(i) - zi(i)
199C
200 xi2(i) = x2(i) - xi(i)
201 yi2(i) = y2(i) - yi(i)
202 zi2(i) = z2(i) - zi(i)
203C
204 xi3(i) = x3(i) - xi(i)
205 yi3(i) = y3(i) - yi(i)
206 zi3(i) = z3(i) - zi(i)
207C
208 xi4(i) = x4(i) - xi(i)
209 yi4(i) = y4(i) - yi(i)
210 zi4(i) = z4(i) - zi(i)
211C
212 sx1 = yi0v(i)*zi1(i) - zi0v(i)*yi1(i)
213 sy1 = zi0v(i)*xi1(i) - xi0v(i)*zi1(i)
214 sz1 = xi0v(i)*yi1(i) - yi0v(i)*xi1(i)
215C
216 sx2 = yi0v(i)*zi2(i) - zi0v(i)*yi2(i)
217 sy2 = zi0v(i)*xi2(i) - xi0v(i)*zi2(i)
218 sz2 = xi0v(i)*yi2(i) - yi0v(i)*xi2(i)
219C
220 sx0 = y01(i)*z02(i) - z01(i)*y02(i)
221 sy0 = z01(i)*x02(i) - x01(i)*z02(i)
222 sz0 = x01(i)*y02(i) - y01(i)*x02(i)
223 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
224C
225 nn=sqrt(s2)
226 n11(i)=sx0*nn
227 n21(i)=sy0*nn
228 n31(i)=sz0*nn
229C
230 area(i)=nn
231C
232 lb1(i) = -(sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
233 lc1(i) = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
234C
235 sx3 = yi0v(i)*zi3(i) - zi0v(i)*yi3(i)
236 sy3 = zi0v(i)*xi3(i) - xi0v(i)*zi3(i)
237 sz3 = xi0v(i)*yi3(i) - yi0v(i)*xi3(i)
238C
239 sx0 = y02(i)*z03(i) - z02(i)*y03(i)
240 sy0 = z02(i)*x03(i) - x02(i)*z03(i)
241 sz0 = x02(i)*y03(i) - y02(i)*x03(i)
242 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
243C
244 nn=sqrt(s2)
245 n12(i)=sx0*nn
246 n22(i)=sy0*nn
247 n32(i)=sz0*nn
248C
249 area(i)=area(i)+nn
250C
251 lb2(i) = -(sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
252 lc2(i) = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
253C
254 sx4 = yi0v(i)*zi4(i) - zi0v(i)*yi4(i)
255 sy4 = zi0v(i)*xi4(i) - xi0v(i)*zi4(i)
256 sz4 = xi0v(i)*yi4(i) - yi0v(i)*xi4(i)
257C
258 sx0 = y03(i)*z04(i) - z03(i)*y04(i)
259 sy0 = z03(i)*x04(i) - x03(i)*z04(i)
260 sz0 = x03(i)*y04(i) - y03(i)*x04(i)
261 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
262C
263 nn=sqrt(s2)
264 n13(i)=sx0*nn
265 n23(i)=sy0*nn
266 n33(i)=sz0*nn
267C
268 area(i)=area(i)+nn
269C
270 lb3(i) = -(sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
271 lc3(i) = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
272C
273 sx0 = y04(i)*z01(i) - z04(i)*y01(i)
274 sy0 = z04(i)*x01(i) - x04(i)*z01(i)
275 sz0 = x04(i)*y01(i) - y04(i)*x01(i)
276 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
277C
278 nn=sqrt(s2)
279 n14(i)=sx0*nn
280 n24(i)=sy0*nn
281 n34(i)=sz0*nn
282C
283 area(i)=area(i)+nn
284 area(i)=half*area(i)
285C
286 lb4(i) = -(sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
287 lc4(i) = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
288 ENDDO
289C---------------------------------------------------------
290C
291 DO i=1,jlt
292 aaa = one/max(em30,x01(i)*x01(i)+y01(i)*y01(i)+z01(i)*z01(i))
293 hlc1(i)= lc1(i)*abs(lc1(i))*aaa
294 hlb4(i)= lb4(i)*abs(lb4(i))*aaa
295 al1(i) = -(xi0v(i)*x01(i)+yi0v(i)*y01(i)+zi0v(i)*z01(i))*aaa
296 al1(i) = max(zero,min(one,al1(i)))
297 aaa = one/max(em30,x02(i)*x02(i)+y02(i)*y02(i)+z02(i)*z02(i))
298 hlc2(i)= lc2(i)*abs(lc2(i))*aaa
299 hlb1(i)= lb1(i)*abs(lb1(i))*aaa
300 al2(i) = -(xi0v(i)*x02(i)+yi0v(i)*y02(i)+zi0v(i)*z02(i))*aaa
301 al2(i) = max(zero,min(one,al2(i)))
302 aaa = one/max(em30,x03(i)*x03(i)+y03(i)*y03(i)+z03(i)*z03(i))
303 hlc3(i)= lc3(i)*abs(lc3(i))*aaa
304 hlb2(i)= lb2(i)*abs(lb2(i))*aaa
305 al3(i) = -(xi0v(i)*x03(i)+yi0v(i)*y03(i)+zi0v(i)*z03(i))*aaa
306 al3(i) = max(zero,min(one,al3(i)))
307 aaa = one/max(em30,x04(i)*x04(i)+y04(i)*y04(i)+z04(i)*z04(i))
308 hlc4(i)= lc4(i)*abs(lc4(i))*aaa
309 hlb3(i)= lb3(i)*abs(lb3(i))*aaa
310 al4(i) = -(xi0v(i)*x04(i)+yi0v(i)*y04(i)+zi0v(i)*z04(i))*aaa
311 al4(i) = max(zero,min(one,al4(i)))
312 ENDDO
313C
314 DO i=1,jlt
315 x12 = x2(i) - x1(i)
316 y12 = y2(i) - y1(i)
317 z12 = z2(i) - z1(i)
318 la = one - lb1(i) - lc1(i)
319C HLA, HLB1, HLC1 necessaires pour triangle angle obtu
320 aaa = one / max(em20,x12*x12+y12*y12+z12*z12)
321 hla= la*abs(la) * aaa
322 inproj(1,i)=0
323 IF(la<zero.AND.
324 + hla<=hlb1(i).AND.hla<=hlc1(i))THEN
325 lb1(i) = (xi2(i)*x12+yi2(i)*y12+zi2(i)*z12) * aaa
326 lb1(i) = max(zero,min(one,lb1(i)))
327 lc1(i) = one - lb1(i)
328 inproj(1,i)=1
329 ELSEIF(lb1(i)<zero.AND.
330 + hlb1(i)<=hlc1(i).AND.hlb1(i)<=hla)THEN
331 lb1(i) = zero
332 lc1(i) = al2(i)
333 inproj(1,i)=1
334 ELSEIF(lc1(i)<zero.AND.
335 + hlc1(i)<=hla.AND.hlc1(i)<=hlb1(i))THEN
336 lc1(i) = zero
337 lb1(i) = al1(i)
338 inproj(1,i)=1
339 ENDIF
340 ENDDO
341C
342 DO i=1,jlt
343 IF(ix3(i)==ix4(i))cycle
344 x23 = x3(i) - x2(i)
345 y23 = y3(i) - y2(i)
346 z23 = z3(i) - z2(i)
347 la = one - lb2(i) - lc2(i)
348C HLA, HLB1, HLC1 necessaires pour triangle angle obtu
349 aaa = one / max(em20,x23*x23+y23*y23+z23*z23)
350 hla= la*abs(la) * aaa
351 inproj(2,i)=0
352 IF(la<zero.AND.
353 + hla<=hlb2(i).AND.hla<=hlc2(i))THEN
354 lb2(i) = (xi3(i)*x23+yi3(i)*y23+zi3(i)*z23)*aaa
355 lb2(i) = max(zero,min(one,lb2(i)))
356 lc2(i) = one - lb2(i)
357 inproj(2,i)=1
358 ELSEIF(lb2(i)<zero.AND.
359 + hlb2(i)<=hlc2(i).AND.hlb2(i)<=hla)THEN
360 lb2(i) = zero
361 lc2(i) = al3(i)
362 inproj(2,i)=1
363 ELSEIF(lc2(i)<zero.AND.
364 + hlc2(i)<=hla.AND.hlc2(i)<=hlb2(i))THEN
365 lc2(i) = zero
366 lb2(i) = al2(i)
367 inproj(2,i)=1
368 ENDIF
369 ENDDO
370C
371 DO i=1,jlt
372 IF(ix3(i)==ix4(i))cycle
373 x34 = x4(i) - x3(i)
374 y34 = y4(i) - y3(i)
375 z34 = z4(i) - z3(i)
376 la = one - lb3(i) - lc3(i)
377C HLA, HLB1, HLC1 necessaires pour triangle angle obtu
378 aaa = one / max(em20,x34*x34+y34*y34+z34*z34)
379 hla= la*abs(la) * aaa
380 inproj(3,i)=0
381 IF(la<zero.AND.
382 + hla<=hlb3(i).AND.hla<=hlc3(i))THEN
383 lb3(i) = (xi4(i)*x34+yi4(i)*y34+zi4(i)*z34)*aaa
384 lb3(i) = max(zero,min(one,lb3(i)))
385 lc3(i) = one - lb3(i)
386 inproj(3,i)=1
387 ELSEIF(lb3(i)<zero.AND.
388 + hlb3(i)<=hlc3(i).AND.hlb3(i)<=hla)THEN
389 lb3(i) = zero
390 lc3(i) = al4(i)
391 inproj(3,i)=1
392 ELSEIF(lc3(i)<zero.AND.
393 + hlc3(i)<=hla.AND.hlc3(i)<=hlb3(i))THEN
394 lc3(i) = zero
395 lb3(i) = al3(i)
396 inproj(3,i)=1
397 ENDIF
398 ENDDO
399C
400 DO i=1,jlt
401 IF(ix3(i)==ix4(i))cycle
402 x41 = x1(i) - x4(i)
403 y41 = y1(i) - y4(i)
404 z41 = z1(i) - z4(i)
405 la = one - lb4(i) - lc4(i)
406C HLA, HLB1, HLC1 necessaires pour triangle angle obtu
407 aaa = one / max(em20,x41*x41+y41*y41+z41*z41)
408 hla= la*abs(la) * aaa
409 inproj(4,i)=0
410 IF(la<zero.AND.
411 + hla<=hlb4(i).AND.hla<=hlc4(i))THEN
412 lb4(i) = (xi1(i)*x41+yi1(i)*y41+zi1(i)*z41)*aaa
413 lb4(i) = max(zero,min(one,lb4(i)))
414 lc4(i) = one - lb4(i)
415 inproj(4,i)=1
416 ELSEIF(lb4(i)<zero.AND.
417 + hlb4(i)<=hlc4(i).AND.hlb4(i)<=hla)THEN
418 lb4(i) = zero
419 lc4(i) = al1(i)
420 inproj(4,i)=1
421 ELSEIF(lc4(i)<zero.AND.
422 + hlc4(i)<=hla.AND.hlc4(i)<=hlb4(i))THEN
423 lc4(i) = zero
424 lb4(i) = al4(i)
425 inproj(4,i)=1
426 ENDIF
427 ENDDO
428C---------------------------------------------------------
429 DO i=1,jlt
430C
431 gap2(i)=gapv(i)*gapv(i)
432C
433 la1 = one - lb1(i) - lc1(i)
434 xp1(i) = lb1(i)*x1(i) + lc1(i)*x2(i) + la1*x0(i)
435 yp1(i) = lb1(i)*y1(i) + lc1(i)*y2(i) + la1*y0(i)
436 zp1(i) = lb1(i)*z1(i) + lc1(i)*z2(i) + la1*z0(i)
437C
438 nx1(i) = xi(i)-xp1(i)
439 ny1(i) = yi(i)-yp1(i)
440 nz1(i) = zi(i)-zp1(i)
441 p1(i) = nx1(i)*nx1(i) + ny1(i)*ny1(i) +nz1(i)*nz1(i)
442C
443 la2 = one - lb2(i) - lc2(i)
444 xp2(i) = lb2(i)*x2(i) + lc2(i)*x3(i) + la2*x0(i)
445 yp2(i) = lb2(i)*y2(i) + lc2(i)*y3(i) + la2*y0(i)
446 zp2(i) = lb2(i)*z2(i) + lc2(i)*z3(i) + la2*z0(i)
447C
448 nx2(i) = xi(i)-xp2(i)
449 ny2(i) = yi(i)-yp2(i)
450 nz2(i) = zi(i)-zp2(i)
451 p2(i) = nx2(i)*nx2(i) + ny2(i)*ny2(i) +nz2(i)*nz2(i)
452C
453 la3 = one - lb3(i) - lc3(i)
454 xp3(i) = lb3(i)*x3(i) + lc3(i)*x4(i) + la3*x0(i)
455 yp3(i) = lb3(i)*y3(i) + lc3(i)*y4(i) + la3*y0(i)
456 zp3(i) = lb3(i)*z3(i) + lc3(i)*z4(i) + la3*z0(i)
457C
458 nx3(i) = xi(i)-xp3(i)
459 ny3(i) = yi(i)-yp3(i)
460 nz3(i) = zi(i)-zp3(i)
461 p3(i) = nx3(i)*nx3(i) + ny3(i)*ny3(i) +nz3(i)*nz3(i)
462C
463 la4 = one - lb4(i) - lc4(i)
464 xp4(i) = lb4(i)*x4(i) + lc4(i)*x1(i) + la4*x0(i)
465 yp4(i) = lb4(i)*y4(i) + lc4(i)*y1(i) + la4*y0(i)
466 zp4(i) = lb4(i)*z4(i) + lc4(i)*z1(i) + la4*z0(i)
467C
468 nx4(i) = xi(i)-xp4(i)
469 ny4(i) = yi(i)-yp4(i)
470 nz4(i) = zi(i)-zp4(i)
471 p4(i) = nx4(i)*nx4(i) + ny4(i)*ny4(i) +nz4(i)*nz4(i)
472C
473 ENDDO
474C
475 DO i=1,jlt
476 END DO
477C
478 DO i=1,jlt
479C
480 IF(ix3(i)/=ix4(i))THEN
481C
482 h =nx1(i)*n11(i) + ny1(i)*n21(i) +nz1(i)*n31(i)
483 ll=p1(i)-h*h
484 IF(inproj(1,i)/=0.AND.ll >= gap2(i))THEN
485 p1(i)=zero
486 ELSE
487 p1(i)=max(zero,gapv(i)-abs(h))
488 END IF
489C
490 h =nx2(i)*n12(i) + ny2(i)*n22(i) +nz2(i)*n32(i)
491 ll=p2(i)-h*h
492 IF(inproj(2,i)/=0.AND.ll >= gap2(i))THEN
493 p2(i)=zero
494 ELSE
495 p2(i)=max(zero,gapv(i)-abs(h))
496 END IF
497C
498 h =nx3(i)*n13(i) + ny3(i)*n23(i) +nz3(i)*n33(i)
499 ll=p3(i)-h*h
500 IF(inproj(3,i)/=0.AND.ll >= gap2(i))THEN
501 p3(i)=zero
502 ELSE
503 p3(i)=max(zero,gapv(i)-abs(h))
504 END IF
505C
506 h =nx4(i)*n14(i) + ny4(i)*n24(i) +nz4(i)*n34(i)
507 ll=p4(i)-h*h
508 IF(inproj(4,i)/=0.AND.ll >= gap2(i))THEN
509 p4(i)=zero
510 ELSE
511 p4(i)=max(zero,gapv(i)-abs(h))
512 END IF
513C
514 pene(i) = max(p1(i),p2(i),p3(i),p4(i))
515C
516 n1(i) = p1(i)*n11(i)+p2(i)*n12(i)+p3(i)*n13(i)+p4(i)*n14(i)
517 n2(i) = p1(i)*n21(i)+p2(i)*n22(i)+p3(i)*n23(i)+p4(i)*n24(i)
518 n3(i) = p1(i)*n31(i)+p2(i)*n32(i)+p3(i)*n33(i)+p4(i)*n34(i)
519 nn=sqrt(n1(i)*n1(i)+n2(i)*n2(i)+n3(i)*n3(i))
520 nn=one/max(em20,nn)
521 n1(i)=n1(i)*nn
522 n2(i)=n2(i)*nn
523 n3(i)=n3(i)*nn
524C
525 la1 = one - lb1(i) - lc1(i)
526 la2 = one - lb2(i) - lc2(i)
527 la3 = one - lb3(i) - lc3(i)
528 la4 = one - lb4(i) - lc4(i)
529C
530 h0 = fourth *
531 . (p1(i)*la1 + p2(i)*la2 + p3(i)*la3 + p4(i)*la4)
532 h1(i) = h0 + p1(i) * lb1(i) + p4(i) * lc4(i)
533 h2(i) = h0 + p2(i) * lb2(i) + p1(i) * lc1(i)
534 h3(i) = h0 + p3(i) * lb3(i) + p2(i) * lc2(i)
535 h4(i) = h0 + p4(i) * lb4(i) + p3(i) * lc3(i)
536 h00 = one/max(em20,h1(i) + h2(i) + h3(i) + h4(i))
537 h1(i) = h1(i) * h00
538 h2(i) = h2(i) * h00
539 h3(i) = h3(i) * h00
540 h4(i) = h4(i) * h00
541C
542 ELSE ! IF(IX3(I)/=IX4(I))THEN
543C
544 h1(i) = lb1(i)
545 h2(i) = lc1(i)
546 h3(i) = one - lb1(i) - lc1(i)
547 h4(i) = zero
548C
549 n1(i) = n11(i)
550 n2(i) = n21(i)
551 n3(i) = n31(i)
552C
553 h =nx1(i)*n11(i) + ny1(i)*n21(i) +nz1(i)*n31(i)
554 ll=p1(i)-h*h
555 IF(inproj(1,i)/=0.AND.ll >= gap2(i))THEN
556 pene(i)=zero
557 ELSE
558 pene(i)=max(zero,gapv(i)-abs(h))
559 END IF
560 END IF
561 ENDDO
562C
563 DO 410 i=1,jlt
564 xp(i)=h1(i)*x1(i)+h2(i)*x2(i)+h3(i)*x3(i)+h4(i)*x4(i)
565 yp(i)=h1(i)*y1(i)+h2(i)*y2(i)+h3(i)*y3(i)+h4(i)*y4(i)
566 zp(i)=h1(i)*z1(i)+h2(i)*z2(i)+h3(i)*z3(i)+h4(i)*z4(i)
567 410 CONTINUE
568C---------------------------------------------------------
569 DO i=1,jlt
570 nx1(i) = xi(i)-xp(i)
571 ny1(i) = yi(i)-yp(i)
572 nz1(i) = zi(i)-zp(i)
573 ENDDO
574C---------------------
575 DO i=1,jlt
576C
577 nvers(1,i) = zero
578 firstimp(1,i) = zero
579C
580 IF(pene(i)==zero)cycle
581C--------------------------------------------------------
582 side=n1(i)*nx1(i)+n2(i)*ny1(i)+n3(i)*nz1(i)
583C
584 IF(ifpen(i)==0)THEN
585 firstimp(1,i) = sign(one,side)
586 IF(firstimp(1,i) < zero)THEN
587 n1(i) = -n1(i)
588 n2(i) = -n2(i)
589 n3(i) = -n3(i)
590 END IF
591 nvers(1,i) = one
592C 1st impact below gap (sorting security)
593 ELSE ! IF(IFPEN(I)==0.OR.TT==ZERO)THEN
594 IF(ifpen(i) < 0)THEN
595 n1(i) = -n1(i)
596 n2(i) = -n2(i)
597 n3(i) = -n3(i)
598C SIDE = -SIDE
599 END IF
600 nvers(1,i)= sign(one,side*ifpen(i))
601 END IF
602C---------------------
603 IF(nvers(1,i)==zero)THEN
604 pene(i)=zero
605 cycle
606 END IF
607C
608C attention a la traversee de la coque => 1E20...,
609C prendre normale = normale au triangle
610 nn=one/
611 . max(em20,sqrt(n1(i)*n1(i)+n2(i)*n2(i)+n3(i)*n3(i)))
612 n1(i)=n1(i)*nn
613 n2(i)=n2(i)*nn
614 n3(i)=n3(i)*nn
615C
616 nx(i)=n1(i)
617 ny(i)=n2(i)
618 nz(i)=n3(i)
619C
620 END DO
621C---------------------
622 DO i=1,jlt
623 IF(pene(i)/=zero)THEN
624 IF(firstimp(1,i) > zero)THEN
625 ifpen(i) = 1
626 ELSEIF(firstimp(1,i) < zero)THEN
627 ifpen(i) =-1
628 END IF
629 END IF
630 ENDDO
631C---------------------
632 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
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