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, J, K, 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 . FAC,FACF, H0, LA1, LA2, LA3, LA4,FACT(MVSIZ),
96 . D1,D2,D3,D4,KN(4,MVSIZ),Q(3,3,MVSIZ),FAC10
97 my_real
98 . Q11,Q12,Q13,Q22,Q23,Q33,H00,
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 1st or no impact
189C------------------------------------
190C It's done in 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 ---- Without friction first ---
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 ----with friction ---
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---------for 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---------temporarily for 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, 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, ISF, 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),facn(mvsiz),facf, fact(mvsiz),fac10
499 my_real
500 . fx,fy,fz,fn,ft1,ft2,fni,fti1,fti2,
501 . t2x(mvsiz), t2y(mvsiz), t2z(mvsiz)
502C-----------------------------------------------
503C
504 fac10 = ten
505 DO i=1,jlt
506 vx(i) = vxi(i) - h1(i)*v(1,ix1(i)) - h2(i)*v(1,ix2(i))
507 . - h3(i)*v(1,ix3(i)) - h4(i)*v(1,ix4(i))
508 vy(i) = vyi(i) - h1(i)*v(2,ix1(i)) - h2(i)*v(2,ix2(i))
509 . - h3(i)*v(2,ix3(i)) - h4(i)*v(2,ix4(i))
510 vz(i) = vzi(i) - h1(i)*v(3,ix1(i)) - h2(i)*v(3,ix2(i))
511 . - h3(i)*v(3,ix3(i)) - h4(i)*v(3,ix4(i))
512 vn(i) = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
513 dx(i) = dxi(i) - h1(i)*d(1,ix1(i)) - h2(i)*d(1,ix2(i))
514 . - h3(i)*d(1,ix3(i)) - h4(i)*d(1,ix4(i))
515 dy(i) = dyi(i) - h1(i)*d(2,ix1(i)) - h2(i)*d(2,ix2(i))
516 . - h3(i)*d(2,ix3(i)) - h4(i)*d(2,ix4(i))
517 dz(i) = dzi(i) - h1(i)*d(3,ix1(i)) - h2(i)*d(3,ix2(i))
518 . - h3(i)*d(3,ix3(i)) - h4(i)*d(3,ix4(i))
519 dn(i) = n1(i)*dx(i) + n2(i)*dy(i) + n3(i)*dz(i)
520 dni(i) = n1(i)*dxi(i) + n2(i)*dyi(i) + n3(i)*dzi(i)
521 ENDDO
522C-------------------------------------------
523 DO i=1,jlt
524 t2x(i) = n2(i)*t1z(i) - n3(i)*t1y(i)
525 t2y(i) = n3(i)*t1x(i) - n1(i)*t1z(i)
526 t2z(i) = n1(i)*t1y(i) - n2(i)*t1x(i)
527 d1t(i) = t1x(i)*dx(i) + t1y(i)*dy(i) + t1z(i)*dz(i)
528 d2t(i) = t2x(i)*dx(i) + t2y(i)*dy(i) + t2z(i)*dz(i)
529 dti1(i) = t1x(i)*dxi(i) + t1y(i)*dyi(i) + t1z(i)*dzi(i)
530 dti2(i) = t2x(i)*dxi(i) + t2y(i)*dyi(i) + t2z(i)*dzi(i)
531 ENDDO
532 IF (scalk<0) THEN
533 isf=1
534 ELSE
535 isf=0
536 ENDIF
537 facf=fac10*abs(scalk)
538 IF (isf==1) THEN
539 DO i=1,jlt
540 IF (vn(i)>zero) THEN
541 facn(i)=stif(i)*facf
542 ELSEIF (vn(i)<zero) THEN
543 facn(i)=stif(i)/facf
544 ELSE
545 facn(i)=stif(i)
546 ENDIF
547 fact(i)=facn(i)*fric
548 ENDDO
549 ELSE
550 DO i=1,jlt
551 facn(i)=stif(i)*facf
552 fact(i)=facn(i)*fric
553 ENDDO
554 ENDIF
555C-------- part nml --------
556 DO i=1,jlt
557 fn = -facn(i)*dni(i)
558 fx=fn*n1(i)
559 fy=fn*n2(i)
560 fz=fn*n3(i)
561 ft1 = -fact(i)*dti1(i)
562 ft2 = -fact(i)*dti2(i)
563 fx = fx + ft1*t1x(i)+ ft2*t2x(i)
564 fy = fy + ft1*t1y(i)+ ft2*t2y(i)
565 fz = fz + ft1*t1z(i)+ ft2*t2z(i)
566 a(1,ix1(i))=a(1,ix1(i))+fx*h1(i)
567 a(1,ix2(i))=a(1,ix2(i))+fx*h2(i)
568 a(1,ix3(i))=a(1,ix3(i))+fx*h3(i)
569 a(1,ix4(i))=a(1,ix4(i))+fx*h4(i)
570 a(2,ix1(i))=a(2,ix1(i))+fy*h1(i)
571 a(2,ix2(i))=a(2,ix2(i))+fy*h2(i)
572 a(2,ix3(i))=a(2,ix3(i))+fy*h3(i)
573 a(2,ix4(i))=a(2,ix4(i))+fy*h4(i)
574 a(3,ix1(i))=a(3,ix1(i))+fz*h1(i)
575 a(3,ix2(i))=a(3,ix2(i))+fz*h2(i)
576 a(3,ix3(i))=a(3,ix3(i))+fz*h3(i)
577 a(3,ix4(i))=a(3,ix4(i))+fz*h4(i)
578 ENDDO
579C-------- part nsl --------
580 DO i=1,jlt
581 fni = facn(i)*dn(i)
582 fx=fni*n1(i)
583 fy=fni*n2(i)
584 fz=fni*n3(i)
585 fti1 = fact(i)*d1t(i)
586 fti2 = fact(i)*d2t(i)
587 fx = fx + fti1*t1x(i)+ fti2*t2x(i)
588 fy = fy + fti1*t1y(i)+ fti2*t2y(i)
589 fz = fz + fti1*t1z(i)+ fti2*t2z(i)
590 ni = index(i)
591 ffi(1,ni)=ffi(1,ni)+fx
592 ffi(2,ni)=ffi(2,ni)+fy
593 ffi(3,ni)=ffi(3,ni)+fz
594 ENDDO
595C
596 RETURN
597 END
598!||====================================================================
599!|| i10kfor3 ../engine/source/interfaces/int10/i10keg3.F
600!||--- called by ------------------------------------------------------
601!|| i10fku3 ../engine/source/interfaces/int10/i10ke3.F
602!||--- uses -----------------------------------------------------
603!|| imp_intm ../engine/share/modules/imp_intm.F
604!||====================================================================
605 SUBROUTINE i10kfor3(JLT ,A ,V ,MS ,CAND_F ,
606 1 STIF ,ITIED ,X1 ,X2 ,X3 ,
607 2 X4 ,Y1 ,Y2 ,Y3 ,Y4 ,
608 3 Z1 ,Z2 ,Z3 ,Z4 ,NSVG ,
609 4 NX1 ,NX2 ,NX3 ,NX4 ,NY1 ,
610 5 NY2 ,NY3 ,NY4 ,NZ1 ,NZ2 ,
611 6 NZ3 ,NZ4 ,LB1 ,LB2 ,LB3 ,
612 7 LB4 ,LC1 ,LC2 ,LC3 ,LC4 ,
613 8 P1 ,P2 ,P3 ,P4 ,NIN ,
614 9 IX1 ,IX2 ,IX3 ,IX4 ,GAPV ,
615 C INDEX ,VXI ,VYI ,VZI ,MSI ,
616 L CN_LOC ,CE_LOC ,XI ,YI ,ZI ,
617 H DXI ,DYI ,DZI ,D ,SCALK )
618C-----------------------------------------------
619C M o d u l e s
620C-----------------------------------------------
621 USE imp_intm
622C-----------------------------------------------
623C I m p l i c i t T y p e s
624C-----------------------------------------------
625#include "implicit_f.inc"
626#include "comlock.inc"
627C-----------------------------------------------
628C G l o b a l P a r a m e t e r s
629C-----------------------------------------------
630#include "mvsiz_p.inc"
631C-----------------------------------------------
632C C o m m o n B l o c k s
633C-----------------------------------------------
634C-----------------------------------------------
635C D u m m y A r g u m e n t s
636C-----------------------------------------------
637 INTEGER JLT,NIN,ITIED
638 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
639 . NSVG(MVSIZ), INDEX(*),CN_LOC(*), CE_LOC(*)
640 my_real
641 . A(3,*), MS(*),X1(*),X2(*),X3(*),X4(*),
642 . Y1(*),Y2(*),Y3(*),Y4(*),Z1(*),Z2(*),Z3(*),Z4(*),
643 . CAND_F(6,*), V(3,*),D(3,*),
644 . VXI(MVSIZ),VYI(MVSIZ),VZI(MVSIZ),MSI(MVSIZ)
645 my_real
646 . NX1(MVSIZ), NX2(MVSIZ), NX3(MVSIZ), NX4(MVSIZ),
647 . NY1(MVSIZ), NY2(MVSIZ), NY3(MVSIZ), NY4(MVSIZ),
648 . NZ1(MVSIZ), NZ2(MVSIZ), NZ3(MVSIZ), NZ4(MVSIZ),
649 . LB1(MVSIZ), LB2(MVSIZ), LB3(MVSIZ), LB4(MVSIZ),
650 . LC1(MVSIZ), LC2(MVSIZ), LC3(MVSIZ), LC4(MVSIZ),
651 . P1(MVSIZ), P2(MVSIZ), P3(MVSIZ), P4(MVSIZ), STIF(MVSIZ),
652 . gapv(mvsiz),
653 . dxi(mvsiz),dyi(mvsiz), dzi(mvsiz),
654 . xi(mvsiz),yi(mvsiz),zi(mvsiz),scalk
655C-----------------------------------------------
656C L o c a l V a r i a b l e s
657C-----------------------------------------------
658 INTEGER I, IG, II, NN
659 INTEGER NS
660 my_real
661 . FXI(MVSIZ), FYI(MVSIZ), FZI(MVSIZ), FNI(MVSIZ),
662 .
663 .
664 .
665 . FT1(MVSIZ), FT2(MVSIZ),
666 . N1(MVSIZ), N2(MVSIZ), N3(MVSIZ), PENE(MVSIZ),
667 . H1(MVSIZ), H2(MVSIZ), H3(MVSIZ), H4(MVSIZ),
668 . vt1(mvsiz), vt2(mvsiz),
669 . nx(mvsiz), ny(mvsiz), nz(mvsiz),
670 . t1x(mvsiz),t1y(mvsiz),t1z(mvsiz),
671 . t2x(mvsiz),t2y(mvsiz),t2z(mvsiz),
672 .
673 . dx(mvsiz), dy(mvsiz), dz(mvsiz), dn(mvsiz),
674 . d1,d2,d3,d4,la1,la2,la3,la4,h0,
675 . norminv,gap2,pene2,fac,fx,fy,fz
676 my_real
677 . fxn(mvsiz), fyn(mvsiz), fzn(mvsiz),
678 . fxt(mvsiz), fyt(mvsiz), fzt(mvsiz)
679C--------------------------------------------------------
680C actualise stif
681C--------------------------------------------------------
682 DO i=1,jlt
683 gap2=gapv(i)*gapv(i)
684C
685 d1 = max(zero, gap2 - p1(i))
686 d2 = max(zero, gap2 - p2(i))
687 d3 = max(zero, gap2 - p3(i))
688 d4 = max(zero, gap2 - p4(i))
689 pene2 = max(d1,d2,d3,d4)
690 IF (pene2<=zero) stif(i) = zero
691 ENDDO
692C--------------------------------------------------------
693C CASE OF MIXED GROUPS
694C--------------------------------------------------------
695 DO i=1,jlt
696 d1 = sqrt(p1(i))
697 p1(i) = max(zero, gapv(i) - d1)
698C
699 d2 = sqrt(p2(i))
700 p2(i) = max(zero, gapv(i) - d2)
701C
702 d3 = sqrt(p3(i))
703 p3(i) = max(zero, gapv(i) - d3)
704C
705 d4 = sqrt(p4(i))
706 p4(i) = max(zero, gapv(i) - d4)
707 ENDDO
708C
709 DO i=1,jlt
710 IF(ix3(i)/=ix4(i))THEN
711 pene(i) = max(p1(i),p2(i),p3(i),p4(i))
712C
713 la1 = one - lb1(i) - lc1(i)
714 la2 = one - lb2(i) - lc2(i)
715 la3 = one - lb3(i) - lc3(i)
716 la4 = one - lb4(i) - lc4(i)
717C
718 h0 = fourth *
719 . (p1(i)*la1 + p2(i)*la2 + p3(i)*la3 + p4(i)*la4)
720 h1(i) = h0 + p1(i) * lb1(i) + p4(i) * lc4(i)
721 h2(i) = h0 + p2(i) * lb2(i) + p1(i) * lc1(i)
722 h3(i) = h0 + p3(i) * lb3(i) + p2(i) * lc2(i)
723 h4(i) = h0 + p4(i) * lb4(i) + p3(i) * lc3(i)
724 h0 = one/max(em20,h1(i) + h2(i) + h3(i) + h4(i))
725 h1(i) = h1(i) * h0
726 h2(i) = h2(i) * h0
727 h3(i) = h3(i) * h0
728 h4(i) = h4(i) * h0
729C
730 ELSE
731 pene(i) = p1(i)
732 n1(i) = nx1(i)
733 n2(i) = ny1(i)
734 n3(i) = nz1(i)
735 h1(i) = lb1(i)
736 h2(i) = lc1(i)
737 h3(i) = one - lb1(i) - lc1(i)
738 h4(i) = zero
739 ENDIF
740 ENDDO
741C
742 DO i=1,jlt
743C correction hourglass
744 IF(ix3(i)/=ix4(i))THEN
745 h0 = -fourth*(h1(i) - h2(i) + h3(i) - h4(i))
746 h0 = min(h0,h2(i),h4(i))
747 h0 = max(h0,-h1(i),-h3(i))
748 h1(i) = h1(i) + h0
749 h2(i) = h2(i) - h0
750 h3(i) = h3(i) + h0
751 h4(i) = h4(i) - h0
752 ENDIF
753 ENDDO
754C
755 DO i=1,jlt
756 dx(i) = dxi(i) - h1(i)*d(1,ix1(i)) - h2(i)*d(1,ix2(i))
757 . - h3(i)*d(1,ix3(i)) - h4(i)*d(1,ix4(i))
758 dy(i) = dyi(i) - h1(i)*d(2,ix1(i)) - h2(i)*d(2,ix2(i))
759 . - h3(i)*d(2,ix3(i)) - h4(i)*d(2,ix4(i))
760 dz(i) = dzi(i) - h1(i)*d(3,ix1(i)) - h2(i)*d(3,ix2(i))
761 . - h3(i)*d(3,ix3(i)) - h4(i)*d(3,ix4(i))
762 ENDDO
763C
764 DO i=1,jlt
765 ii = index(i)
766 IF(cand_f(1,ii)==zero)THEN
767C------------------------------------
768C 1st or no impact
769C------------------------------------
770 ELSE
771C------------------------------------
772C IMPACTS SUIVANTS
773C------------------------------------
774 h1(i) = cand_f(4,ii)
775 h2(i) = cand_f(5,ii)
776 h3(i) = cand_f(6,ii)
777 h4(i) = one - h1(i) - h2(i) - h3(i)
778 ENDIF
779 ENDDO
780C
781 DO i=1,jlt
782 t1x(i) = x3(i) - x1(i)
783 t1y(i) = y3(i) - y1(i)
784 t1z(i) = z3(i) - z1(i)
785 norminv = one/sqrt(t1x(i)**2+t1y(i)**2+t1z(i)**2)
786 t1x(i) = t1x(i)*norminv
787 t1y(i) = t1y(i)*norminv
788 t1z(i) = t1z(i)*norminv
789C
790 t2x(i) = x4(i) - x2(i)
791 t2y(i) = y4(i) - y2(i)
792 t2z(i) = z4(i) - z2(i)
793C
794 nx(i) = t1y(i)*t2z(i) - t1z(i)*t2y(i)
795 ny(i) = t1z(i)*t2x(i) - t1x(i)*t2z(i)
796 nz(i) = t1x(i)*t2y(i) - t1y(i)*t2x(i)
797 norminv = one/sqrt(nx(i)**2+ny(i)**2+nz(i)**2)
798 nx(i) = nx(i)*norminv
799 ny(i) = ny(i)*norminv
800 nz(i) = nz(i)*norminv
801C
802 t2x(i) = ny(i)*t1z(i) - nz(i)*t1y(i)
803 t2y(i) = nz(i)*t1x(i) - nx(i)*t1z(i)
804 t2z(i) = nx(i)*t1y(i) - ny(i)*t1x(i)
805C
806 dn(i) = nx(i)*dx(i) + ny(i)*dy(i) + nz(i)*dz(i)
807 vt1(i) = dx(i)*t1x(i) + dy(i)*t1y(i) + dz(i)*t1z(i)
808 vt2(i) = dx(i)*t2x(i) + dy(i)*t2y(i) + dz(i)*t2z(i)
809 ENDDO
810 fac = abs(scalk)
811 DO i=1,jlt
812 stif(i)=stif(i)*fac
813 ENDDO
814C
815 DO i=1,jlt
816 IF(pene(i)==zero.AND.cand_f(1,index(i))==zero)THEN
817C------------------------------------
818C PAS ENCORE D'IMPACT OU REBOND
819C------------------------------------
820 dn(i) = zero
821 ENDIF
822 ENDDO
823C
824 DO i=1,jlt
825 ii = index(i)
826 fni(i) = cand_f(1,ii) + dn(i) * stif(i)
827 ft1(i) = cand_f(2,ii) + vt1(i) * stif(i)
828 ft2(i) = cand_f(3,ii) + vt2(i) * stif(i)
829 ENDDO
830C
831 DO 100 i=1,jlt
832 IF(itied==0)THEN
833 IF(cand_f(1,index(i))*fni(i)<zero)THEN
834C------------------------------------
835C REBOND
836C------------------------------------
837 fni(i) = zero
838 dn(i) = zero
839 stif(i) = zero
840 ft1(i) = zero
841 ft2(i) = zero
842 ELSE
843C--------
844 ENDIF
845 ENDIF
846C
847 100 CONTINUE
848C-------------------------------------------
849 DO i=1,jlt
850 ii = index(i)
851 fxn(i)= nx(i)*fni(i)
852 fyn(i)= ny(i)*fni(i)
853 fzn(i)= nz(i)*fni(i)
854 fxt(i)= t1x(i)*ft1(i) + t2x(i)*ft2(i)
855 fyt(i)= t1y(i)*ft1(i) + t2y(i)*ft2(i)
856 fzt(i)= t1z(i)*ft1(i) + t2z(i)*ft2(i)
857 fxi(i) = fxn(i) + fxt(i)
858 fyi(i) = fyn(i) + fyt(i)
859 fzi(i) = fzn(i) + fzt(i)
860 ENDDO
861C--------main part-------
862 DO i=1,jlt
863 fx=fxi(i)
864 fy=fyi(i)
865 fz=fzi(i)
866 a(1,ix1(i))=a(1,ix1(i))+fx*h1(i)
867 a(1,ix2(i))=a(1,ix2(i))+fx*h2(i)
868 a(1,ix3(i))=a(1,ix3(i))+fx*h3(i)
869 a(1,ix4(i))=a(1,ix4(i))+fx*h4(i)
870 a(2,ix1(i))=a(2,ix1(i))+fy*h1(i)
871 a(2,ix2(i))=a(2,ix2(i))+fy*h2(i)
872 a(2,ix3(i))=a(2,ix3(i))+fy*h3(i)
873 a(2,ix4(i))=a(2,ix4(i))+fy*h4(i)
874 a(3,ix1(i))=a(3,ix1(i))+fz*h1(i)
875 a(3,ix2(i))=a(3,ix2(i))+fz*h2(i)
876 a(3,ix3(i))=a(3,ix3(i))+fz*h3(i)
877 a(3,ix4(i))=a(3,ix4(i))+fz*h4(i)
878 ENDDO
879C------- Secondary part --------
880 DO i=1,jlt
881 ig=nsvg(i)
882 IF(ig>0)THEN
883 a(1,ig)=a(1,ig)-fxi(i)
884 a(2,ig)=a(2,ig)-fyi(i)
885 a(3,ig)=a(3,ig)-fzi(i)
886 ELSE
887 nn=-ig
888 ns=ind_int(nin)%P(nn)
889 ffi(1,ns)=ffi(1,ns)-fxi(i)
890 ffi(2,ns)=ffi(2,ns)-fyi(i)
891 ffi(3,ns)=ffi(3,ns)-fzi(i)
892 ENDIF
893 ENDDO
894
895C
896 RETURN
897 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:618
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