OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i2for25.F File Reference
#include "implicit_f.inc"
#include "comlock.inc"
#include "mvsiz_p.inc"
#include "com06_c.inc"
#include "com08_c.inc"
#include "scr11_c.inc"
#include "scr14_c.inc"
#include "sms_c.inc"
#include "task_c.inc"
#include "lockon.inc"
#include "lockoff.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine i2for25 (x, v, vr, a, ar, ms, stifn, stifr, weight, irect, nsv, irtl, crst, skew, xini, dx, fini, fsav, fncont, nsn, stfn, stfr, visc, penflag, irot, noint, nodnx_sms, dmint2, sav_for_pena, ms_pena, dt2t, neltst, ityptst, ivisc, h3d_data, fncontp, ftcontp)

Function/Subroutine Documentation

◆ i2for25()

subroutine i2for25 ( x,
v,
vr,
a,
ar,
ms,
stifn,
stifr,
integer, dimension(*) weight,
integer, dimension(4,*) irect,
integer, dimension(*) nsv,
integer, dimension(*) irtl,
crst,
skew,
xini,
dx,
fini,
fsav,
fncont,
integer nsn,
stfn,
stfr,
visc,
integer penflag,
integer irot,
integer noint,
integer, dimension(*) nodnx_sms,
dmint2,
sav_for_pena,
ms_pena,
dt2t,
integer neltst,
integer ityptst,
integer ivisc,
type (h3d_database) h3d_data,
fncontp,
ftcontp )

Definition at line 36 of file i2for25.F.

