OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
rini45.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!|| rini45 ../starter/source/elements/joint/rjoint/rini45.F
25!||--- called by ------------------------------------------------------
26!|| rinit3 ../starter/source/elements/spring/rinit3.F
27!||--- calls -----------------------------------------------------
28!|| ancmsg ../starter/source/output/message/message.F
29!|| get_skew45 ../starter/source/elements/joint/rjoint/rini45.F
30!|| get_u_func ../starter/source/user_interface/uaccess.F
31!|| get_u_func_deri ../starter/source/user_interface/uaccess.F
32!|| get_u_geo ../starter/source/user_interface/uaccess.F
33!|| get_u_pnu ../starter/source/user_interface/uaccess.F
34!|| init_skew45 ../starter/source/elements/joint/rjoint/rini45.F
35!||--- uses -----------------------------------------------------
36!|| message_mod ../starter/share/message_module/message_mod.F
37!||====================================================================
38 SUBROUTINE rini45(NEL ,IOUT ,IPROP , IX, X, XL ,
39 3 MASS ,XINER ,STIFN , STIFR ,VISCM ,
40 4 VISCR,UVAR ,NUVAR,IXR, IXR_KJ,ID,TITR)
41C-----------------------------------------------
42C M o d u l e s
43C-----------------------------------------------
44 USE message_mod
45 USE names_and_titles_mod , ONLY : nchartitle
46 use element_mod , only : nixr
47C-------------------------------------------------------------------------
48C This subroutine initialize springs using user properties.
49C-------------------------------------------------------------------------
50C----------+---------+---+---+--------------------------------------------
51C VAR | SIZE |TYP| RW| DEFINITION
52C----------+---------+---+---+--------------------------------------------
53C IOUT | 1 | I | R | OUTPUT FILE UNIT (L00 file)
54C IPROP | 1 | I | R | PROPERTY NUMBER
55C----------+---------+---+---+--------------------------------------------
56C IX | 4*NEL | I | R | SPRING CONNECTIVITY
57C | IX(1,I) NODE 1 ID
58C | IX(2,I) NODE 2 ID
59C | IX(3,I) OPTIONAL NODE 3 ID
60C | IX(4,I) SPRING ID
61C----------+---------+---+---+--------------------------------------------
62C MASS | NEL | F | W | ELEMENT MASS
63C XINER | NEL | F | W | ELEMENT INERTIA (SPHERICAL)
64C STIFM | NEL | F | W | ELEMENT STIFNESS (TIME STEP)
65C STIFR | NEL | F | W | ELEMENT ROTATION STIFNESS (TIME STEP)
66C VISCM | NEL | F | W | ELEMENT VISCOSITY (TIME STEP)
67C VISCR | NEL | F | W | ELEMENT ROTATION VISCOSITY (TIME STEP)
68C----------+---------+---+---+--------------------------------------------
69C UVAR |NUVAR*NEL| F | W | USER ELEMENT VARIABLES
70C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
71C----------+---------+---+---+--------------------------------------------
72C I m p l i c i t T y p e s
73C-----------------------------------------------
74#include "implicit_f.inc"
75C-----------------------------------------------
76C G l o b a l P a r a m e t e r s
77C-----------------------------------------------
78#include "mvsiz_p.inc"
79C-----------------------------------------------
80C C o m m o n B l o c k s
81C-----------------------------------------------
82#include "vect01_c.inc"
83#include "param_c.inc"
84C----------------------------------------------------------
85C 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
86C----------------------------------------------------------
87 INTEGER NEL,IOUT,IPROP,NUVAR,IX(4,MVSIZ),IXR(NIXR,*),
88 . IXR_KJ(5,*)
89 my_real
90 . mass(nel) ,xiner(nel) ,stifn(nel),xl(mvsiz,3) ,
91 . stifr(nel),viscm(nel) ,viscr(nel),uvar(nuvar,*),
92 . x(3,*)
93 INTEGER ID
94 CHARACTER(LEN=NCHARTITLE) :: TITR
95C-----------------------------------------------
96C L o c a l V a r i a b l e s
97C-----------------------------------------------
98 INTEGER I,IDSK1,IDSK2,JTYP,SKFLG,IFKNX,IFKNY,IFKNZ,
99 . IFKRX,IFKRY,IFKRZ,IFCNX,IFCNY,IFCNZ,IFCRX,IFCRY,IFCRZ,
100 . get_u_pnu,get_skew45,kfunc,il,ifm
101 my_real kxx,kyy,kzz,krx,kry,krz,knn,kr,x1,y1,z1,len2,
102 . k1,k2,k3,k4,k5,k6,c1,c2,c3,c4,c5,c6,ktt,krr,ctt,crr,
103 . cxx,cyy,czz,crx,cry,crz, deri,xf,get_u_func,fm,
104 . ex(lskew),get_u_geo,scf,get_u_func_deri,derimax,kfm,
105 . theta0(3),vect1(lskew),vect2(lskew),stop_angl_min(3),
106 . stop_angl_max(3)
107C-----------------------------------------------
108 EXTERNAL get_u_geo,get_skew45,get_u_func_deri
109 parameter(kfunc=29)
110C=======================================================================
111 jtyp = nint(get_u_geo(1,iprop))
112C
113 kxx = get_u_geo(4,iprop)
114 kyy = get_u_geo(5,iprop)
115 kzz = get_u_geo(6,iprop)
116 krx = get_u_geo(7,iprop)
117 kry = get_u_geo(8,iprop)
118 krz = get_u_geo(9,iprop)
119 knn = get_u_geo(10,iprop)
120 scf = get_u_geo(11,iprop)
121 DO i=1,3
122 stop_angl_min(i) = get_u_geo(35+2*(i-1),iprop)
123 stop_angl_max(i) = get_u_geo(36+2*(i-1),iprop)
124 ENDDO
125 ifknx = get_u_pnu(1,iprop,kfunc)
126 ifkny = get_u_pnu(2,iprop,kfunc)
127 ifknz = get_u_pnu(3,iprop,kfunc)
128 ifkrx = get_u_pnu(4,iprop,kfunc)
129 ifkry = get_u_pnu(5,iprop,kfunc)
130 ifkrz = get_u_pnu(6,iprop,kfunc)
131C----
132 k1 = kxx
133 k2 = kyy
134 k3 = kzz
135 k4 = krx
136 k5 = kry
137 k6 = krz
138 IF (ifknx/=0) THEN
139 xf = get_u_func(ifknx,zero,deri)
140 k1 = max(kxx*deri, em20)
141 ENDIF
142 IF (ifkny/=0) THEN
143 xf = get_u_func(ifkny,zero,deri)
144 k2 = max(kyy*deri, em20)
145 ENDIF
146 IF (ifknz/=0) THEN
147 xf = get_u_func(ifknz,zero,deri)
148 k3 = max(kzz*deri, em20)
149 ENDIF
150 IF (ifkrx/=0) THEN
151 xf = get_u_func(ifkrx,zero,deri)
152 k4 = max(krx*deri, em20)
153 ENDIF
154 IF (ifkry/=0) THEN
155 xf = get_u_func(ifkry,zero,deri)
156 k5 = max(kry*deri, em20)
157 ENDIF
158 IF (ifkrz/=0) THEN
159 xf = get_u_func(ifkrz,zero,deri)
160 k6 = max(krz*deri, em20)
161 ENDIF
162 cxx = get_u_geo(21,iprop)
163 cyy = get_u_geo(22,iprop)
164 czz = get_u_geo(23,iprop)
165 crx = get_u_geo(24,iprop)
166 cry = get_u_geo(25,iprop)
167 crz = get_u_geo(26,iprop)
168C
169 ifcnx = get_u_pnu(7,iprop,kfunc)
170 ifcny = get_u_pnu(8,iprop,kfunc)
171 ifcnz = get_u_pnu(9,iprop,kfunc)
172 ifcrx = get_u_pnu(10,iprop,kfunc)
173 ifcry = get_u_pnu(11,iprop,kfunc)
174 ifcrz = get_u_pnu(12,iprop,kfunc)
175C
176 c1 = cxx
177 c2 = cyy
178 c3 = czz
179 c4 = crx
180 c5 = cry
181 c6 = crz
182 IF (ifcnx/=0) THEN
183 xf = get_u_func(ifcnx,zero,deri)
184 c1 = max(cxx*deri, em20)
185 ENDIF
186 IF (ifcny/=0) THEN
187 xf = get_u_func(ifcny,zero,deri)
188 c2 = max(cyy*deri, em20)
189 ENDIF
190 IF (ifcnz/=0) THEN
191 xf = get_u_func(ifcnz,zero,deri)
192 c3 = max(czz*deri, em20)
193 ENDIF
194 IF (ifcrx/=0) THEN
195 xf = get_u_func(ifcrx,zero,deri)
196 c4 = max(crx*deri, em20)
197 ENDIF
198 IF (ifcry/=0) THEN
199 xf = get_u_func(ifcry,zero,deri)
200 c5 = max(cry*deri, em20)
201 ENDIF
202 IF (ifcrz/=0) THEN
203 xf = get_u_func(ifcrz,zero,deri)
204 c6 = max(crz*deri, em20)
205 ENDIF
206 DO i=1,6
207 ifm = get_u_pnu(12+i,iprop,kfunc)
208 kfm = get_u_geo(40+i,iprop)
209 fm = get_u_geo(46+i,iprop)
210 IF (ifm/=0) THEN
211 derimax = fm*get_u_func_deri(ifm)
212 IF (derimax>kfm) THEN
213 CALL ancmsg(msgid=979,
214 . msgtype=msgwarning,
215 . anmode=aninfo_blind_2,
216 . r1=kfm,
217 . r2=derimax)
218 ENDIF
219 ENDIF
220 END DO
221
222 ktt = max(k1,k2,k3)
223 krr = max(k4,k5,k6)
224 ctt = max(c1,c2,c3)
225 crr = max(c4,c5,c6)
226C--------------------------------------
227C INITIAL ROTATION AND LOCAL FRAME FROM SKEWS
228C--------------------------------------
229 CALL init_skew45(jtyp,iprop,idsk1,idsk2,vect1,vect2,id,titr)
230C
231C--------------------------------------
232C ELEMENT INITIALIZATION
233C--------------------------------------
234 DO i=1,nel
235 mass(i) = zero
236 xiner(i) = zero
237C
238 il = lft+i-1+nft
239C
240C------- computation of local frame
241 ierr=ierr+get_skew45(iout,jtyp,ex,ixr,ixr_kj,il,x,id,titr,
242 . idsk1,idsk2,vect1,vect2,theta0,
243 . stop_angl_min,stop_angl_max)
244C
245 x1 = xl(i,1)
246 y1 = xl(i,2)
247 z1 = xl(i,3)
248 xl(i,1)=ex(1)*x1+ex(2)*y1+ex(3)*z1
249 xl(i,2)=ex(4)*x1+ex(5)*y1+ex(6)*z1
250 xl(i,3)=ex(7)*x1+ex(8)*y1+ex(9)*z1
251C
252 uvar(1,i) = xl(i,1)
253 uvar(2,i) = xl(i,2)
254 uvar(3,i) = xl(i,3)
255 len2=xl(i,1)*xl(i,1)+xl(i,2)*xl(i,2)+xl(i,3)*xl(i,3)
256C
257 IF (idsk2>0) THEN
258 uvar(7,i) = theta0(1)
259 uvar(8,i) = theta0(2)
260 uvar(9,i) = theta0(3)
261 ELSE
262 uvar(7,i) = zero
263 uvar(8,i) = zero
264 uvar(9,i) = zero
265 ENDIF
266C
267 uvar(22,i) = ex(1)
268 uvar(23,i) = ex(2)
269 uvar(24,i) = ex(3)
270 uvar(25,i) = ex(4)
271 uvar(26,i) = ex(5)
272 uvar(27,i) = ex(6)
273 uvar(28,i) = ex(7)
274 uvar(29,i) = ex(8)
275 uvar(30,i) = ex(9)
276C
277C------- storage of stiffness
278 kr = knn*max(scf,len2)
279 uvar(17,i)= knn
280 uvar(18,i)= kr
281 uvar(19,i)= kxx
282 uvar(20,i)= kyy
283 uvar(21,i)= kzz
284C
285 IF(jtyp>=2.AND.jtyp<=4) THEN
286 uvar(31,i)= krx
287 uvar(32,i)= kr
288 uvar(33,i)= kr
289 ELSEIF(jtyp==5) THEN
290 uvar(31,i)= kr
291 uvar(32,i)= kry
292 uvar(33,i)= krz
293 uvar(10,i)= ex(4)
294 uvar(11,i)= ex(5)
295 uvar(12,i)= ex(6)
296 uvar(13,i)= ex(7)
297 uvar(14,i)= ex(8)
298 uvar(15,i)= ex(9)
299 ELSEIF(jtyp>=6.AND.jtyp<=8) THEN
300 uvar(31,i)= kr
301 uvar(32,i)= kr
302 uvar(33,i)= kr
303 ELSE
304 uvar(31,i)= krx
305 uvar(32,i)= kry
306 uvar(33,i)= krz
307 ENDIF
308C
309 uvar(34,i)= zero
310 uvar(35,i)= zero
311 uvar(36,i)= zero
312 uvar(37,i)= zero
313 uvar(38,i)= zero
314 uvar(39,i)= zero
315C
316 stifn(i) = ktt
317 stifr(i) = krr+ktt*len2
318 viscm(i) = ctt
319 viscr(i) = crr
320
321 ENDDO
322C
323 RETURN
324 END
325
326!||====================================================================
327!|| get_skew45 ../starter/source/elements/joint/rjoint/rini45.F
328!||--- called by ------------------------------------------------------
329!|| rini45 ../starter/source/elements/joint/rjoint/rini45.F
330!||--- calls -----------------------------------------------------
331!|| ancmsg ../starter/source/output/message/message.F
332!|| prod_atb ../starter/source/elements/joint/rjoint/rini33.F
333!||--- uses -----------------------------------------------------
334!|| message_mod ../starter/share/message_module/message_mod.F
335!||====================================================================
336 INTEGER FUNCTION get_skew45(IOUT,JTYP,EX,IXR,IXR_KJ,IEL,X,ID,TITR,
337 . IDSK1,IDSK2,VECT1,VECT2,THETA0,
338 . STOP_ANGL_MIN,STOP_ANGL_MAX)
339C-----------------------------------------------
340C M o d u l e s
341C-----------------------------------------------
342 USE message_mod
343 USE names_and_titles_mod , ONLY : nchartitle
344 use element_mod , only : nixr
345C-----------------------------------------------
346C I m p l i c i t T y p e s
347C-----------------------------------------------
348#include "implicit_f.inc"
349C-----------------------------------------------
350C A n a l y s e M o d u l e
351C-----------------------------------------------
352#include "param_c.inc"
353#include "com04_c.inc"
354C----------------------------------------------------------
355C 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
356C----------------------------------------------------------
357 INTEGER iout,ixr(NIXR,*),iel,jtyp,ixr_kj(5,*)
358 my_real ex(lskew),x(3,*)
359 INTEGER id
360 INTEGER, INTENT(INOUT) :: idsk1,idsk2
361 my_real, INTENT(INOUT) :: vect1(lskew),vect2(lskew)
362 my_real, INTENT(INOUT) :: theta0(3)
363 my_real, INTENT(INOUT) :: stop_angl_min(3),stop_angl_max(3)
364 CHARACTER(LEN=NCHARTITLE) :: titr
365C-----------------------------------------------
366C L o c a l V a r i a b l e s
367C-----------------------------------------------
368 INTEGER i,j,k,nnod,n1,n2,n3,n4,ierr1,ielusr,nn(6),hh
369 INTEGER numel_kj,id_kj,nnod2,nnod_req(9),nb_rot
370 my_real pp1,pp2,pp3,pp4,len,scal,exmax
371 my_real vect1_upt(lskew),vect2_upt(lskew),q(lskew),nr,t(3),
372 . cosk,sink,si2,e,ksi,umi,err,scal_sign
373 DATA umi/-1.0/
374C=======================================================================
375C
376 ierr1 = 0
377 nnod = 0
378 nn = 0
379 n1 = 0
380 n2 = 0
381 n3 = 0
382 n4 = 0
383 numel_kj = ixr_kj(1,numelr+1)
384 ielusr = ixr(nixr,iel)
385C---- Min number of nodes per type of joint
386 nnod_req(1) = 2
387 nnod_req(2) = 3
388 nnod_req(3) = 3
389 nnod_req(4) = 3
390 nnod_req(5) = 4
391 nnod_req(6) = 3
392 nnod_req(7) = 3
393 nnod_req(8) = 2
394 nnod_req(9) = 4
395C---- Check of real nodes
396 DO j=1,3
397 IF (ixr(1+j,iel)/=0) THEN
398 nnod = nnod + 1
399 nn(nnod) = ixr(1+j,iel)
400 ENDIF
401 END DO
402C---- Check of additional nodes
403 DO j=1,numel_kj
404 IF (ixr_kj(4,j)==ielusr) id_kj = j
405 END DO
406
407 IF (id_kj>0) THEN
408 DO j=1,3
409 IF (ixr_kj(j,id_kj)/=0) THEN
410 nnod = nnod + 1
411 nn(nnod) = ixr_kj(j,id_kj)
412 ENDIF
413 END DO
414 ENDIF
415
416 nnod2 = nnod
417C
418 len = sqrt((x(1,nn(1))-x(1,nn(2)))**2+(x(2,nn(1))-x(2,nn(2)))**2
419 . +(x(3,nn(1))-x(3,nn(2)))**2)
420 IF ((nnod==2).AND.(len>em10)) nnod2=3
421C
422 IF ((nnod2<nnod_req(jtyp)).AND.(idsk1 == 0)) THEN
423C-----------------------------------------------
424C--- Error - no skew can be defined
425C-----------------------------------------------
426 CALL ancmsg(msgid=936,
427 . msgtype=msgerror,
428 . anmode=aninfo_blind_2,
429 . i1=id,
430 . c1=titr,
431 . i2=ielusr,
432 . i3=jtyp,
433 . i4=nnod_req(jtyp)-nnod2)
434 ELSEIF ((nnod==2).AND.((jtyp==1).OR.(jtyp==8)).AND.(idsk1==0)) THEN
435C-----------------------------------------------
436C--- Global frame is used
437C-----------------------------------------------
438 ex(1)= one
439 ex(2)= zero
440 ex(3)= zero
441 ex(4)= zero
442 ex(5)= one
443 ex(6)= zero
444 ex(7)= zero
445 ex(8)= zero
446 ex(9)= one
447 pp1 = one
448 pp2 = one
449 pp3 = one
450 ELSEIF ((nnod==2).AND.(len>em10)) THEN
451C-----------------------------------------------
452C--- joint with 1 axis - axis defined with 2nd node
453C-----------------------------------------------
454 n1 = nn(1)
455 n2 = nn(2)
456C-- computation of x
457 ex(1)=x(1,n2)-x(1,n1)
458 ex(2)=x(2,n2)-x(2,n1)
459 ex(3)=x(3,n2)-x(3,n1)
460 pp1=sqrt(ex(1)*ex(1)+ex(2)*ex(2)+ex(3)*ex(3))
461 exmax =zero
462 DO j=1,3
463 IF (abs(ex(j))>= exmax) THEN
464 exmax = abs(ex(j))
465 hh = j
466 ENDIF
467 END DO
468C-- computation of y
469 IF (hh<3) THEN
470 ex(4)= -ex(2)
471 ex(5)= ex(1)
472 ex(6)= zero
473 ELSE
474 ex(4)= zero
475 ex(5)= ex(3)
476 ex(6)= -ex(2)
477 ENDIF
478 pp2=sqrt(ex(4)*ex(4)+ex(5)*ex(5)+ex(6)*ex(6))
479C-- computation of z = x ^ y
480 ex(7)=ex(2)*ex(6)-ex(3)*ex(5)
481 ex(8)=ex(3)*ex(4)-ex(1)*ex(6)
482 ex(9)=ex(1)*ex(5)-ex(2)*ex(4)
483 pp3=sqrt(ex(7)*ex(7)+ex(8)*ex(8)+ex(9)*ex(9))
484 ELSEIF (nnod==3) THEN
485C-----------------------------------------------
486C--- joint with 1 axis - axis defined with 3rd node
487C-----------------------------------------------
488 n1 = nn(1)
489 n2 = nn(3)
490C-- computation of x
491 ex(1)=x(1,n2)-x(1,n1)
492 ex(2)=x(2,n2)-x(2,n1)
493 ex(3)=x(3,n2)-x(3,n1)
494 pp1=sqrt(ex(1)*ex(1)+ex(2)*ex(2)+ex(3)*ex(3))
495 exmax =zero
496 DO j=1,3
497 IF (abs(ex(j)) >= exmax) THEN
498 exmax = abs(ex(j))
499 hh = j
500 ENDIF
501 END DO
502C-- computation of y
503 IF (hh<3) THEN
504 ex(4)= -ex(2)
505 ex(5)= ex(1)
506 ex(6)= zero
507 ELSE
508 ex(4)= zero
509 ex(5)= ex(3)
510 ex(6)= -ex(2)
511 ENDIF
512 pp2=sqrt(ex(4)*ex(4)+ex(5)*ex(5)+ex(6)*ex(6))
513C-- computation of z = x ^ y
514 ex(7)=ex(2)*ex(6)-ex(3)*ex(5)
515 ex(8)=ex(3)*ex(4)-ex(1)*ex(6)
516 ex(9)=ex(1)*ex(5)-ex(2)*ex(4)
517 pp3=sqrt(ex(7)*ex(7)+ex(8)*ex(8)+ex(9)*ex(9))
518C-- control of axis's length
519 IF (pp1<=em10) THEN
520 CALL ancmsg(msgid=935,
521 . msgtype=msgerror,
522 . anmode=aninfo_blind_2,
523 . i1=id,
524 . c1=titr,
525 . i2=ielusr)
526 ENDIF
527 ELSEIF ((nnod>=4).AND.(jtyp/=5)) THEN
528C-----------------------------------------------
529C--- General case : skew computed with 3rd and 4th node
530C-----------------------------------------------
531 n1 = nn(1)
532 n2 = nn(3)
533 n3 = nn(4)
534C-- computation of x
535 ex(1)=x(1,n2)-x(1,n1)
536 ex(2)=x(2,n2)-x(2,n1)
537 ex(3)=x(3,n2)-x(3,n1)
538 pp1=sqrt(ex(1)*ex(1)+ex(2)*ex(2)+ex(3)*ex(3))
539C-- computation of y
540 ex(4)=x(1,n3)-x(1,n1)
541 ex(5)=x(2,n3)-x(2,n1)
542 ex(6)=x(3,n3)-x(3,n1)
543 pp2=sqrt(ex(4)*ex(4)+ex(5)*ex(5)+ex(6)*ex(6))
544C-- computation of z = x ^ y
545 ex(7)=ex(2)*ex(6)-ex(3)*ex(5)
546 ex(8)=ex(3)*ex(4)-ex(1)*ex(6)
547 ex(9)=ex(1)*ex(5)-ex(2)*ex(4)
548 pp3=sqrt(ex(7)*ex(7)+ex(8)*ex(8)+ex(9)*ex(9))
549C-- control of orthogonality
550 scal =abs(ex(1)*ex(4)+ex(2)*ex(5)+ex(3)*ex(6))/(pp1*pp2)
551 IF (abs(scal)>=0.98) THEN
552 CALL ancmsg(msgid=1009,
553 . msgtype=msgerror,
554 . anmode=aninfo_blind_2,
555 . i1=id,
556 . c1=titr,
557 . i2=ielusr)
558 ELSE
559 ex(4)=ex(8)*ex(3)-ex(9)*ex(2)
560 ex(5)=ex(9)*ex(1)-ex(7)*ex(3)
561 ex(6)=ex(7)*ex(2)-ex(8)*ex(1)
562 pp2=sqrt(ex(4)*ex(4)+ex(5)*ex(5)+ex(6)*ex(6))
563 ENDIF
564C-- control of N4
565 IF ((n4==n1).OR.(n4==n2).OR.(n4==n3)) THEN
566 CALL ancmsg(msgid=681,
567 . msgtype=msgerror,
568 . anmode=aninfo_blind_2,
569 . i1=ielusr)
570 ENDIF
571C-- control of axis length
572 pp4 = sqrt((x(1,n3)-x(1,n2))**2+(x(2,n3)-x(2,n2))**2
573 . +(x(3,n3)-x(3,n2))**2)
574 IF ((pp1<=em10).OR.(pp2<=em10).OR.(pp4<=em10)) THEN
575 CALL ancmsg(msgid=934,
576 . msgtype=msgerror,
577 . anmode=aninfo_blind_2,
578 . i1=id,
579 . c1=titr,
580 . i2=ielusr)
581 ENDIF
582 ELSEIF ((nnod>=4).AND.(jtyp==5)) THEN
583C-----------------------------------------------
584C--- Universal Joint : local skew
585C-----------------------------------------------
586 n1 = nn(1)
587 n2 = nn(3)
588 n3 = nn(4)
589C-- computation of y
590 ex(4)=x(1,n2)-x(1,n1)
591 ex(5)=x(2,n2)-x(2,n1)
592 ex(6)=x(3,n2)-x(3,n1)
593 pp2=sqrt(ex(4)*ex(4)+ex(5)*ex(5)+ex(6)*ex(6))
594C-- computation of z
595 ex(7)=x(1,n3)-x(1,n1)
596 ex(8)=x(2,n3)-x(2,n1)
597 ex(9)=x(3,n3)-x(3,n1)
598 pp3=sqrt(ex(7)*ex(7)+ex(8)*ex(8)+ex(9)*ex(9))
599C-- computation of x = y ^ z
600 ex(1)=ex(5)*ex(9)-ex(6)*ex(8)
601 ex(2)=ex(6)*ex(7)-ex(4)*ex(9)
602 ex(3)=ex(4)*ex(8)-ex(5)*ex(7)
603 pp1=sqrt(ex(1)*ex(1)+ex(2)*ex(2)+ex(3)*ex(3))
604C-- control of orthogonality
605 scal =abs(ex(7)*ex(4)+ex(8)*ex(5)+ex(9)*ex(6))/(pp1+pp2)
606 IF (abs(scal)>=0.98) THEN
607 CALL ancmsg(msgid=1009,
608 . msgtype=msgerror,
609 . anmode=aninfo_blind_2,
610 . i1=id,
611 . c1=titr,
612 . i2=ielusr)
613 ELSEIF (scal>=1e-4) THEN
614C-- Orthogonalization of the skew
615 CALL ancmsg(msgid=940,
616 . msgtype=msgwarning,
617 . anmode=aninfo_blind_2,
618 . i1=id,
619 . c1=titr,
620 . i2=ielusr)
621C-- computation of z = x ^ y
622 ex(7)=ex(2)*ex(6)-ex(3)*ex(5)
623 ex(8)=ex(3)*ex(4)-ex(1)*ex(6)
624 ex(9)=ex(1)*ex(5)-ex(2)*ex(4)
625 pp3=sqrt(ex(7)*ex(7)+ex(8)*ex(8)+ex(9)*ex(9))
626 ENDIF
627C-- control of axis length
628 pp4 = sqrt((x(1,n3)-x(1,n2))**2+(x(2,n3)-x(2,n2))**2
629 . +(x(3,n3)-x(3,n2))**2)
630 IF ((pp1<=em10).OR.(pp2<=em10).OR.(pp4<=em10)) THEN
631 CALL ancmsg(msgid=934,
632 . msgtype=msgerror,
633 . anmode=aninfo_blind_2,
634 . i1=id,
635 . c1=titr,
636 . i2=ielusr)
637 ENDIF
638 ELSEIF (idsk1 > 0) THEN
639C-----------------------------------------------
640C--- skew1 is used
641C-----------------------------------------------
642 ex(1:9)= vect1(1:9)
643 pp1 = one
644 pp2 = one
645 pp3 = one
646 ELSE
647 CALL ancmsg(msgid=937,
648 . msgtype=msgerror,
649 . anmode=aninfo_blind_2,
650 . i1=id,
651 . c1=titr,
652 . i2=ielusr)
653 ENDIF
654C-----------------------------------------------
655C-- norm
656C-----------------------------------------------
657 ex(1)=ex(1)/pp1
658 ex(2)=ex(2)/pp1
659 ex(3)=ex(3)/pp1
660 ex(4)=ex(4)/pp2
661 ex(5)=ex(5)/pp2
662 ex(6)=ex(6)/pp2
663 ex(7)=ex(7)/pp3
664 ex(8)=ex(8)/pp3
665 ex(9)=ex(9)/pp3
666C-----------------------------------------------
667C-- Initialisation of rotations
668C-----------------------------------------------
669C
670 IF ((idsk1 > 0).AND.(idsk2 > 2)) THEN
671C
672 nb_rot = 3
673 IF ((jtyp==2).OR.(jtyp==3).OR.(jtyp==4)) THEN
674C-- only one init angle allowed on X axis
675 nb_rot = 1
676C check of alignement of skew1
677 scal_sign = sign(one,ex(1)*vect1(1)+ex(2)*vect1(2)+ex(3)*vect1(3))
678 scal = abs(ex(1)*vect1(1)+ex(2)*vect1(2)+ex(3)*vect1(3))
679 IF (scal.LT.0.98) THEN
680C-- misalignement - error
681 CALL ancmsg(msgid=3076,
682 . msgtype=msgerror,
683 . anmode=aninfo_blind_1,
684 . i1=id,
685 . c1=titr,
686 . i2=ielusr,
687 . c2='SKEW1')
688 ELSE
689 IF ((one-scal).GT.em05) THEN
690C-- misalignement - warning
691 CALL ancmsg(msgid=3067,
692 . msgtype=msgwarning,
693 . anmode=aninfo_blind_2,
694 . i1=id,
695 . c1=titr,
696 . i2=ielusr,
697 . c2='SKEW1')
698 ENDIF
699C-- realignement of skew1
700 vect1_upt(1:3) = scal_sign*ex(1:3)
701C-- z1' = x1'^y1
702 vect1_upt(7)=vect1_upt(2)*vect1(6)-vect1_upt(3)*vect1(5)
703 vect1_upt(8)=vect1_upt(3)*vect1(4)-vect1_upt(1)*vect1(6)
704 vect1_upt(9)=vect1_upt(1)*vect1(5)-vect1_upt(2)*vect1(4)
705C-- y1' = z1'^x1'
706 vect1_upt(4)=vect1_upt(8)*vect1_upt(3)-vect1_upt(9)*vect1_upt(2)
707 vect1_upt(5)=vect1_upt(9)*vect1_upt(1)-vect1_upt(7)*vect1_upt(3)
708 vect1_upt(6)=vect1_upt(7)*vect1_upt(2)-vect1_upt(8)*vect1_upt(1)
709 vect1(1:9) = vect1_upt(1:9)
710 ENDIF
711 scal_sign = sign(one,ex(1)*vect2(1)+ex(2)*vect2(2)+ex(3)*vect2(3))
712 scal = abs(ex(1)*vect2(1)+ex(2)*vect2(2)+ex(3)*vect2(3))
713 IF (scal.LT.0.98) THEN
714C-- misalignement - error
715 CALL ancmsg(msgid=3076,
716 . msgtype=msgerror,
717 . anmode=aninfo_blind_1,
718 . i1=id,
719 . c1=titr,
720 . i2=ielusr,
721 . c2='SKEW2')
722 ELSE
723 IF ((one-scal).GT.em05) THEN
724C-- misalignement - warning
725 CALL ancmsg(msgid=3067,
726 . msgtype=msgwarning,
727 . anmode=aninfo_blind_2,
728 . i1=id,
729 . c1=titr,
730 . i2=ielusr,
731 . c2='SKEW2')
732 ENDIF
733C-- realignement of skew2
734 vect2_upt(1:3) = scal_sign*ex(1:3)
735C-- z1' = x1'^y1
736 vect2_upt(7)=vect2_upt(2)*vect2(6)-vect2_upt(3)*vect2(5)
737 vect2_upt(8)=vect2_upt(3)*vect2(4)-vect2_upt(1)*vect2(6)
738 vect2_upt(9)=vect2_upt(1)*vect2(5)-vect2_upt(2)*vect2(4)
739C-- y1' = z1'^x1'
740 vect2_upt(4)=vect2_upt(8)*vect2_upt(3)-vect2_upt(9)*vect2_upt(2)
741 vect2_upt(5)=vect2_upt(9)*vect2_upt(1)-vect2_upt(7)*vect2_upt(3)
742 vect2_upt(6)=vect2_upt(7)*vect2_upt(2)-vect2_upt(8)*vect2_upt(1)
743 vect2(1:9) = vect2_upt(1:9)
744 ENDIF
745 ENDIF
746C
747C-- computation of rotation from skew1 to skew2
748 CALL prod_atb(vect1,vect2,q)
749C
750 e = q(1)+q(5)+q(9)
751 cosk = half * (e - one)
752 cosk = min(cosk,one)
753 cosk = max(cosk,-one)
754 ksi = acos(cosk)
755 sink = sin(ksi)
756 IF(sink==zero) THEN
757 si2 = zero
758 ELSE
759 si2 = half / sink
760 ENDIF
761C
762 t(1) = (q(6) - q(8)) * si2
763 t(2) = (q(7) - q(3)) * si2
764 t(3) = (q(2) - q(4)) * si2
765 nr = sqrt(t(1)*t(1)+t(2)*t(2)+t(3)*t(3))
766 IF (nr/=zero) nr = one/nr
767 t(1) = t(1)*nr
768 t(2) = t(2)*nr
769 t(3) = t(3)*nr
770
771C Vector of initial rotation (in local frame)
772 theta0(1) = t(1)*ksi
773 theta0(2) = t(2)*ksi
774 theta0(3) = t(3)*ksi
775C
776C-- Check of initial angles vs stopping angles
777 DO i=1,nb_rot
778 IF (theta0(i)<0) THEN
779 IF ((theta0(i)<stop_angl_min(i)).AND.(stop_angl_min(i)/=zero)) THEN
780 CALL ancmsg(msgid=1133,
781 . msgtype=msgerror,
782 . anmode=aninfo_blind_1,
783 . i1=id,
784 . c1=titr,
785 . r1=theta0(i),
786 . r2=stop_angl_min(i))
787 ENDIF
788 ELSE
789 IF ((theta0(i)>stop_angl_max(i)).AND.(stop_angl_max(i)/=zero)) THEN
790 CALL ancmsg(msgid=1133,
791 . msgtype=msgerror,
792 . anmode=aninfo_blind_1,
793 . i1=id,
794 . c1=titr,
795 . r1=theta0(i),
796 . r2=stop_angl_max(i))
797 ENDIF
798 ENDIF
799 ENDDO
800C
801 ENDIF
802C-----------------------------------------------
803 get_skew45 = ierr1
804
805 RETURN
806 END
807
808!||====================================================================
809!|| init_skew45 ../starter/source/elements/joint/rjoint/rini45.F
810!||--- called by ------------------------------------------------------
811!|| rini45 ../starter/source/elements/joint/rjoint/rini45.F
812!||--- calls -----------------------------------------------------
813!|| ancmsg ../starter/source/output/message/message.F
814!|| get_u_geo ../starter/source/user_interface/uaccess.F
815!|| get_u_skew ../starter/source/user_interface/uaccess.F
816!||--- uses -----------------------------------------------------
817!|| message_mod ../starter/share/message_module/message_mod.F
818!||====================================================================
819 SUBROUTINE init_skew45(JTYP,IPROP,IDSK1,IDSK2,VECT1,VECT2,ID,TITR)
820C-------------------------------------------------------------------------
821C This subroutine compute the initial angles from skews
822C-------------------------------------------------------------------------
823C-----------------------------------------------
824C M o d u l e s
825C-----------------------------------------------
826 USE message_mod
827 USE names_and_titles_mod , ONLY : nchartitle
828
829C-----------------------------------------------
830C I m p l i c i t T y p e s
831C-----------------------------------------------
832#include "implicit_f.inc"
833C-----------------------------------------------
834C A n a l y s e M o d u l e
835C-----------------------------------------------
836#include "param_c.inc"
837C----------------------------------------------------------
838C 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
839C----------------------------------------------------------
840 INTEGER JTYP,IPROP,ID,IDSK1,IDSK2
841 my_real, INTENT(INOUT) :: VECT1(LSKEW),VECT2(LSKEW)
842 CHARACTER(LEN=NCHARTITLE) :: TITR
843C-----------------------------------------------
844C L o c a l V a r i a b l e s
845C-----------------------------------------------
846 INTEGER I,J,GET_U_SKEW,N1,N2,N3,ISK
847 my_real
848 . Q(LSKEW),GET_U_GEO,NR,T(3),COSK,SINK,
849 . SI2,E,KSI,UMI,ERR,STOP_ANGL
850 DATA umi/-1.0/
851C-----------------------------------------------
852 EXTERNAL get_u_geo,get_u_skew
853C-----------------------------------------------
854
855 idsk1 = nint(get_u_geo(53,iprop))
856 idsk2 = nint(get_u_geo(54,iprop))
857
858C Check skew1 ------
859 IF (idsk1>0) THEN
860 isk = get_u_skew(idsk1,n1,n2,n3,vect1)
861 IF (isk==0) THEN
862 CALL ancmsg(msgid=926,
863 . msgtype=msgerror,
864 . anmode=aninfo_blind_1,
865 . i1=id,
866 . c1=titr,
867 . i2=idsk1)
868 ENDIF
869 ENDIF
870
871C Check skew2 ------
872 IF (idsk2>0) THEN
873 isk = get_u_skew(idsk2,n1,n2,n3,vect2)
874 IF (isk==0) THEN
875 CALL ancmsg(msgid=926,
876 . msgtype=msgerror,
877 . anmode=aninfo_blind_1,
878 . i1=id,
879 . c1=titr,
880 . i2=idsk2)
881 ENDIF
882 IF ((jtyp>4).AND.(jtyp<9)) THEN
883C-- init angle not allowed for joints with no rotations - idsk2 ignored
884 CALL ancmsg(msgid=1134,
885 . msgtype=msgwarning,
886 . anmode=aninfo_blind_1,
887 . i1=id,
888 . c1=titr,
889 . i2=jtyp,
890 . i3=idsk2)
891 idsk2 = 0
892 ENDIF
893 ENDIF
894C
895 RETURN
896 END
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
id
initmumps id
Definition multiplerhs_example.m:4
message_mod
Definition message_mod.F:1257
names_and_titles_mod
Definition names_and_titles_mod.F:1002
names_and_titles_mod::nchartitle
integer, parameter nchartitle
Definition names_and_titles_mod.F:1008
prod_atb
subroutine prod_atb(a, b, x)
Definition rini33.F:553
get_skew45
integer function get_skew45(iout, jtyp, ex, ixr, ixr_kj, iel, x, id, titr, idsk1, idsk2, vect1, vect2, theta0, stop_angl_min, stop_angl_max)
Definition rini45.F:339
init_skew45
subroutine init_skew45(jtyp, iprop, idsk1, idsk2, vect1, vect2, id, titr)
Definition rini45.F:820
rini45
subroutine rini45(nel, iout, iprop, ix, x, xl, mass, xiner, stifn, stifr, viscm, viscr, uvar, nuvar, ixr, ixr_kj, id, titr)
Definition rini45.F:41
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:895
get_u_pnu
integer function get_u_pnu(ivar, ip, k)
Definition uaccess.F:481