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

Go to the source code of this file.

Functions/Subroutines

subroutine rini33 (nel, iout, iprop, ix, xl, mass, xiner, stifn, stifr, viscm, viscr, uvar, nuvar)
integer function get_skew (iout, jtyp, skflg, idsk1, idsk2, u, v, x, a, b)
subroutine qrot33 (skew, t, c, s)
subroutine rot12 (skew, t, c, s)
subroutine prod_abt (a, b, x)
subroutine prod_atb (a, b, x)
subroutine prod_ab (a, b, x)

Function/Subroutine Documentation

◆ get_skew()

integer function get_skew ( integer iout,
integer jtyp,
integer skflg,
integer idsk1,
integer idsk2,
u,
v,
x,
a,
b )

Definition at line 286 of file rini33.F.

287 USE message_mod
288C-----------------------------------------------
289C I m p l i c i t T y p e s
290C-----------------------------------------------
291#include "implicit_f.inc"
292C-----------------------------------------------
293C A n a l y s e M o d u l e
294C-----------------------------------------------
295#include "param_c.inc"
296C----------------------------------------------------------
297C 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
298C----------------------------------------------------------
299 INTEGER IOUT,JTYP,IDSK1,IDSK2,SKFLG
300 my_real u(lskew),v(lskew),x(lskew),a(lskew),b(lskew)
301C-----------------------------------------------
302C L o c a l V a r i a b l e s
303C-----------------------------------------------
304 INTEGER I,J,IERR1,ISK1,ISK2,IP1,IP2, GET_U_SKEW,N1,N2,N3
305 my_real
306 . nx,ny,nz,co,si,ksi,
307 . t(3),q(lskew),q1(lskew),q2(lskew),x1(lskew),x2(lskew)
308C-----------------------------
309 EXTERNAL get_u_skew
310C=======================================================================
311C
312 ierr1 = 0
313 isk1 = get_u_skew(idsk1,n1,n2,n3,u)
314 isk2 = get_u_skew(idsk2,n1,n2,n3,v)
315C
316 IF (jtyp==5) THEN
317C--- universal joint
318 IF ((u(1)*v(1)+u(2)*v(2)+u(3)*v(3))<=em10) THEN
319 x(1) = u(2)*v(3) - u(3)*v(2)
320 x(2) = u(3)*v(1) - u(1)*v(3)
321 x(3) = u(1)*v(2) - u(2)*v(1)
322 nx = sqrt(x(1)*x(1)+x(2)*x(2)+x(3)*x(3))
323 x(1) = x(1) / nx
324 x(2) = x(2) / nx
325 x(3) = x(3) / nx
326 x(4) = u(1)
327 x(5) = u(2)
328 x(6) = u(3)
329 x(7) = v(1)
330 x(8) = v(2)
331 x(9) = v(3)
332C
333 CALL prod_atb(v,u,a)
334 CALL prod_atb(u,v,b)
335 ELSE
336C IERR1 = 1
337C WRITE(*,*)' ** ERROR/NON ORTHOGONAL UNIVERSAL JOINT AXES'
338C WRITE(IOUT,'(//A,A//)')' ** ERROR PROPERTY SET INPUT',
339C . ' NON ORTHOGONAL UNIVERSAL JOINT AXES'
340C IERR1 = IERR1 + 1
341 CALL ancmsg(msgid=389,
342 . msgtype=msgerror,
343 . anmode=aninfo_blind_1)
344c . I1=ID,
345c . C1=TITR)
346 ENDIF
347 ELSE
348C---------------------------
349 IF (skflg==1) THEN
350c---- first skew is used
351 IF (isk2==0) THEN
352 DO i=1,9
353 x(i) = u(i)
354 a(i) = u(i)
355 b(i) = zero
356 ENDDO
357 ELSE
358 DO i=1,9
359 x(i) = u(i)
360 b(i) = zero
361 END DO
362 CALL prod_atb(v,u,a)
363 ENDIF
364C------
365 ELSE
366c---- mean skew is calculated
367 DO i=1,9
368 b(i) = 0.
369 END DO
370 CALL prod_atb(u,v,q)
371 CALL rot12(q, t, co, si)
372 CALL qrot33(q1, t, co, si)
373C
374 CALL prod_atb(v,u,q)
375 CALL rot12(q, t, co, si)
376 CALL qrot33(q2, t, co, si)
377 a(1) = half * (q1(1) + q2(1))
378 a(2) = half * (q1(2) + q2(4))
379 a(3) = half * (q1(3) + q2(7))
380 a(4) = half * (q1(4) + q2(2))
381 a(5) = half * (q1(5) + q2(5))
382 a(6) = half * (q1(6) + q2(8))
383 a(7) = half * (q1(7) + q2(3))
384 a(8) = half * (q1(8) + q2(6))
385 a(9) = half * (q1(9) + q2(9))
386C
387 CALL prod_ab(u,a,x1)
388 CALL prod_abt(v,a,x2)
389 DO i=1,9
390 x(i) = half * (x1(i) + x2(i))
391 END DO
392C-------------------------------------
393 ENDIF
394 ENDIF
395C-----------------------------------------------
396 get_skew = ierr1
397 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
message_mod
Definition message_mod.F:1249
qrot33
subroutine qrot33(skew, t, c, s)
Definition rini33.F:406
prod_atb
subroutine prod_atb(a, b, x)
Definition rini33.F:554
prod_ab
subroutine prod_ab(a, b, x)
Definition rini33.F:576
prod_abt
subroutine prod_abt(a, b, x)
Definition rini33.F:529
get_skew
integer function get_skew(iout, jtyp, skflg, idsk1, idsk2, u, v, x, a, b)
Definition rini33.F:287
rot12
subroutine rot12(skew, t, c, s)
Definition rini33.F:463
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:889
get_u_skew
integer function get_u_skew(idskw, n1, n2, n3, v)
Definition uaccess.F:1128

