OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i2for25.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!|| i2for25 ../engine/source/interfaces/interf/i2for25.F
25!||--- called by ------------------------------------------------------
26!|| intti2f ../engine/source/interfaces/interf/intti2f.F
27!||--- calls -----------------------------------------------------
28!|| i2forces ../engine/source/interfaces/interf/i2forces.F
29!|| i2loceq ../common_source/interf/i2loceq.F
30!|| i2pen_rot ../common_source/interf/i2pen_rot.F
31!|| i2rep ../common_source/interf/i2rep.f
32!|| i2sms25 ../engine/source/interfaces/interf/i2sms25.F
33!||--- uses -----------------------------------------------------
34!|| h3d_mod ../engine/share/modules/h3d_mod.F
35!||====================================================================
36 SUBROUTINE i2for25(X ,V ,VR ,A ,AR ,
37 . MS ,STIFN ,STIFR ,WEIGHT ,IRECT ,
38 . NSV ,IRTL ,CRST ,SKEW ,XINI ,
39 . DX ,FINI ,FSAV ,FNCONT ,NSN ,
40 . STFN ,STFR ,VISC ,PENFLAG ,IROT ,
41 . NOINT ,NODNX_SMS,DMINT2 ,SAV_FOR_PENA ,
42 . MS_PENA ,DT2T ,NELTST ,ITYPTST,IVISC ,
43 . H3D_DATA,FNCONTP ,FTCONTP)
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 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 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 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
481 END SUBROUTINE i2for25
i2for25
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)
Definition i2for25.F:44
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:308
h3d_mod::h3d_database
Definition h3d_mod.F:435