OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
ruser44.F File Reference
#include "implicit_f.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine ruser44 (nel, iout, iprop, uvar, nuvar, fx, fy, fz, xmom, ymom, zmom, e, off, stifm, stifr, viscm, viscr, mass, xiner, dt, xl, vx, ry1, rz1, rx, ry2, rz2, fr_wave_e)

Function/Subroutine Documentation

◆ ruser44()

subroutine ruser44 ( integer nel,
integer iout,
integer iprop,
uvar,
integer nuvar,
fx,
fy,
fz,
xmom,
ymom,
zmom,
e,
off,
stifm,
stifr,
viscm,
viscr,
mass,
xiner,
dt,
xl,
vx,
ry1,
rz1,
rx,
ry2,
rz2,
fr_wave_e )

Definition at line 35 of file ruser44.F.

41C-------------------------------------------------------------------------
42C Crushable frame spring property
43C----------+---------+---+---+--------------------------------------------
44C VAR | SIZE |TYP| RW| DEFINITION
45C----------+---------+---+---+--------------------------------------------
46C IOUT | 1 | I | R | OUTPUT FILE UNIT (L00 file)
47C IPROP | 1 | I | R | PROPERTY NUMBER
48C----------+---------+---+---+--------------------------------------------
49C XL | NEL | F | R | ELEMENT LENGTH
50C----------+---------+---+---+--------------------------------------------
51C UVAR |NUVAR*NEL| F |R/W| USER ELEMENT VARIABLES
52C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
53C-------------------------------------------------------------------------
54C FUNCTION
55C-------------------------------------------------------------------------
56C INTEGER II = GET_U_PNU(I,IP,KK)
57C IFUNCI = GET_U_PNU(I,IP,KFUNC)
58C IPROPI = GET_U_PNU(I,IP,KFUNC)
59C IMATI = GET_U_PNU(I,IP,KMAT)
60C I : VARIABLE INDEX(1 for first variable,...)
61C IP : PROPERTY NUMBER
62C KK : PARAMETER KFUNC,KMAT,KPROP
63C THIS FUNCTION RETURN THE USER STORED FUNCTION(IF KK=KFUNC),
64C MATERIAL(IF KK=KMAT) OR PROPERTY(IF KK=KPROP) NUMBER.
65C SEE LECG29 FOR CORRESPONDING ID STORAGE.
66C-------------------------------------------------------------------------
67C INTEGER IFUNCI = GET_U_MNU(I,IM,KFUNC)
68C I : VARIABLE INDEX(1 for first function)
69C IM : MATERIAL NUMBER
70C KFUNC : ONLY FUNCTION ARE YET AVAILABLE.
71C THIS FUNCTION RETURN THE USER STORED FUNCTION NUMBER(function
72C referred by users materials).
73C SEE LECM29 FOR CORRESPONDING ID STORAGE.
74C-------------------------------------------------------------------------
75C my_real PARAMI = GET_U_GEO(I,IP)
76C I : PARAMETER INDEX(1 for first parameter,...)
77C IP : PROPERTY NUMBER
78C THIS FUNCTION RETURN THE USER GEOMETRY PARAMETERS
79C-------------------------------------------------------------------------
80C my_real PARAMI = GET_U_MAT(I,IM)
81C I : PARAMETER INDEX(1 for first parameter,...)
82C IM : MATERIAL NUMBER
83C THIS FUNCTION RETURN THE USER MATERIAL PARAMETERS
84C NOTE: GET_U_MAT(0,IMAT) RETURN THE DENSITY
85C-------------------------------------------------------------------------
86C INTEGER PID = GET_U_PID(IP)
87C IP : PROPERTY NUMBER
88C THIS FUNCTION RETURN THE USER PROPERTY ID CORRESPONDING TO
89C USER PROPERTY NUMBER IP.
90C-------------------------------------------------------------------------
91C INTEGER MID = GET_U_MID(IM)
92C IM : MATERIAL NUMBER
93C THIS FUNCTION RETURN THE USER MATERIAL ID CORRESPONDING TO
94C USER MATERIAL NUMBER IM.
95C-------------------------------------------------------------------------
96C my_real Y = GET_U_FUNC(IFUNC,X,DYDX)
97C IFUNC : function number obtained by
98C IFUNC = GET_U_MNU(I,IM,KFUNC) or IFUNC = GET_U_PNU(I,IP,KFUNC)
99C X : X value
100C DYDX : slope dY/dX
101C THIS FUNCTION RETURN Y(X)
102C-------------------------------------------------------------------------
103C I m p l i c i t T y p e s
104C-----------------------------------------------
105#include "implicit_f.inc"
106C----------------------------------------------------------
107C D u m m y A r g u m e n t s a n d F u n c t i o n
108C----------------------------------------------------------
109 INTEGER IOUT,NEL,NUVAR,IPROP,ICO,
110 . GET_U_PNU,GET_U_PID,GET_U_MID,GET_U_MNU,
111 . KFUNC,KMAT,KPROP
112 my_real
113 . uvar(nuvar,*),dt ,
114 . fx(*), fy(*), fz(*), e(*), vx(*),mass(*) ,xiner(*),
115 . ry1(*), rz1(*), off(*), xmom(*), ymom(*),
116 . zmom(*), rx(*), ry2(*), rz2(*),xl(*),
117 . stifm(*) ,stifr(*) , viscm(*) ,viscr(*) ,fr_wave_e(*) ,
118 . get_u_mat, get_u_geo, get_u_func
119 EXTERNAL get_u_mnu,get_u_pnu,get_u_mid,get_u_pid,
120 . get_u_mat,get_u_geo, get_u_func
121 parameter(kfunc=29)
122 parameter(kmat=31)
123 parameter(kprop=47)
124C-----------------------------------------------
125C L o c a l V a r i a b l e s
126C-----------------------------------------------
127 INTEGER I,ID1,ID2,ID10,ID20,
128 . IFUN_XP,IFUN_XMI,IFUN_XXPI,IFUN_XXMI,IFUN_YY1PI,
129 . IFUN_YY1MI,IFUN_YY2PI,IFUN_YY2MI,IFUN_ZZ1PI,
130 . IFUN_ZZ1MI,IFUN_ZZ2PI,IFUN_ZZ2MI,
131 . IFUN_XMR,IFUN_XXPR,IFUN_XXMR,IFUN_YY1PR,
132 . IFUN_YY1MR,IFUN_YY2PR,IFUN_YY2MR,IFUN_ZZ1PR,
133 . IFUN_ZZ1MR,IFUN_ZZ2PR,IFUN_ZZ2MR,
134 . IFUN_DAMP_XI,IFUN_DAMP_YI,IFUN_DAMP_ZI,
135 . IFUN_DAMP_XXI,IFUN_DAMP_YYI,IFUN_DAMP_ZZI
136 INTEGER IFUNXM(NEL),IFUNXXM(NEL),IFUNYY1M(NEL),
137 . IFUNZZ1M(NEL),IFUNYY2M(NEL),IFUNZZ2M(NEL),
138 . IFUNXP(NEL),IFUNXXP(NEL),IFUNYY1P(NEL),
139 . IFUNZZ1P(NEL),IFUNYY2P(NEL),IFUNZZ2P(NEL),
140 . JPOSXM(NEL),JPOSXXM(NEL),JPOSYY1M(NEL),
141 . JPOSZZ1M(NEL),JPOSYY2M(NEL),JPOSZZ2M(NEL),
142 . JPOSXP(NEL),JPOSXXP(NEL),JPOSYY1P(NEL),
143 . JPOSZZ1P(NEL),JPOSYY2P(NEL),JPOSZZ2P(NEL),
144 . JPOS_DAMP_X(NEL),JPOS_DAMP_Y(NEL),JPOS_DAMP_Z(NEL),
145 . JPOS_DAMP_XX(NEL),JPOS_DAMP_YY(NEL),JPOS_DAMP_ZZ(NEL),
146 . IFUN_DAMP_X(NEL),IFUN_DAMP_Y(NEL),IFUN_DAMP_Z(NEL),
147 . IFUN_DAMP_XX(NEL),IFUN_DAMP_YY(NEL),IFUN_DAMP_ZZ(NEL)
148 my_real ::
149 . xlimg,xlim,xxlim,yy1lim,yy2lim,zz1lim,zz2lim,
150 . k11,k44,k55,k66,k5b,k6c,leng2,leng,uleng0,r,ax,ay,az,
151 . fscal_x,fscal_rx,fscal_ry1,fscal_ry2,fscal_rz1,fscal_rz2,
152 . a,b,d,ff1,ff2,r1,mm,ff,fscal_damp_x,fscal_damp_y,
153 . fscal_damp_z,fscal_damp_xx,fscal_damp_yy,
154 . fscal_damp_zz,f_x,f_y,f_z,f_xx,f_yy,f_zz,alpha,ncf,
155 . idamping
156 my_real , DIMENSION(NEL) ::
157 . dydx,x,fxm,fxp,xx,fxxm,fxxp,yy1,fyy1m,fyy1p,
158 . zz1,fzz1m,fzz1p,yy2,fyy2m,fyy2p,zz2,fzz2m,fzz2p,
159 . uleng,my1,my2,mz1,mz2,fx_damp,fy_damp,fz_damp,
160 . fxx_damp,fyy_damp,fzz_damp,xdot,xxdot,yy1dot,zz1dot,yy2dot,zz2dot,
161 . xdot_old,xxdot_old,yy1dot_old,zz1dot_old,yy2dot_old,zz2dot_old,
162 . elodotx,elodoty,elodotz,elodotxx,elodotyy,elodotzz
163C=======================================================================
164 xlimg = get_u_geo(1,iprop)
165 xlim = get_u_geo(2,iprop)
166 xxlim = get_u_geo(3,iprop)
167 yy1lim = get_u_geo(4,iprop)
168 zz1lim = get_u_geo(5,iprop)
169 yy2lim = get_u_geo(6,iprop)
170 zz2lim = get_u_geo(7,iprop)
171 ico = nint(get_u_geo(16,iprop))
172 fscal_x = get_u_geo(17,iprop)
173 fscal_rx = get_u_geo(18,iprop)
174 fscal_ry1 = get_u_geo(19,iprop)
175 fscal_ry2 = get_u_geo(20,iprop)
176 fscal_rz1 = get_u_geo(21,iprop)
177 fscal_rz2 = get_u_geo(22,iprop)
178!------------------------------------------
179 ifun_xmi = get_u_pnu(1,iprop,kfunc)
180 ifun_xxmi = get_u_pnu(2,iprop,kfunc)
181 ifun_yy1mi= get_u_pnu(3,iprop,kfunc)
182 ifun_zz1mi= get_u_pnu(4,iprop,kfunc)
183 ifun_yy2mi= get_u_pnu(5,iprop,kfunc)
184 ifun_zz2mi= get_u_pnu(6,iprop,kfunc)
185 ifun_xp = get_u_pnu(7,iprop,kfunc)
186 ifun_xxpi = get_u_pnu(8,iprop,kfunc)
187 ifun_yy1pi= get_u_pnu(9,iprop,kfunc)
188 ifun_zz1pi= get_u_pnu(10,iprop,kfunc)
189 ifun_yy2pi= get_u_pnu(11,iprop,kfunc)
190 ifun_zz2pi= get_u_pnu(12,iprop,kfunc)
191 ifun_xmr = get_u_pnu(13,iprop,kfunc)
192 ifun_xxmr = get_u_pnu(14,iprop,kfunc)
193 ifun_yy1mr= get_u_pnu(15,iprop,kfunc)
194 ifun_zz1mr= get_u_pnu(16,iprop,kfunc)
195 ifun_yy2mr= get_u_pnu(17,iprop,kfunc)
196 ifun_zz2mr= get_u_pnu(18,iprop,kfunc)
197 ifun_xxpr = get_u_pnu(19,iprop,kfunc)
198 ifun_yy1pr= get_u_pnu(20,iprop,kfunc)
199 ifun_zz1pr= get_u_pnu(21,iprop,kfunc)
200 ifun_yy2pr= get_u_pnu(22,iprop,kfunc)
201 ifun_zz2pr= get_u_pnu(23,iprop,kfunc)
202!------------------------------------------
203! strain rate filtering number of cycles:
204 ncf = get_u_geo(35,iprop)
205! flag for damping activation
206 idamping = get_u_geo(36,iprop)
207!
208! add damping:
209!
210 IF (idamping > zero) THEN
211 ifun_damp_xi = get_u_pnu(24,iprop,kfunc)
212 ifun_damp_yi = get_u_pnu(25,iprop,kfunc)
213 ifun_damp_zi = get_u_pnu(26,iprop,kfunc)
214 ifun_damp_xxi = get_u_pnu(27,iprop,kfunc)
215 ifun_damp_yyi = get_u_pnu(28,iprop,kfunc)
216 ifun_damp_zzi = get_u_pnu(29,iprop,kfunc)
217!
218 fscal_damp_x = get_u_geo(23,iprop)
219 fscal_damp_y = get_u_geo(24,iprop)
220 fscal_damp_z = get_u_geo(25,iprop)
221 fscal_damp_xx = get_u_geo(26,iprop)
222 fscal_damp_yy = get_u_geo(27,iprop)
223 fscal_damp_zz = get_u_geo(28,iprop)
224! coefficients for function damping:
225 f_x = get_u_geo(29,iprop)
226 f_y = get_u_geo(30,iprop)
227 f_z = get_u_geo(31,iprop)
228 f_xx = get_u_geo(32,iprop)
229 f_yy = get_u_geo(33,iprop)
230 f_zz = get_u_geo(34,iprop)
231 ENDIF ! IF (IDAMPING > ZERO)
232!
233! add filtering:
234!
235 IF (ncf > zero) THEN
236 alpha = two/(ncf+one)
237!
238 DO i=1,nel
239 xdot_old(i) = uvar(37,i) ! old strain rate
240 xxdot_old(i) = uvar(38,i)
241 yy1dot_old(i) = uvar(39,i)
242 zz1dot_old(i) = uvar(40,i)
243 yy2dot_old(i) = uvar(41,i)
244 zz2dot_old(i) = uvar(42,i)
245 ENDDO
246!
247 DO i=1,nel
248 uleng0 = uvar(30,i)
249 xdot(i) = vx(i) * uleng0 ! current strain rate
250 xxdot(i) = rx(i)
251 yy1dot(i) = ry1(i)
252 zz1dot(i) = rz1(i)
253 yy2dot(i) = ry2(i)
254 zz2dot(i) = rz2(i)
255 ENDDO
256! filter current strain rate
257 DO i=1,nel
258 xdot(i) = alpha * xdot(i) + (one-alpha) * xdot_old(i)
259 xxdot(i) = alpha * xxdot(i) + (one-alpha) * xxdot_old(i)
260 yy1dot(i) = alpha * yy1dot(i) + (one-alpha) * yy1dot_old(i)
261 zz1dot(i) = alpha * zz1dot(i) + (one-alpha) * zz1dot_old(i)
262 yy2dot(i) = alpha * yy2dot(i) + (one-alpha) * yy2dot_old(i)
263 zz2dot(i) = alpha * zz2dot(i) + (one-alpha) * zz2dot_old(i)
264 ENDDO
265! save new filtered strain rate:
266 DO i=1,nel
267 uvar(37,i) = xdot(i) ! new filtered strain rate
268 uvar(38,i) = xxdot(i)
269 uvar(39,i) = yy1dot(i)
270 uvar(40,i) = zz1dot(i)
271 uvar(41,i) = yy2dot(i)
272 uvar(42,i) = zz2dot(i)
273 ENDDO
274! total filtered strain :
275 DO i=1,nel
276!! ULENG0 = UVAR(30,I)
277!! X(I) = UVAR(1,I) + DT * XDOT(I) * ULENG0 ! filtered strain
278 x(i) = uvar(1,i) + dt * xdot(i) ! filtered strain
279 xx(i) = uvar(2,i) + dt * xxdot(i)
280 yy1(i) = uvar(3,i) + dt * yy1dot(i)
281 zz1(i) = uvar(4,i) + dt * zz1dot(i)
282 yy2(i) = uvar(5,i) + dt * yy2dot(i)
283 zz2(i) = uvar(6,i) + dt * zz2dot(i)
284 ENDDO
285 ENDIF ! IF (NCF > ZERO)
286C=======================================================================
287C ELASTIQUE
288C=======================================================================
289 DO i=1,nel
290 my1(i) = -ymom(i)+half*xl(i)*fz(i)
291 my2(i) = ymom(i)+half*xl(i)*fz(i)
292 mz1(i) = -zmom(i)-half*xl(i)*fy(i)
293 mz2(i) = zmom(i)-half*xl(i)*fy(i)
294C
295 leng = max(xl(i),em20)
296 leng2 = leng*leng
297 uleng(i) = one/leng
298C
299 k11 = uvar(19,i)
300 k44 = uvar(20,i)
301 k55 = uvar(21,i)*leng2
302 k66 = uvar(22,i)*leng2
303 k5b = uvar(23,i)*leng2
304 k6c = uvar(24,i)*leng2
305C
306 fx(i) = fx(i) + dt * k11 * vx(i)
307 xmom(i)= xmom(i)+ dt * k44 * rx(i)
308 my1(i) = my1(i) + dt * (k55 * ry1(i) + k5b * ry2(i))
309 my2(i) = my2(i) + dt * (k5b * ry1(i) + k55 * ry2(i))
310 mz1(i) = mz1(i) + dt * (k66 * rz1(i) + k6c * rz2(i))
311 mz2(i) = mz2(i) + dt * (k6c * rz1(i) + k66 * rz2(i))
312C
313 uleng0 = uvar(30,i)
314C
315 IF (ncf == zero) THEN ! no filtering
316 x(i) = uvar(1,i) + dt * vx(i) * uleng0 ! strain
317 xx(i) = uvar(2,i) + dt * rx(i)
318 yy1(i) = uvar(3,i) + dt * ry1(i)
319 zz1(i) = uvar(4,i) + dt * rz1(i)
320 yy2(i) = uvar(5,i) + dt * ry2(i)
321 zz2(i) = uvar(6,i) + dt * rz2(i)
322 ENDIF ! IF (NCF == 0)
323C
324 jposxm(i) = nint(uvar(7,i))
325 jposxxm(i) = nint(uvar(8,i))
326 jposyy1m(i) = nint(uvar(9,i))
327 jposzz1m(i) = nint(uvar(10,i))
328 jposyy2m(i) = nint(uvar(11,i))
329 jposzz2m(i) = nint(uvar(12,i))
330 jposxp(i) = nint(uvar(13,i))
331 jposxxp(i) = nint(uvar(14,i))
332 jposyy1p(i) = nint(uvar(15,i))
333 jposzz1p(i) = nint(uvar(16,i))
334 jposyy2p(i) = nint(uvar(17,i))
335 jposzz2p(i) = nint(uvar(18,i))
336C
337 ifunxp(i) = ifun_xp
338 ifunxm(i) = ifun_xmi
339 ifunxxm(i) = ifun_xxmi
340 ifunxxp(i) = ifun_xxpi
341 ifunyy1m(i) = ifun_yy1mi
342 ifunyy1p(i) = ifun_yy1pi
343 ifunyy2m(i) = ifun_yy2mi
344 ifunyy2p(i) = ifun_yy2pi
345 ifunzz1m(i) = ifun_zz1mi
346 ifunzz1p(i) = ifun_zz1pi
347 ifunzz2m(i) = ifun_zz2mi
348 ifunzz2p(i) = ifun_zz2pi
349!
350! add function damping:
351!
352 IF (idamping > zero) THEN
353 ifun_damp_x(i) = ifun_damp_xi
354 ifun_damp_y(i) = ifun_damp_yi
355 ifun_damp_z(i) = ifun_damp_zi
356 ifun_damp_xx(i) = ifun_damp_xxi
357 ifun_damp_yy(i) = ifun_damp_yyi
358 ifun_damp_zz(i) = ifun_damp_zzi
359!
360 jpos_damp_x(i) = nint(uvar(31,i))
361 jpos_damp_y(i) = nint(uvar(32,i))
362 jpos_damp_z(i) = nint(uvar(33,i))
363 jpos_damp_xx(i) = nint(uvar(34,i))
364 jpos_damp_yy(i) = nint(uvar(35,i))
365 jpos_damp_zz(i) = nint(uvar(36,i))
366 ENDIF ! IF (IDAMPING > ZERO)
367 ENDDO
368!
369 IF (idamping > zero) THEN
370! elongation rate for damping (linear + function):
371 DO i=1,nel
372 elodotx(i) = vx(i) ! FX
373 elodoty(i) = - half * (rz2(i) + rz1(i)) * xl(i) ! FY
374 elodotz(i) = half * (ry2(i) + ry1(i)) * xl(i) ! FZ
375 elodotxx(i) = rx(i) ! XMOM
376 elodotyy(i) = ry2(i) - ry1(i) ! YMOM
377 elodotzz(i) = rz2(i) - rz1(i) ! ZMOM
378 ENDDO
379 ENDIF ! IF (IDAMPING > ZERO)
380C--------------------------------
381 IF (ico == 0) THEN ! => ELASTO PLASTIQUE - Classique
382C--------------------------------
383 DO i=1,nel
384 id1 = nint(uvar(29,i))/2
385 id2 = nint(uvar(29,i)) - 2*id1
386 id10=id1
387 id20=id2
388 IF (fr_wave_e(i) == one) THEN
389 jposxm(i) = 0
390 ifunxm(i) = ifun_xmr
391 ENDIF
392 IF (x(i) < xlimg) THEN !!! collapsed element
393 fr_wave_e(i)=one
394 ELSE
395 fr_wave_e(i)=zero
396 ENDIF
397C
398 IF (x(i) < xlim .OR. abs(xx(i)) > xxlim) THEN
399 id1 = 1
400 id2 = 1
401 ENDIF
402 IF (abs(yy1(i)) > yy1lim .OR. abs(zz1(i)) > zz1lim) id1 = 1
403 IF (abs(yy2(i)) > yy2lim .OR. abs(zz2(i)) > zz2lim) id2 = 1
404C
405 IF (id1 == 1) THEN
406 ifunxm(i) = ifun_xmr
407 ifunxxm(i) = ifun_xxmr
408 ifunxxp(i) = ifun_xxpr
409 ifunyy1m(i) = ifun_yy1mr
410 ifunyy1p(i) = ifun_yy1pr
411 ifunzz1m(i) = ifun_zz1mr
412 ifunzz1p(i) = ifun_zz1pr
413 IF (id10 == 0) THEN
414 jposxm(i) = 0
415 jposxxm(i) = 0
416 jposxxp(i) = 0
417 jposyy1m(i) = 0
418 jposzz1m(i) = 0
419 jposyy1p(i) = 0
420 jposzz1p(i) = 0
421 ENDIF
422 ENDIF
423 IF (id2 == 1) THEN
424 ifunxm(i) = ifun_xmr
425 ifunxxm(i) = ifun_xxmr
426 ifunxxp(i) = ifun_xxpr
427 ifunyy2m(i) = ifun_yy2mr
428 ifunyy2p(i) = ifun_yy2pr
429 ifunzz2m(i) = ifun_zz2mr
430 ifunzz2p(i) = ifun_zz2pr
431 IF (id20 == 0) THEN
432 jposxm(i) = 0
433 jposxxm(i) = 0
434 jposxxp(i) = 0
435 jposyy2m(i) = 0
436 jposzz2m(i) = 0
437 jposyy2p(i) = 0
438 jposzz2p(i) = 0
439 ENDIF
440 ENDIF
441 uvar(29,i) = 2*id1 + id2
442 ENDDO
443!
444 CALL get_v_func(ifunxm,nel,x,dydx,fxm,jposxm)
445 CALL get_v_func(ifunxp,nel,x,dydx,fxp,jposxp)
446 CALL get_v_func(ifunxxm,nel,xx,dydx,fxxm,jposxxm)
447 CALL get_v_func(ifunxxp,nel,xx,dydx,fxxp,jposxxp)
448 CALL get_v_func(ifunyy1m,nel,yy1,dydx,fyy1m,jposyy1m)
449 CALL get_v_func(ifunyy1p,nel,yy1,dydx,fyy1p,jposyy1p)
450 CALL get_v_func(ifunzz1m,nel,zz1,dydx,fzz1m,jposzz1m)
451 CALL get_v_func(ifunzz1p,nel,zz1,dydx,fzz1p,jposzz1p)
452 CALL get_v_func(ifunyy2m,nel,yy2,dydx,fyy2m,jposyy2m)
453 CALL get_v_func(ifunyy2p,nel,yy2,dydx,fyy2p,jposyy2p)
454 CALL get_v_func(ifunzz2m,nel,zz2,dydx,fzz2m,jposzz2m)
455 CALL get_v_func(ifunzz2p,nel,zz2,dydx,fzz2p,jposzz2p)
456!-------------------------------------------------------------
457!
458! add damping function interpolation:
459!
460 IF (idamping > zero) THEN
461 IF (ifun_damp_xi > 0) THEN
462 DO i=1,nel
463 elodotx(i) = elodotx(i) / f_x
464 ENDDO
465 CALL get_v_func(ifun_damp_x,nel,elodotx,dydx,fx_damp,jpos_damp_x)
466 ENDIF
467!
468 IF (ifun_damp_yi > 0) THEN
469 DO i=1,nel
470 elodoty(i) = elodoty(i) / f_y
471 ENDDO
472 CALL get_v_func(ifun_damp_y,nel,elodoty,dydx,fy_damp,jpos_damp_y)
473 ENDIF
474!
475 IF (ifun_damp_zi > 0) THEN
476 DO i=1,nel
477 elodotz(i) = elodotz(i) / f_z
478 ENDDO
479 CALL get_v_func(ifun_damp_z,nel,elodotz,dydx,fz_damp,jpos_damp_z)
480 ENDIF
481!
482 IF (ifun_damp_xxi > 0) THEN
483 DO i=1,nel
484 elodotxx(i) = elodotxx(i) / f_xx
485 ENDDO
486 CALL get_v_func(ifun_damp_xx,nel,elodotxx,dydx,fxx_damp,jpos_damp_xx)
487 ENDIF
488!
489 IF (ifun_damp_yyi > 0) THEN
490 DO i=1,nel
491 elodotyy(i) = elodotyy(i) / f_yy
492 ENDDO
493 CALL get_v_func(ifun_damp_yy,nel,elodotyy,dydx,fyy_damp,jpos_damp_yy)
494 ENDIF
495!
496 IF (ifun_damp_zzi > 0) THEN
497 DO i=1,nel
498 elodotzz(i) = elodotzz(i) / f_zz
499 ENDDO
500 CALL get_v_func(ifun_damp_zz,nel,elodotzz,dydx,fzz_damp,jpos_damp_zz)
501 ENDIF
502 ENDIF ! IF (IDAMPING > ZERO)
503!
504 DO i=1,nel
505 fxp(i) = fxp(i) * fscal_x
506 fxm(i) = fxm(i) * fscal_x
507 fxxp(i) = fxxp(i) * fscal_rx
508 fxxm(i) = fxxm(i) * fscal_rx
509 fyy1m(i) = fyy1m(i)* fscal_ry1
510 fyy1p(i) = fyy1p(i)* fscal_ry1
511 fzz1m(i) = fzz1m(i)* fscal_rz1
512 fzz1p(i) = fzz1p(i)* fscal_rz1
513 fyy2m(i) = fyy2m(i)* fscal_ry2
514 fyy2p(i) = fyy2p(i)* fscal_ry2
515 fzz2m(i) = fzz2m(i)* fscal_rz2
516 fzz2p(i) = fzz2p(i)* fscal_rz2
517 ENDDO
518!
519 DO i=1,nel
520 uvar(1,i) = x(i)
521 uvar(2,i) = xx(i)
522 uvar(3,i) = yy1(i)
523 uvar(4,i) = zz1(i)
524 uvar(5,i) = yy2(i)
525 uvar(6,i) = zz2(i)
526 uvar(7,i) = jposxm(i)
527 uvar(8,i) = jposxxm(i)
528 uvar(9,i) = jposyy1m(i)
529 uvar(10,i) = jposzz1m(i)
530 uvar(11,i) = jposyy2m(i)
531 uvar(12,i) = jposzz2m(i)
532 uvar(13,i) = jposxp(i)
533 uvar(14,i) = jposxxp(i)
534 uvar(15,i) = jposyy1p(i)
535 uvar(16,i) = jposzz1p(i)
536 uvar(17,i) = jposyy2p(i)
537 uvar(18,i) = jposzz2p(i)
538!---
539! add damping:
540!---
541 IF (idamping > zero) THEN
542 uvar(31,i) = jpos_damp_x(i)
543 uvar(32,i) = jpos_damp_y(i)
544 uvar(33,i) = jpos_damp_z(i)
545 uvar(34,i) = jpos_damp_xx(i)
546 uvar(35,i) = jpos_damp_yy(i)
547 uvar(36,i) = jpos_damp_zz(i)
548 ENDIF ! IF IDAMPING > ZERO)
549!---
550 fx(i) = max(min(fx(i) ,fxp(i)),fxm(i) )
551 xmom(i) = max(min(xmom(i),fxxp(i)),fxxm(i))
552 my1(i) = max(min(my1(i),fyy1p(i)),fyy1m(i))
553 mz1(i) = max(min(mz1(i),fzz1p(i)),fzz1m(i))
554 my2(i) = max(min(my2(i),fyy2p(i)),fyy2m(i))
555 mz2(i) = max(min(mz2(i),fzz2p(i)),fzz2m(i))
556 fz(i) = (my1(i)+my2(i)) * uleng(i)
557 fy(i) = -(mz1(i)+mz2(i)) * uleng(i)
558 ymom(i) = -half*(my1(i)-my2(i))
559 zmom(i) = -half*(mz1(i)-mz2(i))
560C
561 viscm(i) = zero
562 viscr(i) = zero
563 stifm(i) = uvar(25,i)
564 stifr(i) = uvar(26,i)
565 mass(i) = uvar(27,i)
566 xiner(i) = uvar(28,i)
567!
568! add damping term (linear or function):
569!
570 IF (idamping > zero) THEN
571 IF (ifun_damp_x(i) == 0) THEN ! ! linear damping
572 fx(i) = fx(i) + fscal_damp_x * elodotx(i)
573 ELSE
574 fx(i) = fx(i) + fscal_damp_x * fx_damp(i)
575 ENDIF
576!
577 IF (ifun_damp_y(i) == 0) THEN ! ! linear damping
578 fy(i) = fy(i) + fscal_damp_y * elodoty(i)
579 ELSE
580 fy(i) = fy(i) + fscal_damp_y * fy_damp(i)
581 ENDIF
582!
583 IF (ifun_damp_z(i) == 0) THEN ! ! linear damping
584 fz(i) = fz(i) + fscal_damp_z * elodotz(i)
585 ELSE
586 fz(i) = fz(i) + fscal_damp_z * fz_damp(i)
587 ENDIF
588!
589 IF (ifun_damp_xx(i) == 0) THEN ! ! linear damping
590 xmom(i) = xmom(i) + fscal_damp_xx * elodotxx(i)
591 ELSE
592 xmom(i) = xmom(i) + fscal_damp_xx * fxx_damp(i)
593 ENDIF
594!
595 IF (ifun_damp_yy(i) == 0) THEN ! ! linear damping
596 ymom(i) = ymom(i) + fscal_damp_yy * elodotyy(i)
597 ELSE
598 ymom(i) = ymom(i) + fscal_damp_yy * fyy_damp(i)
599 ENDIF
600!
601 IF (ifun_damp_zz(i) == 0) THEN ! ! linear damping
602 zmom(i) = zmom(i) + fscal_damp_zz * elodotzz(i)
603 ELSE
604 zmom(i) = zmom(i) + fscal_damp_zz * fzz_damp(i)
605 ENDIF
606 ENDIF ! IF (IDAMPING > ZERO)
607!
608 ENDDO
609C--------------------------------
610 ELSE ! => ELASTO PLASTIQUE - Directions couplees
611C--------------------------------
612 DO i=1,nel
613 id1 = nint(uvar(29,i))/2
614 id2 = nint(uvar(29,i)) - 2*id1
615 id10=id1
616 id20=id2
617 IF (fr_wave_e(i) == one) THEN !!! collapsed element
618 jposxm(i) = 0
619 ifunxm(i) = ifun_xmr
620 ENDIF
621 IF (x(i) < xlimg) THEN
622 fr_wave_e(i)=one
623 ELSE
624 fr_wave_e(i)=zero
625 ENDIF
626C
627 IF (x(i) < xlim .OR. abs(xx(i)) > xxlim) THEN
628 id1 = 1
629 id2 = 1
630 ENDIF
631 IF (abs(yy1(i)) > yy1lim .OR. abs(zz1(i)) > zz1lim) id1 = 1
632 IF (abs(yy2(i)) > yy2lim .OR. abs(zz2(i)) > zz2lim) id2 = 1
633C
634 IF (id1 == 1) THEN
635 ifunxm(i) = ifun_xmr
636 ifunxxm(i) = ifun_xxmr
637 ifunxxp(i) = ifun_xxpr
638 ifunyy1m(i) = ifun_yy1mr
639 ifunyy1p(i) = ifun_yy1pr
640 ifunzz1m(i) = ifun_zz1mr
641 ifunzz1p(i) = ifun_zz1pr
642 IF (id10 == 0) THEN
643 jposxm(i) = 0
644 jposxxm(i) = 0
645 jposxxp(i) = 0
646 jposyy1m(i) = 0
647 jposzz1m(i) = 0
648 jposyy1p(i) = 0
649 jposzz1p(i) = 0
650 ENDIF
651 ENDIF
652 IF (id2 == 1) THEN
653 ifunxm(i) = ifun_xmr
654 ifunxxm(i) = ifun_xxmr
655 ifunxxp(i) = ifun_xxpr
656 ifunyy2m(i) = ifun_yy2mr
657 ifunyy2p(i) = ifun_yy2pr
658 ifunzz2m(i) = ifun_zz2mr
659 ifunzz2p(i) = ifun_zz2pr
660 IF (id20 == 0) THEN
661 jposxm(i) = 0
662 jposxxm(i) = 0
663 jposxxp(i) = 0
664 jposyy2m(i) = 0
665 jposzz2m(i) = 0
666 jposyy2p(i) = 0
667 jposzz2p(i) = 0
668 ENDIF
669 ENDIF
670 uvar(29,i) = 2*id1 + id2
671 ENDDO
672C
673 CALL get_v_func(ifunxm,nel,x,dydx,fxm,jposxm)
674 CALL get_v_func(ifunxp,nel,x,dydx,fxp,jposxp)
675 CALL get_v_func(ifunxxm,nel,xx,dydx,fxxm,jposxxm)
676 CALL get_v_func(ifunxxp,nel,xx,dydx,fxxp,jposxxp)
677 CALL get_v_func(ifunyy1m,nel,yy1,dydx,fyy1m,jposyy1m)
678 CALL get_v_func(ifunyy1p,nel,yy1,dydx,fyy1p,jposyy1p)
679 CALL get_v_func(ifunzz1m,nel,zz1,dydx,fzz1m,jposzz1m)
680 CALL get_v_func(ifunzz1p,nel,zz1,dydx,fzz1p,jposzz1p)
681 CALL get_v_func(ifunyy2m,nel,yy2,dydx,fyy2m,jposyy2m)
682 CALL get_v_func(ifunyy2p,nel,yy2,dydx,fyy2p,jposyy2p)
683 CALL get_v_func(ifunzz2m,nel,zz2,dydx,fzz2m,jposzz2m)
684 CALL get_v_func(ifunzz2p,nel,zz2,dydx,fzz2p,jposzz2p)
685!
686! add damping function interpolation:
687!
688 IF (idamping > zero) THEN
689 IF (ifun_damp_xi > 0) THEN
690 DO i=1,nel
691 elodotx(i) = elodotx(i) / f_x
692 ENDDO
693 CALL get_v_func(ifun_damp_x,nel,elodotx,dydx,fx_damp,jpos_damp_x)
694 ENDIF
695!
696 IF (ifun_damp_yi > 0) THEN
697 DO i=1,nel
698 elodoty(i) = elodoty(i) / f_y
699 ENDDO
700 CALL get_v_func(ifun_damp_y,nel,elodoty,dydx,fy_damp,jpos_damp_y)
701 ENDIF
702!
703 IF (ifun_damp_zi > 0) THEN
704 DO i=1,nel
705 elodotz(i) = elodotz(i) / f_z
706 ENDDO
707 CALL get_v_func(ifun_damp_z,nel,elodotz,dydx,fz_damp,jpos_damp_z)
708 ENDIF
709!
710 IF (ifun_damp_xxi > 0) THEN
711 DO i=1,nel
712 elodotxx(i) = elodotxx(i) / f_xx
713 ENDDO
714 CALL get_v_func(ifun_damp_xx,nel,elodotxx,dydx,fxx_damp,jpos_damp_xx)
715 ENDIF
716!
717 IF (ifun_damp_yyi > 0) THEN
718 DO i=1,nel
719 elodotyy(i) = elodotyy(i) / f_yy
720 ENDDO
721 CALL get_v_func(ifun_damp_yy,nel,elodotyy,dydx,fyy_damp,jpos_damp_yy)
722 ENDIF
723!
724 IF (ifun_damp_zzi > 0) THEN
725 DO i=1,nel
726 elodotzz(i) = elodotzz(i) / f_zz
727 ENDDO
728 CALL get_v_func(ifun_damp_zz,nel,elodotzz,dydx,fzz_damp,jpos_damp_zz)
729 ENDIF
730 ENDIF ! IF (IDAMPING > ZERO)
731!
732 DO i=1,nel
733 fxp(i) = fxp(i) * fscal_x
734 fxm(i) = fxm(i) * fscal_x
735 fxxp(i) = fxxp(i) * fscal_rx
736 fxxm(i) = fxxm(i) * fscal_rx
737 fyy1m(i) = fyy1m(i)* fscal_ry1
738 fyy1p(i) = fyy1p(i)* fscal_ry1
739 fzz1m(i) = fzz1m(i)* fscal_rz1
740 fzz1p(i) = fzz1p(i)* fscal_rz1
741 fyy2m(i) = fyy2m(i)* fscal_ry2
742 fyy2p(i) = fyy2p(i)* fscal_ry2
743 fzz2m(i) = fzz2m(i)* fscal_rz2
744 fzz2p(i) = fzz2p(i)* fscal_rz2
745 ENDDO
746C
747 DO i=1,nel
748 uvar(1,i) = x(i)
749 uvar(2,i) = xx(i)
750 uvar(3,i) = yy1(i)
751 uvar(4,i) = zz1(i)
752 uvar(5,i) = yy2(i)
753 uvar(6,i) = zz2(i)
754 uvar(7,i) = jposxm(i)
755 uvar(8,i) = jposxxm(i)
756 uvar(9,i) = jposyy1m(i)
757 uvar(10,i) = jposzz1m(i)
758 uvar(11,i) = jposyy2m(i)
759 uvar(12,i) = jposzz2m(i)
760 uvar(13,i) = jposxp(i)
761 uvar(14,i) = jposxxp(i)
762 uvar(15,i) = jposyy1p(i)
763 uvar(16,i) = jposzz1p(i)
764 uvar(17,i) = jposyy2p(i)
765 uvar(18,i) = jposzz2p(i)
766!---
767! add damping:
768!---
769 IF (idamping > zero) THEN
770 uvar(31,i) = jpos_damp_x(i)
771 uvar(32,i) = jpos_damp_y(i)
772 uvar(33,i) = jpos_damp_z(i)
773 uvar(34,i) = jpos_damp_xx(i)
774 uvar(35,i) = jpos_damp_yy(i)
775 uvar(36,i) = jpos_damp_zz(i)
776 ENDIF ! IF IDAMPING > ZERO)
777!
778 my1(i) = max(min(my1(i),fyy1p(i)),fyy1m(i))
779 mz1(i) = max(min(mz1(i),fzz1p(i)),fzz1m(i))
780 my2(i) = max(min(my2(i),fyy2p(i)),fyy2m(i))
781 mz2(i) = max(min(mz2(i),fzz2p(i)),fzz2m(i))
782!
783 xmom(i) = max(min(xmom(i),fxxp(i)),fxxm(i))
784!
785 fyy1p(i) = max(em20, fyy1p(i))
786 fyy1m(i) = -max(em20,-fyy1m(i))
787 fyy2p(i) = max(em20, fyy2p(i))
788 fyy2m(i) = -max(em20,-fyy2m(i))
789!
790 fzz1p(i) = max(em20, fzz1p(i))
791 fzz1m(i) = -max(em20,-fzz1m(i))
792 fzz2p(i) = max(em20, fzz2p(i))
793 fzz2m(i) = -max(em20,-fzz2m(i))
794!
795 fxxp(i) = max(em20, fxxp(i))
796 fxxm(i) = -max(em20,-fxxm(i))
797!
798C---------------------------------------------------
799C Mo = R M
800C Mn = Mo [l^2 Fo^2 - l^2Fo F + Mo M ]/[l^2Fo^2+Mo^2]
801C Fn = Fo [ Mo^2 - Mo M + l^2Fo F ]/[l^2Fo^2+Mo^2]
802C Fn = Fo [ 1 - Mn/Mo ]
803C Fn = Fo [ 1 - [ l^2Fo^2 - l^2Fo F + Mo M ]/[l^2Fo^2+Mo^2]]
804C Fn = Fo [ 1 - [ l^2Fo^2 - l^2Fo F + R M^2 ]/[l^2Fo^2+R^2M^2]]
805C Fn = Fo [ 1 - [ l^2Fo(Fo - F) + R M^2 ]/[l^2Fo^2+R^2M^2]]
806C Mn = Mo [ 1 - Fn/Fo ]
807C Rn = Mn/M = Mo/M [ 1 - Fn/Fo ]
808C Rn = R [ 1 - Fn/Fo ]
809C si l = Mo/Fo
810C Fn = Fo [ 1/2 + F/2Fo - 1/2R ]
811C Rn = R [ 1 - Fn/Fo ] ! twice written
812C---------------------------------------------------
813!
814! - NODE 1 -
815!
816 ax = max(xmom(i)/fxxp(i),xmom(i)/fxxm(i))
817 ay = max(my1(i)/fyy1p(i),my1(i)/fyy1m(i))
818 az = max(mz1(i)/fzz1p(i),mz1(i)/fzz1m(i))
819!
820 mm = sqrt(ax*ax+ay*ay+az*az)
821!
822 ff = fx(i)/min(-em20,fxm(i))
823 r = (one+mm-ff)/max(em20,two*mm)
824 IF (r > zero .AND. r < one) THEN
825 ff1 = half*(one-mm+ff)*fxm(i)
826 ELSE
827 ff1 = fxm(i)
828 r = max(zero,min(one,r))
829 ENDIF
830 r1 = r
831!
832! - NODE 2 -
833!
834 ay = max(my2(i)/fyy2p(i),my2(i)/fyy2m(i))
835 az = max(mz2(i)/fzz2p(i),mz2(i)/fzz2m(i))
836 mm = sqrt(ax*ax+ay*ay+az*az)
837 r = (one+mm-ff)/max(em20,two*mm)
838 IF (r > zero .AND. r < one) THEN
839 ff2 = half*(one-mm+ff)*fxm(i)
840 ELSE
841 ff2 = fxm(i)
842 r = max(zero,min(one,r))
843 ENDIF
844 my1(i) = r1*my1(i)
845 mz1(i) = r1*mz1(i)
846 my2(i) = r*my2(i)
847 mz2(i) = r*mz2(i)
848C
849 xmom(i) = min(r,r1)*xmom(i)
850C
851 fx(i) = max(min(fx(i),fxp(i)), ff1, ff2)
852C
853 fz(i) = (my1(i)+my2(i)) * uleng(i)
854 fy(i) = -(mz1(i)+mz2(i)) * uleng(i)
855 ymom(i) = -half*(my1(i)-my2(i))
856 zmom(i) = -half*(mz1(i)-mz2(i))
857C
858 viscm(i) = zero
859 viscr(i) = zero
860 stifm(i) = uvar(25,i)
861 stifr(i) = uvar(26,i)
862 mass(i) = uvar(27,i)
863 xiner(i) = uvar(28,i)
864!
865! add damping term (linear or function):
866!
867 IF (idamping > zero) THEN
868 IF (ifun_damp_x(i) == 0) THEN ! ! linear damping
869 fx(i) = fx(i) + fscal_damp_x * elodotx(i)
870 ELSE
871 fx(i) = fx(i) + fscal_damp_x * fx_damp(i)
872 ENDIF
873!
874 IF (ifun_damp_y(i) == 0) THEN ! ! linear damping
875 fy(i) = fy(i) + fscal_damp_y * elodoty(i)
876 ELSE
877 fy(i) = fy(i) + fscal_damp_y * fy_damp(i)
878 ENDIF
879!
880 IF (ifun_damp_z(i) == 0) THEN ! ! linear damping
881 fz(i) = fz(i) + fscal_damp_z * elodotz(i)
882 ELSE
883 fz(i) = fz(i) + fscal_damp_z * fz_damp(i)
884 ENDIF
885!
886 IF (ifun_damp_xx(i) == 0) THEN ! ! linear damping
887 xmom(i) = xmom(i) + fscal_damp_xx * elodotxx(i)
888 ELSE
889 xmom(i) = xmom(i) + fscal_damp_xx * fxx_damp(i)
890 ENDIF
891!
892 IF (ifun_damp_yy(i) == 0) THEN ! ! linear damping
893 ymom(i) = ymom(i) + fscal_damp_yy * elodotyy(i)
894 ELSE
895 ymom(i) = ymom(i) + fscal_damp_yy * fyy_damp(i)
896 ENDIF
897!
898 IF (ifun_damp_zz(i) == 0) THEN ! ! linear damping
899 zmom(i) = zmom(i) + fscal_damp_zz * elodotzz(i)
900 ELSE
901 zmom(i) = zmom(i) + fscal_damp_zz * fzz_damp(i)
902 ENDIF
903 ENDIF ! IF (IDAMPING > ZERO)
904!
905 ENDDO
906C-------------------------------
907 ENDIF
908C-------------------------------
909 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
get_u_pid
integer function get_u_pid(ip)
Definition uaccess.F:626
get_u_pnu
integer function get_u_pnu(ivar, ip, k)
Definition uaccess.F:482
get_u_mid
integer function get_u_mid(im)
Definition uaccess.F:668
get_u_mnu
integer function get_u_mnu(ivar, im, k)
Definition uaccess.F:565
get_v_func
subroutine get_v_func(ifunc, llt, xx, dydx, yy, jpos)
Definition ufunc.F:230