◆ prod_ab()

subroutine prod_ab ( a,
b,
x )

Definition at line 575 of file rini33.F.

576#include "implicit_f.inc"
577#include "param_c.inc"
578 my_real a(lskew),b(lskew),x(lskew)
579C
580 x(1)=a(1)*b(1)+a(4)*b(2)+a(7)*b(3)
581 x(2)=a(2)*b(1)+a(5)*b(2)+a(8)*b(3)
582 x(3)=a(3)*b(1)+a(6)*b(2)+a(9)*b(3)
583 x(4)=a(1)*b(4)+a(4)*b(5)+a(7)*b(6)
584 x(5)=a(2)*b(4)+a(5)*b(5)+a(8)*b(6)
585 x(6)=a(3)*b(4)+a(6)*b(5)+a(9)*b(6)
586 x(7)=a(1)*b(7)+a(4)*b(8)+a(7)*b(9)
587 x(8)=a(2)*b(7)+a(5)*b(8)+a(8)*b(9)
588 x(9)=a(3)*b(7)+a(6)*b(8)+a(9)*b(9)
589 RETURN

◆ prod_abt()

subroutine prod_abt ( a,
b,
x )

Definition at line 528 of file rini33.F.

529#include "implicit_f.inc"
530#include "param_c.inc"
531 INTEGER I,J
532 my_real a(lskew),b(lskew),x(lskew)
533C
534 x(1)=a(1)*b(1)+a(4)*b(4)+a(7)*b(7)
535 x(4)=a(1)*b(2)+a(4)*b(5)+a(7)*b(8)
536 x(7)=a(1)*b(3)+a(4)*b(6)+a(7)*b(9)
537 x(2)=a(2)*b(1)+a(5)*b(4)+a(8)*b(7)
538 x(5)=a(2)*b(2)+a(5)*b(5)+a(8)*b(8)
539 x(8)=a(2)*b(3)+a(5)*b(6)+a(8)*b(9)
540 x(3)=a(3)*b(1)+a(6)*b(4)+a(9)*b(7)
541 x(6)=a(3)*b(2)+a(6)*b(5)+a(9)*b(8)
542 x(9)=a(3)*b(3)+a(6)*b(6)+a(9)*b(9)
543 RETURN

