OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i10keg3.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!|| i10keg3 ../engine/source/interfaces/int10/i10keg3.F
25!||--- called by ------------------------------------------------------
26!|| i10ke3 ../engine/source/interfaces/int10/i10ke3.F
27!||--- uses -----------------------------------------------------
28!|| imp_intm ../engine/share/modules/imp_intm.F
29!||====================================================================
30 SUBROUTINE i10keg3(JLT ,A ,V ,MS ,FRIC ,
31 1 NX1 ,NX2 ,NX3 ,NX4 ,NY1 ,
32 2 NY2 ,NY3 ,NY4 ,NZ1 ,NZ2 ,
33 3 NZ3 ,NZ4 ,LB1 ,LB2 ,LB3 ,
34 4 LB4 ,LC1 ,LC2 ,LC3 ,LC4 ,
35 5 P1 ,P2 ,P3 ,P4 ,NIN ,
36 6 IX1 ,IX2 ,IX3 ,IX4 ,NSVG ,
37 7 GAPV ,ITIED ,CAND_F ,INDEX ,STIF ,
38 8 VXI ,VYI ,VZI ,MSI ,X1 ,
39 9 X2 ,X3 ,X4 ,Y1 ,Y2 ,
40 A Y3 ,Y4 ,Z1 ,Z2 ,Z3 ,
41 B Z4 ,KI11 ,KI12 ,KJ11 ,KJ12 ,
42 C KK11 ,KK12 ,KL11 ,KL12 ,OFF ,
43 D SCALK ,LREM )
44C-----------------------------------------------
45C M o d u l e s
46C-----------------------------------------------
47 USE imp_intm
48C-----------------------------------------------
49C I m p l i c i t T y p e s
50C-----------------------------------------------
51#include "implicit_f.inc"
52C-----------------------------------------------
53C G l o b a l P a r a m e t e r s
54C-----------------------------------------------
55#include "mvsiz_p.inc"
56C-----------------------------------------------
57C C o m m o n B l o c k s
58C-----------------------------------------------
59#include "com01_c.inc"
60#include "com08_c.inc"
61#include "impl1_c.inc"
62C-----------------------------------------------
63C D u m m y A r g u m e n t s
64C-----------------------------------------------
65 INTEGER JLT, ITIED,NIN,LREM
66 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
67 . NSVG(MVSIZ), INDEX(MVSIZ)
68 my_real
69 . A(3,*), MS(*), V(3,*),X1(*),X2(*),X3(*),X4(*),
70 . Y1(*),Y2(*),Y3(*),Y4(*),Z1(*),Z2(*),Z3(*),Z4(*),
71 . CAND_F(6,*),FRIC,OFF(*),SCALK,
72 . VXI(MVSIZ),VYI(MVSIZ),VZI(MVSIZ),MSI(MVSIZ)
73 my_real
74 . NX1(MVSIZ), NX2(MVSIZ), NX3(MVSIZ), NX4(MVSIZ),
75 . NY1(MVSIZ), NY2(MVSIZ), NY3(MVSIZ), NY4(MVSIZ),
76 . NZ1(MVSIZ), NZ2(MVSIZ), NZ3(MVSIZ), NZ4(MVSIZ),
77 . LB1(MVSIZ), LB2(MVSIZ), LB3(MVSIZ), LB4(MVSIZ),
78 . lc1(mvsiz), lc2(mvsiz), lc3(mvsiz), lc4(mvsiz),
79 . p1(mvsiz), p2(mvsiz), p3(mvsiz), p4(mvsiz), stif(mvsiz),
80 . gapv(mvsiz),ki11(3,3,mvsiz),kj11(3,3,mvsiz),
81 . kk11(3,3,mvsiz),kl11(3,3,mvsiz),ki12(3,3,mvsiz),
82 . kj12(3,3,mvsiz),kk12(3,3,mvsiz),kl12(3,3,mvsiz)
83C-----------------------------------------------
84C L o c a l V a r i a b l e s
85C-----------------------------------------------
86 INTEGER I, J1, J, K,IG,ISF,NN,NS,JLTF,NE,II
87 my_real
88 . N1(MVSIZ), N2(MVSIZ), N3(MVSIZ), PENE(MVSIZ),
89 . H1(MVSIZ), H2(MVSIZ), H3(MVSIZ), H4(MVSIZ),
90 . VT1(MVSIZ), VT2(MVSIZ),FNI(MVSIZ),
91 . NX(MVSIZ), NY(MVSIZ), NZ(MVSIZ),
92 . T1X(MVSIZ),T1Y(MVSIZ),T1Z(MVSIZ),
93 . T2X(MVSIZ),T2Y(MVSIZ),T2Z(MVSIZ),NORMINV,
94 . VX(MVSIZ), VY(MVSIZ), VZ(MVSIZ), VN(MVSIZ),
95 . S2,FAC,FACF, H0, LA1, LA2, LA3, LA4,FACT(MVSIZ),
96 . D1,D2,D3,D4,A1,A2,A3,A4,KN(4,MVSIZ),Q(3,3,MVSIZ),FAC10
97 my_real
98 . PREC,Q11,Q12,Q13,Q22,Q23,Q33,H00,VTX,VTY,VTZ,VT,
99 . kt1,kt2,kt3,kt4,q1,q2
100C-----------------------------------------------
101 fric =one
102 fac10 = ten
103 IF (imp_int7==3) THEN
104 DO i=1,jlt
105 d1 = sqrt(p1(i))
106 p1(i) = fourth*gapv(i)
107 d2 = sqrt(p2(i))
108 p2(i) = fourth*gapv(i)
109 d3 = sqrt(p3(i))
110 p3(i) = fourth*gapv(i)
111 d4 = sqrt(p4(i))
112 p4(i) = fourth*gapv(i)
113 ENDDO
114 ELSE
115 DO i=1,jlt
116C
117 d1 = sqrt(p1(i))
118 p1(i) = max(zero, gapv(i) - d1)
119C
120 d2 = sqrt(p2(i))
121 p2(i) = max(zero, gapv(i) - d2)
122C
123 d3 = sqrt(p3(i))
124 p3(i) = max(zero, gapv(i) - d3)
125C
126 d4 = sqrt(p4(i))
127 p4(i) = max(zero, gapv(i) - d4)
128 ENDDO
129 ENDIF !(IMP_INT7==3)
130C
131 DO i=1,jlt
132 IF(ix3(i)/=ix4(i))THEN
133 pene(i) = max(p1(i),p2(i),p3(i),p4(i))
134C
135 la1 = one - lb1(i) - lc1(i)
136 la2 = one - lb2(i) - lc2(i)
137 la3 = one - lb3(i) - lc3(i)
138 la4 = one - lb4(i) - lc4(i)
139C
140 h0 = fourth *
141 . (p1(i)*la1 + p2(i)*la2 + p3(i)*la3 + p4(i)*la4)
142 h1(i) = h0 + p1(i) * lb1(i) + p4(i) * lc4(i)
143 h2(i) = h0 + p2(i) * lb2(i) + p1(i) * lc1(i)
144 h3(i) = h0 + p3(i) * lb3(i) + p2(i) * lc2(i)
145 h4(i) = h0 + p4(i) * lb4(i) + p3(i) * lc3(i)
146 h00 = one/max(em20,h1(i) + h2(i) + h3(i) + h4(i))
147 h1(i) = h1(i) * h00
148 h2(i) = h2(i) * h00
149 h3(i) = h3(i) * h00
150 h4(i) = h4(i) * h00
151C
152 ELSE
153 pene(i) = p1(i)
154 n1(i) = nx1(i)
155 n2(i) = ny1(i)
156 n3(i) = nz1(i)
157 h1(i) = lb1(i)
158 h2(i) = lc1(i)
159 h3(i) = one - lb1(i) - lc1(i)
160 h4(i) = zero
161 ENDIF
162 ENDDO
163C
164C DO I=1,JLT
165C S2 = ONE/MAX(EM30,SQRT(N1(I)**2 + N2(I)**2 + N3(I)**2))
166C N1(I) = N1(I)*S2
167C N2(I) = N2(I)*S2
168C N3(I) = N3(I)*S2
169C ENDDO
170C
171 DO i=1,jlt
172C correction hourglass
173 IF(ix3(i)/=ix4(i))THEN
174 h0 = -fourth*(h1(i) - h2(i) + h3(i) - h4(i))
175 h0 = min(h0,h2(i),h4(i))
176 h0 = max(h0,-h1(i),-h3(i))
177 h1(i) = h1(i) + h0
178 h2(i) = h2(i) - h0
179 h3(i) = h3(i) + h0
180 h4(i) = h4(i) - h0
181 ENDIF
182 ENDDO
183C-------------------------------------------
184 DO i=1,jlt
185 ii = index(i)
186 IF(cand_f(1,ii)==zero)THEN
187C------------------------------------
188C 1ER IMPACT ou PAS d'IMPACT
189C------------------------------------
190C c'est fait dans i10for3.F CAND_F(4,II) = H1(I)
191C CAND_F(5,II) = H2(I)
192C CAND_F(6,II) = H3(I)
193 ELSE
194C------------------------------------
195C IMPACTS SUIVANTS
196C------------------------------------
197 h1(i) = cand_f(4,ii)
198 h2(i) = cand_f(5,ii)
199 h3(i) = cand_f(6,ii)
200 h4(i) = one - h1(i) - h2(i) - h3(i)
201 ENDIF
202 ENDDO
203C
204 DO i=1,jlt
205 vx(i) = vxi(i) - h1(i)*v(1,ix1(i)) - h2(i)*v(1,ix2(i))
206 . - h3(i)*v(1,ix3(i)) - h4(i)*v(1,ix4(i))
207 vy(i) = vyi(i) - h1(i)*v(2,ix1(i)) - h2(i)*v(2,ix2(i))
208 . - h3(i)*v(2,ix3(i)) - h4(i)*v(2,ix4(i))
209 vz(i) = vzi(i) - h1(i)*v(3,ix1(i)) - h2(i)*v(3,ix2(i))
210 . - h3(i)*v(3,ix3(i)) - h4(i)*v(3,ix4(i))
211 ENDDO
212C
213 DO i=1,jlt
214 t1x(i) = x3(i) - x1(i)
215 t1y(i) = y3(i) - y1(i)
216 t1z(i) = z3(i) - z1(i)
217 norminv = one/sqrt(t1x(i)**2+t1y(i)**2+t1z(i)**2)
218 t1x(i) = t1x(i)*norminv
219 t1y(i) = t1y(i)*norminv
220 t1z(i) = t1z(i)*norminv
221C
222 t2x(i) = x4(i) - x2(i)
223 t2y(i) = y4(i) - y2(i)
224 t2z(i) = z4(i) - z2(i)
225C
226 nx(i) = t1y(i)*t2z(i) - t1z(i)*t2y(i)
227 ny(i) = t1z(i)*t2x(i) - t1x(i)*t2z(i)
228 nz(i) = t1x(i)*t2y(i) - t1y(i)*t2x(i)
229 norminv = one/sqrt(nx(i)**2+ny(i)**2+nz(i)**2)
230 nx(i) = nx(i)*norminv
231 ny(i) = ny(i)*norminv
232 nz(i) = nz(i)*norminv
233C
234 t2x(i) = ny(i)*t1z(i) - nz(i)*t1y(i)
235 t2y(i) = nz(i)*t1x(i) - nx(i)*t1z(i)
236 t2z(i) = nx(i)*t1y(i) - ny(i)*t1x(i)
237C
238 vn(i) = vx(i)*nx(i) + vy(i)*ny(i) + vz(i)*nz(i)
239 vt1(i) = vx(i)*t1x(i) + vy(i)*t1y(i) + vz(i)*t1z(i)
240 vt2(i) = vx(i)*t2x(i) + vy(i)*t2y(i) + vz(i)*t2z(i)
241 ENDDO
242C
243 DO i=1,jlt
244 IF(pene(i)==zero.AND.cand_f(1,index(i))==zero)THEN
245C------------------------------------
246C PAS ENCORE D'IMPACT OU REBOND
247C------------------------------------
248 vn(i) = zero
249 vt1(i) = zero
250 vt2(i) = zero
251 ENDIF
252 ENDDO
253C
254 DO i=1,jlt
255 ii = index(i)
256 fni(i) = cand_f(1,ii) + vn(i) * dt1 * stif(i)
257 ENDDO
258C
259 DO 100 i=1,jlt
260 IF(itied==0)THEN
261 IF(cand_f(1,index(i))*fni(i)<zero)THEN
262C------------------------------------
263C REBOND
264C------------------------------------
265 fni(i) = zero
266 vn(i) = zero
267 vt1(i) = zero
268 vt2(i) = zero
269 stif(i) = zero
270 ELSE
271C--------
272 ENDIF
273 ELSE
274 stif(i) = stif(i) * abs(vn(i)) * dt1/max(pene(i),em10)
275 ENDIF
276C
277 100 CONTINUE
278C
279C---------------------------------
280C ----sans frottement d'abord---
281 DO i=1,jlt
282 IF (abs(vt1(i))>zero.OR.abs(vt2(i))>zero) THEN
283 q(1,1,i)=t1x(i)
284 q(1,2,i)=t1y(i)
285 q(1,3,i)=t1z(i)
286 q(3,1,i)=nx(i)
287 q(3,2,i)=ny(i)
288 q(3,3,i)=nz(i)
289 q(2,1,i)=t2x(i)
290 q(2,2,i)=t2y(i)
291 q(2,3,i)=t2z(i)
292 fact(i)=fric
293 ELSE
294 fact(i)=zero
295 ENDIF
296 ENDDO
297 IF (scalk<0) THEN
298 isf=1
299 ELSE
300 isf=0
301 ENDIF
302 facf=fac10*abs(scalk)
303 IF (isf==1) THEN
304 DO i=1,jlt
305 IF (vn(i)>zero) THEN
306 fac=stif(i)*facf
307 ELSEIF (vn(i)<zero) THEN
308 fac=stif(i)/facf
309 ELSE
310 fac=stif(i)
311 ENDIF
312 kn(1,i)=fac*h1(i)
313 kn(2,i)=fac*h2(i)
314 kn(3,i)=fac*h3(i)
315 kn(4,i)=fac*h4(i)
316 fact(i)=fac*fact(i)
317 ENDDO
318 ELSE
319 DO i=1,jlt
320 fac=stif(i)*facf
321 kn(1,i)=fac*h1(i)
322 kn(2,i)=fac*h2(i)
323 kn(3,i)=fac*h3(i)
324 kn(4,i)=fac*h4(i)
325 fact(i)=fac*fact(i)
326 ENDDO
327 ENDIF
328 DO i=1,jlt
329 q11=nx(i)*nx(i)
330 q12=nx(i)*ny(i)
331 q13=nx(i)*nz(i)
332 q22=ny(i)*ny(i)
333 q23=ny(i)*nz(i)
334 q33=nz(i)*nz(i)
335 ki11(1,1,i)=kn(1,i)*q11
336 ki11(1,2,i)=kn(1,i)*q12
337 ki11(1,3,i)=kn(1,i)*q13
338 ki11(2,2,i)=kn(1,i)*q22
339 ki11(2,3,i)=kn(1,i)*q23
340 ki11(3,3,i)=kn(1,i)*q33
341 kj11(1,1,i)=kn(2,i)*q11
342 kj11(1,2,i)=kn(2,i)*q12
343 kj11(1,3,i)=kn(2,i)*q13
344 kj11(2,2,i)=kn(2,i)*q22
345 kj11(2,3,i)=kn(2,i)*q23
346 kj11(3,3,i)=kn(2,i)*q33
347 kk11(1,1,i)=kn(3,i)*q11
348 kk11(1,2,i)=kn(3,i)*q12
349 kk11(1,3,i)=kn(3,i)*q13
350 kk11(2,2,i)=kn(3,i)*q22
351 kk11(2,3,i)=kn(3,i)*q23
352 kk11(3,3,i)=kn(3,i)*q33
353 kl11(1,1,i)=kn(4,i)*q11
354 kl11(1,2,i)=kn(4,i)*q12
355 kl11(1,3,i)=kn(4,i)*q13
356 kl11(2,2,i)=kn(4,i)*q22
357 kl11(2,3,i)=kn(4,i)*q23
358 kl11(3,3,i)=kn(4,i)*q33
359 ENDDO
360C ----avec frottement ---
361 DO j=1,3
362 DO k=j,3
363 DO i=1,jlt
364 IF (fact(i)>zero) THEN
365 q1 =q(1,j,i)*q(1,k,i)
366 q2 =q(2,j,i)*q(2,k,i)
367 fac=fact(i)*(q1+q2)
368 kt1=fac*h1(i)
369 ki11(j,k,i)=ki11(j,k,i)+kt1
370 kt2=fac*h2(i)
371 kj11(j,k,i)=kj11(j,k,i)+kt2
372 kt3=fac*h3(i)
373 kk11(j,k,i)=kk11(j,k,i)+kt3
374 kt4=fac*h4(i)
375 kl11(j,k,i)=kl11(j,k,i)+kt4
376 ENDIF
377 ENDDO
378 ENDDO
379 ENDDO
380C
381 DO j=1,3
382 DO k=j,3
383 DO i=1,jlt
384 ki12(j,k,i)=-ki11(j,k,i)
385 kj12(j,k,i)=-kj11(j,k,i)
386 kk12(j,k,i)=-kk11(j,k,i)
387 kl12(j,k,i)=-kl11(j,k,i)
388 ENDDO
389 ENDDO
390 ENDDO
391 DO j=1,3
392 DO k=j+1,3
393 DO i=1,jlt
394 ki12(k,j,i)=-ki11(j,k,i)
395 kj12(k,j,i)=-kj11(j,k,i)
396 kk12(k,j,i)=-kk11(j,k,i)
397 kl12(k,j,i)=-kl11(j,k,i)
398 ENDDO
399 ENDDO
400 ENDDO
401C
402 DO i=1,jlt
403 off(i)=one
404 ENDDO
405 IF (nspmd>1) THEN
406 IF ((intp_d)>0) THEN
407 DO i=1,jlt
408 IF(nsvg(i)<0) THEN
409 nn=-nsvg(i)
410 ns=ind_int(nin)%P(nn)
411C---------pour diag_ss---
412 ffi(1,ns)=zero
413 ffi(2,ns)=zero
414 ffi(3,ns)=zero
415 dfi(1,ns)=zero
416 dfi(2,ns)=zero
417 dfi(3,ns)=zero
418 ENDIF
419 ENDDO
420 ELSE
421 jltf = 0
422 DO i=1,jlt
423 IF(nsvg(i)<0) THEN
424 nn=-nsvg(i)
425 jltf = jltf + 1
426 ne=shf_int(nin) + jltf +lrem
427 ns=ind_int(nin)%P(nn)
428 stifs(ne)=stif(i)
429 h_e(1,ne)=h1(i)
430 h_e(2,ne)=h2(i)
431 h_e(3,ne)=h3(i)
432 h_e(4,ne)=h4(i)
433 n_e(1,ne)=nx(i)
434 n_e(2,ne)=ny(i)
435 n_e(3,ne)=nz(i)
436C---------pour temporairement diag_ss---
437 ffi(1,ns)=zero
438 ffi(2,ns)=zero
439 ffi(3,ns)=zero
440 dfi(1,ns)=zero
441 dfi(2,ns)=zero
442 dfi(3,ns)=zero
443 ENDIF
444 ENDDO
445 ENDIF
446 ENDIF
447C
448 RETURN
449 END
450!||====================================================================
451!|| i10frf3 ../engine/source/interfaces/int10/i10keg3.F
452!||--- called by ------------------------------------------------------
453!|| i10forcf3 ../engine/source/interfaces/int10/i10ke3.F
454!||--- uses -----------------------------------------------------
455!|| imp_intm ../engine/share/modules/imp_intm.f
456!||====================================================================
457 SUBROUTINE i10frf3(JLT ,A ,V ,MS ,FRIC ,
458 1 N1 ,N2 ,N3 ,T1X ,T1Y ,
459 2 T1Z ,H1 ,H2 ,H3 ,H4 ,
460 3 IX1 ,IX2 ,IX3 ,IX4 ,INDEX ,
461 4 VXI ,VYI ,VZI ,MSI ,DXI ,
462 5 DYI ,DZI ,STIF ,NIN ,D ,
463 6 SCALK )
464C-----------------------------------------------
465C M o d u l e s
466C-----------------------------------------------
467 USE imp_intm
468C-----------------------------------------------
469C I m p l i c i t T y p e s
470C-----------------------------------------------
471#include "implicit_f.inc"
472C-----------------------------------------------
473C G l o b a l P a r a m e t e r s
474C-----------------------------------------------
475#include "mvsiz_p.inc"
476C-----------------------------------------------
477C D u m m y A r g u m e n t s
478C-----------------------------------------------
479 INTEGER JLT, INACTI,NIN
480 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
481 . INDEX(MVSIZ)
482 my_real
483 . A(3,*), MS(*), V(3,*),D(3,*),
484 . FRIC,SCALK,DXI(MVSIZ),DYI(MVSIZ),DZI(MVSIZ),
485 . H1(MVSIZ), H2(MVSIZ), H3(MVSIZ), H4(MVSIZ),
486 . vxi(mvsiz),vyi(mvsiz),vzi(mvsiz),msi(mvsiz)
487 my_real
488 . n1(mvsiz), n2(mvsiz), n3(mvsiz),stif(mvsiz),
489 . t1x(mvsiz), t1y(mvsiz), t1z(mvsiz)
490C-----------------------------------------------
491C L o c a l V a r i a b l e s
492C-----------------------------------------------
493 INTEGER I, J1, J, K,IG,ISF,NN,NS,NI
494 my_real
495 . VX(MVSIZ), VY(MVSIZ), VZ(MVSIZ), VN(MVSIZ),
496 . DX(MVSIZ), DY(MVSIZ), DZ(MVSIZ), DN(MVSIZ),
497 . DNI(MVSIZ),D1T(MVSIZ),D2T(MVSIZ), DTI1(MVSIZ),
498 . dti2(mvsiz),s2,facn(mvsiz),facf, fact(mvsiz),fac10
499 my_real
500 . fx,fy,fz,fn,ft1,ft2,fni,fti1,fti2,vtx,vty,vtz,vt,
501 . t2x(mvsiz), t2y(mvsiz), t2z(mvsiz),
502 . kt1,kt2,kt3,kt4,q1,q2
503
504C-----------------------------------------------
505C
506 fac10 = ten
507 DO i=1,jlt
508 vx(i) = vxi(i) - h1(i)*v(1,ix1(i)) - h2(i)*v(1,ix2(i))
509 . - h3(i)*v(1,ix3(i)) - h4(i)*v(1,ix4(i))
510 vy(i) = vyi(i) - h1(i)*v(2,ix1(i)) - h2(i)*v(2,ix2(i))
511 . - h3(i)*v(2,ix3(i)) - h4(i)*v(2,ix4(i))
512 vz(i) = vzi(i) - h1(i)*v(3,ix1(i)) - h2(i)*v(3,ix2(i))
513 . - h3(i)*v(3,ix3(i)) - h4(i)*v(3,ix4(i))
514 vn(i) = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
515 dx(i) = dxi(i) - h1(i)*d(1,ix1(i)) - h2(i)*d(1,ix2(i))
516 . - h3(i)*d(1,ix3(i)) - h4(i)*d(1,ix4(i))
517 dy(i) = dyi(i) - h1(i)*d(2,ix1(i)) - h2(i)*d(2,ix2(i))
518 . - h3(i)*d(2,ix3(i)) - h4(i)*d(2,ix4(i))
519 dz(i) = dzi(i) - h1(i)*d(3,ix1(i)) - h2(i)*d(3,ix2(i))
520 . - h3(i)*d(3,ix3(i)) - h4(i)*d(3,ix4(i))
521 dn(i) = n1(i)*dx(i) + n2(i)*dy(i) + n3(i)*dz(i)
522 dni(i) = n1(i)*dxi(i) + n2(i)*dyi(i) + n3(i)*dzi(i)
523 ENDDO
524C-------------------------------------------
525 DO i=1,jlt
526 t2x(i) = n2(i)*t1z(i) - n3(i)*t1y(i)
527 t2y(i) = n3(i)*t1x(i) - n1(i)*t1z(i)
528 t2z(i) = n1(i)*t1y(i) - n2(i)*t1x(i)
529 d1t(i) = t1x(i)*dx(i) + t1y(i)*dy(i) + t1z(i)*dz(i)
530 d2t(i) = t2x(i)*dx(i) + t2y(i)*dy(i) + t2z(i)*dz(i)
531 dti1(i) = t1x(i)*dxi(i) + t1y(i)*dyi(i) + t1z(i)*dzi(i)
532 dti2(i) = t2x(i)*dxi(i) + t2y(i)*dyi(i) + t2z(i)*dzi(i)
533 ENDDO
534 IF (scalk<0) THEN
535 isf=1
536 ELSE
537 isf=0
538 ENDIF
539 facf=fac10*abs(scalk)
540 IF (isf==1) THEN
541 DO i=1,jlt
542 IF (vn(i)>zero) THEN
543 facn(i)=stif(i)*facf
544 ELSEIF (vn(i)<zero) THEN
545 facn(i)=stif(i)/facf
546 ELSE
547 facn(i)=stif(i)
548 ENDIF
549 fact(i)=facn(i)*fric
550 ENDDO
551 ELSE
552 DO i=1,jlt
553 facn(i)=stif(i)*facf
554 fact(i)=facn(i)*fric
555 ENDDO
556 ENDIF
557C--------partie NML-------
558 DO i=1,jlt
559 fn = -facn(i)*dni(i)
560 fx=fn*n1(i)
561 fy=fn*n2(i)
562 fz=fn*n3(i)
563 ft1 = -fact(i)*dti1(i)
564 ft2 = -fact(i)*dti2(i)
565 fx = fx + ft1*t1x(i)+ ft2*t2x(i)
566 fy = fy + ft1*t1y(i)+ ft2*t2y(i)
567 fz = fz + ft1*t1z(i)+ ft2*t2z(i)
568 a(1,ix1(i))=a(1,ix1(i))+fx*h1(i)
569 a(1,ix2(i))=a(1,ix2(i))+fx*h2(i)
570 a(1,ix3(i))=a(1,ix3(i))+fx*h3(i)
571 a(1,ix4(i))=a(1,ix4(i))+fx*h4(i)
572 a(2,ix1(i))=a(2,ix1(i))+fy*h1(i)
573 a(2,ix2(i))=a(2,ix2(i))+fy*h2(i)
574 a(2,ix3(i))=a(2,ix3(i))+fy*h3(i)
575 a(2,ix4(i))=a(2,ix4(i))+fy*h4(i)
576 a(3,ix1(i))=a(3,ix1(i))+fz*h1(i)
577 a(3,ix2(i))=a(3,ix2(i))+fz*h2(i)
578 a(3,ix3(i))=a(3,ix3(i))+fz*h3(i)
579 a(3,ix4(i))=a(3,ix4(i))+fz*h4(i)
580 ENDDO
581C--------partie NSL-------
582 DO i=1,jlt
583 fni = facn(i)*dn(i)
584 fx=fni*n1(i)
585 fy=fni*n2(i)
586 fz=fni*n3(i)
587 fti1 = fact(i)*d1t(i)
588 fti2 = fact(i)*d2t(i)
589 fx = fx + fti1*t1x(i)+ fti2*t2x(i)
590 fy = fy + fti1*t1y(i)+ fti2*t2y(i)
591 fz = fz + fti1*t1z(i)+ fti2*t2z(i)
592 ni = index(i)
593 ffi(1,ni)=ffi(1,ni)+fx
594 ffi(2,ni)=ffi(2,ni)+fy
595 ffi(3,ni)=ffi(3,ni)+fz
596 ENDDO
597C
598 RETURN
599 END
600!||====================================================================
601!|| i10kfor3 ../engine/source/interfaces/int10/i10keg3.F
602!||--- called by ------------------------------------------------------
603!|| i10fku3 ../engine/source/interfaces/int10/i10ke3.F
604!||--- uses -----------------------------------------------------
605!|| imp_intm ../engine/share/modules/imp_intm.F
606!||====================================================================
607 SUBROUTINE i10kfor3(JLT ,A ,V ,MS ,CAND_F ,
608 1 STIF ,ITIED ,X1 ,X2 ,X3 ,
609 2 X4 ,Y1 ,Y2 ,Y3 ,Y4 ,
610 3 Z1 ,Z2 ,Z3 ,Z4 ,NSVG ,
611 4 NX1 ,NX2 ,NX3 ,NX4 ,NY1 ,
612 5 NY2 ,NY3 ,NY4 ,NZ1 ,NZ2 ,
613 6 NZ3 ,NZ4 ,LB1 ,LB2 ,LB3 ,
614 7 LB4 ,LC1 ,LC2 ,LC3 ,LC4 ,
615 8 P1 ,P2 ,P3 ,P4 ,NIN ,
616 9 IX1 ,IX2 ,IX3 ,IX4 ,GAPV ,
617 C INDEX ,VXI ,VYI ,VZI ,MSI ,
618 L CN_LOC ,CE_LOC ,XI ,YI ,ZI ,
619 H DXI ,DYI ,DZI ,D ,SCALK )
620C-----------------------------------------------
621C M o d u l e s
622C-----------------------------------------------
623 USE imp_intm
624C-----------------------------------------------
625C I m p l i c i t T y p e s
626C-----------------------------------------------
627#include "implicit_f.inc"
628#include "comlock.inc"
629C-----------------------------------------------
630C G l o b a l P a r a m e t e r s
631C-----------------------------------------------
632#include "mvsiz_p.inc"
633C-----------------------------------------------
634C C o m m o n B l o c k s
635C-----------------------------------------------
636C-----------------------------------------------
637C D u m m y A r g u m e n t s
638C-----------------------------------------------
639 INTEGER JLT,NIN,ITIED
640 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
641 . NSVG(MVSIZ), INDEX(*),CN_LOC(*), CE_LOC(*)
642 my_real
643 . A(3,*), MS(*),X1(*),X2(*),X3(*),X4(*),
644 . Y1(*),Y2(*),Y3(*),Y4(*),Z1(*),Z2(*),Z3(*),Z4(*),
645 . CAND_F(6,*), V(3,*),D(3,*),
646 . VXI(MVSIZ),VYI(MVSIZ),VZI(MVSIZ),MSI(MVSIZ)
647 my_real
648 . NX1(MVSIZ), NX2(MVSIZ), NX3(MVSIZ), NX4(MVSIZ),
649 . NY1(MVSIZ), NY2(MVSIZ), NY3(MVSIZ), NY4(MVSIZ),
650 . NZ1(MVSIZ), NZ2(MVSIZ), NZ3(MVSIZ), NZ4(MVSIZ),
651 . LB1(MVSIZ), LB2(MVSIZ), LB3(MVSIZ), LB4(MVSIZ),
652 . LC1(MVSIZ), LC2(MVSIZ), LC3(MVSIZ), LC4(MVSIZ),
653 . P1(MVSIZ), P2(MVSIZ), P3(MVSIZ), P4(MVSIZ), STIF(MVSIZ),
654 . gapv(mvsiz),
655 . dxi(mvsiz),dyi(mvsiz), dzi(mvsiz),
656 . xi(mvsiz),yi(mvsiz),zi(mvsiz),scalk
657C-----------------------------------------------
658C L o c a l V a r i a b l e s
659C-----------------------------------------------
660 INTEGER I, J1, IG, II , K0,NBINTER,K1S,K,J,NN,JG
661 INTEGER JJ,KK,IN,IE,NSUB,IBID,NI,NS
662 my_real
663 . FXI(MVSIZ), FYI(MVSIZ), FZI(MVSIZ), FNI(MVSIZ),
664 . FX1(MVSIZ), FX2(MVSIZ), FX3(MVSIZ), FX4(MVSIZ),
665 . FY1(MVSIZ), FY2(MVSIZ), FY3(MVSIZ), FY4(MVSIZ),
666 . FZ1(MVSIZ), FZ2(MVSIZ), FZ3(MVSIZ), FZ4(MVSIZ),
667 . FT1(MVSIZ), FT2(MVSIZ),
668 . N1(MVSIZ), N2(MVSIZ), N3(MVSIZ), PENE(MVSIZ),
669 . H1(MVSIZ), H2(MVSIZ), H3(MVSIZ), H4(MVSIZ),
670 . vt1(mvsiz), vt2(mvsiz),
671 . nx(mvsiz), ny(mvsiz), nz(mvsiz),
672 . t1x(mvsiz),t1y(mvsiz),t1z(mvsiz),
673 . t2x(mvsiz),t2y(mvsiz),t2z(mvsiz),
674 . vx(mvsiz), vy(mvsiz), vz(mvsiz), vn(mvsiz),
675 . dx(mvsiz), dy(mvsiz), dz(mvsiz), dn(mvsiz),
676 . s2,d1,d2,d3,d4,a1,a2,a3,a4,la1,la2,la3,la4,h0,
677 . norminv,bid,gap2,pene2,fac,fx,fy,fz
678 my_real
679 . fxn(mvsiz), fyn(mvsiz), fzn(mvsiz),
680 . fxt(mvsiz), fyt(mvsiz), fzt(mvsiz)
681C--------------------------------------------------------
682C actualise stif
683C--------------------------------------------------------
684 DO i=1,jlt
685 gap2=gapv(i)*gapv(i)
686C
687 d1 = max(zero, gap2 - p1(i))
688 d2 = max(zero, gap2 - p2(i))
689 d3 = max(zero, gap2 - p3(i))
690 d4 = max(zero, gap2 - p4(i))
691 pene2 = max(d1,d2,d3,d4)
692 IF (pene2<=zero) stif(i) = zero
693 ENDDO
694C--------------------------------------------------------
695C CAS DES PAQUETS MIXTES
696C--------------------------------------------------------
697 DO i=1,jlt
698 d1 = sqrt(p1(i))
699 p1(i) = max(zero, gapv(i) - d1)
700C
701 d2 = sqrt(p2(i))
702 p2(i) = max(zero, gapv(i) - d2)
703C
704 d3 = sqrt(p3(i))
705 p3(i) = max(zero, gapv(i) - d3)
706C
707 d4 = sqrt(p4(i))
708 p4(i) = max(zero, gapv(i) - d4)
709 ENDDO
710C
711 DO i=1,jlt
712 IF(ix3(i)/=ix4(i))THEN
713 pene(i) = max(p1(i),p2(i),p3(i),p4(i))
714C
715 la1 = one - lb1(i) - lc1(i)
716 la2 = one - lb2(i) - lc2(i)
717 la3 = one - lb3(i) - lc3(i)
718 la4 = one - lb4(i) - lc4(i)
719C
720 h0 = fourth *
721 . (p1(i)*la1 + p2(i)*la2 + p3(i)*la3 + p4(i)*la4)
722 h1(i) = h0 + p1(i) * lb1(i) + p4(i) * lc4(i)
723 h2(i) = h0 + p2(i) * lb2(i) + p1(i) * lc1(i)
724 h3(i) = h0 + p3(i) * lb3(i) + p2(i) * lc2(i)
725 h4(i) = h0 + p4(i) * lb4(i) + p3(i) * lc3(i)
726 h0 = one/max(em20,h1(i) + h2(i) + h3(i) + h4(i))
727 h1(i) = h1(i) * h0
728 h2(i) = h2(i) * h0
729 h3(i) = h3(i) * h0
730 h4(i) = h4(i) * h0
731C
732 ELSE
733 pene(i) = p1(i)
734 n1(i) = nx1(i)
735 n2(i) = ny1(i)
736 n3(i) = nz1(i)
737 h1(i) = lb1(i)
738 h2(i) = lc1(i)
739 h3(i) = one - lb1(i) - lc1(i)
740 h4(i) = zero
741 ENDIF
742 ENDDO
743C
744 DO i=1,jlt
745C correction hourglass
746 IF(ix3(i)/=ix4(i))THEN
747 h0 = -fourth*(h1(i) - h2(i) + h3(i) - h4(i))
748 h0 = min(h0,h2(i),h4(i))
749 h0 = max(h0,-h1(i),-h3(i))
750 h1(i) = h1(i) + h0
751 h2(i) = h2(i) - h0
752 h3(i) = h3(i) + h0
753 h4(i) = h4(i) - h0
754 ENDIF
755 ENDDO
756C
757 DO i=1,jlt
758 dx(i) = dxi(i) - h1(i)*d(1,ix1(i)) - h2(i)*d(1,ix2(i))
759 . - h3(i)*d(1,ix3(i)) - h4(i)*d(1,ix4(i))
760 dy(i) = dyi(i) - h1(i)*d(2,ix1(i)) - h2(i)*d(2,ix2(i))
761 . - h3(i)*d(2,ix3(i)) - h4(i)*d(2,ix4(i))
762 dz(i) = dzi(i) - h1(i)*d(3,ix1(i)) - h2(i)*d(3,ix2(i))
763 . - h3(i)*d(3,ix3(i)) - h4(i)*d(3,ix4(i))
764 ENDDO
765C
766 DO i=1,jlt
767 ii = index(i)
768 IF(cand_f(1,ii)==zero)THEN
769C------------------------------------
770C 1ER IMPACT ou PAS d'IMPACT
771C------------------------------------
772 ELSE
773C------------------------------------
774C IMPACTS SUIVANTS
775C------------------------------------
776 h1(i) = cand_f(4,ii)
777 h2(i) = cand_f(5,ii)
778 h3(i) = cand_f(6,ii)
779 h4(i) = one - h1(i) - h2(i) - h3(i)
780 ENDIF
781 ENDDO
782C
783 DO i=1,jlt
784 t1x(i) = x3(i) - x1(i)
785 t1y(i) = y3(i) - y1(i)
786 t1z(i) = z3(i) - z1(i)
787 norminv = one/sqrt(t1x(i)**2+t1y(i)**2+t1z(i)**2)
788 t1x(i) = t1x(i)*norminv
789 t1y(i) = t1y(i)*norminv
790 t1z(i) = t1z(i)*norminv
791C
792 t2x(i) = x4(i) - x2(i)
793 t2y(i) = y4(i) - y2(i)
794 t2z(i) = z4(i) - z2(i)
795C
796 nx(i) = t1y(i)*t2z(i) - t1z(i)*t2y(i)
797 ny(i) = t1z(i)*t2x(i) - t1x(i)*t2z(i)
798 nz(i) = t1x(i)*t2y(i) - t1y(i)*t2x(i)
799 norminv = one/sqrt(nx(i)**2+ny(i)**2+nz(i)**2)
800 nx(i) = nx(i)*norminv
801 ny(i) = ny(i)*norminv
802 nz(i) = nz(i)*norminv
803C
804 t2x(i) = ny(i)*t1z(i) - nz(i)*t1y(i)
805 t2y(i) = nz(i)*t1x(i) - nx(i)*t1z(i)
806 t2z(i) = nx(i)*t1y(i) - ny(i)*t1x(i)
807C
808 dn(i) = nx(i)*dx(i) + ny(i)*dy(i) + nz(i)*dz(i)
809 vt1(i) = dx(i)*t1x(i) + dy(i)*t1y(i) + dz(i)*t1z(i)
810 vt2(i) = dx(i)*t2x(i) + dy(i)*t2y(i) + dz(i)*t2z(i)
811 ENDDO
812 fac = abs(scalk)
813 DO i=1,jlt
814 stif(i)=stif(i)*fac
815 ENDDO
816C
817 DO i=1,jlt
818 IF(pene(i)==zero.AND.cand_f(1,index(i))==zero)THEN
819C------------------------------------
820C PAS ENCORE D'IMPACT OU REBOND
821C------------------------------------
822 dn(i) = zero
823 ENDIF
824 ENDDO
825C
826 DO i=1,jlt
827 ii = index(i)
828 fni(i) = cand_f(1,ii) + dn(i) * stif(i)
829 ft1(i) = cand_f(2,ii) + vt1(i) * stif(i)
830 ft2(i) = cand_f(3,ii) + vt2(i) * stif(i)
831 ENDDO
832C
833 DO 100 i=1,jlt
834 IF(itied==0)THEN
835 IF(cand_f(1,index(i))*fni(i)<zero)THEN
836C------------------------------------
837C REBOND
838C------------------------------------
839 fni(i) = zero
840 dn(i) = zero
841 stif(i) = zero
842 ft1(i) = zero
843 ft2(i) = zero
844 ELSE
845C--------
846 ENDIF
847 ENDIF
848C
849 100 CONTINUE
850C-------------------------------------------
851 DO i=1,jlt
852 ii = index(i)
853 fxn(i)= nx(i)*fni(i)
854 fyn(i)= ny(i)*fni(i)
855 fzn(i)= nz(i)*fni(i)
856 fxt(i)= t1x(i)*ft1(i) + t2x(i)*ft2(i)
857 fyt(i)= t1y(i)*ft1(i) + t2y(i)*ft2(i)
858 fzt(i)= t1z(i)*ft1(i) + t2z(i)*ft2(i)
859 fxi(i) = fxn(i) + fxt(i)
860 fyi(i) = fyn(i) + fyt(i)
861 fzi(i) = fzn(i) + fzt(i)
862 ENDDO
863C--------main part-------
864 DO i=1,jlt
865 fx=fxi(i)
866 fy=fyi(i)
867 fz=fzi(i)
868 a(1,ix1(i))=a(1,ix1(i))+fx*h1(i)
869 a(1,ix2(i))=a(1,ix2(i))+fx*h2(i)
870 a(1,ix3(i))=a(1,ix3(i))+fx*h3(i)
871 a(1,ix4(i))=a(1,ix4(i))+fx*h4(i)
872 a(2,ix1(i))=a(2,ix1(i))+fy*h1(i)
873 a(2,ix2(i))=a(2,ix2(i))+fy*h2(i)
874 a(2,ix3(i))=a(2,ix3(i))+fy*h3(i)
875 a(2,ix4(i))=a(2,ix4(i))+fy*h4(i)
876 a(3,ix1(i))=a(3,ix1(i))+fz*h1(i)
877 a(3,ix2(i))=a(3,ix2(i))+fz*h2(i)
878 a(3,ix3(i))=a(3,ix3(i))+fz*h3(i)
879 a(3,ix4(i))=a(3,ix4(i))+fz*h4(i)
880 ENDDO
881C--------secnd part-------
882 DO i=1,jlt
883 ig=nsvg(i)
884 IF(ig>0)THEN
885 a(1,ig)=a(1,ig)-fxi(i)
886 a(2,ig)=a(2,ig)-fyi(i)
887 a(3,ig)=a(3,ig)-fzi(i)
888 ELSE
889 nn=-ig
890 ns=ind_int(nin)%P(nn)
891 ffi(1,ns)=ffi(1,ns)-fxi(i)
892 ffi(2,ns)=ffi(2,ns)-fyi(i)
893 ffi(3,ns)=ffi(3,ns)-fzi(i)
894 ENDIF
895 ENDDO
896
897C
898 RETURN
899 END
i10kfor3
subroutine i10kfor3(jlt, a, v, ms, cand_f, stif, itied, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, nsvg, nx1, nx2, nx3, nx4, ny1, ny2, ny3, ny4, nz1, nz2, nz3, nz4, lb1, lb2, lb3, lb4, lc1, lc2, lc3, lc4, p1, p2, p3, p4, nin, ix1, ix2, ix3, ix4, gapv, index, vxi, vyi, vzi, msi, cn_loc, ce_loc, xi, yi, zi, dxi, dyi, dzi, d, scalk)
Definition i10keg3.F:620
i10keg3
subroutine i10keg3(jlt, a, v, ms, fric, nx1, nx2, nx3, nx4, ny1, ny2, ny3, ny4, nz1, nz2, nz3, nz4, lb1, lb2, lb3, lb4, lc1, lc2, lc3, lc4, p1, p2, p3, p4, nin, ix1, ix2, ix3, ix4, nsvg, gapv, itied, cand_f, index, stif, vxi, vyi, vzi, msi, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, ki11, ki12, kj11, kj12, kk11, kk12, kl11, kl12, off, scalk, lrem)
Definition i10keg3.F:44
i10frf3
subroutine i10frf3(jlt, a, v, ms, fric, n1, n2, n3, t1x, t1y, t1z, h1, h2, h3, h4, ix1, ix2, ix3, ix4, index, vxi, vyi, vzi, msi, dxi, dyi, dzi, stif, nin, d, scalk)
Definition i10keg3.F:464
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
imp_intm
Definition imp_intm.F:118
imp_intm::shf_int
integer, dimension(:), allocatable shf_int
Definition imp_intm.F:136
imp_intm::intp_d
integer intp_d
Definition imp_intm.F:173
imp_intm::ind_int
type(int_pointer2), dimension(:), allocatable ind_int
Definition imp_intm.F:133