OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
r3def3.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!|| r3def3 ../engine/source/elements/spring/r3def3.F
25!||--- called by ------------------------------------------------------
26!|| rforc3 ../engine/source/elements/spring/rforc3.F
27!||--- calls -----------------------------------------------------
28!|| redef3 ../engine/source/elements/spring/redef3.F90
29!|| table_interp ../engine/source/tools/curve/table_tools.F
30!||--- uses -----------------------------------------------------
31!|| element_mod ../common_source/modules/elements/element_mod.F90
32!|| interface_table_mod ../engine/share/modules/table_mod.F
33!|| python_funct_mod ../common_source/modules/python_mod.F90
34!|| redef3_mod ../engine/source/elements/spring/redef3.F90
35!|| table_mod ../engine/share/modules/table_mod.F
36!||====================================================================
37 SUBROUTINE r3def3( PYTHON,
38 1 GEO, F, AL0, E,
39 2 DL, NPF, TF, OFF,
40 3 DPL, FEP, DPX2, DFS,
41 4 V, IXR, DF, ANIM,
42 5 IPOS, IGEO, AL0_ERR,
43 6 X1DP, X2DP, X3DP, YIELD,
44 7 TABLE, INIFRIC, NGL, MGN,
45 8 EX, EY, EZ, XK,
46 9 XM, XC, AK, EX2,
47 A EY2, EZ2, NC1, NC2,
48 B NC3, NUVAR, UVAR, DL0,
49 C NEL, NFT, STF, SANIN,
50 D IRESP, SNPC)
51C-----------------------------------------------
52C M o d u l e s
53C-----------------------------------------------
54 USE redef3_mod
55 USE table_mod
56 USE interface_table_mod
57 USE python_funct_mod
58 use element_mod , only : nixr
59C-----------------------------------------------
60C I m p l i c i t T y p e s
61C-----------------------------------------------
62#include "implicit_f.inc"
63#include "comlock.inc"
64C-----------------------------------------------
65C G l o b a l P a r a m e t e r s
66C-----------------------------------------------
67#include "mvsiz_p.inc"
68C-----------------------------------------------
69C C o m m o n B l o c k s
70C-----------------------------------------------
71#include "units_c.inc"
72#include "scr17_c.inc"
73#include "param_c.inc"
74#include "com04_c.inc"
75#include "com08_c.inc"
76#include "scr14_c.inc"
77#include "com01_c.inc"
78#include "impl1_c.inc"
79C-----------------------------------------------
80C D u m m y A r g u m e n t s
81C-----------------------------------------------
82 type(python_), intent(inout) :: PYTHON
83 INTEGER, INTENT(IN) :: SANIN
84 INTEGER, INTENT(IN) :: STF !< Size of TF
85 INTEGER, INTENT(IN) :: NFT
86 INTEGER, INTENT(IN) :: NEL
87 INTEGER, INTENT(IN) :: IRESP !< Single Presision Flag
88 INTEGER, INTENT(IN) :: SNPC !< Size of NPF
89 INTEGER NPF(SNPC),IXR(NIXR,*),IGEO(NPROPGI,*),NGL(*),MGN(*),
90 . NC1(*),NC2(*),NC3(*),NUVAR
91C REAL
92 my_real
93 . GEO(NPROPG,*), F(*), AL0(*), E(*), DL(*), TF(STF), OFF(*),
94 . DPL(*), FEP(*), DPX2(*),DFS(*),V(3,*),DF(*),ANIM(*),
95 . ipos(*),al0_err(*),yield(*),inifric(*),
96 . ex(mvsiz),ey(mvsiz),ez(mvsiz),xk(mvsiz),
97 . xm(mvsiz),xc(mvsiz),ak(mvsiz),ex2(mvsiz),ey2(mvsiz),ez2(mvsiz),
98 . uvar(nuvar,*),dl0(*)
99 DOUBLE PRECISION X1DP(3,*),X2DP(3,*),X3DP(3,*)
100 TYPE(ttable) TABLE(*)
101 TARGET :: uvar
102C-----------------------------------------------
103c L o c a l V a r i a b l e s
104C-----------------------------------------------
105 INTEGER IFUNC2(MVSIZ),
106 . IECROU(MVSIZ), IFUNC(MVSIZ), IFV(MVSIZ),
107 . NINDX, INDX(MVSIZ),I,ILENG,IPID,J,
108 . IFUNC3(MVSIZ),ITABF(MVSIZ),IFRIC(MVSIZ)
109C REAL
110 my_real
111 . XL0(MVSIZ),DMN(MVSIZ),DMX(MVSIZ),
112 . dlold(mvsiz),b(mvsiz), d(mvsiz), dv(mvsiz),
113 . ff(mvsiz),lscale(mvsiz),ee(mvsiz),
114 . gf3(mvsiz),fric(mvsiz),xscalf(mvsiz),yscalf(mvsiz),
115 . f_min(mvsiz),f_max(mvsiz),epla(mvsiz)
116 my_real
117 . beta, mu,fmax, ddx ,vx21,vy21,vz21,vl21,
118 . vx32,vy32,vz32,vl32,not_used,not_used2(2)
119 my_real,
120 . DIMENSION(1) :: xx
121 DOUBLE PRECISION EXDP(MVSIZ) ,EYDP(MVSIZ) ,EZDP(MVSIZ),
122 . EX2DP(MVSIZ),EY2DP(MVSIZ),EZ2DP(MVSIZ),
123 . al1dp(mvsiz),al2dp(mvsiz),aldp(mvsiz),
124 . al0dp(mvsiz)
125 my_real :: max_slope(mvsiz)
126 my_real ,DIMENSION(:), POINTER :: xx_old
127 TARGET :: not_used2
128C-----------------------------------------------
129C
130 not_used = zero
131 not_used2 = zero
132C
133 DO i=1,nel
134 epla(i) = zero
135 ipid = mgn(i)
136 xm(i) =geo(1,ipid)
137 xk(i) =geo(2,ipid) ! K?
138 xc(i) =geo(3,ipid) ! C?
139 iecrou(i)=nint(geo(7,ipid))
140 ak(i) = geo(10,ipid)
141 b(i) = geo(11,ipid)
142 d(i) = geo(13,ipid)
143 ee(i) = geo(40,ipid)
144 gf3(i) = geo(132,ipid)
145 ff(i) = geo(18,ipid)
146 dmn(i) = geo(15,ipid)
147 dmx(i) = geo(16,ipid)
148 fric(i) = geo(17,ipid)
149 xscalf(i)= geo(20,ipid)
150 lscale(i)= geo(39,ipid)
151 ifunc(i) = igeo(101,ipid)
152 ifv(i) = igeo(102,ipid)
153 ifunc2(i)= igeo(103,ipid)
154 ifunc3(i)= igeo(119,ipid)
155 itabf(i) = igeo(201,ipid)
156 ifric(i) = igeo(202,ipid)
157 f_min(i) = geo(138,ipid)
158 f_max(i) = geo(139,ipid)
159 yscalf(i)= geo(140,ipid)
160 max_slope(i) = geo(141,ipid)
161 ENDDO
162C
163 DO i=1,nel
164 exdp(i)=x2dp(1,i)-x1dp(1,i)
165 eydp(i)=x2dp(2,i)-x1dp(2,i)
166 ezdp(i)=x2dp(3,i)-x1dp(3,i)
167 ex2dp(i)=x2dp(1,i)-x3dp(1,i)
168 ey2dp(i)=x2dp(2,i)-x3dp(2,i)
169 ez2dp(i)=x2dp(3,i)-x3dp(3,i)
170 dlold(i)=dl(i)
171 ENDDO
172C
173 IF (inispri /= 0 .and. tt == zero) THEN
174 DO i=1,nel
175 dlold(i)=dl0(i)
176 ENDDO
177 ENDIF
178C
179 DO i=1,nel
180 al1dp(i)=sqrt(exdp(i)*exdp(i)+eydp(i)*eydp(i)+ezdp(i)*ezdp(i))
181 al2dp(i)=sqrt(ex2dp(i)*ex2dp(i)+ey2dp(i)*ey2dp(i)+
182 . ez2dp(i)*ez2dp(i))
183 aldp(i) =al1dp(i) + al2dp(i)
184 al1dp(i)= max(al1dp(i),em15)
185 exdp(i) = exdp(i)/al1dp(i)
186 eydp(i) = eydp(i)/al1dp(i)
187 ezdp(i) = ezdp(i)/al1dp(i)
188 al2dp(i)= max(al2dp(i),em15)
189 ex2dp(i)= ex2dp(i)/al2dp(i)
190 ey2dp(i)= ey2dp(i)/al2dp(i)
191 ez2dp(i)= ez2dp(i)/al2dp(i)
192 ENDDO
193C
194 DO i=1,nel
195 ex(i) = exdp(i)
196 ey(i) = eydp(i)
197 ez(i) = ezdp(i)
198 ex2(i)= ex2dp(i)
199 ey2(i)= ey2dp(i)
200 ez2(i)= ez2dp(i)
201 ENDDO
202!
203 IF (inispri /= 0 .and. tt == zero) THEN
204 DO i=1,nel
205 xl0(i)= al0(i)
206! if not initialized length
207 IF (xl0(i) == zero) xl0(i) = aldp(i)
208 ENDDO
209 ENDIF
210!
211 IF (tt == zero) THEN
212 DO i=1,nel
213 al0(i)= aldp(i) ! cast double vers My_real
214 ENDDO
215 ENDIF
216C
217 IF (scodver >= 101) THEN
218 IF (tt == zero) THEN
219 DO i=1,nel
220 al0_err(i)=aldp(i)-al0(i) ! difference between double and my_real
221 ENDDO
222 ENDIF
223 ENDIF
224!
225 IF ( inispri /= 0 .and. tt == zero) THEN
226 DO i=1,nel
227 al0(i)= xl0(i)
228 ENDDO
229 ENDIF
230!
231 DO i=1,nel
232 al0dp(i)= al0(i) ! cast My_real en double
233 ENDDO
234!
235 IF (scodver >= 101) THEN
236 DO i=1,nel
237 al0dp(i)= al0dp(i)+al0_err(i) ! Al_dp must be recalculated as well in order to guarantee the absolute coh rence between al0_dp and al_dp
238 ENDDO
239 ENDIF
240C
241 IF (ismdisp > 0) THEN
242 DO i=1,nel
243 vx21= v(1,nc2(i)) - v(1,nc1(i))
244 vy21= v(2,nc2(i)) - v(2,nc1(i))
245 vz21= v(3,nc2(i)) - v(3,nc1(i))
246 vl21 = vx21*ex(i)+vy21*ey(i)+vz21*ez(i)
247 vx32= v(1,nc2(i)) - v(1,nc3(i))
248 vy32= v(2,nc2(i)) - v(2,nc3(i))
249 vz32= v(3,nc2(i)) - v(3,nc3(i))
250 vl32 = vx32*ex2(i)+vy32*ey2(i)+vz32*ez2(i)
251 dl(i)= dlold(i)+(vl21+vl32)*dt1
252 ENDDO
253 ELSE
254 DO i=1,nel
255 dl(i)= (aldp(i)-al0dp(i))
256 ENDDO
257 ENDIF !(ISMDISP>0) THEN
258C
259 DO i=1,nel
260 ileng=nint(geo(93,mgn(i)))
261 IF (ileng /= 0) THEN
262 xl0(i)=al0dp(i)
263 ELSE
264 xl0(i) = one
265 ENDIF
266 ENDDO
267C
268 IF (nuvar > 0) THEN
269 xx_old => uvar(1,1:nel)
270 ELSE
271 xx_old => not_used2
272 ENDIF
273 CALL redef3(python,
274 1 f, xk, dl, fep,
275 2 dlold, dpl, tf, npf,
276 3 xc, off, e, dpx2,
277 4 anim, anim_fe(11),ipos,
278 5 xl0, dmn, dmx, dv,
279 6 ff, lscale, ee, gf3,
280 7 ifunc3, yield, aldp, ak,
281 8 b, d, iecrou, ifunc,
282 9 ifv, ifunc2, epla, xx_old,
283 a nel, nft, stf, sanin, dt1,
284 b iresp, impl_s, idyna, snpc,
285 c max_slope=max_slope)
286 nindx = 0
287 DO i=1,nel
288 IF (off(i) == one .AND. dmx(i) /= zero .AND. dmn(i) /= zero) THEN
289 IF (dl(i) > dmx(i) .OR. dl(i) < dmn(i)) THEN
290 off(i)=zero
291 nindx = nindx + 1
292 indx(nindx) = i
293 idel7nok = 1
294 ENDIF
295 ENDIF
296 ENDDO
297C
298 DO j=1,nindx
299 i = indx(j)
300#include "lockon.inc"
301 WRITE(iout, 1000) ngl(i)
302 WRITE(istdo,1100) ngl(i),tt
303#include "lockoff.inc"
304 ENDDO
305C
306 DO i=1,nel
307 xm(i)=xm(i)*xl0(i)
308 xk(i)=xk(i)/xl0(i)
309C--- for time step compute adding +max slope of h
310C Should it be replaced by the
311C XC(I)=(XC(I)+GEO(141,MGN(I)))/XL0(I)
312 xc(i) = (xc(i) + max_slope(i)) / xl0(i)
313 ENDDO
314C
315 DO i=1,nel
316 IF (ifric(i) == 0) THEN ! old approach
317 IF (itabf(i) > zero) THEN
318 ddx = exdp(i) * (v(1,nc2(i)) - v(1,nc1(i)))
319 . + eydp(i) * (v(2,nc2(i)) - v(2,nc1(i)))
320 . + ezdp(i) * (v(3,nc2(i)) - v(3,nc1(i)))
321 . + ex2dp(i)* (v(1,nc3(i)) - v(1,nc2(i)))
322 . + ey2dp(i)* (v(2,nc3(i)) - v(2,nc2(i)))
323 . + ez2dp(i)* (v(3,nc3(i)) - v(3,nc2(i)))
324 dfs(i)= dfs(i) + ddx * dt1 * xk(i)
325 beta = pi - acos(exdp(i)*ex2dp(i)+eydp(i)*ey2dp(i)+
326 . ezdp(i)*ez2dp(i))
327 xx(1) = abs(dfs(i)*xscalf(i))
328 CALL table_interp(table(itabf(i)),xx,mu)
329 mu = mu*yscalf(i)
330 fmax = max(zero,f(i) * tanh(half*mu*beta))
331 fmax = min(fmax,abs(dfs(i)))
332 dfs(i)= sign(fmax,dfs(i))
333 df(i) = dfs(i)
334 ELSEIF (fric(i) > zero) THEN
335 ddx = exdp(i) * (v(1,nc2(i)) - v(1,nc1(i)))
336 . + eydp(i) * (v(2,nc2(i)) - v(2,nc1(i)))
337 . + ezdp(i) * (v(3,nc2(i)) - v(3,nc1(i)))
338 . + ex2dp(i)* (v(1,nc3(i)) - v(1,nc2(i)))
339 . + ey2dp(i)* (v(2,nc3(i)) - v(2,nc2(i)))
340 . + ez2dp(i)* (v(3,nc3(i)) - v(3,nc2(i)))
341 dfs(i)= dfs(i) + ddx * dt1 * xk(i)
342 beta = pi - acos(exdp(i)*ex2dp(i)+eydp(i)*ey2dp(i)+
343 . ezdp(i)*ez2dp(i))
344 fmax = max(zero,f(i) * tanh(half*fric(i)*beta))
345 fmax = min(fmax,abs(dfs(i)))
346 dfs(i)= sign(fmax,dfs(i))
347 df(i) = dfs(i)
348C df=dfs for restart compatibility if fric=0
349C and DFS = EV (NB9) not existing before 23F and old 23F
350 ELSE
351 df(i)= zero
352 ENDIF ! IF (ITABF(I) > ZERO)
353 ELSEIF (ifric(i) > 0) THEN ! new approach
354 IF (itabf(i) > zero) THEN
355 ddx = exdp(i) * (v(1,nc2(i)) - v(1,nc1(i)))
356 . + eydp(i) * (v(2,nc2(i)) - v(2,nc1(i)))
357 . + ezdp(i) * (v(3,nc2(i)) - v(3,nc1(i)))
358 . + ex2dp(i)* (v(1,nc3(i)) - v(1,nc2(i)))
359 . + ey2dp(i)* (v(2,nc3(i)) - v(2,nc2(i)))
360 . + ez2dp(i)* (v(3,nc3(i)) - v(3,nc2(i)))
361 dfs(i)= dfs(i) + ddx * dt1 * xk(i)
362 beta = pi - acos(exdp(i)*ex2dp(i)+eydp(i)*ey2dp(i)+
363 . ezdp(i)*ez2dp(i))
364cc XX(1) = ABS(DFS(I)*XSCALF(I)) --> wrong only positive (abscisa) side function used
365 xx(1) = dfs(i)*xscalf(i)
366 CALL table_interp(table(itabf(i)),xx,mu)
367 mu = mu*yscalf(i)
368C---
369C force limit reached (positive and negative)
370C---
371cc IF (F_MIN(I) /= ZERO .OR. F_MAX(I) /= 0) THEN
372cc IF (DFS(I) <= F_MIN(I) .OR. DFS(I) >= F_MAX(I)) MU = FRIC(I)
373cc ELSE ! no limit force defined
374cc ENDIF
375 IF (dfs(i) <= f_min(i) .OR. dfs(i) >= f_max(i)) inifric(i) = one
376 IF (inifric(i) == one) mu = fric(i) ! reset friction to initial
377C---
378 fmax = max(zero,f(i) * tanh(half*mu*beta))
379 fmax = min(fmax,abs(dfs(i)))
380 dfs(i)= sign(fmax,dfs(i))
381 df(i) = dfs(i)
382 ELSEIF (fric(i) > zero .AND. itabf(i) == 0) THEN
383 ddx = exdp(i) * (v(1,nc2(i)) - v(1,nc1(i)))
384 . + eydp(i) * (v(2,nc2(i)) - v(2,nc1(i)))
385 . + ezdp(i) * (v(3,nc2(i)) - v(3,nc1(i)))
386 . + ex2dp(i)* (v(1,nc3(i)) - v(1,nc2(i)))
387 . + ey2dp(i)* (v(2,nc3(i)) - v(2,nc2(i)))
388 . + ez2dp(i)* (v(3,nc3(i)) - v(3,nc2(i)))
389 dfs(i)= dfs(i) + ddx * dt1 * xk(i)
390 beta = pi - acos(exdp(i)*ex2dp(i)+eydp(i)*ey2dp(i)+
391 . ezdp(i)*ez2dp(i))
392 fmax = max(zero,f(i) * tanh(half*fric(i)*beta))
393 fmax = min(fmax,abs(dfs(i)))
394 dfs(i)= sign(fmax,dfs(i))
395 df(i) = dfs(i)
396 ELSE
397 df(i)= zero
398 ENDIF ! IF (ITABF(I) > ZERO)
399 ENDIF ! IF (IFRIC(I) == 0)
400 ENDDO ! DO I=1,NEL
401C---
402 1000 FORMAT(1x,'-- RUPTURE OF SPRING ELEMENT NUMBER ',i10)
403 1100 FORMAT(1x,'-- RUPTURE OF SPRING ELEMENT :',i10,' AT TIME :',g11.4)
404C---
405 RETURN
406 END
interface_table_mod::table_interp
Definition table_mod.F:243
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
interface_table_mod
Definition table_mod.F:239
table_mod
Definition table_mod.F:113
r3def3
subroutine r3def3(python, geo, f, al0, e, dl, npf, tf, off, dpl, fep, dpx2, dfs, v, ixr, df, anim, ipos, igeo, al0_err, x1dp, x2dp, x3dp, yield, table, inifric, ngl, mgn, ex, ey, ez, xk, xm, xc, ak, ex2, ey2, ez2, nc1, nc2, nc3, nuvar, uvar, dl0, nel, nft, stf, sanin, iresp, snpc)
Definition r3def3.F:51
table_mod::ttable
Definition table_mod.F:126