◆ prod_atb()

subroutine prod_atb ( a,
b,
x )

Definition at line 553 of file rini33.F.

554#include "implicit_f.inc"
555#include "param_c.inc"
556 my_real a(lskew),b(lskew),x(lskew)
557C
558 x(1)=a(1)*b(1)+a(2)*b(2)+a(3)*b(3)
559 x(2)=a(4)*b(1)+a(5)*b(2)+a(6)*b(3)
560 x(3)=a(7)*b(1)+a(8)*b(2)+a(9)*b(3)
561 x(4)=a(1)*b(4)+a(2)*b(5)+a(3)*b(6)
562 x(5)=a(4)*b(4)+a(5)*b(5)+a(6)*b(6)
563 x(6)=a(7)*b(4)+a(8)*b(5)+a(9)*b(6)
564 x(7)=a(1)*b(7)+a(2)*b(8)+a(3)*b(9)
565 x(8)=a(4)*b(7)+a(5)*b(8)+a(6)*b(9)
566 x(9)=a(7)*b(7)+a(8)*b(8)+a(9)*b(9)
567 RETURN

◆ qrot33()

subroutine qrot33 ( skew,
t,
c,
s )

Definition at line 405 of file rini33.F.

406C-------------------------------------------------------------------------
407C I m p l i c i t T y p e s
408C-----------------------------------------------
409#include "implicit_f.inc"
410C-----------------------------------------------
411C C o m m o n B l o c k s
412C-----------------------------------------------
413#include "param_c.inc"
414C----------------------------------------------------------
415C 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
416C----------------------------------------------------------
417 my_real t(3), skew(lskew)
418C-----------------------------------------------
419C L o c a l V a r i a b l e s
420C-----------------------------------------------
421 INTEGER I
422 my_real e11,e22,e33,e12,e21,e13,e31,e23,e32,
423 . u1,u2,u3, s,c,ci,si2
424C=======================================================================
425 ci = one - c
426 u1 = t(1)
427 u2 = t(2)
428 u3 = t(3)
429
430 e11 = u1 * u1 *ci + c
431 e22 = u2 * u2 *ci + c
432 e33 = u3 * u3 *ci + c
433
434 e12 = u1 * u2 * ci
435 e21 = e12 + u3 * s
436 e12 = e12 - u3 * s
437
438 e13 = u1 * u3 * ci
439 e31 = e13 - u2 * s
440 e13 = e13 + u2 * s
441
442 e23 = u2 * u3 * ci
443 e32 = e23 + u1 * s
444 e23 = e23 - u1 * s
445C
446 skew(1) = e11
447 skew(4) = e12
448 skew(7) = e13
449 skew(2) = e21
450 skew(5) = e22
451 skew(8) = e23
452 skew(3) = e31
453 skew(6 ) = e32
454 skew(9) = e33
455 RETURN

◆ rini33()

subroutine rini33 ( integer nel,
integer iout,
integer iprop,
integer, dimension(4,mvsiz) ix,
xl,
mass,
xiner,
stifn,
stifr,
viscm,
viscr,
uvar,
integer nuvar )

Definition at line 35 of file rini33.F.

