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