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,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 . flx,fly,flz,fs(3),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,stbrk,dti, dxt
95 my_real
96 . h(4,mvsiz),fx(4),fy(4),fz(4),fmx(4),fmy(4),fmz(4),
97 . rx(4),ry(4),rz(4),rm(3),rs(3),vs(3),vm(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 velocities
317 dlx = v1*dt1
318 dly = v2*dt1
319 dlz = v3*dt1
320C-------------- DX == Relative displacement
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 ELSEIF(ms_pena(i)==zero.OR.ms_pena(ix1(k))==zero.OR.
336 . ms_pena(ix2(k))==zero.OR.
337 . ms_pena(ix3(k))==zero.OR.
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 Second 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 n/t components of the nodal forces -> 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 IF(idtmins==2.OR.idtmins_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
i2forces
subroutine i2forces(x, fs, fx, fy, fz, irect, nir, fsav, fncont, fncontp, ftcontp, weight, h3d_data, nsl, h)
Definition i2forces.F:52
i2loceq
subroutine i2loceq(nir, rs, rx, ry, rz, fmx, fmy, fmz, h, stifm)
Definition i2loceq.F:40
i2pen_rot
subroutine i2pen_rot(skew, tt, dt1, stif, rs, rm, v1, v2, v3, rx, ry, rz, va, vb, vc, vd)
Definition i2pen_rot.F:34
i2rep
subroutine i2rep(x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z, nir)
Definition i2rep.F:48
i2sms25
subroutine i2sms25(jlt, ix1, ix2, ix3, ix4, nsvg, h, stif, noint, dmint2, nodnx_sms, vis, stifm, dti)
Definition i2sms25.F:34
min
#define min(a, b)
Definition macros.h:20
h3d_mod
Definition h3d_mod.F:311