38 USE message_mod
39C-------------------------------------------------------------------------
40C This subroutine initialize springs using user properties.
41C-------------------------------------------------------------------------
42C----------+---------+---+---+--------------------------------------------
43C VAR | SIZE |TYP| RW| DEFINITION
44C----------+---------+---+---+--------------------------------------------
45C IOUT | 1 | I | R | OUTPUT FILE UNIT (L00 file)
46C IPROP | 1 | I | R | PROPERTY NUMBER
47C----------+---------+---+---+--------------------------------------------
48C IX | 4*NEL | I | R | SPRING CONNECTIVITY
49C | IX(1,I) NODE 1 ID
50C | IX(2,I) NODE 2 ID
51C | IX(3,I) OPTIONAL NODE 3 ID
52C | IX(4,I) SPRING ID
53C----------+---------+---+---+--------------------------------------------
54C MASS | NEL | F | W | ELEMENT MASS
55C XINER | NEL | F | W | ELEMENT INERTIA (SPHERICAL)
56C STIFM | NEL | F | W | ELEMENT STIFNESS (TIME STEP)
57C STIFR | NEL | F | W | ELEMENT ROTATION STIFNESS (TIME STEP)
58C VISCM | NEL | F | W | ELEMENT VISCOSITY (TIME STEP)
59C VISCR | NEL | F | W | ELEMENT ROTATION VISCOSITY (TIME STEP)
60C----------+---------+---+---+--------------------------------------------
61C UVAR |NUVAR*NEL| F | W | USER ELEMENT VARIABLES
62C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
63C----------+---------+---+---+--------------------------------------------
64C I m p l i c i t T y p e s
65C-----------------------------------------------
66#include "implicit_f.inc"
67C-----------------------------------------------
68C G l o b a l P a r a m e t e r s
69C-----------------------------------------------
70#include "mvsiz_p.inc"
71C-----------------------------------------------
72C C o m m o n B l o c k s
73C-----------------------------------------------
74#include "param_c.inc"
75C----------------------------------------------------------
76C 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
77C----------------------------------------------------------
78 INTEGER NEL,IOUT,IPROP,NUVAR,IX(4,MVSIZ)
79 my_real
80 . mass(nel) ,xiner(nel) ,stifn(nel),xl(mvsiz,3) ,
81 . stifr(nel),viscm(nel) ,viscr(nel),uvar(nuvar,*)
82C-----------------------------------------------
83C L o c a l V a r i a b l e s
84C-----------------------------------------------
85 INTEGER I,IDSK1,IDSK2,JTYP,SKFLG,IFKNX,IFKNY,IFKNZ,
86 . IFKRX,IFKRY,IFKRZ,IFCNX,IFCNY,IFCNZ,IFCRX,IFCRY,IFCRZ,
87 . GET_U_PNU,GET_SKEW,KFUNC,KMAT,KPROP
88 my_real kxx,kyy,kzz,krx,kry,krz,knn,kr,x1,y1,z1,len2,
89 . k1,k2,k3,k4,k5,k6,c1,c2,c3,c4,c5,c6,ktt,krr,ctt,crr,
90 . cxx,cyy,czz,crx,cry,crz, deri,xf,get_u_func,
91 . u(lskew),v(lskew),a(lskew),b(lskew),ex(lskew),get_u_geo
92C-----------------------------------------------
93 EXTERNAL get_u_geo,get_skew
94 parameter(kfunc=29)
95 parameter(kmat=31)
96 parameter(kprop=33)
97C=======================================================================
98 jtyp = nint(get_u_geo(1,iprop))
99 idsk1= nint(get_u_geo(2,iprop))
100 idsk2= nint(get_u_geo(3,iprop))
101 skflg= nint(get_u_geo(14,iprop))
102 kxx = get_u_geo(4,iprop)
103 kyy = get_u_geo(5,iprop)
104 kzz = get_u_geo(6,iprop)
105 krx = get_u_geo(7,iprop)
106 kry = get_u_geo(8,iprop)
107 krz = get_u_geo(9,iprop)
108 knn = get_u_geo(10,iprop)
109 ifknx = get_u_pnu(1,iprop,kfunc)
110 ifkny = get_u_pnu(2,iprop,kfunc)
111 ifknz = get_u_pnu(3,iprop,kfunc)
112 ifkrx = get_u_pnu(4,iprop,kfunc)
113 ifkry = get_u_pnu(5,iprop,kfunc)
114 ifkrz = get_u_pnu(6,iprop,kfunc)
115C----
116 k1 = kxx
117 k2 = kyy
118 k3 = kzz
119 k4 = krx
120 k5 = kry
121 k6 = krz
122 IF (ifknx/=0) THEN
123 xf = get_u_func(ifknx,zero,deri)
124 k1 = max(kxx*deri, em20)
125 ENDIF
126 IF (ifkny/=0) THEN
127 xf = get_u_func(ifkny,zero,deri)
128 k2 = max(kyy*deri, em20)
129 ENDIF
130 IF (ifknz/=0) THEN
131 xf = get_u_func(ifknz,zero,deri)
132 k3 = max(kzz*deri, em20)
133 ENDIF
134 IF (ifkrx/=0) THEN
135 xf = get_u_func(ifkrx,zero,deri)
136 k4 = max(krx*deri, em20)
137 ENDIF
138 IF (ifkry/=0) THEN
139 xf = get_u_func(ifkry,zero,deri)
140 k5 = max(kry*deri, em20)
141 ENDIF
142 IF (ifkrz/=0) THEN
143 xf = get_u_func(ifkrz,zero,deri)
144 k6 = max(krz*deri, em20)
145 ENDIF
146 cxx = get_u_geo(21,iprop)
147 cyy = get_u_geo(22,iprop)
148 czz = get_u_geo(23,iprop)
149 crx = get_u_geo(24,iprop)
150 cry = get_u_geo(25,iprop)
151 crz = get_u_geo(26,iprop)
152C
153 ifcnx = get_u_pnu(7,iprop,kfunc)
154 ifcny = get_u_pnu(8,iprop,kfunc)
155 ifcnz = get_u_pnu(9,iprop,kfunc)
156 ifcrx = get_u_pnu(10,iprop,kfunc)
157 ifcry = get_u_pnu(11,iprop,kfunc)
158 ifcrz = get_u_pnu(12,iprop,kfunc)
159C
160 c1 = cxx
161 c2 = cyy
162 c3 = czz
163 c4 = crx
164 c5 = cry
165 c6 = crz
166 IF (ifcnx/=0) THEN
167 xf = get_u_func(ifcnx,zero,deri)
168 c1 = max(cxx*deri, em20)
169 ENDIF
170 IF (ifcny/=0) THEN
171 xf = get_u_func(ifcny,zero,deri)
172 c2 = max(cyy*deri, em20)
173 ENDIF
174 IF (ifcnz/=0) THEN
175 xf = get_u_func(ifcnz,zero,deri)
176 c3 = max(czz*deri, em20)
177 ENDIF
178 IF (ifcrx/=0) THEN
179 xf = get_u_func(ifcrx,zero,deri)
180 c4 = max(crx*deri, em20)
181 ENDIF
182 IF (ifcry/=0) THEN
183 xf = get_u_func(ifcry,zero,deri)
184 c5 = max(cry*deri, em20)
185 ENDIF
186 IF (ifcrz/=0) THEN
187 xf = get_u_func(ifcrz,zero,deri)
188 c6 = max(crz*deri, em20)
189 ENDIF
190 ktt = max(k1,k2,k3)
191 krr = max(k4,k5,k6)
192 ctt = max(c1,c2,c3)
193 crr = max(c4,c5,c6)
194C------- local frame
195 ierr=ierr+get_skew(iout,jtyp,skflg,idsk1,idsk2,u,v,ex,a,b)
196 DO i=1,nel
197 x1 = xl(i,1)
198 y1 = xl(i,2)
199 z1 = xl(i,3)
200 xl(i,1)=ex(1)*x1+ex(2)*y1+ex(3)*z1
201 xl(i,2)=ex(4)*x1+ex(5)*y1+ex(6)*z1
202 xl(i,3)=ex(7)*x1+ex(8)*y1+ex(9)*z1
203 ENDDO
204C--------------------------------------
205C ELEMENT INITIALIZATION
206C--------------------------------------
207 DO i=1,nel
208 mass(i) = zero
209 xiner(i) = zero
210 uvar(1,i) = xl(i,1)
211 uvar(2,i) = xl(i,2)
212 uvar(3,i) = xl(i,3)
213 len2=xl(i,1)*xl(i,1)+xl(i,2)*xl(i,2)+xl(i,3)*xl(i,3)
214 uvar(4,i) = a(1)
215 uvar(5,i) = a(2)
216 uvar(6,i) = a(3)
217 uvar(7,i) = a(4)
218 uvar(8,i) = a(5)
219 uvar(9,i) = a(6)
220 uvar(10,i)= a(7)
221 uvar(11,i)= a(8)
222 uvar(12,i)= a(9)
223 uvar(22,i)= ex(1)
224 uvar(23,i)= ex(2)
225 uvar(24,i)= ex(3)
226 uvar(25,i)= ex(4)
227 uvar(26,i)= ex(5)
228 uvar(27,i)= ex(6)
229 uvar(28,i)= ex(7)
230 uvar(29,i)= ex(8)
231 uvar(30,i)= ex(9)
232C
233 kr = knn*max(one,len2)
234 uvar(19,i)= kxx
235 uvar(20,i)= kyy
236 uvar(21,i)= kzz
237
238 IF(jtyp>=2.AND.jtyp<=4) THEN
239 uvar(31,i)= krx
240 uvar(32,i)= kr
241 uvar(33,i)= kr
242 ELSEIF(jtyp==5) THEN
243 uvar(31,i)= kr
244 uvar(32,i)= kry
245 uvar(33,i)= krz
246 ELSEIF(jtyp>=6.AND.jtyp<=8) THEN
247 uvar(31,i)= kr
248 uvar(32,i)= kr
249 uvar(33,i)= kr
250 ELSE
251 uvar(31,i)= krx
252 uvar(32,i)= kry
253 uvar(33,i)= krz
254 ENDIF
255C
256 uvar(34,i)= zero
257 uvar(35,i)= zero
258 uvar(36,i)= zero
259 uvar(37,i)= zero
260 uvar(38,i)= zero
261 uvar(39,i)= zero
262C
263 stifn(i) = ktt
264 stifr(i) = krr+ktt*len2
265 viscm(i) = ctt
266 viscr(i) = crr
267 ENDDO
268C
269 RETURN
max
#define max(a, b)
Definition macros.h:21
get_u_pnu
integer function get_u_pnu(ivar, ip, k)
Definition uaccess.F:482

