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