OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
rgbodfp.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!|| rgbodfp ../engine/source/constraints/general/rbody/rgbodfp.F
25!||--- called by ------------------------------------------------------
26!|| rbyfor ../engine/source/constraints/general/rbody/rbyfor.F
27!||--- calls -----------------------------------------------------
28!|| rotbmr ../engine/source/tools/skew/rotbmr.F
29!|| sum_6_float ../engine/source/system/parit.F
30!||--- uses -----------------------------------------------------
31!|| my_alloc_mod ../common_source/tools/memory/my_alloc.F90
32!||====================================================================
33 SUBROUTINE rgbodfp(AF ,AM ,X ,FS ,RBY ,
34 2 NOD ,M ,IN ,VR ,STIFN,
35 3 STIFR,FOPTA ,WEIGHT ,MS ,V ,
36 4 IFLAG,ICODR ,ISKEW ,SKEW ,RBF6 ,
37 5 NSN ,NATIV_SMS,FBSAV6,IPARSENS,
38 6 NODREAC,FTHREAC,CPTREAC,ISPHER)
39C-----------------------------------------------
40C M o d u l e s
41C-----------------------------------------------
42 USE my_alloc_mod
43C-----------------------------------------------
44C I m p l i c i t T y p e s
45C-----------------------------------------------
46#include "implicit_f.inc"
47C-----------------------------------------------
48C C o m m o n B l o c k s
49C-----------------------------------------------
50#include "com01_c.inc"
51#include "com08_c.inc"
52#include "impl1_c.inc"
53#include "param_c.inc"
54#include "parit_c.inc"
55#include "sms_c.inc"
56C-----------------------------------------------
57C D u m m y A r g u m e n t s
58C-----------------------------------------------
59 INTEGER , INTENT(IN ) :: ISPHER
60 INTEGER NOD(*), WEIGHT(*), ICODR(*), ISKEW(*),
61 . IFLAG,NSN,M, NATIV_SMS(*),IPARSENS,
62 . NODREAC(*),CPTREAC
63 my_real
64 . af(3,*), am(3,*), x(3,*), fs(*), rby(*), vr(3,*), skew(lskew,*),
65 . stifn(*),stifr(*),fopta(6),in(*),ms(*), v(3,*),
66 . fthreac(6,*)
67 DOUBLE PRECISION RBF6(8,6)
68 DOUBLE PRECISION
69 . FBSAV6(12,6)
70C-----------------------------------------------
71C L o c a l V a r i a b l e s
72C-----------------------------------------------
73 INTEGER I, N, LCOD, ISK, K
74C REAL
75 my_real
76 . WA1, WA2, WA3, DD, VI(3),II1,II2,II3,II4,II5,II6,II7,II8,II9,
77 . ENROTT,ENCINT,XMASST,XMOMTT,YMOMTT,ZMOMTT,MAS,AFM1,AFM2,AFM3,
78 . arm1,arm2,arm3,stfm,stfrm,
79 . det, il1,il2,il3,il4,il5,il6,il7,il8,il9
80 my_real,DIMENSION(:), ALLOCATABLE :: f1
81 my_real,DIMENSION(:), ALLOCATABLE :: f2
82 my_real,DIMENSION(:), ALLOCATABLE :: f3
83 my_real,DIMENSION(:), ALLOCATABLE :: f4
84 my_real,DIMENSION(:), ALLOCATABLE :: f5
85 my_real,DIMENSION(:), ALLOCATABLE :: f6
86 my_real,DIMENSION(:), ALLOCATABLE :: f7
87 my_real,DIMENSION(:), ALLOCATABLE :: f8
88C-----------------------------------------------
89 CALL my_alloc(f1,nsn)
90 CALL my_alloc(f2,nsn)
91 CALL my_alloc(f3,nsn)
92 CALL my_alloc(f4,nsn)
93 CALL my_alloc(f5,nsn)
94 CALL my_alloc(f6,nsn)
95 CALL my_alloc(f7,nsn)
96 CALL my_alloc(f8,nsn)
97C-----------------------------------------------
98C
99c M =NBY(1)
100C rajout optimisation pour spmd
101 IF (m<0) RETURN
102C premiere passe : calcul de force sur tous les r.b.
103 IF (iflag==1) THEN
104
105C CORRECTION DE L'INERTIE DU RIGID BODY POUR DT NODAL
106 rby(10) = max(rby(10),in(m))
107 rby(11) = max(rby(11),in(m))
108 rby(12) = max(rby(12),in(m))
109C
110c NSN =NBY(2)
111C
112 IF(idtmins/=2)THEN
113 DO i=1,nsn
114 n = nod(i)
115 IF(weight(n) == 1)THEN
116 dd = (x(1,n)-x(1,m))**2+(x(2,n)-x(2,m))**2+(x(3,n)-x(3,m))**2
117 f1(i) = stifn(n)
118 f2(i) = stifr(n)+dd*stifn(n)
119 f3(i) = af(1,n)
120 f4(i) = af(2,n)
121 f5(i) = af(3,n)
122 f6(i) = am(1,n)
123 . +(x(2,n)-x(2,m))*af(3,n)-(x(3,n)-x(3,m))*af(2,n)
124 f7(i) = am(2,n)
125 . +(x(3,n)-x(3,m))*af(1,n)-(x(1,n)-x(1,m))*af(3,n)
126 f8(i) = am(3,n)
127 . +(x(1,n)-x(1,m))*af(2,n)-(x(2,n)-x(2,m))*af(1,n)
128 ELSE
129 f1(i) = zero
130 f2(i) = zero
131 f3(i) = zero
132 f4(i) = zero
133 f5(i) = zero
134 f6(i) = zero
135 f7(i) = zero
136 f8(i) = zero
137 END IF
138 ENDDO
139 ELSE ! IF(IDTMINS/=2)THEN
140 DO i=1,nsn
141 n = nod(i)
142 IF(weight(n) == 1)THEN
143 dd = (x(1,n)-x(1,m))**2+(x(2,n)-x(2,m))**2+(x(3,n)-x(3,m))**2
144 f1(i) = stifn(n)
145 f2(i) = stifr(n)+dd*stifn(n)
146C
147 f3(i) = af(1,n)
148 f4(i) = af(2,n)
149 f5(i) = af(3,n)
150 f6(i) = am(1,n)
151 . +(x(2,n)-x(2,m))*af(3,n)-(x(3,n)-x(3,m))*af(2,n)
152 f7(i) = am(2,n)
153 . +(x(3,n)-x(3,m))*af(1,n)-(x(1,n)-x(1,m))*af(3,n)
154 f8(i) = am(3,n)
155 . +(x(1,n)-x(1,m))*af(2,n)-(x(2,n)-x(2,m))*af(1,n)
156 ELSE
157 f1(i) = zero
158 f2(i) = zero
159 f3(i) = zero
160 f4(i) = zero
161 f5(i) = zero
162 f6(i) = zero
163 f7(i) = zero
164 f8(i) = zero
165 END IF
166 ENDDO
167 END IF
168C
169C Traitement Parith/ON avant echange
170C
171 DO k = 1, 6
172 rbf6(1,k) = zero
173 rbf6(2,k) = zero
174 rbf6(3,k) = zero
175 rbf6(4,k) = zero
176 rbf6(5,k) = zero
177 rbf6(6,k) = zero
178 rbf6(7,k) = zero
179 rbf6(8,k) = zero
180 END DO
181 IF(iparit > 0)THEN
182 CALL sum_6_float(1 ,nsn ,f1, rbf6(1,1),8)
183 CALL sum_6_float(1 ,nsn ,f2, rbf6(2,1),8)
184 CALL sum_6_float(1 ,nsn ,f3, rbf6(3,1),8)
185 CALL sum_6_float(1 ,nsn ,f4, rbf6(4,1),8)
186 CALL sum_6_float(1 ,nsn ,f5, rbf6(5,1),8)
187 CALL sum_6_float(1 ,nsn ,f6, rbf6(6,1),8)
188 CALL sum_6_float(1 ,nsn ,f7, rbf6(7,1),8)
189 CALL sum_6_float(1 ,nsn ,f8, rbf6(8,1),8)
190 ELSE
191 DO i=1,nsn
192 k = 1
193 rbf6(1,k) = rbf6(1,k) + f1(i)
194 rbf6(2,k) = rbf6(2,k) + f2(i)
195 rbf6(3,k) = rbf6(3,k) + f3(i)
196 rbf6(4,k) = rbf6(4,k) + f4(i)
197 rbf6(5,k) = rbf6(5,k) + f5(i)
198 rbf6(6,k) = rbf6(6,k) + f6(i)
199 rbf6(7,k) = rbf6(7,k) + f7(i)
200 rbf6(8,k) = rbf6(8,k) + f8(i)
201 END DO
202 END IF
203c
204C
205C phase 2 : sauvegarde anim + rby
206C
207 ELSEIF (iflag==2) THEN
208c NSN =NBY(2)
209C
210C Traitement Parith/ON apres echange
211C
212 stifn(m) = stifn(m)+
213 + rbf6(1,1)+rbf6(1,2)+rbf6(1,3)+
214 + rbf6(1,4)+rbf6(1,5)+rbf6(1,6)
215 stifr(m) = stifr(m)+
216 + rbf6(2,1)+rbf6(2,2)+rbf6(2,3)+
217 + rbf6(2,4)+rbf6(2,5)+rbf6(2,6)
218c
219 afm1 = rbf6(3,1)+rbf6(3,2)+rbf6(3,3)+
220 + rbf6(3,4)+rbf6(3,5)+rbf6(3,6)
221 afm2 = rbf6(4,1)+rbf6(4,2)+rbf6(4,3)+
222 + rbf6(4,4)+rbf6(4,5)+rbf6(4,6)
223 afm3 = rbf6(5,1)+rbf6(5,2)+rbf6(5,3)+
224 + rbf6(5,4)+rbf6(5,5)+rbf6(5,6)
225 arm1 = rbf6(6,1)+rbf6(6,2)+rbf6(6,3)+
226 + rbf6(6,4)+rbf6(6,5)+rbf6(6,6)
227 arm2 = rbf6(7,1)+rbf6(7,2)+rbf6(7,3)+
228 + rbf6(7,4)+rbf6(7,5)+rbf6(7,6)
229 arm3 = rbf6(8,1)+rbf6(8,2)+rbf6(8,3)+
230 + rbf6(8,4)+rbf6(8,5)+rbf6(8,6)
231c
232 af(1,m) = af(1,m)+
233 + rbf6(3,1)+rbf6(3,2)+rbf6(3,3)+
234 + rbf6(3,4)+rbf6(3,5)+rbf6(3,6)
235 af(2,m) = af(2,m)+
236 + rbf6(4,1)+rbf6(4,2)+rbf6(4,3)+
237 + rbf6(4,4)+rbf6(4,5)+rbf6(4,6)
238 af(3,m) = af(3,m)+
239 + rbf6(5,1)+rbf6(5,2)+rbf6(5,3)+
240 + rbf6(5,4)+rbf6(5,5)+rbf6(5,6)
241 am(1,m) = am(1,m)+
242 + rbf6(6,1)+rbf6(6,2)+rbf6(6,3)+
243 + rbf6(6,4)+rbf6(6,5)+rbf6(6,6)
244 am(2,m) = am(2,m)+
245 + rbf6(7,1)+rbf6(7,2)+rbf6(7,3)+
246 + rbf6(7,4)+rbf6(7,5)+rbf6(7,6)
247 am(3,m) = am(3,m)+
248 + rbf6(8,1)+rbf6(8,2)+rbf6(8,3)+
249 + rbf6(8,4)+rbf6(8,5)+rbf6(8,6)
250C
251 IF(imconv==1)THEN
252 !---
253 ! reaction forces/moments in rbody TH/RBODY
254 !---
255!! FS(1)=FS(1)+AF(1,M)*DT1*WEIGHT(M)
256!! FS(2)=FS(2)+AF(2,M)*DT1*WEIGHT(M)
257!! FS(3)=FS(3)+AF(3,M)*DT1*WEIGHT(M)
258!! FS(4)=FS(4)+AM(1,M)*DT1*WEIGHT(M)
259!! FS(5)=FS(5)+AM(2,M)*DT1*WEIGHT(M)
260!! FS(6)=FS(6)+AM(3,M)*DT1*WEIGHT(M)
261 fs(1)=fs(1)+afm1*dt1*weight(m)
262 fs(2)=fs(2)+afm2*dt1*weight(m)
263 fs(3)=fs(3)+afm3*dt1*weight(m)
264 fs(4)=fs(4)+arm1*dt1*weight(m)
265 fs(5)=fs(5)+arm2*dt1*weight(m)
266 fs(6)=fs(6)+arm3*dt1*weight(m)
267 !---
268 ! reaction forces/moments in "main node" of RBODY for TH/NODE
269 !---
270 IF (cptreac > 0) THEN
271 IF(nodreac(m) > 0) THEN
272 fthreac(4,nodreac(m)) = fthreac(4,nodreac(m)) - am(1,m)*dt1*weight(m)
273 fthreac(5,nodreac(m)) = fthreac(5,nodreac(m)) - am(2,m)*dt1*weight(m)
274 fthreac(6,nodreac(m)) = fthreac(6,nodreac(m)) - am(3,m)*dt1*weight(m)
275 ENDIF
276 ENDIF ! IF (CPTREAC . and. NODREAC(M) > 0)
277 ENDIF
278C
279C
280 !---
281 ! reaction forces/moments in RBODY /ANIM/VECT/FOPT
282 !---
283!! FOPTA(1)=AF(1,M)*WEIGHT(M)
284!! FOPTA(2)=AF(2,M)*WEIGHT(M)
285!! FOPTA(3)=AF(3,M)*WEIGHT(M)
286!! FOPTA(4)=AM(1,M)*WEIGHT(M)
287!! FOPTA(5)=AM(2,M)*WEIGHT(M)
288!! FOPTA(6)=AM(3,M)*WEIGHT(M)
289 fopta(1)=afm1*weight(m)
290 fopta(2)=afm2*weight(m)
291 fopta(3)=afm3*weight(m)
292 fopta(4)=arm1*weight(m)
293 fopta(5)=arm2*weight(m)
294 fopta(6)=arm3*weight(m)
295C
296 IF (iparsens /= 0)THEN
297 DO i=1,6
298 fbsav6(1,i) = rbf6(3,i)*weight(m)
299 fbsav6(2,i) = rbf6(4,i)*weight(m)
300 fbsav6(3,i) = rbf6(5,i)*weight(m)
301 fbsav6(4,i) = rbf6(6,i)*weight(m)
302 fbsav6(5,i) = rbf6(7,i)*weight(m)
303 fbsav6(6,i) = rbf6(8,i)*weight(m)
304 ENDDO
305 fbsav6(7,1) = af(1,m)*weight(m)
306 fbsav6(8,1) = af(2,m)*weight(m)
307 fbsav6(9,1) = af(3,m)*weight(m)
308 fbsav6(10,1) = am(1,m)*weight(m)
309 fbsav6(11,1) = am(2,m)*weight(m)
310 fbsav6(12,1) = am(3,m)*weight(m)
311 ENDIF
312c
313 isk =iskew(m)
314 lcod=icodr(m)
315 IF(lcod/=0.AND.imconv==1)THEN
316C
317C rotation de la matrice d'orientation (directions principales)
318 IF(n2d==0) THEN
319 vi(1)=rby(1)*vr(1,m)+rby(2)*vr(2,m)+rby(3)*vr(3,m)
320 vi(2)=rby(4)*vr(1,m)+rby(5)*vr(2,m)+rby(6)*vr(3,m)
321 vi(3)=rby(7)*vr(1,m)+rby(8)*vr(2,m)+rby(9)*vr(3,m)
322 ELSEIF(n2d==1) THEN
323 vi(1)=zero
324 vi(2)=zero
325 vi(3)=rby(9)*vr(3,m)
326 ELSE
327 vi(1)=rby(1)*vr(1,m)
328 vi(2)=zero
329 vi(3)=zero
330 ENDIF
331
332 CALL rotbmr(vi,rby,dt1)
333C
334C matrice d'inertie en repere global
335 IF(n2d==0) THEN
336C ANALYSE 3D
337 ii1=rby(10)*rby(1)
338 ii2=rby(10)*rby(2)
339 ii3=rby(10)*rby(3)
340 ii4=rby(11)*rby(4)
341 ii5=rby(11)*rby(5)
342 ii6=rby(11)*rby(6)
343 ii7=rby(12)*rby(7)
344 ii8=rby(12)*rby(8)
345 ii9=rby(12)*rby(9)
346C
347 rby(17)=rby(1)*ii1 + rby(4)*ii4 + rby(7)*ii7
348 rby(18)=rby(1)*ii2 + rby(4)*ii5 + rby(7)*ii8
349 rby(19)=rby(1)*ii3 + rby(4)*ii6 + rby(7)*ii9
350 rby(20)=rby(2)*ii1 + rby(5)*ii4 + rby(8)*ii7
351 rby(21)=rby(2)*ii2 + rby(5)*ii5 + rby(8)*ii8
352 rby(22)=rby(2)*ii3 + rby(5)*ii6 + rby(8)*ii9
353 rby(23)=rby(3)*ii1 + rby(6)*ii4 + rby(9)*ii7
354 rby(24)=rby(3)*ii2 + rby(6)*ii5 + rby(9)*ii8
355 rby(25)=rby(3)*ii3 + rby(6)*ii6 + rby(9)*ii9
356C
357C ajout des termes [Iglobal] vr ^ vr
358 wa1=rby(17)*vr(1,m)+rby(18)*vr(2,m)+rby(19)*vr(3,m)
359 wa2=rby(20)*vr(1,m)+rby(21)*vr(2,m)+rby(22)*vr(3,m)
360 wa3=rby(23)*vr(1,m)+rby(24)*vr(2,m)+rby(25)*vr(3,m)
361C
362 ELSEIF(n2d==1) THEN
363C ANALYSE 2D : Axisymmetry
364C I= A 0 0
365C 0 A 0
366C 0 0 B
367C
368 rby(17)=rby(10)
369 rby(21)=rby(11)
370 rby(25)=rby(12)
371C ajout des termes [Iglobal] vr ^ vr
372 wa1=zero
373 wa2=zero
374 wa3=rby(12)*vr(3,m)
375C
376 ELSE
377C ANALYSE 2D : Plane strain Inertia matrix
378C I= A 0 0
379C 0 B D
380C 0 D C
381 ii1=rby(10)*rby(1)
382 ii5=rby(11)*rby(5)
383 ii6=rby(11)*rby(6)
384 ii8=rby(12)*rby(8)
385 ii9=rby(12)*rby(9)
386C
387 rby(17)=rby(1)*ii1
388 rby(21)=rby(5)*ii5 + rby(8)*ii8
389 rby(22)=rby(5)*ii6 + rby(8)*ii9
390 rby(23)=zero
391 rby(24)=rby(6)*ii5 + rby(9)*ii8
392 rby(25)=rby(6)*ii6 + rby(9)*ii9
393C
394C ajout des termes [Iglobal] vr ^ vr
395 wa1=rby(17)*vr(1,m)
396 wa2=zero
397 wa3=zero
398 ENDIF
399C
400 IF (impl_s==0) THEN
401 am(1,m)=am(1,m)+(wa2*vr(3,m)-wa3*vr(2,m))
402 am(2,m)=am(2,m)+(wa3*vr(1,m)-wa1*vr(3,m))
403 am(3,m)=am(3,m)+(wa1*vr(2,m)-wa2*vr(1,m))
404 END IF
405
406C
407 IF(isk==1)THEN
408C------------------
409C REPERE GLOBAL :
410C Resolution [Iglobal] gama = M, compte-tenu des conditions aux limites
411C Ex : gamaz=0
412C | Ixx Ixy Ixz | { gamax } { Mx }
413C | Iyx Iyy Iyz | { gamay } = { My }
414C | Izx Izy Izz | { 0 } { Mz + DMz } DMz inconnue
415C equivaut a
416C | Ixx Ixy | { gamax } { Mx }
417C | Iyx Iyy | { gamay } = { My }
418C et gamaz=0
419C------------------
420 IF(lcod==1)THEN
421 det=one/(rby(17)*rby(21)-rby(18)*rby(20))
422 wa1=am(1,m)
423 wa2=am(2,m)
424 am(1,m)=( rby(21)*wa1-rby(20)*wa2)*det
425 am(2,m)=(-rby(18)*wa1+rby(17)*wa2)*det
426 am(3,m)=zero
427 ELSEIF(lcod==2)THEN
428 det=one/(rby(17)*rby(25)-rby(19)*rby(23))
429 wa1=am(1,m)
430 wa2=am(3,m)
431 am(1,m)=( rby(25)*wa1-rby(23)*wa2)*det
432 am(2,m)=zero
433 am(3,m)=(-rby(19)*wa1+rby(17)*wa2)*det
434 ELSEIF(lcod==3)THEN
435 am(1,m)=am(1,m)/rby(17)
436 am(2,m)=zero
437 am(3,m)=zero
438 ELSEIF(lcod==4)THEN
439 det=one/(rby(21)*rby(25)-rby(22)*rby(24))
440 wa1=am(2,m)
441 wa2=am(3,m)
442 am(1,m)=zero
443 am(2,m)=( rby(25)*wa1-rby(24)*wa2)*det
444 am(3,m)=(-rby(22)*wa1+rby(21)*wa2)*det
445 ELSEIF(lcod==5)THEN
446 am(1,m)=zero
447 am(2,m)=am(2,m)/rby(21)
448 am(3,m)=zero
449 ELSEIF(lcod==6)THEN
450 am(1,m)=zero
451 am(2,m)=zero
452 am(3,m)=am(3,m)/rby(25)
453 ELSEIF(lcod==7)THEN
454 am(1,m)=zero
455 am(2,m)=zero
456 am(3,m)=zero
457 ENDIF
458 ELSE
459C-------------------
460C REPERE OBLIQUE
461C------------------
462C
463C Passage dans le skew : (vitesse), moments, matrice d'inertie.
464C
465C WA1=VR(1,M)
466C WA2=VR(2,M)
467C WA3=VR(3,M)
468C VL(1,M)=SKEW(1,ISK)*WA1+SKEW(4,ISK)*WA2+SKEW(7,ISK)*WA3
469C VL(2,M)=SKEW(2,ISK)*WA1+SKEW(5,ISK)*WA2+SKEW(8,ISK)*WA3
470C VL(3,M)=SKEW(3,ISK)*WA1+SKEW(6,ISK)*WA2+SKEW(9,ISK)*WA3
471C
472 wa1=am(1,m)
473 wa2=am(2,m)
474 wa3=am(3,m)
475 am(1,m)=skew(1,isk)*wa1+skew(2,isk)*wa2+skew(3,isk)*wa3
476 am(2,m)=skew(4,isk)*wa1+skew(5,isk)*wa2+skew(6,isk)*wa3
477 am(3,m)=skew(7,isk)*wa1+skew(8,isk)*wa2+skew(9,isk)*wa3
478C
479C Resolution ds le repere local, compte-tenu des conditions aux limites
480C Ex : v3+gama3*dt12=0
481C | IL1 IL2 IL3 | { gama1 } { M1 }
482C | IL4 IL5 IL6 | { gama2 } = { M2 }
483C | IL7 IL8 IL9 | { -v3/dt12 } { M3 + DM3 } DM3 inconnue
484C equivaut a
485C | IL1 IL2 IL3 | { gama1 } { M1 } | IL1 IL2 IL3 | { 0 }
486C | IL4 IL5 IL6 | { gama2 } = { M2 } + | IL4 IL5 IL6 | { 0 }
487C | IL7 IL8 IL9 | { 0 } { M3 + DM3 } | IL7 IL8 IL9 | { v3/dt12 }
488C
489C pas de solution => gama3=0, v'3=0 (reimpose dans la condition limite)
490C
491C | IL1 IL2 IL3 | { gama1 } { M1 } | IL1 IL2 IL3 | { gama1 } { M1 }
492C | IL4 IL5 IL6 | { gama2 } = { M2 } <==> | IL4 IL5 IL6 | { gama2 } = { M2 }
493C | IL7 IL8 IL9 | { 0 } { M3 + DM3 }
494C
495 IF(lcod==1)THEN
496C12345678901234567890123456789012345678901234567890123456789012345678901
497 ii1=rby(17)*skew(1,isk)+rby(18)*skew(2,isk)+rby(19)*skew(3,isk)
498 ii2=rby(17)*skew(4,isk)+rby(18)*skew(5,isk)+rby(19)*skew(6,isk)
499 ii4=rby(20)*skew(1,isk)+rby(21)*skew(2,isk)+rby(22)*skew(3,isk)
500 ii5=rby(20)*skew(4,isk)+rby(21)*skew(5,isk)+rby(22)*skew(6,isk)
501 ii7=rby(23)*skew(1,isk)+rby(24)*skew(2,isk)+rby(25)*skew(3,isk)
502 ii8=rby(23)*skew(4,isk)+rby(24)*skew(5,isk)+rby(25)*skew(6,isk)
503 il1=skew(1,isk)*ii1+skew(2,isk)*ii4+skew(3,isk)*ii7
504 il2=skew(1,isk)*ii2+skew(2,isk)*ii5+skew(3,isk)*ii8
505 il4=skew(4,isk)*ii1+skew(5,isk)*ii4+skew(6,isk)*ii7
506 il5=skew(4,isk)*ii2+skew(5,isk)*ii5+skew(6,isk)*ii8
507C
508 det=one/(il1*il5-il2*il4)
509 wa1=am(1,m)
510 wa2=am(2,m)
511 am(1,m)=( il5*wa1-il4*wa2)*det
512 am(2,m)=(-il2*wa1+il1*wa2)*det
513 am(3,m)=zero
514C VL(3)=ZERO
515 ELSEIF(lcod==2)THEN
516 ii1=rby(17)*skew(1,isk)+rby(18)*skew(2,isk)+rby(19)*skew(3,isk)
517 ii3=rby(17)*skew(7,isk)+rby(18)*skew(8,isk)+rby(19)*skew(9,isk)
518 ii4=rby(20)*skew(1,isk)+rby(21)*skew(2,isk)+rby(22)*skew(3,isk)
519 ii6=rby(20)*skew(7,isk)+rby(21)*skew(8,isk)+rby(22)*skew(9,isk)
520 ii7=rby(23)*skew(1,isk)+rby(24)*skew(2,isk)+rby(25)*skew(3,isk)
521 ii9=rby(23)*skew(7,isk)+rby(24)*skew(8,isk)+rby(25)*skew(9,isk)
522 il1=skew(1,isk)*ii1+skew(2,isk)*ii4+skew(3,isk)*ii7
523 il3=skew(1,isk)*ii3+skew(2,isk)*ii6+skew(3,isk)*ii9
524 il7=skew(7,isk)*ii1+skew(8,isk)*ii4+skew(9,isk)*ii7
525 il9=skew(7,isk)*ii3+skew(8,isk)*ii6+skew(9,isk)*ii9
526C
527 det=one/(il1*il9-il3*il7)
528 wa1=am(1,m)
529 wa2=am(3,m)
530 am(1,m)=( il9*wa1-il7*wa2)*det
531 am(2,m)=zero
532 am(3,m)=(-il3*wa1+il1*wa2)*det
533C VL(2)=ZERO
534 ELSEIF(lcod==3)THEN
535 ii1=rby(17)*skew(1,isk)+rby(18)*skew(2,isk)+rby(19)*skew(3,isk)
536 ii4=rby(20)*skew(1,isk)+rby(21)*skew(2,isk)+rby(22)*skew(3,isk)
537 ii7=rby(23)*skew(1,isk)+rby(24)*skew(2,isk)+rby(25)*skew(3,isk)
538 il1=skew(1,isk)*ii1+skew(2,isk)*ii4+skew(3,isk)*ii7
539C
540 am(1,m)=am(1,m)/il1
541 am(2,m)=zero
542 am(3,m)=zero
543C VL(2)=ZERO
544C VL(3)=ZERO
545 ELSEIF(lcod==4)THEN
546 ii2=rby(17)*skew(4,isk)+rby(18)*skew(5,isk)+rby(19)*skew(6,isk)
547 ii3=rby(17)*skew(7,isk)+rby(18)*skew(8,isk)+rby(19)*skew(9,isk)
548 ii5=rby(20)*skew(4,isk)+rby(21)*skew(5,isk)+rby(22)*skew(6,isk)
549 ii6=rby(20)*skew(7,isk)+rby(21)*skew(8,isk)+rby(22)*skew(9,isk)
550 ii8=rby(23)*skew(4,isk)+rby(24)*skew(5,isk)+rby(25)*skew(6,isk)
551 ii9=rby(23)*skew(7,isk)+rby(24)*skew(8,isk)+rby(25)*skew(9,isk)
552 il5=skew(4,isk)*ii2+skew(5,isk)*ii5+skew(6,isk)*ii8
553 il6=skew(4,isk)*ii3+skew(5,isk)*ii6+skew(6,isk)*ii9
554 il8=skew(7,isk)*ii2+skew(8,isk)*ii5+skew(9,isk)*ii8
555 il9=skew(7,isk)*ii3+skew(8,isk)*ii6+skew(9,isk)*ii9
556C
557 det=one/(il5*il9-il6*il8)
558 wa1=am(2,m)
559 wa2=am(3,m)
560 am(1,m)=zero
561 am(2,m)=( il9*wa1-il8*wa2)*det
562 am(3,m)=(-il6*wa1+il5*wa2)*det
563C VL(1)=ZERO
564 ELSEIF(lcod==5)THEN
565 ii2=rby(17)*skew(4,isk)+rby(18)*skew(5,isk)+rby(19)*skew(6,isk)
566 ii5=rby(20)*skew(4,isk)+rby(21)*skew(5,isk)+rby(22)*skew(6,isk)
567 ii8=rby(23)*skew(4,isk)+rby(24)*skew(5,isk)+rby(25)*skew(6,isk)
568 il5=skew(4,isk)*ii2+skew(5,isk)*ii5+skew(6,isk)*ii8
569C
570 am(1,m)=zero
571 am(2,m)=am(2,m)/il5
572 am(3,m)=zero
573C VL(1)=ZERO
574C VL(3)=ZERO
575 ELSEIF(lcod==6)THEN
576 ii3=rby(17)*skew(7,isk)+rby(18)*skew(8,isk)+rby(19)*skew(9,isk)
577 ii6=rby(20)*skew(7,isk)+rby(21)*skew(8,isk)+rby(22)*skew(9,isk)
578 ii9=rby(23)*skew(7,isk)+rby(24)*skew(8,isk)+rby(25)*skew(9,isk)
579 il9=skew(7,isk)*ii3+skew(8,isk)*ii6+skew(9,isk)*ii9
580C
581 am(1,m)=zero
582 am(2,m)=zero
583 am(3,m)=am(3,m)/il9
584C VL(1)=ZERO
585C VL(2)=ZERO
586 ELSEIF(lcod==7)THEN
587 am(1,m)=zero
588 am(2,m)=zero
589 am(3,m)=zero
590C VL(1)=ZERO
591C VL(2)=ZERO
592C VL(3)=ZERO
593 ENDIF
594 wa1=am(1,m)
595 wa2=am(2,m)
596 wa3=am(3,m)
597 am(1,m)=skew(1,isk)*wa1+skew(4,isk)*wa2+skew(7,isk)*wa3
598 am(2,m)=skew(2,isk)*wa1+skew(5,isk)*wa2+skew(8,isk)*wa3
599 am(3,m)=skew(3,isk)*wa1+skew(6,isk)*wa2+skew(9,isk)*wa3
600 ENDIF
601C
602C CALCUL D'ONE PSEUDO MOMENT:
603C EN FAIT L'ACCELERATION DE ROTATION * INERTIE DU NOEUD MAIN (IMIN DU RB)
604 am(1,m)=in(m)*am(1,m)
605 am(2,m)=in(m)*am(2,m)
606 am(3,m)=in(m)*am(3,m)
607C
608 ELSEIF (imconv==1) THEN
609 IF(n2d==0) THEN
610 IF(ispher==1) THEN
611! do it just if this principal base is used somewhere else
612 vi(1)=rby(1)*vr(1,m)+rby(2)*vr(2,m)+rby(3)*vr(3,m)
613 vi(2)=rby(4)*vr(1,m)+rby(5)*vr(2,m)+rby(6)*vr(3,m)
614 vi(3)=rby(7)*vr(1,m)+rby(8)*vr(2,m)+rby(9)*vr(3,m)
615 CALL rotbmr(vi,rby,dt1)
616 ELSE
617 wa1=am(1,m)
618 wa2=am(2,m)
619 wa3=am(3,m)
620C les contributions des vr ne sont ajoutees que sur le processeur main
621 vi(1)=rby(1)*vr(1,m)+rby(2)*vr(2,m)+rby(3)*vr(3,m)
622 vi(2)=rby(4)*vr(1,m)+rby(5)*vr(2,m)+rby(6)*vr(3,m)
623 vi(3)=rby(7)*vr(1,m)+rby(8)*vr(2,m)+rby(9)*vr(3,m)
624C Update of skew of principal inerties - VI is the rotation vector so it is unchanged
625 CALL rotbmr(vi,rby,dt1)
626C repere globale -> repere d'inertie principale
627 am(1,m)=rby(1)*wa1+rby(2)*wa2+rby(3)*wa3
628 am(2,m)=rby(4)*wa1+rby(5)*wa2+rby(6)*wa3
629 am(3,m)=rby(7)*wa1+rby(8)*wa2+rby(9)*wa3
630 IF (impl_s==0) THEN
631 am(1,m) = am(1,m) + (rby(11)-rby(12))*vi(2)*vi(3)
632 am(2,m) = am(2,m) + (rby(12)-rby(10))*vi(3)*vi(1)
633 am(3,m) = am(3,m) + (rby(10)-rby(11))*vi(1)*vi(2)
634 END IF
635C
636C CALCUL D'ONE PSEUDO MOMENT:
637C EN FAIT L'ACCELERATION DE ROTATION * INERTIE DU NOEUD MAIN (IMIN DU RB)
638C
639 IF (impl_s==0) THEN
640 wa1 = am(1,m)*in(m)/rby(10)
641 wa2 = am(2,m)*in(m)/rby(11)
642 wa3 = am(3,m)*in(m)/rby(12)
643 ELSE
644 wa1 = am(1,m)
645 wa2 = am(2,m)
646 wa3 = am(3,m)
647 END IF
648C repere d'inertie principale -> repere globale
649 am(1,m)=rby(1)*wa1+rby(4)*wa2+rby(7)*wa3
650 am(2,m)=rby(2)*wa1+rby(5)*wa2+rby(8)*wa3
651 am(3,m)=rby(3)*wa1+rby(6)*wa2+rby(9)*wa3
652C MATRICE d'inertie -> repere globale
653 ii1=rby(10)*rby(1)
654 ii2=rby(10)*rby(2)
655 ii3=rby(10)*rby(3)
656 ii4=rby(11)*rby(4)
657 ii5=rby(11)*rby(5)
658 ii6=rby(11)*rby(6)
659 ii7=rby(12)*rby(7)
660 ii8=rby(12)*rby(8)
661 ii9=rby(12)*rby(9)
662C
663 rby(17)=rby(1)*ii1 + rby(4)*ii4 + rby(7)*ii7
664 rby(18)=rby(1)*ii2 + rby(4)*ii5 + rby(7)*ii8
665 rby(19)=rby(1)*ii3 + rby(4)*ii6 + rby(7)*ii9
666 rby(20)=rby(2)*ii1 + rby(5)*ii4 + rby(8)*ii7
667 rby(21)=rby(2)*ii2 + rby(5)*ii5 + rby(8)*ii8
668 rby(22)=rby(2)*ii3 + rby(5)*ii6 + rby(8)*ii9
669 rby(23)=rby(3)*ii1 + rby(6)*ii4 + rby(9)*ii7
670 rby(24)=rby(3)*ii2 + rby(6)*ii5 + rby(9)*ii8
671 rby(25)=rby(3)*ii3 + rby(6)*ii6 + rby(9)*ii9
672 END IF !(ISPHER==1) THEN
673C
674 ELSEIF(n2d==1) THEN
675C
676C CALCUL D'ONE PSEUDO MOMENT:
677C EN FAIT L'ACCELERATION DE ROTATION * INERTIE DU NOEUD MAIN (IMIN DU RB)
678C
679 IF (impl_s==0) THEN
680 wa1 = zero
681 wa2 = zero
682 wa3 = am(3,m)*in(m)/rby(12)
683 ELSE
684 wa1 = zero
685 wa2 = zero
686 wa3 = am(3,m)
687 END IF
688C MATRICE d'inertie -> repere globale
689C
690 rby(17)=rby(10)
691 rby(21)=rby(11)
692 rby(25)=rby(12)
693C
694 ELSEIF(n2d==2) THEN
695 wa1=am(1,m)
696 wa2=zero
697 wa3=zero
698C repere globale -> repere d'inertie principale
699 am(1,m)=rby(1)*wa1
700 am(2,m)=zero
701 am(3,m)=zero
702C les contributions des vr ne sont ajoutees que sur le processeur main
703 vi(1)=rby(1)*vr(1,m)
704 vi(2)=zero
705 vi(3)=zero
706 CALL rotbmr(vi,rby,dt1)
707C
708C CALCUL D'ONE PSEUDO MOMENT:
709C EN FAIT L'ACCELERATION DE ROTATION * INERTIE DU NOEUD MAIN (IMIN DU RB)
710C
711 IF (impl_s==0) THEN
712 wa1 = am(1,m)*in(m)/rby(10)
713 wa2 = zero
714 wa3 = zero
715 ELSE
716 wa1 = am(1,m)
717 wa2 = am(2,m)
718 wa3 = am(3,m)
719 END IF
720C repere d'inertie principale -> repere globale
721 am(1,m)=rby(1)*wa1
722 am(2,m)=rby(2)*wa1
723 am(3,m)=rby(3)*wa1
724C MATRICE d'inertie -> repere globale
725 ii1=rby(10)*rby(1)
726 ii5=rby(11)*rby(5)
727 ii6=rby(11)*rby(6)
728 ii8=rby(12)*rby(8)
729 ii9=rby(12)*rby(9)
730C
731 rby(17)=rby(1)*ii1
732 rby(21)=rby(5)*ii5 + rby(8)*ii8
733 rby(22)=rby(5)*ii6 + rby(8)*ii9
734 rby(23)=zero
735 rby(24)=rby(6)*ii5 + rby(9)*ii8
736 rby(25)=rby(6)*ii6 + rby(9)*ii9
737C
738 ENDIF
739 ENDIF
740
741 DO i=1,nsn
742 n = nod(i)
743C
744C AF, AM does not need to be reset and it will still be used for computing reaction forces...
745C AF(1,N)= ZERO
746C AF(2,N)= ZERO
747C AF(3,N)= ZERO
748C
749C AM(1,N)= ZERO
750C AM(2,N)= ZERO
751C AM(3,N)= ZERO
752C
753 stifr(n)= em20
754 stifn(n)= em20
755 ENDDO
756 ENDIF
757C
758 DEALLOCATE(f1)
759 DEALLOCATE(f2)
760 DEALLOCATE(f3)
761 DEALLOCATE(f4)
762 DEALLOCATE(f5)
763 DEALLOCATE(f6)
764 DEALLOCATE(f7)
765 DEALLOCATE(f8)
766
767 RETURN
768 END
769
770C--------------------------------------------------------------------------------
max
#define max(a, b)
Definition macros.h:21
sum_6_float
subroutine sum_6_float(jft, jlt, f, f6, n)
Definition parit.F:64
rgbodfp
subroutine rgbodfp(af, am, x, fs, rby, nod, m, in, vr, stifn, stifr, fopta, weight, ms, v, iflag, icodr, iskew, skew, rbf6, nsn, nativ_sms, fbsav6, iparsens, nodreac, fthreac, cptreac, ispher)
Definition rgbodfp.F:39
rotbmr
subroutine rotbmr(vr, rby, dt)
Definition rotbmr.F:35