◆ rot12()

subroutine rot12 ( skew,
t,
c,
s )

Definition at line 462 of file rini33.F.

463C-------------------------------------------------------------------------
464C I m p l i c i t T y p e s
465C-----------------------------------------------
466#include "implicit_f.inc"
467C-----------------------------------------------
468C C o m m o n B l o c k s
469C-----------------------------------------------
470#include "param_c.inc"
471C----------------------------------------------------------
472C 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
473C----------------------------------------------------------
474 my_real rot(3), t(3), skew(lskew), ksi, s, c
475C-----------------------------------------------
476C L o c a l V a r i a b l e s
477C-----------------------------------------------
478 INTEGER I
479 my_real e,e11,e22,e33,e12,e21,e13,e31,e23,e32,si2,nr,
480 . umi
481 DATA umi/-1.0/
482C=======================================================================
483C
484 pi = 2.*atan2(one,zero)
485C--- first skew rotation
486 e11 = skew(1)
487 e12 = skew(4)
488 e13 = skew(7)
489 e21 = skew(2)
490 e22 = skew(5)
491 e23 = skew(8)
492 e31 = skew(3)
493 e32 = skew(6)
494 e33 = skew(9)
495 e = e11+e22+e33
496 c = half * (e - 1.)
497 c = min(c,one)
498 c = max(c,umi)
499 ksi = acos(c)
500 s = sin(ksi)
501 IF(s==zero) THEN
502 si2 = zero
503 ELSE
504 si2 = half / s
505 ENDIF
506 t(1) = (e32 - e23) * si2
507 t(2) = (e13 - e31) * si2
508 t(3) = (e21 - e12) * si2
509 nr = sqrt(t(1)*t(1)+t(2)*t(2)+t(3)*t(3))
510 IF (nr/=zero) nr = one/nr
511 t(1) = t(1)*nr
512 t(2) = t(2)*nr
513 t(3) = t(3)*nr
514C
515 c = half*(c+ one)
516 s = sqrt(one-c)
517 c = sqrt(c)
518C
519 RETURN
min
#define min(a, b)
Definition macros.h:20