OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i2for28p_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!|| i2for28p_pen ../engine/source/interfaces/interf/i2for28p_pen.F
25!||--- called by ------------------------------------------------------
26!|| i2for28p ../engine/source/interfaces/interf/i2for28p.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 i2for28p_pen(
36 . X ,V ,VR ,A ,AR ,
37 . MS ,IN ,STIFN ,STIFR ,WEIGHT ,
38 . NSV ,IRTL ,CRST ,SKEW ,DX ,
39 . DR ,FINI ,FSAV ,FNCONT ,NSN ,
40 . I0 ,I2SIZE ,IADI2 ,FSKYI2 ,STFN ,
41 . STFR ,VISC ,PENFLAG ,IROT ,NOINT ,
42 . NODNX_SMS,DMINT2 ,DT2T ,NELTST ,ITYPTST ,
43 . IRECT ,INDXP ,IADX ,
44 . H3D_DATA ,FNCONTP ,FTCONTP)
45C-----------------------------------------------
46C M o d u l e s
47C-----------------------------------------------
48 USE h3d_mod
49C-----------------------------------------------
50C I m p l i c i t T y p e s
51C-----------------------------------------------
52#include "implicit_f.inc"
53#include "comlock.inc"
54C-----------------------------------------------
55C G l o b a l P a r a m e t e r s
56C-----------------------------------------------
57#include "mvsiz_p.inc"
58C-----------------------------------------------
59C D u m m y A r g u m e n t s
60C-----------------------------------------------
61 INTEGER NSN, I0,I2SIZE,PENFLAG,IROT, NOINT,NELTST,ITYPTST
62 INTEGER IRECT(4,*),NSV(*),IRTL(*),WEIGHT(*),IADI2(4,*),
63 . NODNX_SMS(*),IADX(*),INDXP(*)
64C REAL
65 my_real
66 . VISC,DT2T
67 my_real
68 . X(3,*),VR(3,*),V(3,*),A(3,*),AR(3,*),MS(*),CRST(2,*),
69 . SKEW(9,*),DX(3,*),DR(3,*),FINI(6,*),FSAV(*),FNCONT(3,*),
70 . stifn(*),stifr(*),stfn(*),stfr(*),fskyi2(i2size,*),
71 . dmint2(4,*),in(*),fncontp(3,*) ,ftcontp(3,*)
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),I0BASE
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 i0base = i0
113C
114 nsvg(1:mvsiz) = 0
115C----------------
116 DO kk=1,nsn,mvsiz
117C
118 llt=min(nsn-kk+1,mvsiz)
119 DO k=1,llt
120C
121 ii = indxp(kk+k-1)
122 IF (ii == 0) cycle
123 i = nsv(ii)
124C
125 IF (i > 0) THEN
126 nsvg(k) = i
127 w = weight(i)
128 s = crst(1,ii)
129 t = crst(2,ii)
130 l = irtl(ii)
131C
132 ix1(k) = irect(1,l)
133 ix2(k) = irect(2,l)
134 ix3(k) = irect(3,l)
135 ix4(k) = irect(4,l)
136C
137 nir= 4
138 sp = one + s
139 sm = one - s
140 tp = fourth*(one + t)
141 tm = fourth*(one - t)
142C
143 h(1,k)=fourth
144 h(2,k)=fourth
145 h(3,k)=fourth
146 h(4,k)=fourth
147C
148 h2(1,k)=tm*sm
149 h2(2,k)=tm*sp
150 h2(3,k)=tp*sp
151 h2(4,k)=tp*sm
152C
153 IF (ix3(k) == ix4(k)) THEN
154 nir = 3
155 h(1,k)=third
156 h(2,k)=third
157 h(3,k)=third
158 h(4,k) = zero
159 ENDIF
160C------------------------------------------------
161C rep local facette main
162C------------------------------------------------
163 x1 = x(1,ix1(k))
164 y1 = x(2,ix1(k))
165 z1 = x(3,ix1(k))
166 x2 = x(1,ix2(k))
167 y2 = x(2,ix2(k))
168 z2 = x(3,ix2(k))
169 x3 = x(1,ix3(k))
170 y3 = x(2,ix3(k))
171 z3 = x(3,ix3(k))
172 x4 = x(1,ix4(k))
173 y4 = x(2,ix4(k))
174 z4 = x(3,ix4(k))
175 xs = x(1,i)
176 ys = x(2,i)
177 zs = x(3,i)
178 vsx = v(1,i)
179 vsy = v(2,i)
180 vsz = v(3,i)
181 vx1 = v(1,ix1(k))
182 vy1 = v(2,ix1(k))
183 vz1 = v(3,ix1(k))
184 vx2 = v(1,ix2(k))
185 vy2 = v(2,ix2(k))
186 vz2 = v(3,ix2(k))
187 vx3 = v(1,ix3(k))
188 vy3 = v(2,ix3(k))
189 vz3 = v(3,ix3(k))
190 vx4 = v(1,ix4(k))
191 vy4 = v(2,ix4(k))
192 vz4 = v(3,ix4(k))
193 IF (iroddl == 1 .AND. in(i) > zero) THEN
194 vrsx = vr(1,i)
195 vrsy = vr(2,i)
196 vrsz = vr(3,i)
197 vrx1 = vr(1,ix1(k))
198 vry1 = vr(2,ix1(k))
199 vrz1 = vr(3,ix1(k))
200 vrx2 = vr(1,ix2(k))
201 vry2 = vr(2,ix2(k))
202 vrz2 = vr(3,ix2(k))
203 vrx3 = vr(1,ix3(k))
204 vry3 = vr(2,ix3(k))
205 vrz3 = vr(3,ix3(k))
206 vrx4 = vr(1,ix4(k))
207 vry4 = vr(2,ix4(k))
208 vrz4 = vr(3,ix4(k))
209 ENDIF
210C---------------------
211 CALL i2rep(x1 ,x2 ,x3 ,x4 ,
212 . y1 ,y2 ,y3 ,y4 ,
213 . z1 ,z2 ,z3 ,z4 ,
214 . e1x ,e1y ,e1z ,
215 . e2x ,e2y ,e2z ,
216 . e3x ,e3y ,e3z ,nir )
217C------------------------------------------------
218 IF (nir == 4) THEN
219 xm = x1*h2(1,k) + x2*h2(2,k) + x3*h2(3,k) + x4*h2(4,k)
220 ym = y1*h2(1,k) + y2*h2(2,k) + y3*h2(3,k) + y4*h2(4,k)
221 zm = z1*h2(1,k) + z2*h2(2,k) + z3*h2(3,k) + z4*h2(4,k)
222 x0 = (x1 + x2 + x3 + x4)/nir
223 y0 = (y1 + y2 + y3 + y4)/nir
224 z0 = (z1 + z2 + z3 + z4)/nir
225
226 xm = xm - x0
227 ym = ym - y0
228 zm = zm - z0
229 xs = xs - x0
230 ys = ys - y0
231 zs = zs - z0
232 xsm = xs - xm
233 ysm = ys - ym
234 zsm = zs - zm
235c
236 vx0 = (vx1 + vx2 + vx3 + vx4)/nir
237 vy0 = (vy1 + vy2 + vy3 + vy4)/nir
238 vz0 = (vz1 + vz2 + vz3 + vz4)/nir
239 vmx = vx1*h(1,k) + vx2*h(2,k) + vx3*h(3,k) + vx4*h(4,k) - vx0
240 vmy = vy1*h(1,k) + vy2*h(2,k) + vy3*h(3,k) + vy4*h(4,k) - vy0
241 vmz = vz1*h(1,k) + vz2*h(2,k) + vz3*h(3,k) + vz4*h(4,k) - vz0
242C
243 ELSE
244 x0 = (x1 + x2 + x3)/nir
245 y0 = (y1 + y2 + y3)/nir
246 z0 = (z1 + z2 + z3)/nir
247
248 xm = x1*h(1,k) + x2*h(2,k) + x3*h(3,k)
249 ym = y1*h(1,k) + y2*h(2,k) + y3*h(3,k)
250 zm = z1*h(1,k) + z2*h(2,k) + z3*h(3,k)
251
252 xm = xm - x0
253 ym = ym - y0
254 zm = zm - z0
255 xs = xs - x0
256 ys = ys - y0
257 zs = zs - z0
258 xsm = xs - xm
259 ysm = ys - ym
260 zsm = zs - zm
261
262 vx0 = (vx1 + vx2 + vx3)/nir
263 vy0 = (vy1 + vy2 + vy3)/nir
264 vz0 = (vz1 + vz2 + vz3)/nir
265 vmx = vx1*h(1,k) + vx2*h(2,k) + vx3*h(3,k) - vx0
266 vmy = vy1*h(1,k) + vy2*h(2,k) + vy3*h(3,k) - vy0
267 vmz = vz1*h(1,k) + vz2*h(2,k) + vz3*h(3,k) - vz0
268 ENDIF
269 x1 = x1 - x0
270 y1 = y1 - y0
271 z1 = z1 - z0
272 x2 = x2 - x0
273 y2 = y2 - y0
274 z2 = z2 - z0
275 x3 = x3 - x0
276 y3 = y3 - y0
277 z3 = z3 - z0
278 x4 = x4 - x0
279 y4 = y4 - y0
280 z4 = z4 - z0
281 vsx = vsx - vx0
282 vsy = vsy - vy0
283 vsz = vsz - vz0
284C
285c global -> local
286c
287 rs(1) = xs*e1x + ys*e1y + zs*e1z
288 rs(2) = xs*e2x + ys*e2y + zs*e2z
289 rs(3) = xs*e3x + ys*e3y + zs*e3z
290 rm(1) = xm*e1x + ym*e1y + zm*e1z
291 rm(2) = xm*e2x + ym*e2y + zm*e2z
292 rm(3) = xm*e3x + ym*e3y + zm*e3z
293c
294 rx(1) = e1x*x1 + e1y*y1 + e1z*z1
295 ry(1) = e2x*x1 + e2y*y1 + e2z*z1
296 rz(1) = e3x*x1 + e3y*y1 + e3z*z1
297 rx(2) = e1x*x2 + e1y*y2 + e1z*z2
298 ry(2) = e2x*x2 + e2y*y2 + e2z*z2
299 rz(2) = e3x*x2 + e3y*y2 + e3z*z2
300 rx(3) = e1x*x3 + e1y*y3 + e1z*z3
301 ry(3) = e2x*x3 + e2y*y3 + e2z*z3
302 rz(3) = e3x*x3 + e3y*y3 + e3z*z3
303 rx(4) = e1x*x4 + e1y*y4 + e1z*z4
304 ry(4) = e2x*x4 + e2y*y4 + e2z*z4
305 rz(4) = e3x*x4 + e3y*y4 + e3z*z4
306C
307 IF (nir==3) THEN
308 rx(4)=zero
309 ry(4)=zero
310 rz(4)=zero
311 END IF
312C
313 vs(1) = vsx*e1x + vsy*e1y + vsz*e1z
314 vs(2) = vsx*e2x + vsy*e2y + vsz*e2z
315 vs(3) = vsx*e3x + vsy*e3y + vsz*e3z
316 IF (iroddl == 1 .AND. in(i) > zero) THEN
317 vrs(1) = vrsx*e1x + vrsy*e1y + vrsz*e1z
318 vrs(2) = vrsx*e2x + vrsy*e2y + vrsz*e2z
319 vrs(3) = vrsx*e3x + vrsy*e3y + vrsz*e3z
320 ENDIF
321
322 vm(1) = vmx*e1x + vmy*e1y + vmz*e1z
323 vm(2) = vmx*e2x + vmy*e2y + vmz*e2z
324 vm(3) = vmx*e3x + vmy*e3y + vmz*e3z
325
326 va(1) = vx1*e1x + vy1*e1y + vz1*e1z
327 va(2) = vx1*e2x + vy1*e2y + vz1*e2z
328 va(3) = vx1*e3x + vy1*e3y + vz1*e3z
329
330 vb(1) = vx2*e1x + vy2*e1y + vz2*e1z
331 vb(2) = vx2*e2x + vy2*e2y + vz2*e2z
332 vb(3) = vx2*e3x + vy2*e3y + vz2*e3z
333
334 vc(1) = vx3*e1x + vy3*e1y + vz3*e1z
335 vc(2) = vx3*e2x + vy3*e2y + vz3*e2z
336 vc(3) = vx3*e3x + vy3*e3y + vz3*e3z
337
338 vd(1) = vx4*e1x + vy4*e1y + vz4*e1z
339 vd(2) = vx4*e2x + vy4*e2y + vz4*e2z
340 vd(3) = vx4*e3x + vy4*e3y + vz4*e3z
341C
342C--------- Local displacement
343 IF (tt == zero) THEN
344 dx(1,ii) = zero
345 dx(2,ii) = zero
346 dx(3,ii) = zero
347 fini(1,ii) = zero
348 fini(2,ii) = zero
349 fini(3,ii) = zero
350 dr(1,ii) = zero
351 dr(2,ii) = zero
352 dr(3,ii) = zero
353 fini(4,ii) = zero
354 fini(5,ii) = zero
355 fini(6,ii) = zero
356 ENDIF
357C
358 vx = vs(1) - vm(1)
359 vy = vs(2) - vm(2)
360 vz = vs(3) - vm(3)
361
362C--------- Vi = Vi -VR ^ MS
363 CALL i2pen_rot28(
364 . skew(1,ii) ,tt ,dt1 ,stbrk,
365 . rs ,rm ,vx ,vy ,vz ,
366 . rx ,ry ,rz ,va ,vb ,
367 . vc ,vd ,vrm ,vrs ,det ,
368 . b1 ,b2 ,b3 ,c1 ,c2 ,
369 . c3 ,in(i))
370C
371 vrx = vrs(1) - vrm(1)
372 vry = vrs(2) - vrm(2)
373 vrz = vrs(3) - vrm(3)
374
375C------------- vers increm en vitesses
376 dlx = vx*dt1
377 dly = vy*dt1
378 dlz = vz*dt1
379 drx = vrx*dt1
380 dry = vry*dt1
381 drz = vrz*dt1
382C------------- DX == deplacement relatif
383 dx(1,ii) = dx(1,ii) + dlx
384 dx(2,ii) = dx(2,ii) + dly
385 dx(3,ii) = dx(3,ii) + dlz
386 dr(1,ii) = dr(1,ii) + drx
387 dr(2,ii) = dr(2,ii) + dry
388 dr(3,ii) = dr(3,ii) + drz
389C
390C------------------------------------------------
391C Calcul de la force
392C------------------------------------------------
393C
394 flx = dx(1,ii) * stfn(ii)
395 fly = dx(2,ii) * stfn(ii)
396 flz = dx(3,ii) * stfn(ii)
397C
398 IF(ms(i)==zero.OR.ms(ix1(k))==zero.OR.
399 . ms(ix2(k))==zero.OR.
400 . ms(ix3(k))==zero.OR.
401 . ms(ix4(k))==zero)THEN
402 mharm = zero
403 ELSEIF(nir==4)THEN
404 mharm = one/ms(i) +
405 . one/ms(ix1(k)) + one/ms(ix2(k)) + one/ms(ix3(k)) + one/ms(ix4(k))
406 mharm = one/mharm
407 ELSE
408 mharm = one/ms(i) +
409 . one/ms(ix1(k)) + one/ms(ix2(k)) + one/ms(ix3(k))
410 mharm = one/mharm
411 END IF
412 dkm = two*stfn(ii)*mharm
413 vis(k) = visc*sqrt(dkm)
414C
415 fxv = vis(k) * vx
416 fyv = vis(k) * vy
417 fzv = vis(k) * vz
418c
419 dxt = dx(1,ii)**2 + dx(2,ii)**2 + dx(3,ii)**2
420 econtt = econtt + half*stfn(ii)*dxt*w
421
422 econvt = econvt + (fxv*vx + fyv*vy + fzv*vz)*dt1*w
423c
424 flx = flx + fxv
425 fly = fly + fyv
426 flz = flz + fzv
427C
428 fs(1) = e1x*flx + e2x*fly + e3x*flz
429 fs(2) = e1y*flx + e2y*fly + e3y*flz
430 fs(3) = e1z*flx + e2z*fly + e3z*flz
431C
432C------------------------------------------------
433C Calcul du Moment
434C------------------------------------------------
435C
436 IF (iroddl == 1 .AND. in(i) > zero) THEN
437C
438C-- Secnd node shell of spring
439C
440 IF(in(i)==zero.OR.in(ix1(k))==zero.OR.
441 . in(ix2(k))==zero.OR.
442 . in(ix3(k))==zero.OR.
443 . in(ix4(k))==zero)THEN
444 inharm = zero
445 ELSEIF(nir==4)THEN
446 inharm = one/in(i) +
447 . one/in(ix1(k)) + one/in(ix2(k)) + one/in(ix3(k)) + one/in(ix4(k))
448 inharm = one/inharm
449 ELSE
450 inharm = one/in(i) +
451 . one/in(ix1(k)) + one/in(ix2(k)) + one/in(ix3(k))
452 inharm = one/inharm
453 END IF
454C
455 dki = two*stfr(ii)*inharm
456 visr(k) = visc*sqrt(dki)
457C
458 mlx = dr(1,ii) * stfr(ii)
459 mly = dr(2,ii) * stfr(ii)
460 mlz = dr(3,ii) * stfr(ii)
461C
462 mvx = visr(k) * vrx
463 mvy = visr(k) * vry
464 mvz = visr(k) * vrz
465C
466 dxt = dr(1,ii)**2 + dr(2,ii)**2 + dr(3,ii)**2
467 econtt = econtt + half*stfr(ii)*dxt*w
468
469 econvt = econvt + (mvx*vrx
470 . + mvy*vry
471 . + mvz*vrz)*dt1*w
472C
473 mlx = mlx + mvx
474 mly = mly + mvy
475 mlz = mlz + mvz
476C
477 mgx = e1x*mlx + e2x*mly + e3x*mlz
478 mgy = e1y*mlx + e2y*mly + e3y*mlz
479 mgz = e1z*mlx + e2z*mly + e3z*mlz
480C
481 mrx = half*(ysm*fs(3) - zsm*fs(2))
482 mry = half*(zsm*fs(1) - xsm*fs(3))
483 mrz = half*(xsm*fs(2) - ysm*fs(1))
484C
485 ELSE
486C
487C-- Secnd node of solids
488C
489 mgx = zero
490 mgy = zero
491 mgz = zero
492C
493 mrx = ysm*fs(3) - zsm*fs(2)
494 mry = zsm*fs(1) - xsm*fs(3)
495 mrz = xsm*fs(2) - ysm*fs(1)
496C
497 ENDIF
498C
499C------------------------------------------------
500C Computation of stiffness for nodal time step
501C------------------------------------------------
502C
503 stf = stfn(ii)*(visc + sqrt(visc**2 + (one+stbrk)))**2
504C
505 len2 = xsm**2+ysm**2+zsm**2
506 str = (stfr(ii)+stfn(ii)*len2)*(visc + sqrt(visc**2 + one))**2
507C
508C----------------------------------------------------
509C Secnd forces/moments -> global coordinates
510C----------------------------------------------------
511C
512 a(1,i) = a(1,i) - fs(1)
513 a(2,i) = a(2,i) - fs(2)
514 a(3,i) = a(3,i) - fs(3)
515 stifn(i) = stifn(i) + stf
516C
517C for SMS ::
518 stif(k) = (one+stbrk)*stfn(ii)
519C
520 IF (iroddl == 1) THEN
521 IF (in(i)>zero) THEN
522 ar(1,i) = ar(1,i) - mgx + mrx
523 ar(2,i) = ar(2,i) - mgy + mry
524 ar(3,i) = ar(3,i) - mgz + mrz
525 stifr(i) = stifr(i) + str
526 ENDIF
527 ENDIF
528C
529C----------------------------------------------------
530C Main forces/moments
531C----------------------------------------------------
532C
533C---- Transfer or moments in forces
534C
535 mttx=e1x*(mgx+mrx) + e1y*(mgy+mry) + e1z*(mgz+mrz) + rm(2)*flz - rm(3)*fly
536 mtty=e2x*(mgx+mrx) + e2y*(mgy+mry) + e2z*(mgz+mrz) + rm(3)*flx - rm(1)*flz
537 mttz=e3x*(mgx+mrx) + e3y*(mgy+mry) + e3z*(mgz+mrz) + rm(1)*fly - rm(2)*flx
538C
539 a1=det*(mttx*b1+mtty*c3+mttz*c2)
540 a2=det*(mtty*b2+mttz*c1+mttx*c3)
541 a3=det*(mttz*b3+mttx*c2+mtty*c1)
542C
543 derx = (b1+c3+c2)
544 dery = (b2+c1+c3)
545 derz = (b3+c2+c1)
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
561 i0 = i0base + iadx(ii)
562 jj = 1
563 nn = iadi2(jj,i0)
564c
565 fskyi2(1,nn) = fx(jj)
566 fskyi2(2,nn) = fy(jj)
567 fskyi2(3,nn) = fz(jj)
568 fskyi2(4,nn) = zero
569 fskyi2(5,nn) = abs(stf*h(jj,k))+stf_mom(k)
570 IF (iroddl == 1) THEN
571 fskyi2(6,nn) = zero
572 fskyi2(7,nn) = zero
573 fskyi2(8,nn) = zero
574 fskyi2(9,nn) = zero
575 fskyi2(10,nn)= zero
576 ENDIF
577c
578 jj = 2
579 nn = iadi2(jj,i0)
580 fskyi2(1,nn) = fx(jj)
581 fskyi2(2,nn) = fy(jj)
582 fskyi2(3,nn) = fz(jj)
583 fskyi2(4,nn) = zero
584 fskyi2(5,nn) = abs(stf*h(jj,k))+stf_mom(k)
585 IF (iroddl == 1) THEN
586 fskyi2(6,nn) = zero
587 fskyi2(7,nn) = zero
588 fskyi2(8,nn) = zero
589 fskyi2(9,nn) = zero
590 fskyi2(10,nn)= zero
591 ENDIF
592c
593 jj = 3
594 nn = iadi2(jj,i0)
595 fskyi2(1,nn) = fx(jj)
596 fskyi2(2,nn) = fy(jj)
597 fskyi2(3,nn) = fz(jj)
598 fskyi2(4,nn) = zero
599 fskyi2(5,nn) = abs(stf*h(jj,k))+stf_mom(k)
600 IF (iroddl == 1) THEN
601 fskyi2(6,nn) = zero
602 fskyi2(7,nn) = zero
603 fskyi2(8,nn) = zero
604 fskyi2(9,nn) = zero
605 fskyi2(10,nn)= zero
606 ENDIF
607c
608 IF (nir==4) THEN
609 jj = 4
610 nn = iadi2(jj,i0)
611 fskyi2(1,nn) = fx(jj)
612 fskyi2(2,nn) = fy(jj)
613 fskyi2(3,nn) = fz(jj)
614 fskyi2(4,nn) = zero
615 fskyi2(5,nn) = abs(stf*h(jj,k))+stf_mom(k)
616 IF (iroddl == 1) THEN
617 fskyi2(6,nn) = zero
618 fskyi2(7,nn) = zero
619 fskyi2(8,nn) = zero
620 fskyi2(9,nn) = zero
621 fskyi2(10,nn)= zero
622 ENDIF
623 ELSE
624 jj = 4
625 nn = iadi2(jj,i0)
626 fskyi2(1,nn) = zero
627 fskyi2(2,nn) = zero
628 fskyi2(3,nn) = zero
629 fskyi2(4,nn) = zero
630 fskyi2(5,nn) = zero
631 IF (iroddl == 1) THEN
632 fskyi2(6,nn) = zero
633 fskyi2(7,nn) = zero
634 fskyi2(8,nn) = zero
635 fskyi2(9,nn) = zero
636 fskyi2(10,nn)= zero
637 ENDIF
638 ENDIF
639C
640 ENDIF
641C------------------------------------------------
642 fini(1,ii) = flx
643 fini(2,ii) = fly
644 fini(3,ii) = flz
645 IF (iroddl == 1 .AND. in(i) > zero) THEN
646 fini(4,ii) = mlx
647 fini(5,ii) = mly
648 fini(6,ii) = mlz
649 ENDIF
650C------------------------------------------------
651C composantes N/T de la forces nodale -> output
652C------------------------------------------------
653 hl(1:4) = h(1:4,k)
654 CALL i2forces(x ,fs ,fx ,fy ,fz ,
655 . irect(1,l),nir ,fsav ,fncont ,fncontp,
656 . ftcontp ,weight ,h3d_data,i ,hl)
657
658C----------
659 ELSE ! desactivated secnd node
660 nsvg(k)= -i
661 l = irtl(ii)
662C
663 ix1(k) = irect(1,l)
664 ix2(k) = irect(2,l)
665 ix3(k) = irect(3,l)
666 ix4(k) = irect(4,l)
667 stif(k)= zero
668 vis(k) = zero
669 IF (weight(-i) == 1) THEN ! stokage ZERO pour noeuds delete par idel2
670 i0 = i0base + iadx(ii)
671 jj = 1
672 nn = iadi2(jj,i0)
673 fskyi2(1,nn) = zero
674 fskyi2(2,nn) = zero
675 fskyi2(3,nn) = zero
676 fskyi2(4,nn) = zero
677 fskyi2(5,nn) = zero
678 IF (iroddl == 1) THEN
679 fskyi2(6,nn) = zero
680 fskyi2(7,nn) = zero
681 fskyi2(8,nn) = zero
682 fskyi2(9,nn) = zero
683 fskyi2(10,nn)= zero
684 ENDIF
685 jj = 2
686 nn = iadi2(jj,i0)
687 fskyi2(1,nn) = zero
688 fskyi2(2,nn) = zero
689 fskyi2(3,nn) = zero
690 fskyi2(4,nn) = zero
691 fskyi2(5,nn) = zero
692 IF (iroddl == 1) THEN
693 fskyi2(6,nn) = zero
694 fskyi2(7,nn) = zero
695 fskyi2(8,nn) = zero
696 fskyi2(9,nn) = zero
697 fskyi2(10,nn)= zero
698 ENDIF
699 jj = 3
700 nn = iadi2(jj,i0)
701 fskyi2(1,nn) = zero
702 fskyi2(2,nn) = zero
703 fskyi2(3,nn) = zero
704 fskyi2(4,nn) = zero
705 fskyi2(5,nn) = zero
706 IF (iroddl == 1) THEN
707 fskyi2(6,nn) = zero
708 fskyi2(7,nn) = zero
709 fskyi2(8,nn) = zero
710 fskyi2(9,nn) = zero
711 fskyi2(10,nn)= zero
712 ENDIF
713 jj = 4
714 nn = iadi2(jj,i0)
715 fskyi2(1,nn) = zero
716 fskyi2(2,nn) = zero
717 fskyi2(3,nn) = zero
718 fskyi2(4,nn) = zero
719 fskyi2(5,nn) = zero
720 IF (iroddl == 1) THEN
721 fskyi2(6,nn) = zero
722 fskyi2(7,nn) = zero
723 fskyi2(8,nn) = zero
724 fskyi2(9,nn) = zero
725 fskyi2(10,nn)= zero
726 ENDIF
727 ENDIF
728 ENDIF
729 ENDDO
730c
731 IF(idtmins==2.OR.idtmins_int/=0)THEN
732 dti=dt2t
733 CALL i2sms28(llt ,ix1 ,ix2 ,ix3 ,ix4 ,
734 2 nsvg ,h ,stif ,noint ,
735 3 dmint2(1,kk),nodnx_sms ,vis ,
736 4 stf_mom ,dti)
737 IF(dti<dt2t)THEN
738 dt2t = dti
739 neltst = noint
740 ityptst = 10
741 ENDIF
742 END IF
743 ENDDO
744C----------
745#include "lockon.inc"
746 econt = econt + econtt ! Elastic energy
747 econtd = econtd + econvt ! Damping Elastic energy
748 fsav(26) = fsav(26) + econtt !
749 fsav(28) = fsav(28) + econvt
750#include "lockoff.inc"
751C-----------
752 RETURN
753 END SUBROUTINE i2for28p_pen
i2for28p_pen
subroutine i2for28p_pen(x, v, vr, a, ar, ms, in, stifn, stifr, weight, nsv, irtl, crst, skew, dx, dr, fini, fsav, fncont, nsn, i0, i2size, iadi2, fskyi2, stfn, stfr, visc, penflag, irot, noint, nodnx_sms, dmint2, dt2t, neltst, ityptst, irect, indxp, iadx, h3d_data, fncontp, ftcontp)
Definition i2for28p_pen.F:45
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