OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
rmatpon.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!|| rmatpon ../engine/source/materials/mat/mat013/rmatpon.F
25!||--- called by ------------------------------------------------------
26!|| rmatforp ../engine/source/materials/mat/mat013/rmatforp.F
27!||--- calls -----------------------------------------------------
28!|| rotbmr ../engine/source/tools/skew/rotbmr.F
29!|| sum_6_float ../engine/source/system/parit.F
30!||====================================================================
31 SUBROUTINE rmatpon(AF ,AM ,X ,RBY ,NOD ,
32 2 NBY ,STIFN,STIFR,WEIGHT,NSL ,
33 3 RBF6 ,ICOD ,ARBY ,VRBY ,ARRBY,
34 4 VRRBY,IFLAG )
35C-----------------------------------------------
36C I m p l i c i t T y p e s
37C-----------------------------------------------
38#include "implicit_f.inc"
39C-----------------------------------------------
40C C o m m o n B l o c k s
41C-----------------------------------------------
42#include "com08_c.inc"
43C-----------------------------------------------
44C D u m m y A r g u m e n t s
45C-----------------------------------------------
46 INTEGER NOD(*), NBY(*), WEIGHT(*), ICOD(2,*),
47 . NSL, IFLAG
48C REAL
49 my_real
50 . af(3,*), am(3,*), x(3,*), rby(*),
51 . stifn(*),stifr(*),
52 . vrrby(3,*),arby(3,*),arrby(3,*),vrby(3,*)
53 DOUBLE PRECISION RBF6(6,6)
54C-----------------------------------------------
55C L o c a l V a r i a b l e s
56C-----------------------------------------------
57 INTEGER M, I, N, LCOD, ISK, K
58C REAL
59 my_real
60 . vi(3),xg,yg,zg,ii1,ii2,ii3,ii4,ii5,ii6,ii7,ii8,
61 . ii9,wa1,wa2,wa3,det,in,msrby,rx,ry,rz,rb(12),
62 . f1(nsl), f2(nsl), f3(nsl), f4(nsl),
63 . f5(nsl), f6(nsl)
64C-----------------------------------------------
65C premiere passe : calcul de force sur tous les r.b.
66 m = nby(1)
67 IF( m < 0) RETURN
68C
69 IF (iflag == 1) THEN
70C
71c NSL = NBY(2)
72 xg = rby(2)
73 yg = rby(3)
74 zg = rby(4)
75 DO i=1,nsl
76 n = nod(i)
77 IF (weight(n) == 1) THEN
78 f1(i) = af(1,n)
79 f2(i) = af(2,n)
80 f3(i) = af(3,n)
81 rx = x(1,n) - xg
82 ry = x(2,n) - yg
83 rz = x(3,n) - zg
84 f4(i) = am(1,n) + ry*af(3,n) - rz*af(2,n)
85 f5(i) = am(2,n) + rz*af(1,n) - rx*af(3,n)
86 f6(i) = am(3,n) + rx*af(2,n) - ry*af(1,n)
87 ELSE
88 f1(i) = zero
89 f2(i) = zero
90 f3(i) = zero
91 f4(i) = zero
92 f5(i) = zero
93 f6(i) = zero
94 END IF
95 ENDDO
96C
97C Traitement Parith/ON avant echange
98C
99 DO k = 1, 6
100 rbf6(1,k) = zero
101 rbf6(2,k) = zero
102 rbf6(3,k) = zero
103 rbf6(4,k) = zero
104 rbf6(5,k) = zero
105 rbf6(6,k) = zero
106 END DO
107 CALL sum_6_float(1 ,nsl ,f1, rbf6(1,1), 6)
108 CALL sum_6_float(1 ,nsl ,f2, rbf6(2,1), 6)
109 CALL sum_6_float(1 ,nsl ,f3, rbf6(3,1), 6)
110 CALL sum_6_float(1 ,nsl ,f4, rbf6(4,1), 6)
111 CALL sum_6_float(1 ,nsl ,f5, rbf6(5,1), 6)
112 CALL sum_6_float(1 ,nsl ,f6, rbf6(6,1), 6)
113C
114C phase 2 :
115C
116 ELSEIF (iflag==2) THEN
117C
118C Traitement Parith/ON apres echange
119C
120 arby(1,m) = rbf6(1,1)+rbf6(1,2)+rbf6(1,3)+
121 + rbf6(1,4)+rbf6(1,5)+rbf6(1,6)
122 arby(2,m) = rbf6(2,1)+rbf6(2,2)+rbf6(2,3)+
123 + rbf6(2,4)+rbf6(2,5)+rbf6(2,6)
124 arby(3,m) = rbf6(3,1)+rbf6(3,2)+rbf6(3,3)+
125 + rbf6(3,4)+rbf6(3,5)+rbf6(3,6)
126 arrby(1,m) = rbf6(4,1)+rbf6(4,2)+rbf6(4,3)+
127 + rbf6(4,4)+rbf6(4,5)+rbf6(4,6)
128 arrby(2,m) = rbf6(5,1)+rbf6(5,2)+rbf6(5,3)+
129 + rbf6(5,4)+rbf6(5,5)+rbf6(5,6)
130 arrby(3,m) = rbf6(6,1)+rbf6(6,2)+rbf6(6,3)+
131 + rbf6(6,4)+rbf6(6,5)+rbf6(6,6)
132C
133 lcod=icod(1,m)
134 rb(1)= rby(5)
135 rb(2)= rby(6)
136 rb(3)= rby(7)
137 rb(4)= rby(8)
138 rb(5)= rby(9)
139 rb(6)= rby(10)
140 rb(7)= rby(11)
141 rb(8)= rby(12)
142 rb(9)= rby(13)
143 rb(10)= rby(14)
144 rb(11)= rby(15)
145 rb(12)= rby(16)
146 in = rby(17)
147 IF(lcod > 0)THEN
148C rotation de la matrice d'orientation (directions principales)
149 vi(1)=rb(1)*vrrby(1,m) + rb(2)*vrrby(2,m)
150 . + rb(3)*vrrby(3,m)
151 vi(2)=rb(4)*vrrby(1,m) + rb(5)*vrrby(2,m)
152 . + rb(6)*vrrby(3,m)
153 vi(3)=rb(7)*vrrby(1,m) + rb(8)*vrrby(2,m)
154 . + rb(9)*vrrby(3,m)
155 CALL rotbmr(vi,rb,dt1)
156C
157C matrice d'inertie en repere global
158 ii1=rb(10)*rb(1)
159 ii2=rb(10)*rb(2)
160 ii3=rb(10)*rb(3)
161 ii4=rb(11)*rb(4)
162 ii5=rb(11)*rb(5)
163 ii6=rb(11)*rb(6)
164 ii7=rb(12)*rb(7)
165 ii8=rb(12)*rb(8)
166 ii9=rb(12)*rb(9)
167C
168 rby(18)=rb(1)*ii1 + rb(4)*ii4 + rb(7)*ii7
169 rby(19)=rb(1)*ii2 + rb(4)*ii5 + rb(7)*ii8
170 rby(20)=rb(1)*ii3 + rb(4)*ii6 + rb(7)*ii9
171 rby(21)=rb(2)*ii1 + rb(5)*ii4 + rb(8)*ii7
172 rby(22)=rb(2)*ii2 + rb(5)*ii5 + rb(8)*ii8
173 rby(23)=rb(2)*ii3 + rb(5)*ii6 + rb(8)*ii9
174 rby(24)=rb(3)*ii1 + rb(6)*ii4 + rb(9)*ii7
175 rby(25)=rb(3)*ii2 + rb(6)*ii5 + rb(9)*ii8
176 rby(26)=rb(3)*ii3 + rb(6)*ii6 + rb(9)*ii9
177C
178C ajout des termes [Iglobal] vr ^ vr
179 wa1=rby(18)*vrrby(1,m)+rby(19)*vrrby(2,m)+rby(20)*vrrby(3,m)
180 wa2=rby(21)*vrrby(1,m)+rby(22)*vrrby(2,m)+rby(23)*vrrby(3,m)
181 wa3=rby(24)*vrrby(1,m)+rby(25)*vrrby(2,m)+rby(26)*vrrby(3,m)
182C
183 arrby(1,m)=arrby(1,m) + (wa2*vrrby(3,m)-wa3*vrrby(2,m))
184 arrby(2,m)=arrby(2,m) + (wa3*vrrby(1,m)-wa1*vrrby(3,m))
185 arrby(3,m)=arrby(3,m) + (wa1*vrrby(2,m)-wa2*vrrby(1,m))
186C------------------
187C REPERE GLOBAL :
188C Resolution [Iglobal] gama = M, compte-tenu des conditions aux limites
189C Ex : gamaz=0
190C | Ixx Ixy Ixz | { gamax } { Mx }
191C | Iyx Iyy Iyz | { gamay } = { My }
192C | Izx Izy Izz | { 0 } { Mz + DMz } DMz inconnue
193C equivaut a
194C | Ixx Ixy | { gamax } { Mx }
195C | Iyx Iyy | { gamay } = { My }
196C et gamaz=0
197C------------------
198 IF(lcod == 1)THEN
199 det=one/(rby(18)*rby(22) - rby(19)*rby(21))
200 wa1=arrby(1,m)
201 wa2=arrby(2,m)
202 arrby(1,m)=( rby(22)*wa1 - rby(21)*wa2)*det
203 arrby(2,m)=(-rby(19)*wa1 + rby(18)*wa2)*det
204 arrby(3,m)= zero
205 vrrby(3,m) = zero
206 ELSEIF(lcod == 2)THEN
207 det=one/(rby(18)*rby(26) - rby(20)*rby(24))
208 wa1=arrby(1,m)
209 wa2=arrby(3,m)
210 arrby(1,m)=( rby(26)*wa1 - rby(24)*wa2)*det
211 arrby(2,m)= zero
212 arrby(3,m)=(-rby(20)*wa1 + rby(18)*wa2)*det
213 vrrby(2,m) = zero
214 ELSEIF(lcod == 3)THEN
215 arrby(1,m)=arrby(1,m)/rby(18)
216 arrby(2,m)=zero
217 arrby(3,m)=zero
218 vrrby(2,m) = zero
219 vrrby(3,m) = zero
220 ELSEIF(lcod == 4)THEN
221 det=one/(rby(22)*rby(26) - rby(23)*rby(25))
222 wa1=arrby(2,m)
223 wa2=arrby(3,m)
224 arrby(1,m)=zero
225 arrby(2,m)=( rby(26)*wa1 - rby(25)*wa2)*det
226 arrby(3,m)=(-rby(23)*wa1 + rby(22)*wa2)*det
227 vrrby(1,m) = zero
228 ELSEIF(lcod == 5)THEN
229 arrby(1,m) =zero
230 arrby(2,m) =arrby(2,m)/rby(22)
231 arrby(3,m) = zero
232 vrrby(1,m) = zero
233 vrrby(3,m) = zero
234 ELSEIF(lcod == 6)THEN
235 arrby(1,m)=zero
236 arrby(2,m)=zero
237 arrby(3,m)=arrby(3,m)/rby(26)
238 vrrby(1,m) = zero
239 vrrby(2,m) = zero
240 ELSEIF(lcod == 7)THEN
241 arrby(1,m) = zero
242 arrby(2,m) = zero
243 arrby(3,m) = zero
244 vrrby(1,m) = zero
245 vrrby(2,m) = zero
246 vrrby(3,m) = zero
247 ENDIF
248C
249 ELSE
250C------------------------
251C calcul des accelerations dans le r p re globale du main
252C----------------------------
253C correction des momemt M = m - VRRBY X I.VRRBY
254C
255 wa1=arrby(1,m)
256 wa2=arrby(2,m)
257 wa3=arrby(3,m)
258C repere globale -> repere d'inertie principale
259 arrby(1,m) = rb(1)*wa1 + rb(2)*wa2 + rb(3)*wa3
260 arrby(2,m) = rb(4)*wa1 + rb(5)*wa2 + rb(6)*wa3
261 arrby(3,m) = rb(7)*wa1 + rb(8)*wa2 + rb(9)*wa3
262C les contributions des VRRBY ne sont ajoutees que sur le processeur main
263 vi(1) = rb(1)*vrrby(1,m)+rb(2)*vrrby(2,m)+rb(3)*vrrby(3,m)
264 vi(2) = rb(4)*vrrby(1,m)+rb(5)*vrrby(2,m)+rb(6)*vrrby(3,m)
265 vi(3) = rb(7)*vrrby(1,m)+rb(8)*vrrby(2,m)+rb(9)*vrrby(3,m)
266C
267 CALL rotbmr(vi,rb,dt1)
268 arrby(1,m) = arrby(1,m)
269 . + ((rb(11)-rb(12))*vi(2)*vi(3))
270 arrby(2,m) = arrby(2,m)
271 . + ((rb(12)-rb(10))*vi(3)*vi(1))
272 arrby(3,m) = arrby(3,m)
273 . + ((rb(10)-rb(11))*vi(1)*vi(2))
274C CALCUL D'ONE PSEUDO MOMENT:
275C EN FAIT L'ACCELERATION DE ROTATION * INERTIE DU NOEUD MAIN (IMIN DU RB)
276C
277 wa1 = arrby(1,m)*in/rb(10)
278 wa2 = arrby(2,m)*in/rb(11)
279 wa3 = arrby(3,m)*in/rb(12)
280C
281 arrby(1,m) = rb(1)*wa1 + rb(4)*wa2 + rb(7)*wa3
282 arrby(2,m) = rb(2)*wa1 + rb(5)*wa2 + rb(8)*wa3
283 arrby(3,m) = rb(3)*wa1 + rb(6)*wa2 + rb(9)*wa3
284C
285C MATRICE d'inertie -> repere globale
286 ii1=rb(10)*rb(1)
287 ii2=rb(10)*rb(2)
288 ii3=rb(10)*rb(3)
289 ii4=rb(11)*rb(4)
290 ii5=rb(11)*rb(5)
291 ii6=rb(11)*rb(6)
292 ii7=rb(12)*rb(7)
293 ii8=rb(12)*rb(8)
294 ii9=rb(12)*rb(9)
295C
296 rby(18)=rb(1)*ii1 + rb(4)*ii4 + rb(7)*ii7
297 rby(19)=rb(1)*ii2 + rb(4)*ii5 + rb(7)*ii8
298 rby(20)=rb(1)*ii3 + rb(4)*ii6 + rb(7)*ii9
299 rby(21)=rb(2)*ii1 + rb(5)*ii4 + rb(8)*ii7
300 rby(22)=rb(2)*ii2 + rb(5)*ii5 + rb(8)*ii8
301 rby(23)=rb(2)*ii3 + rb(5)*ii6 + rb(8)*ii9
302 rby(24)=rb(3)*ii1 + rb(6)*ii4 + rb(9)*ii7
303 rby(25)=rb(3)*ii2 + rb(6)*ii5 + rb(9)*ii8
304 rby(26)=rb(3)*ii3 + rb(6)*ii6 + rb(9)*ii9
305 ENDIF
306C
307c NSL =NBY(2)
308 msrby = rby(1)
309 in = rby(17)
310 IF(msrby > 0) THEN
311 arby(1,m) = arby(1,m) / msrby
312 arby(2,m) = arby(2,m) / msrby
313 arby(3,m) = arby(3,m) / msrby
314 ELSE
315 arby(1,m) = zero
316 arby(2,m) = zero
317 arby(3,m) = zero
318 ENDIF
319
320 IF(in > 0) THEN
321 arrby(1,m) = arrby(1,m) / in
322 arrby(2,m) = arrby(2,m) / in
323 arrby(3,m) = arrby(3,m) / in
324 ELSE
325 arrby(1,m) = zero
326 arrby(2,m) = zero
327 arrby(3,m) = zero
328 ENDIF
329C
330 lcod = icod(2,m)
331 IF(lcod == 1) THEN
332 arby(3,m) = zero
333 vrby(3,m) = zero
334 ELSEIF(lcod == 2) THEN
335 arby(2,m) = zero
336 vrby(2,m) = zero
337 ELSEIF(lcod == 3) THEN
338 arby(2,m) = zero
339 arby(3,m) = zero
340 vrby(2,m) = zero
341 vrby(3,m) = zero
342 ELSEIF(lcod == 4) THEN
343 arby(1,m) = zero
344 vrby(1,m) = zero
345 ELSEIF(lcod == 5) THEN
346 arby(1,m) = zero
347 arby(3,m) = zero
348 vrby(1,m) = zero
349 vrby(3,m) = zero
350 ELSEIF(lcod == 6) THEN
351 arby(1,m) = zero
352 arby(2,m) = zero
353 vrby(1,m) = zero
354 vrby(2,m) = zero
355 ELSEIF(lcod == 7) THEN
356 arby(1,m) = zero
357 arby(2,m) = zero
358 arby(3,m) = zero
359 vrby(1,m) = zero
360 vrby(2,m) = zero
361 vrby(3,m) = zero
362 ENDIF
363C
364 DO i=1,nsl
365 n = nod(i)
366C
367 af(1,n)= zero
368 af(2,n)= zero
369 af(3,n)= zero
370C
371 am(1,n)= zero
372 am(2,n)= zero
373 am(3,n)= zero
374C
375 stifr(n)= em20
376 stifn(n)= em20
377 ENDDO
378 ENDIF
379C
380 RETURN
381 END
my_real
#define my_real
Definition cppsort.cpp:32
sum_6_float
subroutine sum_6_float(jft, jlt, f, f6, n)
Definition parit.F:64
rmatpon
subroutine rmatpon(af, am, x, rby, nod, nby, stifn, stifr, weight, nsl, rbf6, icod, arby, vrby, arrby, vrrby, iflag)
Definition rmatpon.F:35
rotbmr
subroutine rotbmr(vr, rby, dt)
Definition rotbmr.F:35