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