OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i2loceq.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!|| i2loceq ../common_source/interf/i2loceq.F
25!||--- called by ------------------------------------------------------
26!|| i2_dtn_25 ../starter/source/interfaces/inter3d1/i2_dtn.F
27!|| i2_dtn_27_pen ../starter/source/interfaces/inter3d1/i2_dtn_27.F
28!|| i2for25 ../engine/source/interfaces/interf/i2for25.F
29!|| i2for25p ../engine/source/interfaces/interf/i2for25p.F
30!|| i2for26 ../engine/source/interfaces/interf/i2for26.F
31!|| i2for26p ../engine/source/interfaces/interf/i2for26p.F
32!|| i2for27_pen ../engine/source/interfaces/interf/i2for27_pen.F
33!|| i2for27p_pen ../engine/source/interfaces/interf/i2for27p_pen.F
34!||--- calls -----------------------------------------------------
35!|| inv3 ../engine/source/elements/joint/rskew33.F
36!||====================================================================
37 SUBROUTINE i2loceq(
38 . NIR ,RS ,RX ,RY ,RZ ,
39 . FMX ,FMY ,FMZ ,H ,STIFM )
40C-----------------------------------------------
41C I m p l i c i t T y p e s
42C-----------------------------------------------
43#include "implicit_f.inc"
44C-----------------------------------------------
45C D u m m y A r g u m e n t s
46C-----------------------------------------------
47 INTEGER NIR
48C REAL
49 my_real
50 . rs(3),rx(4),ry(4),rz(4),fmx(4),fmy(4),fmz(4),h(4),stifm
51C-----------------------------------------------
52C C o m m o n B l o c k s
53C-----------------------------------------------
54#include "units_c.inc"
55#include "param_c.inc"
56#include "com01_c.inc"
57C-----------------------------------------------
58C L o c a l V a r i a b l e s
59C-----------------------------------------------
60 INTEGER J
61 my_real MX,MY,MZ,D1,D2,DD,RAX,RBX,RAY,RBY, MXY,MYX,MYZ,MZY,MXZ,MZX,MXX,MYY,MZZ,RXX,RYY,FXX,FYY,FA,FB
62 my_real rsx(4),rsy(4),rsz(4),df(nir),mat(3,3),mati(3,3),bx,by,bz
63C=======================================================================
64C Moment equilibrium in LOCAL coordinates
65C-------------------------------------------------
66 rsx(1) = rx(1) - rs(1)
67 rsy(1) = ry(1) - rs(2)
68 rsz(1) = rz(1) - rs(3)
69 rsx(2) = rx(2) - rs(1)
70 rsy(2) = ry(2) - rs(2)
71 rsz(2) = rz(2) - rs(3)
72 rsx(3) = rx(3) - rs(1)
73 rsy(3) = ry(3) - rs(2)
74 rsz(3) = rz(3) - rs(3)
75 rsx(4) = rx(4) - rs(1)
76 rsy(4) = ry(4) - rs(2)
77 rsz(4) = rz(4) - rs(3)
78c
79 mx = zero
80 my = zero
81 mz = zero
82 DO j=1,nir
83 mx = mx + rsy(j)*fmz(j) - rsz(j)*fmy(j)
84 my = my + rsz(j)*fmx(j) - rsx(j)*fmz(j)
85 mz = mz + rsx(j)*fmy(j) - rsy(j)*fmx(j)
86 ENDDO
87C
88 bx= rsx(1)*h(1)+rsx(2)*h(2)+rsx(3)*h(3)+rsx(4)*h(4)
89 by= rsy(1)*h(1)+rsy(2)*h(2)+rsy(3)*h(3)+rsy(4)*h(4)
90 bz= rsz(1)*h(1)+rsz(2)*h(2)+rsz(3)*h(3)+rsz(4)*h(4)
91C-------------------------------------------------
92 IF (nir == 4) THEN
93 rxx = rx(2)+rx(3)-rx(1)-rx(4)
94 ryy = ry(3)+ry(4)-ry(1)-ry(2)
95c RXX = ABS(RXX)
96c RYY = ABS(RYY)
97C--- Moment Z
98c
99 mxy = rsx(1)*fmy(1)+rsx(2)*fmy(2)+rsx(3)*fmy(3)+rsx(4)*fmy(4)
100 myx = rsy(1)*fmx(1)+rsy(2)*fmx(2)+rsy(3)*fmx(3)+rsy(4)*fmx(4)
101 mzz = mxy - myx
102 fyy = mxy/rxx
103 fxx = -myx/ryy
104C
105 stifm=max(abs(bx) / abs(rxx),abs(by) / abs(ryy))
106C
107 fmx(1) = fmx(1) - fxx
108 fmx(2) = fmx(2) - fxx
109 fmx(3) = fmx(3) + fxx
110 fmx(4) = fmx(4) + fxx
111
112 fmy(1) = fmy(1) + fyy
113 fmy(2) = fmy(2) - fyy
114 fmy(3) = fmy(3) - fyy
115 fmy(4) = fmy(4) + fyy
116c
117c--- Moments Mx, MY
118c
119 myz = rsy(1)*fmz(1)+rsy(2)*fmz(2)+rsy(3)*fmz(3)+rsy(4)*fmz(4)
120 mzy = rsz(1)*fmy(1)+rsz(2)*fmy(2)+rsz(3)*fmy(3)+rsz(4)*fmy(4)
121 mxx = myz - mzy
122 mzx = rsz(1)*fmx(1)+rsz(2)*fmx(2)+rsz(3)*fmx(3)+rsz(4)*fmx(4)
123 mxz = rsx(1)*fmz(1)+rsx(2)*fmz(2)+rsx(3)*fmz(3)+rsx(4)*fmz(4)
124 myy = mzx - mxz
125c
126 rax = rx(1) - rx(3)
127 rbx = rx(4) - rx(2)
128 ray = ry(1) - ry(3)
129 rby = ry(4) - ry(2)
130 d1 = -mxx*rbx - myy*rby
131 d2 = mxx*rax + myy*ray
132 dd = ray*rbx - rax*rby
133 fa = d1 / dd
134 fb = d2 / dd
135C
136 stifm=max(stifm,
137 . sqrt((by*by+bz*bz)*rbx*rbx+(bz*bz+bx*bx)*rby*rby)/abs(dd),
138 . sqrt((by*by+bz*bz)*rax*rax+(bz*bz+bx*bx)*ray*ray)/abs(dd))
139C
140 fmz(1) = fmz(1) + fa
141 fmz(2) = fmz(2) - fb
142 fmz(3) = fmz(3) - fa
143 fmz(4) = fmz(4) + fb
144c
145
146 ELSEIF (nir == 3) THEN
147 rxx = rx(2)+rx(3)-two*rx(1)
148 ryy = two*ry(3)-ry(1)-ry(2)
149C--- Moment Z
150C
151 mxy = rsx(1)*fmy(1)+rsx(2)*fmy(2)+rsx(3)*fmy(3)
152 myx = rsy(1)*fmx(1)+rsy(2)*fmx(2)+rsy(3)*fmx(3)
153 mzz = mxy - myx
154
155 fyy = -mxy/rxx
156 fxx = -myx/ryy
157C
158 stifm=max(abs(bx) / abs(rxx),abs(by) / abs(ryy))
159C
160 fmx(1) = fmx(1) - fxx
161 fmx(2) = fmx(2) - fxx
162 fmx(3) = fmx(3) + two * fxx
163
164 fmy(1) = fmy(1) - two * fyy
165 fmy(2) = fmy(2) + fyy
166 fmy(3) = fmy(3) + fyy
167c
168 mx=mx+(two*rsz(1)-rsz(2)-rsz(3))*fyy
169 my=my+(two*rsz(3)-rsz(1)-rsz(2))*fxx
170c--- Moments Mx, MY
171c
172 mat(1,1) = rsy(1)
173 mat(1,2) = rsy(2)
174 mat(1,3) = rsy(3)
175 mat(2,1) = rsx(1)
176 mat(2,2) = rsx(2)
177 mat(2,3) = rsx(3)
178 mat(3,1) = one
179 mat(3,2) = one
180 mat(3,3) = one
181 CALL inv3(mat,mati)
182 df(1) = -mati(1,1)*mx + mati(1,2)*my
183 df(2) = -mati(2,1)*mx + mati(2,2)*my
184 df(3) = -mati(3,1)*mx + mati(3,2)*my
185 DO j=1,nir
186 fmz(j) = fmz(j) + df(j)
187 ENDDO
188C
189 bx= rsx(1)*h(1)+rsx(2)*h(2)+rsx(3)*h(3)
190 by= rsy(1)*h(1)+rsy(2)*h(2)+rsy(3)*h(3)
191 bz= rsz(1)*h(1)+rsz(2)*h(2)+rsz(3)*h(3)
192 stifm= max(stifm,sqrt(bx*bx+by*by+bz*bz)*sqrt(max(
193 . mati(1,1)*mati(1,1)+mati(2,1)*mati(2,1)+mati(3,1)*mati(3,1),
194 . mati(1,2)*mati(1,2)+mati(2,2)*mati(2,2)+mati(3,2)*mati(3,2),
195 . mati(1,3)*mati(1,3)+mati(2,3)*mati(2,3)+mati(3,3)*mati(3,3))))
196c MX = ZERO
197c MY = ZERO
198c MZ = ZERO
199c DO J=1,NIR
200c MX = MX + RSY(J)*FMZ(J) - RSZ(J)*FMY(J)
201c MY = MY + RSZ(J)*FMX(J) - RSX(J)*FMZ(J)
202c MZ = MZ + RSX(J)*FMY(J) - RSY(J)*FMX(J)
203c ENDDO
204C print *,mx,my,mz
205C-----
206 ENDIF
207C-----------
208 RETURN
209 END
210!||====================================================================
211!|| i2loceq_27 ../common_source/interf/i2loceq.F
212!||--- called by ------------------------------------------------------
213!|| i2_dtn_27_cin ../starter/source/interfaces/inter3d1/i2_dtn_27.F
214!|| i2for27_cin ../engine/source/interfaces/interf/i2for27_cin.F
215!|| i2for27p_cin ../engine/source/interfaces/interf/i2for27p_cin.F
216!||--- calls -----------------------------------------------------
217!|| inv3 ../engine/source/elements/joint/rskew33.F
218!||====================================================================
219 SUBROUTINE i2loceq_27(
220 . NIR ,RS ,RX ,RY ,RZ ,
221 . FMX ,FMY ,FMZ ,H ,STIFM ,
222 . MXS ,MYS ,MZS ,STIFMR ,BETAX ,
223 . BETAY)
224C-----------------------------------------------
225C I m p l i c i t T y p e s
226C-----------------------------------------------
227#include "implicit_f.inc"
228C-----------------------------------------------
229C D u m m y A r g u m e n t s
230C-----------------------------------------------
231 INTEGER NIR
232C REAL
233 my_real
234 . RS(3),RX(4),RY(4),RZ(4),FMX(4),FMY(4),FMZ(4),H(4),STIFM,MXS,MYS,MZS,STIFMR,
235 . BETAX,BETAY
236C-----------------------------------------------
237C C o m m o n B l o c k s
238C-----------------------------------------------
239#include "units_c.inc"
240#include "param_c.inc"
241#include "com01_c.inc"
242C-----------------------------------------------
243C L o c a l V a r i a b l e s
244C-----------------------------------------------
245 INTEGER J
246 my_real MX,MY,MZ,D1,D2,DD,RAX,RBX,RAY,RBY,MXY,MYX,MYZ,MZY,MXZ,MZX,MXX,MYY,MZZ,RXX,RYY,FXX,FYY,FA,FB
247 my_real RSX(4),RSY(4),RSZ(4),DF(NIR),MAT(3,3),MATI(3,3),BX,BY,BZ,ALPHA
248C=======================================================================
249C Moment equilibrium in LOCAL coordinates
250C-------------------------------------------------
251 rsx(1) = rx(1) - rs(1)
252 rsy(1) = ry(1) - rs(2)
253 rsz(1) = rz(1) - rs(3)
254 rsx(2) = rx(2) - rs(1)
255 rsy(2) = ry(2) - rs(2)
256 rsz(2) = rz(2) - rs(3)
257 rsx(3) = rx(3) - rs(1)
258 rsy(3) = ry(3) - rs(2)
259 rsz(3) = rz(3) - rs(3)
260 rsx(4) = rx(4) - rs(1)
261 rsy(4) = ry(4) - rs(2)
262 rsz(4) = rz(4) - rs(3)
263c
264 mx = zero
265 my = zero
266 mz = zero
267 DO j=1,nir
268 mx = mx + rsy(j)*fmz(j) - rsz(j)*fmy(j)
269 my = my + rsz(j)*fmx(j) - rsx(j)*fmz(j)
270 mz = mz + rsx(j)*fmy(j) - rsy(j)*fmx(j)
271 ENDDO
272C
273 bx= rsx(1)*h(1)+rsx(2)*h(2)+rsx(3)*h(3)+rsx(4)*h(4)
274 by= rsy(1)*h(1)+rsy(2)*h(2)+rsy(3)*h(3)+rsy(4)*h(4)
275 bz= rsz(1)*h(1)+rsz(2)*h(2)+rsz(3)*h(3)+rsz(4)*h(4)
276C
277C-------------------------------------------------
278 IF (nir == 4) THEN
279 rxx = rx(2)+rx(3)-rx(1)-rx(4)
280 ryy = ry(3)+ry(4)-ry(1)-ry(2)
281c RXX = ABS(RXX)
282c RYY = ABS(RYY)
283C--- Moment Z
284c
285 alpha = abs(rxx)/(abs(rxx)+abs(ryy))
286C
287C RQ : ALPHA is a coeff used for a optimised repartition of moment MZS on FXX and FYY
288C
289 mxy = rsx(1)*fmy(1)+rsx(2)*fmy(2)+rsx(3)*fmy(3)+rsx(4)*fmy(4)-alpha*mzs
290 myx = rsy(1)*fmx(1)+rsy(2)*fmx(2)+rsy(3)*fmx(3)+rsy(4)*fmx(4)+(one-alpha)*mzs
291 mzz = mxy - myx
292C
293C RQ : BETAX = BETAY = 1 in standard situation - moment mz must be transferred either only by FXX or FXX if switch to 1D formulation
294C standard situation -> FXX = -MYX/RYY and FYY = MXY/RXX
295C 1D formulation in X -> FXX = 0 and FYY = MZ / RXX
296C 1D formulation in Y -> FYY = 0 and FXX = MZ / RYY
297C
298 fyy = (betay*mxy-(one-betax)*myx)/rxx
299 fxx = ((one-betay)*mxy-betax*myx)/ryy
300C
301 stifm=max(sqrt((one-betax)*by*by+ betay*bx*bx)/abs(rxx),sqrt(betax*by*by+ (one-betay)*bx*bx)/abs(ryy))
302 stifmr=one/(abs(rxx)+abs(ryy))
303C
304 fmx(1) = fmx(1) - fxx
305 fmx(2) = fmx(2) - fxx
306 fmx(3) = fmx(3) + fxx
307 fmx(4) = fmx(4) + fxx
308
309 fmy(1) = fmy(1) + fyy
310 fmy(2) = fmy(2) - fyy
311 fmy(3) = fmy(3) - fyy
312 fmy(4) = fmy(4) + fyy
313c
314c--- Moments Mx, MY
315C
316C RQ : BETAX = BETAY = 1 in standard situation - used to impose that mx=0 or my=0 if switch to 1D formulation
317c
318 myz = rsy(1)*fmz(1)+rsy(2)*fmz(2)+rsy(3)*fmz(3)+rsy(4)*fmz(4)
319 mzy = rsz(1)*fmy(1)+rsz(2)*fmy(2)+rsz(3)*fmy(3)+rsz(4)*fmy(4)
320 mxx = betax*(myz - mzy - mxs)
321 mzx = rsz(1)*fmx(1)+rsz(2)*fmx(2)+rsz(3)*fmx(3)+rsz(4)*fmx(4)
322 mxz = rsx(1)*fmz(1)+rsx(2)*fmz(2)+rsx(3)*fmz(3)+rsx(4)*fmz(4)
323 myy = betay*(mzx - mxz - mys)
324c
325 rax = rx(1) - rx(3)
326 rbx = rx(4) - rx(2)
327 ray = ry(1) - ry(3)
328 rby = ry(4) - ry(2)
329 d1 = -mxx*rbx - myy*rby
330 d2 = mxx*rax + myy*ray
331 dd = ray*rbx - rax*rby
332 fa = d1 / dd
333 fb = d2 / dd
334C
335 stifm=max(stifm,
336 . sqrt((by*by+bz*bz)*rbx*rbx*betax+(bz*bz+bx*bx)*rby*rby*betay)/abs(dd),
337 . sqrt((by*by+bz*bz)*rax*rax*betax+(bz*bz+bx*bx)*ray*ray*betay)/abs(dd))
338C
339 stifmr=max(stifmr,(abs(rbx*betax)+abs(rby*betay))/abs(dd),(abs(rax*betax)+abs(ray*betay))/abs(dd))
340C
341 fmz(1) = fmz(1) + fa
342 fmz(2) = fmz(2) - fb
343 fmz(3) = fmz(3) - fa
344 fmz(4) = fmz(4) + fb
345c
346
347 ELSEIF (nir == 3) THEN
348 rxx = rx(2)+rx(3)-two*rx(1)
349 ryy = two*ry(3)-ry(1)-ry(2)
350C
351C RQ : ALPHA is a coeff used for a optimised repartition of moment MZS on FXX and FYY
352C
353 alpha = abs(rxx)/(abs(rxx)+abs(ryy))
354C
355C--- Moment Z
356C
357C
358 mxy = rsx(1)*fmy(1)+rsx(2)*fmy(2)+rsx(3)*fmy(3)-alpha*mzs
359 myx = rsy(1)*fmx(1)+rsy(2)*fmx(2)+rsy(3)*fmx(3)+(one-alpha)*mzs
360 mzz = mxy - myx
361C
362C RQ : BETAX = BETAY = 1 in standard situation - moment mz must be transferred either only by FXX or FXX if switch to 1D formulation
363C standard situation -> FXX = -MYX/RYY and FYY = MXY/RXX
364C 1D formulation in X -> FXX = 0 and FYY = MZ / RXX
365C 1D formulation in Y -> FYY = 0 and FXX = MZ / RYY
366C
367 fyy = -(betay*mxy-(one-betax)*myx)/rxx
368 fxx = ((one-betay)*mxy-betax*myx)/ryy
369C
370 stifm=max(sqrt((one-betax)*by*by+ betay*bx*bx)/abs(rxx),sqrt(betax*by*by+ (one-betay)*bx*bx)/abs(ryy))
371 stifmr=one/(abs(rxx)+abs(ryy))
372C
373 fmx(1) = fmx(1) - fxx
374 fmx(2) = fmx(2) - fxx
375 fmx(3) = fmx(3) + two * fxx
376
377 fmy(1) = fmy(1) - two * fyy
378 fmy(2) = fmy(2) + fyy
379 fmy(3) = fmy(3) + fyy
380c
381c--- Moments Mx, MY
382C
383C RQ : BETAX = BETAY = 1 in standard situation - used to impose that mx=0 or my=0 if switch to 1D formulation
384c
385 mx=betax*(mx+(two*rsz(1)-rsz(2)-rsz(3))*fyy - mxs)
386 my=betay*(my+(two*rsz(3)-rsz(1)-rsz(2))*fxx - mys)
387c
388 mat(1,1) = rsy(1)
389 mat(1,2) = rsy(2)
390 mat(1,3) = rsy(3)
391 mat(2,1) = rsx(1)
392 mat(2,2) = rsx(2)
393 mat(2,3) = rsx(3)
394 mat(3,1) = one
395 mat(3,2) = one
396 mat(3,3) = one
397 CALL inv3(mat,mati)
398 df(1) = -mati(1,1)*mx + mati(1,2)*my
399 df(2) = -mati(2,1)*mx + mati(2,2)*my
400 df(3) = -mati(3,1)*mx + mati(3,2)*my
401 DO j=1,nir
402 fmz(j) = fmz(j) + df(j)
403 ENDDO
404C
405 bx= rsx(1)*h(1)+rsx(2)*h(2)+rsx(3)*h(3)
406 by= rsy(1)*h(1)+rsy(2)*h(2)+rsy(3)*h(3)
407 bz= rsz(1)*h(1)+rsz(2)*h(2)+rsz(3)*h(3)
408 stifm= max(stifm,sqrt(bx*bx+by*by+bz*bz)*sqrt(max(
409 . betax*(mati(1,1)*mati(1,1)+mati(2,1)*mati(2,1)+mati(3,1)*mati(3,1)),
410 . betay*(mati(1,2)*mati(1,2)+mati(2,2)*mati(2,2)+mati(3,2)*mati(3,2)),
411 . mati(1,3)*mati(1,3)+mati(2,3)*mati(2,3)+mati(3,3)*mati(3,3))))
412C
413 stifmr=max(stifmr,abs(mati(1,1))+abs(mati(1,2)),abs(mati(2,1))+abs(mati(2,2)),abs(mati(3,1))+abs(mati(3,2)))
414C
415C-----
416 ENDIF
417C-----------
418 RETURN
419 END
420
my_real
#define my_real
Definition cppsort.cpp:32
i2loceq
subroutine i2loceq(nir, rs, rx, ry, rz, fmx, fmy, fmz, h, stifm)
Definition i2loceq.F:40
i2loceq_27
subroutine i2loceq_27(nir, rs, rx, ry, rz, fmx, fmy, fmz, h, stifm, mxs, mys, mzs, stifmr, betax, betay)
Definition i2loceq.F:224
inv3
subroutine inv3(a, b)
Definition inv3.F:29
max
#define max(a, b)
Definition macros.h:21