OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i2for28_pen.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!|| i2for28_pen ../engine/source/interfaces/interf/i2for28_pen.F
25!||--- called by ------------------------------------------------------
26!|| i2for28 ../engine/source/interfaces/interf/i2for28.F
27!||--- calls -----------------------------------------------------
28!|| i2forces ../engine/source/interfaces/interf/i2forces.F
29!|| i2pen_rot28 ../common_source/interf/i2pen_rot.F
30!|| i2rep ../common_source/interf/i2rep.F
31!|| i2sms28 ../engine/source/interfaces/interf/i2sms28.F
32!||--- uses -----------------------------------------------------
33!|| h3d_mod ../engine/share/modules/h3d_mod.F
34!||====================================================================
35 SUBROUTINE i2for28_pen(X ,V ,VR ,A ,AR ,
36 . MS_PENA ,IN ,STIFN ,STIFR ,WEIGHT ,
37 . NSV ,IRTL ,CRST ,SKEW ,DX ,
38 . DR ,FINI ,FSAV ,FNCONT ,NSN ,
39 . STFN ,STFR ,VISC ,PENFLAG ,IROT ,
40 . NOINT ,NODNX_SMS,DMINT2 ,SAV_FOR_PENA,IRECT,
41 . DT2T ,NELTST ,ITYPTST ,INDXP ,SAV_INER_POFF,
42 . H3D_DATA,FNCONTP ,FTCONTP)
43C-----------------------------------------------
44C M o d u l e s
45C-----------------------------------------------
46 USE h3d_mod
47C-----------------------------------------------
48C I m p l i c i t T y p e s
49C-----------------------------------------------
50#include "implicit_f.inc"
51#include "comlock.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 D u m m y A r g u m e n t s
58C-----------------------------------------------
59 INTEGER NSN,PENFLAG,IROT, NOINT,NELTST,ITYPTST
60 INTEGER IRECT(4,*),NSV(*),IRTL(*),WEIGHT(*),INDXP(NSN),
61 . NODNX_SMS(*)
62C REAL
63 my_real
64 . VISC,DT2T
65 my_real
66 . x(3,*),vr(3,*),v(3,*),a(3,*),ar(3,*),dr(3,*),skew(9,*),
67 . dx(3,*),fini(6,*),ms_pena(*),in(*),stifn(*),stifr(*),stfn(*),stfr(*),
68 . crst(2,*),fsav(*),fncont(3,*),
69 . dmint2(4,*),sav_for_pena(8,*),fncontp(3,*) ,ftcontp(3,*)
70 my_real
71 . sav_iner_poff(*)
72 TYPE (H3D_DATABASE) :: H3D_DATA
73C-----------------------------------------------
74C C o m m o n B l o c k s
75C-----------------------------------------------
76#include "com01_c.inc"
77#include "com06_c.inc"
78#include "com08_c.inc"
79#include "scr11_c.inc"
80#include "scr14_c.inc"
81#include "sms_c.inc"
82#include "task_c.inc"
83C-----------------------------------------------
84C L o c a l V a r i a b l e s
85C-----------------------------------------------
86 INTEGER NIR,I,J,II,JJ,L,W,NN,KK,K,LLT,
87 . IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
88 . NSVG(MVSIZ)
89C REAL
90 my_real
91 . S,T,SP,SM,TP,TM,ECONTT,ECONVT,E1X,E1Y,E1Z,E2X,E2Y,E2Z,E3X,E3Y,E3Z,
92 . FNORM,FLX,FLY,FLZ,FS(3),XSM,YSM,ZSM,XM,YM,ZM,
93 . x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,x0,y0,z0,xs,ys,zs,stf_mom(mvsiz),
94 . vx1,vx2,vx3,vx4,vy1,vy2,vy3,vy4,vz1,vz2,vz3,vz4,dlx,dly,dlz,
95 . vx0,vy0,vz0,vsx,vsy,vsz,vmx,vmy,vmz,vx,vy,vz,dtinv,stf,
96 . fxv,fyv,fzv,drx,dry,drz,stbrk,dti,mharm,dkm,det,b1,b2,b3,c1,c2,c3,
97 . a1,a2,a3,mttx,mtty,mttz,derx,dery,derz, dxt
98 my_real
99 . h(4,mvsiz),fn(3),ft(3),fx(4),fy(4),fz(4),fmx(4),fmy(4),fmz(4),
100 . rx(4),ry(4),rz(4),rm(3),rs(3),v0(3),vs(3),vm(3),
101 . stif(mvsiz), vis(mvsiz), va(3),vb(3),vc(3),vd(3),h2(4,mvsiz)
102 my_real
103 . vrm(3),vrs(3),
104 . vrx0,vrx1,vrx2,vrx3,vrx4,vry0,vry1,vry2,vry3,vry4,vrz0,vrz1,vrz2,vrz3,vrz4,
105 . vrsx,vrsy,vrsz,vrx,vry,vrz,mlx,mly,mlz,mx(4),my(4),mz(4),mrx,mry,mrz,
106 . mgx,mgy,mgz,msx,msy,msz,mvx,mvy,mvz,str,visr(mvsiz),dki,inharm,len2,
107 . hl(4)
108C=======================================================================
109 i7kglo = 1
110 econtt = zero
111 econvt = zero
112 nsvg(1:mvsiz) = 0
113C----------------
114 DO kk=1,nsn,mvsiz
115C
116 llt=min(nsn-kk+1,mvsiz)
117 DO k=1,llt
118C
119 ii = indxp(kk+k-1)
120 IF (ii == 0) cycle
121 i = nsv(ii)
122C
123 IF (i > 0) THEN
124 nsvg(k) = i
125 w = weight(i)
126 s = crst(1,ii)
127 t = crst(2,ii)
128 l = irtl(ii)
129C
130 ix1(k) = irect(1,l)
131 ix2(k) = irect(2,l)
132 ix3(k) = irect(3,l)
133 ix4(k) = irect(4,l)
134C
135 nir= 4
136 sp = one + s
137 sm = one - s
138 tp = fourth*(one + t)
139 tm = fourth*(one - t)
140C
141 h(1,k)=fourth
142 h(2,k)=fourth
143 h(3,k)=fourth
144 h(4,k)=fourth
145C
146 h2(1,k)=tm*sm
147 h2(2,k)=tm*sp
148 h2(3,k)=tp*sp
149 h2(4,k)=tp*sm
150C
151 IF (ix3(k) == ix4(k)) THEN
152 nir = 3
153 h(1,k)=third
154 h(2,k)=third
155 h(3,k)=third
156 h(4,k) = zero
157 ENDIF
158C------------------------------------------------
159C rep local facette main
160C------------------------------------------------
161 x1 = x(1,ix1(k))
162 y1 = x(2,ix1(k))
163 z1 = x(3,ix1(k))
164 x2 = x(1,ix2(k))
165 y2 = x(2,ix2(k))
166 z2 = x(3,ix2(k))
167 x3 = x(1,ix3(k))
168 y3 = x(2,ix3(k))
169 z3 = x(3,ix3(k))
170 x4 = x(1,ix4(k))
171 y4 = x(2,ix4(k))
172 z4 = x(3,ix4(k))
173 xs = x(1,i)
174 ys = x(2,i)
175 zs = x(3,i)
176 vsx = v(1,i)
177 vsy = v(2,i)
178 vsz = v(3,i)
179 vx1 = v(1,ix1(k))
180 vy1 = v(2,ix1(k))
181 vz1 = v(3,ix1(k))
182 vx2 = v(1,ix2(k))
183 vy2 = v(2,ix2(k))
184 vz2 = v(3,ix2(k))
185 vx3 = v(1,ix3(k))
186 vy3 = v(2,ix3(k))
187 vz3 = v(3,ix3(k))
188 vx4 = v(1,ix4(k))
189 vy4 = v(2,ix4(k))
190 vz4 = v(3,ix4(k))
191 IF (iroddl == 1 .AND. sav_iner_poff(i) > zero) THEN
192 vrsx = vr(1,i)
193 vrsy = vr(2,i)
194 vrsz = vr(3,i)
195 vrx1 = vr(1,ix1(k))
196 vry1 = vr(2,ix1(k))
197 vrz1 = vr(3,ix1(k))
198 vrx2 = vr(1,ix2(k))
199 vry2 = vr(2,ix2(k))
200 vrz2 = vr(3,ix2(k))
201 vrx3 = vr(1,ix3(k))
202 vry3 = vr(2,ix3(k))
203 vrz3 = vr(3,ix3(k))
204 vrx4 = vr(1,ix4(k))
205 vry4 = vr(2,ix4(k))
206 vrz4 = vr(3,ix4(k))
207 ENDIF
208C---------------------
209 CALL i2rep(x1 ,x2 ,x3 ,x4 ,
210 . y1 ,y2 ,y3 ,y4 ,
211 . z1 ,z2 ,z3 ,z4 ,
212 . e1x ,e1y ,e1z ,
213 . e2x ,e2y ,e2z ,
214 . e3x ,e3y ,e3z ,nir )
215C------------------------------------------------
216 IF (nir == 4) THEN
217 xm = x1*h2(1,k) + x2*h2(2,k) + x3*h2(3,k) + x4*h2(4,k)
218 ym = y1*h2(1,k) + y2*h2(2,k) + y3*h2(3,k) + y4*h2(4,k)
219 zm = z1*h2(1,k) + z2*h2(2,k) + z3*h2(3,k) + z4*h2(4,k)
220 x0 = (x1 + x2 + x3 + x4)/nir
221 y0 = (y1 + y2 + y3 + y4)/nir
222 z0 = (z1 + z2 + z3 + z4)/nir
223
224 xm = xm - x0
225 ym = ym - y0
226 zm = zm - z0
227 xs = xs - x0
228 ys = ys - y0
229 zs = zs - z0
230 xsm = xs - xm
231 ysm = ys - ym
232 zsm = zs - zm
233c
234 vx0 = (vx1 + vx2 + vx3 + vx4)/nir
235 vy0 = (vy1 + vy2 + vy3 + vy4)/nir
236 vz0 = (vz1 + vz2 + vz3 + vz4)/nir
237 vmx = vx1*h(1,k) + vx2*h(2,k) + vx3*h(3,k) + vx4*h(4,k) - vx0
238 vmy = vy1*h(1,k) + vy2*h(2,k) + vy3*h(3,k) + vy4*h(4,k) - vy0
239 vmz = vz1*h(1,k) + vz2*h(2,k) + vz3*h(3,k) + vz4*h(4,k) - vz0
240C
241 ELSE
242 x0 = (x1 + x2 + x3)/nir
243 y0 = (y1 + y2 + y3)/nir
244 z0 = (z1 + z2 + z3)/nir
245
246 xm = x1*h(1,k) + x2*h(2,k) + x3*h(3,k)
247 ym = y1*h(1,k) + y2*h(2,k) + y3*h(3,k)
248 zm = z1*h(1,k) + z2*h(2,k) + z3*h(3,k)
249
250 xm = xm - x0
251 ym = ym - y0
252 zm = zm - z0
253 xs = xs - x0
254 ys = ys - y0
255 zs = zs - z0
256 xsm = xs - xm
257 ysm = ys - ym
258 zsm = zs - zm
259
260 vx0 = (vx1 + vx2 + vx3)/nir
261 vy0 = (vy1 + vy2 + vy3)/nir
262 vz0 = (vz1 + vz2 + vz3)/nir
263 vmx = vx1*h(1,k) + vx2*h(2,k) + vx3*h(3,k) - vx0
264 vmy = vy1*h(1,k) + vy2*h(2,k) + vy3*h(3,k) - vy0
265 vmz = vz1*h(1,k) + vz2*h(2,k) + vz3*h(3,k) - vz0
266 ENDIF
267 x1 = x1 - x0
268 y1 = y1 - y0
269 z1 = z1 - z0
270 x2 = x2 - x0
271 y2 = y2 - y0
272 z2 = z2 - z0
273 x3 = x3 - x0
274 y3 = y3 - y0
275 z3 = z3 - z0
276 x4 = x4 - x0
277 y4 = y4 - y0
278 z4 = z4 - z0
279 vsx = vsx - vx0
280 vsy = vsy - vy0
281 vsz = vsz - vz0
282C
283c global -> local
284c
285 rs(1) = xs*e1x + ys*e1y + zs*e1z
286 rs(2) = xs*e2x + ys*e2y + zs*e2z
287 rs(3) = xs*e3x + ys*e3y + zs*e3z
288 rm(1) = xm*e1x + ym*e1y + zm*e1z
289 rm(2) = xm*e2x + ym*e2y + zm*e2z
290 rm(3) = xm*e3x + ym*e3y + zm*e3z
291c
292 rx(1) = e1x*x1 + e1y*y1 + e1z*z1
293 ry(1) = e2x*x1 + e2y*y1 + e2z*z1
294 rz(1) = e3x*x1 + e3y*y1 + e3z*z1
295 rx(2) = e1x*x2 + e1y*y2 + e1z*z2
296 ry(2) = e2x*x2 + e2y*y2 + e2z*z2
297 rz(2) = e3x*x2 + e3y*y2 + e3z*z2
298 rx(3) = e1x*x3 + e1y*y3 + e1z*z3
299 ry(3) = e2x*x3 + e2y*y3 + e2z*z3
300 rz(3) = e3x*x3 + e3y*y3 + e3z*z3
301 rx(4) = e1x*x4 + e1y*y4 + e1z*z4
302 ry(4) = e2x*x4 + e2y*y4 + e2z*z4
303 rz(4) = e3x*x4 + e3y*y4 + e3z*z4
304C
305 IF (nir==3) THEN
306 rx(4)=zero
307 ry(4)=zero
308 rz(4)=zero
309 END IF
310C
311 vs(1) = vsx*e1x + vsy*e1y + vsz*e1z
312 vs(2) = vsx*e2x + vsy*e2y + vsz*e2z
313 vs(3) = vsx*e3x + vsy*e3y + vsz*e3z
314 IF (iroddl == 1 .AND. sav_iner_poff(i) > zero) THEN
315 vrs(1) = vrsx*e1x + vrsy*e1y + vrsz*e1z
316 vrs(2) = vrsx*e2x + vrsy*e2y + vrsz*e2z
317 vrs(3) = vrsx*e3x + vrsy*e3y + vrsz*e3z
318 ENDIF
319
320 vm(1) = vmx*e1x + vmy*e1y + vmz*e1z
321 vm(2) = vmx*e2x + vmy*e2y + vmz*e2z
322 vm(3) = vmx*e3x + vmy*e3y + vmz*e3z
323
324 va(1) = vx1*e1x + vy1*e1y + vz1*e1z
325 va(2) = vx1*e2x + vy1*e2y + vz1*e2z
326 va(3) = vx1*e3x + vy1*e3y + vz1*e3z
327
328 vb(1) = vx2*e1x + vy2*e1y + vz2*e1z
329 vb(2) = vx2*e2x + vy2*e2y + vz2*e2z
330 vb(3) = vx2*e3x + vy2*e3y + vz2*e3z
331
332 vc(1) = vx3*e1x + vy3*e1y + vz3*e1z
333 vc(2) = vx3*e2x + vy3*e2y + vz3*e2z
334 vc(3) = vx3*e3x + vy3*e3y + vz3*e3z
335
336 vd(1) = vx4*e1x + vy4*e1y + vz4*e1z
337 vd(2) = vx4*e2x + vy4*e2y + vz4*e2z
338 vd(3) = vx4*e3x + vy4*e3y + vz4*e3z
339C
340C--------- Local displacement
341 IF (tt == zero) THEN
342 dx(1,ii) = zero
343 dx(2,ii) = zero
344 dx(3,ii) = zero
345 fini(1,ii) = zero
346 fini(2,ii) = zero
347 fini(3,ii) = zero
348 dr(1,ii) = zero
349 dr(2,ii) = zero
350 dr(3,ii) = zero
351 fini(4,ii) = zero
352 fini(5,ii) = zero
353 fini(6,ii) = zero
354 ENDIF
355C
356 vx = vs(1) - vm(1)
357 vy = vs(2) - vm(2)
358 vz = vs(3) - vm(3)
359
360C--------- Vi = Vi -VR ^ MS
361 CALL i2pen_rot28(
362 . skew(1,ii) ,tt ,dt1 ,stbrk,
363 . rs ,rm ,vx ,vy ,vz ,
364 . rx ,ry ,rz ,va ,vb ,
365 . vc ,vd ,vrm ,vrs ,det ,
366 . b1 ,b2 ,b3 ,c1 ,c2 ,
367 . c3 ,sav_iner_poff(i))
368C
369 vrx = vrs(1) - vrm(1)
370 vry = vrs(2) - vrm(2)
371 vrz = vrs(3) - vrm(3)
372
373C------------- vers increm en vitesses
374 dlx = vx*dt1
375 dly = vy*dt1
376 dlz = vz*dt1
377 drx = vrx*dt1
378 dry = vry*dt1
379 drz = vrz*dt1
380C------------- DX == deplacement relatif
381 dx(1,ii) = dx(1,ii) + dlx
382 dx(2,ii) = dx(2,ii) + dly
383 dx(3,ii) = dx(3,ii) + dlz
384 dr(1,ii) = dr(1,ii) + drx
385 dr(2,ii) = dr(2,ii) + dry
386 dr(3,ii) = dr(3,ii) + drz
387C
388C------------------------------------------------
389C Calcul de la force
390C------------------------------------------------
391C
392 flx = dx(1,ii) * stfn(ii)
393 fly = dx(2,ii) * stfn(ii)
394 flz = dx(3,ii) * stfn(ii)
395C
396 IF(ms_pena(i)==zero.OR.ms_pena(ix1(k))==zero.OR.
397 . ms_pena(ix2(k))==zero.OR.
398 . ms_pena(ix3(k))==zero.OR.
399 . ms_pena(ix4(k))==zero)THEN
400 mharm = zero
401 ELSEIF(nir==4)THEN
402 mharm = one/ms_pena(i) +
403 . one/ms_pena(ix1(k)) + one/ms_pena(ix2(k)) + one/ms_pena(ix3(k)) + one/ms_pena(ix4(k))
404 mharm = one/mharm
405 ELSE
406 mharm = one/ms_pena(i) +
407 . one/ms_pena(ix1(k)) + one/ms_pena(ix2(k)) + one/ms_pena(ix3(k))
408 mharm = one/mharm
409 END IF
410 dkm = two*stfn(ii)*mharm
411 vis(k) = visc*sqrt(dkm)
412C
413 fxv = vis(k) * vx
414 fyv = vis(k) * vy
415 fzv = vis(k) * vz
416c
417 dxt = dx(1,ii)**2 + dx(2,ii)**2 + dx(3,ii)**2
418 econtt = econtt + half*stfn(ii)*dxt*w
419
420 econvt = econvt + (fxv*vx + fyv*vy + fzv*vz)*dt1*w
421c
422 flx = flx + fxv
423 fly = fly + fyv
424 flz = flz + fzv
425C
426 fs(1) = e1x*flx + e2x*fly + e3x*flz
427 fs(2) = e1y*flx + e2y*fly + e3y*flz
428 fs(3) = e1z*flx + e2z*fly + e3z*flz
429C
430C------------------------------------------------
431C Calcul du Moment
432C------------------------------------------------
433C
434 IF (iroddl == 1 .AND. sav_iner_poff(i) > zero) THEN
435C
436C-- Secnd node shell of spring
437C
438 IF(sav_iner_poff(i)==zero.OR.sav_iner_poff(ix1(k))==zero.OR.
439 . sav_iner_poff(ix2(k))==zero.OR.
440 . sav_iner_poff(ix3(k))==zero.OR.
441 . sav_iner_poff(ix4(k))==zero)THEN
442 inharm = zero
443 ELSEIF(nir==4)THEN
444 inharm = one/sav_iner_poff(i) +
445 . one/sav_iner_poff(ix1(k)) + one/sav_iner_poff(ix2(k)) +
446 . one/sav_iner_poff(ix3(k)) + one/sav_iner_poff(ix4(k))
447 inharm = one/inharm
448 ELSE
449 inharm = one/sav_iner_poff(i) +
450 . one/sav_iner_poff(ix1(k)) + one/sav_iner_poff(ix2(k)) + one/sav_iner_poff(ix3(k))
451 inharm = one/inharm
452 END IF
453C
454 dki = two*stfr(ii)*inharm
455 visr(k) = visc*sqrt(dki)
456C
457 mlx = dr(1,ii) * stfr(ii)
458 mly = dr(2,ii) * stfr(ii)
459 mlz = dr(3,ii) * stfr(ii)
460C
461 mvx = visr(k) * vrx
462 mvy = visr(k) * vry
463 mvz = visr(k) * vrz
464C
465 dxt = dr(1,ii)**2 + dr(2,ii)**2 + dr(3,ii)**2
466 econtt = econtt + half*stfr(ii)*dxt*w
467
468 econvt = econvt + (mvx*vrx
469 . + mvy*vry
470 . + mvz*vrz)*dt1*w
471C
472 mlx = mlx + mvx
473 mly = mly + mvy
474 mlz = mlz + mvz
475C
476 mgx = e1x*mlx + e2x*mly + e3x*mlz
477 mgy = e1y*mlx + e2y*mly + e3y*mlz
478 mgz = e1z*mlx + e2z*mly + e3z*mlz
479C
480 mrx = half*(ysm*fs(3) - zsm*fs(2))
481 mry = half*(zsm*fs(1) - xsm*fs(3))
482 mrz = half*(xsm*fs(2) - ysm*fs(1))
483C
484 ELSE
485C
486C-- Secnd node of solids
487C
488 mgx = zero
489 mgy = zero
490 mgz = zero
491C
492 mrx = ysm*fs(3) - zsm*fs(2)
493 mry = zsm*fs(1) - xsm*fs(3)
494 mrz = xsm*fs(2) - ysm*fs(1)
495C
496 ENDIF
497
498C------------------------------------------------
499C Computation of stiffness for nodal time step
500C------------------------------------------------
501C
502 stf = stfn(ii)*(visc + sqrt(visc**2 + (one+stbrk)))**2
503C
504 len2 = xsm**2+ysm**2+zsm**2
505 str = (stfr(ii)+stfn(ii)*len2)*(visc + sqrt(visc**2 + one))**2
506C
507C----------------------------------------------------
508C Secnd forces/moments -> global coordinates
509C----------------------------------------------------
510C
511 sav_for_pena(1,i) = sav_for_pena(1,i) - fs(1)
512 sav_for_pena(2,i) = sav_for_pena(2,i) - fs(2)
513 sav_for_pena(3,i) = sav_for_pena(3,i) - fs(3)
514 sav_for_pena(4,i) = sav_for_pena(4,i) + stf
515C
516C for SMS ::
517 stif(k) = (one+stbrk)*stfn(ii)
518C
519 IF (iroddl == 1) THEN
520 IF (sav_iner_poff(i)>zero) THEN
521 sav_for_pena(5,i) = sav_for_pena(5,i) - mgx + mrx
522 sav_for_pena(6,i) = sav_for_pena(6,i) - mgy + mry
523 sav_for_pena(7,i) = sav_for_pena(7,i) - mgz + mrz
524 sav_for_pena(8,i) = sav_for_pena(8,i) + str
525 ENDIF
526 ENDIF
527C
528C----------------------------------------------------
529C Main forces/moments
530C----------------------------------------------------
531C
532C---- Transfer or moments in forces
533C
534 mttx=e1x*(mgx+mrx) + e1y*(mgy+mry) + e1z*(mgz+mrz) + rm(2)*flz - rm(3)*fly
535 mtty=e2x*(mgx+mrx) + e2y*(mgy+mry) + e2z*(mgz+mrz) + rm(3)*flx - rm(1)*flz
536 mttz=e3x*(mgx+mrx) + e3y*(mgy+mry) + e3z*(mgz+mrz) + rm(1)*fly - rm(2)*flx
537C
538 a1=det*(mttx*b1+mtty*c3+mttz*c2)
539 a2=det*(mtty*b2+mttz*c1+mttx*c3)
540 a3=det*(mttz*b3+mttx*c2+mtty*c1)
541C
542 derx = (b1+c3+c2)
543 dery = (b2+c1+c3)
544 derz = (b3+c2+c1)
545C for SMS ::
546 stf_mom(k) = det*max(derx,dery,derz)*(str+stf*(xm*xm+ym*ym+zm*zm))
547C
548 DO j=1,4
549 fmx(j) = h(j,k)*flx + a2*rz(j) - a3*ry(j)
550 fmy(j) = h(j,k)*fly + a3*rx(j) - a1*rz(j)
551 fmz(j) = h(j,k)*flz + a1*ry(j) - a2*rx(j)
552 ENDDO
553C
554 DO j=1,4
555 fx(j) = e1x*fmx(j) + e2x*fmy(j) + e3x*fmz(j)
556 fy(j) = e1y*fmx(j) + e2y*fmy(j) + e3y*fmz(j)
557 fz(j) = e1z*fmx(j) + e2z*fmy(j) + e3z*fmz(j)
558 ENDDO
559C
560 IF (w == 1) THEN
561C
562 a(1,ix1(k)) = a(1,ix1(k)) + fx(1)
563 a(2,ix1(k)) = a(2,ix1(k)) + fy(1)
564 a(3,ix1(k)) = a(3,ix1(k)) + fz(1)
565 stifn(ix1(k)) = stifn(ix1(k))+abs(stf*h(1,k))+stf_mom(k)
566c
567 a(1,ix2(k)) = a(1,ix2(k)) + fx(2)
568 a(2,ix2(k)) = a(2,ix2(k)) + fy(2)
569 a(3,ix2(k)) = a(3,ix2(k)) + fz(2)
570 stifn(ix2(k)) = stifn(ix2(k))+abs(stf*h(2,k))+stf_mom(k)
571c
572 a(1,ix3(k)) = a(1,ix3(k)) + fx(3)
573 a(2,ix3(k)) = a(2,ix3(k)) + fy(3)
574 a(3,ix3(k)) = a(3,ix3(k)) + fz(3)
575 stifn(ix3(k)) = stifn(ix3(k))+abs(stf*h(3,k))+stf_mom(k)
576c
577 IF (nir==4) THEN
578 a(1,ix4(k)) = a(1,ix4(k)) + fx(4)
579 a(2,ix4(k)) = a(2,ix4(k)) + fy(4)
580 a(3,ix4(k)) = a(3,ix4(k)) + fz(4)
581 stifn(ix4(k)) = stifn(ix4(k))+abs(stf*h(4,k))+stf_mom(k)
582 ENDIF
583 ENDIF
584C
585
586C------------------------------------------------
587 fini(1,ii) = flx
588 fini(2,ii) = fly
589 fini(3,ii) = flz
590 IF (iroddl == 1 .AND. sav_iner_poff(i) > zero) THEN
591 fini(4,ii) = mlx
592 fini(5,ii) = mly
593 fini(6,ii) = mlz
594 ENDIF
595C------------------------------------------------
596C composantes N/T de la forces nodale -> output
597C------------------------------------------------
598 hl(1:4) = h(1:4,k)
599 CALL i2forces(x ,fs ,fx ,fy ,fz ,
600 . irect(1,l),nir ,fsav ,fncont ,fncontp,
601 . ftcontp ,weight ,h3d_data,i ,hl)
602C----------
603 ELSE ! desactivated secnd node
604 nsvg(k)= -i
605 l = irtl(ii)
606C
607 ix1(k) = irect(1,l)
608 ix2(k) = irect(2,l)
609 ix3(k) = irect(3,l)
610 ix4(k) = irect(4,l)
611 stif(k)= zero
612 vis(k) = zero
613 ENDIF
614 ENDDO
615c
616 IF(idtmins==2.OR.idtmins_int/=0)THEN
617 dti=dt2t
618 CALL i2sms28(llt ,ix1 ,ix2 ,ix3 ,ix4 ,
619 2 nsvg ,h ,stif ,noint ,
620 3 dmint2(1,kk),nodnx_sms ,vis ,
621 4 stf_mom ,dti )
622 IF(dti<dt2t)THEN
623 dt2t = dti
624 neltst = noint
625 ityptst = 10
626 ENDIF
627 END IF
628 ENDDO
629C----------
630#include "lockon.inc"
631 econt = econt + econtt ! Elastic energy
632 econtd = econtd + econvt ! Damping Elastic energy
633 fsav(26) = fsav(26) + econtt
634 fsav(28) = fsav(28) + econvt
635#include "lockoff.inc"
636C-----------
637 RETURN
638 END SUBROUTINE i2for28_pen
i2for28_pen
subroutine i2for28_pen(x, v, vr, a, ar, ms_pena, in, stifn, stifr, weight, nsv, irtl, crst, skew, dx, dr, fini, fsav, fncont, nsn, stfn, stfr, visc, penflag, irot, noint, nodnx_sms, dmint2, sav_for_pena, irect, dt2t, neltst, ityptst, indxp, sav_iner_poff, h3d_data, fncontp, ftcontp)
Definition i2for28_pen.F:43
i2forces
subroutine i2forces(x, fs, fx, fy, fz, irect, nir, fsav, fncont, fncontp, ftcontp, weight, h3d_data, nsl, h)
Definition i2forces.F:52
i2pen_rot28
subroutine i2pen_rot28(skew, tt, dt1, stif, rs, rm, vx, vy, vz, rx, ry, rz, va, vb, vc, vd, vrm, vrs, det, b1, b2, b3, c1, c2, c3, in_secnd)
Definition i2pen_rot.F:403
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
i2sms28
subroutine i2sms28(jlt, ix1, ix2, ix3, ix4, nsvg, h, stif, noint, dmint2, nodnx_sms, vis, stf_mom, dti)
Definition i2sms28.F:34
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
h3d_mod
Definition h3d_mod.F:308