44C-----------------------------------------------
45C M o d u l e s
46C-----------------------------------------------
47 USE h3d_mod
48C-----------------------------------------------
49C I m p l i c i t T y p e s
50C-----------------------------------------------
51#include "implicit_f.inc"
52#include "comlock.inc"
53C-----------------------------------------------
54C G l o b a l P a r a m e t e r s
55C-----------------------------------------------
56#include "mvsiz_p.inc"
57C-----------------------------------------------
58C D u m m y A r g u m e n t s
59C-----------------------------------------------
60 INTEGER NSN,PENFLAG,IROT, NOINT,NELTST,ITYPTST
61 INTEGER IRECT(4,*),NSV(*),IRTL(*),WEIGHT(*),
62 . NODNX_SMS(*),IVISC
63C REAL
64 my_real
65 . visc,dt2t
66 my_real
67 . x(3,*),vr(3,*),v(3,*),a(3,*),ar(3,*),xini(3,*),skew(9,*),
68 . dx(3,*),fini(3,*),ms(*),stifn(*),stifr(*),stfn(*),stfr(*),
69 . crst(2,*),fsav(*),fncont(3,*),fncontp(3,*) ,ftcontp(3,*),
70 . dmint2(4,*),sav_for_pena(4,*),ms_pena(*)
71 TYPE (H3D_DATABASE) :: H3D_DATA
72C-----------------------------------------------
73C C o m m o n B l o c k s
74C-----------------------------------------------
75#include "com06_c.inc"
76#include "com08_c.inc"
77#include "scr11_c.inc"
78#include "scr14_c.inc"
79#include "sms_c.inc"
80#include "task_c.inc"
81C-----------------------------------------------
82C L o c a l V a r i a b l e s
83C-----------------------------------------------
84 INTEGER NIR,I,J,II,JJ,L,W,KK,K,LLT,
85 . IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
86 . NSVG(MVSIZ)
87C REAL
88 my_real
89 . S,T,SP,SM,TP,TM,ECONTT,ECONVT,E1X,E1Y,E1Z,E2X,E2Y,E2Z,E3X,E3Y,E3Z,
90 . FNORM,FLX,FLY,FLZ,FS(3),DDX,DDY,DDZ,XSM,YSM,ZSM,XM,YM,ZM,
91 . X1,X2,X3,X4,Y1,Y2,Y3,Y4,Z1,Z2,Z3,Z4,X0,Y0,Z0,XS,YS,ZS,STF,
92 . VX1,VX2,VX3,VX4,VY1,VY2,VY3,VY4,VZ1,VZ2,VZ3,VZ4,DLX,DLY,DLZ,
93 . VX0,VY0,VZ0,VSX,VSY,VSZ,VMX,VMY,VMZ,V1,V2,V3,DTINV,
94 . FXV,FYV,FZV,DRX,DRY,DRZ,STBRK,DTI, DXT
95 my_real
96 . H(4,MVSIZ),FN(3),FT(3),FX(4),FY(4),FZ(4),FMX(4),FMY(4),FMZ(4),
97 . RX(4),RY(4),RZ(4),RM(3),RS(3),V0(3),VS(3),VM(3),DXOLD(3),
98 . STIF(MVSIZ), VIS(MVSIZ), STIFM(MVSIZ),HL(4)
99 my_real
100 . MHARM,DKM,VA(3),VB(3),VC(3),VD(3),VISCL
101
102C=======================================================================
103 I7KGLO = 1
104 ECONTT = ZERO
105 ECONVT = ZERO
106C----------------
107 DO KK=1,NSN,MVSIZ
108C
109 LLT=MIN(NSN-KK+1,MVSIZ)
110 DO K=1,LLT
111C
112 II= KK+K-1
113 I = NSV(II)
114C
115 IF (I > 0) THEN
116 NSVG(K) = I
117 W = WEIGHT(I)
118 S = CRST(1,II)
119 T = CRST(2,II)
120 L = IRTL(II)
121C
122 IX1(K) = IRECT(1,L)
123 IX2(K) = IRECT(2,L)
124 IX3(K) = IRECT(3,L)
125 IX4(K) = IRECT(4,L)
126 NIR= 4
127 SP = ONE + S
128 SM = ONE - S
129 TP = FOURTH*(ONE + T)
130 TM = FOURTH*(ONE - T)
131 H(1,K)=TM*SM
132 H(2,K)=TM*SP
133 H(3,K)=TP*SP
134 H(4,K)=TP*SM
135 IF (IX3(K) == IX4(K)) THEN
136 NIR = 3
137 H(3,K) = H(3,K) + H(4,K)
138 H(4,K) = ZERO
139 ENDIF
140C------------------------------------------------
141C rep local facette main
142C------------------------------------------------
143 X1 = X(1,IX1(K))
144 Y1 = X(2,IX1(K))
145 Z1 = X(3,IX1(K))
146 X2 = X(1,IX2(K))
147 Y2 = X(2,IX2(K))
148 Z2 = X(3,IX2(K))
149 X3 = X(1,IX3(K))
150 Y3 = X(2,IX3(K))
151 Z3 = X(3,IX3(K))
152 X4 = X(1,IX4(K))
153 Y4 = X(2,IX4(K))
154 Z4 = X(3,IX4(K))
155 XS = X(1,I)
156 YS = X(2,I)
157 ZS = X(3,I)
158 VSX = V(1,I)
159 VSY = V(2,I)
160 VSZ = V(3,I)
161 VX1 = V(1,IX1(K))
162 VY1 = V(2,IX1(K))
163 VZ1 = V(3,IX1(K))
164 VX2 = V(1,IX2(K))
165 VY2 = V(2,IX2(K))
166 VZ2 = V(3,IX2(K))
167 VX3 = V(1,IX3(K))
168 VY3 = V(2,IX3(K))
169 VZ3 = V(3,IX3(K))
170 VX4 = V(1,IX4(K))
171 VY4 = V(2,IX4(K))
172 VZ4 = V(3,IX4(K))
173C---------------------
174 CALL I2REP(X1 ,X2 ,X3 ,X4 ,
175 . Y1 ,Y2 ,Y3 ,Y4 ,
176 . Z1 ,Z2 ,Z3 ,Z4 ,
177 . E1X ,E1Y ,E1Z ,
178 . E2X ,E2Y ,E2Z ,
179 . E3X ,E3Y ,E3Z ,NIR )
180C------------------------------------------------
181 IF (NIR == 4) THEN
182 XM = X1*H(1,K) + X2*H(2,K) + X3*H(3,K) + X4*H(4,K)
183 YM = Y1*H(1,K) + Y2*H(2,K) + Y3*H(3,K) + Y4*H(4,K)
184 ZM = Z1*H(1,K) + Z2*H(2,K) + Z3*H(3,K) + Z4*H(4,K)
185 X0 = (X1 + X2 + X3 + X4)/NIR
186 Y0 = (Y1 + Y2 + Y3 + Y4)/NIR
187 Z0 = (Z1 + Z2 + Z3 + Z4)/NIR
188
189 XM = XM - X0
190 YM = YM - Y0
191 ZM = ZM - Z0
192 XS = XS - X0
193 YS = YS - Y0
194 ZS = ZS - Z0
195 XSM = XS - XM
196 YSM = YS - YM
197 ZSM = ZS - ZM
198c
199 VX0 = (VX1 + VX2 + VX3 + VX4)/NIR
200 VY0 = (VY1 + VY2 + VY3 + VY4)/NIR
201 VZ0 = (VZ1 + VZ2 + VZ3 + VZ4)/NIR
202 VMX = VX1*H(1,K) + VX2*H(2,K) + VX3*H(3,K) + VX4*H(4,K) - VX0
203 VMY = VY1*H(1,K) + VY2*H(2,K) + VY3*H(3,K) + VY4*H(4,K) - VY0
204 VMZ = VZ1*H(1,K) + VZ2*H(2,K) + VZ3*H(3,K) + VZ4*H(4,K) - VZ0
205 ELSE
206 X0 = (X1 + X2 + X3)/NIR
207 Y0 = (Y1 + Y2 + Y3)/NIR
208 Z0 = (Z1 + Z2 + Z3)/NIR
209
210 XM = X1*H(1,K) + X2*H(2,K) + X3*H(3,K)
211 YM = Y1*H(1,K) + Y2*H(2,K) + Y3*H(3,K)
212 ZM = Z1*H(1,K) + Z2*H(2,K) + Z3*H(3,K)
213
214 XM = XM - X0
215 YM = YM - Y0
216 ZM = ZM - Z0
217 XS = XS - X0
218 YS = YS - Y0
219 ZS = ZS - Z0
220 XSM = XS - XM
221 YSM = YS - YM
222 ZSM = ZS - ZM
223
224 VX0 = (VX1 + VX2 + VX3)/NIR
225 VY0 = (VY1 + VY2 + VY3)/NIR
226 VZ0 = (VZ1 + VZ2 + VZ3)/NIR
227 VMX = VX1*H(1,K) + VX2*H(2,K) + VX3*H(3,K) - VX0
228 VMY = VY1*H(1,K) + VY2*H(2,K) + VY3*H(3,K) - VY0
229 VMZ = VZ1*H(1,K) + VZ2*H(2,K) + VZ3*H(3,K) - VZ0
230 ENDIF
231 X1 = X1 - X0
232 Y1 = Y1 - Y0
233 Z1 = Z1 - Z0
234 X2 = X2 - X0
235 Y2 = Y2 - Y0
236 Z2 = Z2 - Z0
237 X3 = X3 - X0
238 Y3 = Y3 - Y0
239 Z3 = Z3 - Z0
240 X4 = X4 - X0
241 Y4 = Y4 - Y0
242 Z4 = Z4 - Z0
243 VSX = VSX - VX0
244 VSY = VSY - VY0
245 VSZ = VSZ - VZ0
246C
247c global -> local
248c
249 RS(1) = XS*E1X + YS*E1Y + ZS*E1Z
250 RS(2) = XS*E2X + YS*E2Y + ZS*E2Z
251 RS(3) = XS*E3X + YS*E3Y + ZS*E3Z
252 RM(1) = XM*E1X + YM*E1Y + ZM*E1Z
253 RM(2) = XM*E2X + YM*E2Y + ZM*E2Z
254 RM(3) = XM*E3X + YM*E3Y + ZM*E3Z
255c
256 RX(1) = E1X*X1 + E1Y*Y1 + E1Z*Z1
257 RY(1) = E2X*X1 + E2Y*Y1 + E2Z*Z1
258 RZ(1) = E3X*X1 + E3Y*Y1 + E3Z*Z1
259 RX(2) = E1X*X2 + E1Y*Y2 + E1Z*Z2
260 RY(2) = E2X*X2 + E2Y*Y2 + E2Z*Z2
261 RZ(2) = E3X*X2 + E3Y*Y2 + E3Z*Z2
262 RX(3) = E1X*X3 + E1Y*Y3 + E1Z*Z3
263 RY(3) = E2X*X3 + E2Y*Y3 + E2Z*Z3
264 RZ(3) = E3X*X3 + E3Y*Y3 + E3Z*Z3
265 RX(4) = E1X*X4 + E1Y*Y4 + E1Z*Z4
266 RY(4) = E2X*X4 + E2Y*Y4 + E2Z*Z4
267 RZ(4) = E3X*X4 + E3Y*Y4 + E3Z*Z4
268C
269 IF(NIR==3)THEN
270 RX(4)=ZERO
271 RY(4)=ZERO
272 RZ(4)=ZERO
273 END IF
274C
275 VS(1) = VSX*E1X + VSY*E1Y + VSZ*E1Z
276 VS(2) = VSX*E2X + VSY*E2Y + VSZ*E2Z
277 VS(3) = VSX*E3X + VSY*E3Y + VSZ*E3Z
278 VM(1) = VMX*E1X + VMY*E1Y + VMZ*E1Z
279 VM(2) = VMX*E2X + VMY*E2Y + VMZ*E2Z
280 VM(3) = VMX*E3X + VMY*E3Y + VMZ*E3Z
281C
282 VA(1) = VX1*E1X + VY1*E1Y + VZ1*E1Z
283 VA(2) = VX1*E2X + VY1*E2Y + VZ1*E2Z
284 VA(3) = VX1*E3X + VY1*E3Y + VZ1*E3Z
285
286 VB(1) = VX2*E1X + VY2*E1Y + VZ2*E1Z
287 VB(2) = VX2*E2X + VY2*E2Y + VZ2*E2Z
288 VB(3) = VX2*E3X + VY2*E3Y + VZ2*E3Z
289
290 VC(1) = VX3*E1X + VY3*E1Y + VZ3*E1Z
291 VC(2) = VX3*E2X + VY3*E2Y + VZ3*E2Z
292 VC(3) = VX3*E3X + VY3*E3Y + VZ3*E3Z
293
294 VD(1) = VX4*E1X + VY4*E1Y + VZ4*E1Z
295 VD(2) = VX4*E2X + VY4*E2Y + VZ4*E2Z
296 VD(3) = VX4*E3X + VY4*E3Y + VZ4*E3Z
297C
298 V1 = VS(1) - VM(1)
299 V2 = VS(2) - VM(2)
300 V3 = VS(3) - VM(3)
301C
302C--------- Local displacement
303 IF (TT == ZERO) THEN
304 DX(1,II) = ZERO
305 DX(2,II) = ZERO
306 DX(3,II) = ZERO
307 FINI(1,II) = ZERO
308 FINI(2,II) = ZERO
309 FINI(3,II) = ZERO
310 ENDIF
311C--------- Vi = Vi -VR ^ MS
312 CALL I2PEN_ROT(SKEW(1,II) ,TT ,DT1 ,STBRK,
313 . RS ,RM ,V1 ,V2 ,V3 ,
314 . RX ,RY ,RZ ,VA ,VB ,
315 . VC ,VD)
316C------------- vers increm en vitesses
317 DLX = V1*DT1
318 DLY = V2*DT1
319 DLZ = V3*DT1
320C------------- DX == deplacement relatif
321 DX(1,II) = DX(1,II) + DLX
322 DX(2,II) = DX(2,II) + DLY
323 DX(3,II) = DX(3,II) + DLZ
324C------------------------------------------------
325C Total force
326C------------------------------------------------
327 FLX = DX(1,II) * STFN(II)
328 FLY = DX(2,II) * STFN(II)
329 FLZ = DX(3,II) * STFN(II)
330 VISCL = VISC
331C
332 IF (IVISC==1) THEN
333C-- Old visc formulation from Radioss V14 --
334 MHARM = MS_PENA(I)
335.OR..OR. ELSEIF(MS_PENA(I)==ZEROMS_PENA(IX1(K))==ZERO
336.OR. . MS_PENA(IX2(K))==ZERO
337.OR. . MS_PENA(IX3(K))==ZERO
338 . MS_PENA(IX4(K))==ZERO)THEN
339 MHARM = ZERO
340 VISCL = ZERO
341 ELSEIF(NIR==4)THEN
342 MHARM = ONE/MS_PENA(I) +
343 . ONE/MS_PENA(IX1(K)) + ONE/MS_PENA(IX2(K)) + ONE/MS_PENA(IX3(K)) + ONE/MS_PENA(IX4(K))
344 MHARM = ONE/MHARM
345 ELSE
346 MHARM = ONE/MS_PENA(I) +
347 . ONE/MS_PENA(IX1(K)) + ONE/MS_PENA(IX2(K)) + ONE/MS_PENA(IX3(K))
348 MHARM = ONE/MHARM
349 END IF
350 DKM = TWO*STFN(II)*MHARM
351 VIS(K) = VISC*SQRT(DKM)
352C
353 FXV = VIS(K) * V1
354 FYV = VIS(K) * V2
355 FZV = VIS(K) * V3
356c
357 DXT = DX(1,II)*DX(1,II) + DX(2,II)*DX(2,II) + DX(3,II)*DX(3,II)
358 ECONTT = ECONTT + HALF*STFN(II)*DXT*W
359
360 ECONVT = ECONVT + (FXV*V1 + FYV*V2 + FZV*V3)*DT1*W
361c
362 FLX = FLX + FXV
363 FLY = FLY + FYV
364 FLZ = FLZ + FZV
365C
366 DO J=1,4
367 FMX(J) = H(J,K)*FLX
368 FMY(J) = H(J,K)*FLY
369 FMZ(J) = H(J,K)*FLZ
370 ENDDO
371C----------------------------------------------------
372c update main forces (moment balance)
373 CALL I2LOCEQ( NIR ,RS ,RX ,RY ,RZ ,
374 . FMX ,FMY ,FMZ ,H(1,K) ,STIFM(K))
375C----------------------------------------------------
376C Secnd forces -> global coordinates
377C----------------------------------------------------
378 DO J=1,4
379 FX(J) = E1X*FMX(J) + E2X*FMY(J) + E3X*FMZ(J)
380 FY(J) = E1Y*FMX(J) + E2Y*FMY(J) + E3Y*FMZ(J)
381 FZ(J) = E1Z*FMX(J) + E2Z*FMY(J) + E3Z*FMZ(J)
382 ENDDO
383 FS(1) = ZERO
384 FS(2) = ZERO
385 FS(3) = ZERO
386 DO J=1,NIR
387 FS(1) = FS(1) + FX(J)
388 FS(2) = FS(2) + FY(J)
389 FS(3) = FS(3) + FZ(J)
390 ENDDO
391c A(1,I) = A(1,I) - FS(1)
392c A(2,I) = A(2,I) - FS(2)
393c A(3,I) = A(3,I) - FS(3)
394 SAV_FOR_PENA(1,I)= SAV_FOR_PENA(1,I)-FS(1)
395 SAV_FOR_PENA(2,I)= SAV_FOR_PENA(2,I)-FS(2)
396 SAV_FOR_PENA(3,I)= SAV_FOR_PENA(3,I)-FS(3)
397C
398 IF (IVISC==1) THEN
399C-- Old visc formulation from Radioss V14 --
400 DTINV = ZERO
401 IF (DT1 > ZERO) DTINV=ONE/DT1
402 STF = (ONE+STBRK)*STFN(II) + TWO*VIS(K)*DTINV
403 ELSE
404 STF = STFN(II)*(VISCL + SQRT(VISCL**2 + (ONE+STBRK)))**2
405 ENDIF
406C
407 SAV_FOR_PENA(4,I)= SAV_FOR_PENA(4,I)+STF
408 STIF(K) = (ONE+STBRK)*STFN(II)
409C----------------------------------------------------
410C Main forces
411C----------------------------------------------------
412 IF (W == 1) THEN
413 A(1,IX1(K)) = A(1,IX1(K)) + FX(1)
414 A(2,IX1(K)) = A(2,IX1(K)) + FY(1)
415 A(3,IX1(K)) = A(3,IX1(K)) + FZ(1)
416 STIFN(IX1(K)) = STIFN(IX1(K))+ABS(STF*H(1,K))+STIFM(K)*STF
417c
418 A(1,IX2(K)) = A(1,IX2(K)) + FX(2)
419 A(2,IX2(K)) = A(2,IX2(K)) + FY(2)
420 A(3,IX2(K)) = A(3,IX2(K)) + FZ(2)
421 STIFN(IX2(K)) = STIFN(IX2(K))+ABS(STF*H(2,K))+STIFM(K)*STF
422c
423 A(1,IX3(K)) = A(1,IX3(K)) + FX(3)
424 A(2,IX3(K)) = A(2,IX3(K)) + FY(3)
425 A(3,IX3(K)) = A(3,IX3(K)) + FZ(3)
426 STIFN(IX3(K)) = STIFN(IX3(K))+ABS(STF*H(3,K))+STIFM(K)*STF
427c
428 A(1,IX4(K)) = A(1,IX4(K)) + FX(4)
429 A(2,IX4(K)) = A(2,IX4(K)) + FY(4)
430 A(3,IX4(K)) = A(3,IX4(K)) + FZ(4)
431 STIFN(IX4(K)) = STIFN(IX4(K))+ABS(STF*H(4,K))
432 . +STIFM(K)*STF*SIGN(ONE,ABS(H(4,K)))
433 ENDIF
434
435C------------------------------------------------
436 FINI(1,II) = FLX
437 FINI(2,II) = FLY
438 FINI(3,II) = FLZ
439C------------------------------------------------
440C composantes N/T de la forces nodale -> output
441C------------------------------------------------
442 HL(1:4) = H(1:4,K)
443 CALL I2FORCES(X ,FS ,FX ,FY ,FZ ,
444 . IRECT(1,L),NIR ,FSAV ,FNCONT ,FNCONTP,
445 . FTCONTP ,WEIGHT ,H3D_DATA,I ,HL)
446
447C----------
448 ELSE
449 NSVG(K)= -I
450 L = IRTL(II)
451C
452 IX1(K) = IRECT(1,L)
453 IX2(K) = IRECT(2,L)
454 IX3(K) = IRECT(3,L)
455 IX4(K) = IRECT(4,L)
456 STIF(K)= ZERO
457 VIS(K) = ZERO
458 ENDIF
459 ENDDO
460.OR. IF(IDTMINS==2IDTMINS_INT/=0)THEN
461 DTI=DT2T
462 CALL I2SMS25(LLT ,IX1 ,IX2 ,IX3 ,IX4 ,
463 2 NSVG ,H ,STIF ,NOINT,DMINT2(1,KK),
464 3 NODNX_SMS ,VIS ,STIFM ,DTI)
465 IF(DTI<DT2T)THEN
466 DT2T = DTI
467 NELTST = NOINT
468 ITYPTST = 10
469 ENDIF
470 END IF
471 ENDDO
472C----------
473#include "lockon.inc"
474 ECONT = ECONT + ECONTT ! Elastic energy
475 ECONTD = ECONTD + ECONVT ! Damping Elastic energy
476 FSAV(26) = FSAV(26) + ECONTT
477 FSAV(28) = FSAV(28) + ECONVT
478#include "lockoff.inc"
479C-----------
480 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
h3d_mod
Definition h3d_mod.F:308