OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mass-fluid_tg.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!|| mass_fluid_tg ../starter/source/fluid/mass-fluid_tg.F
25!||--- called by ------------------------------------------------------
26!|| hm_read_bem ../starter/source/loads/bem/hm_read_bem.f
27!||--- calls -----------------------------------------------------
28!|| ancmsg ../starter/source/output/message/message.F
29!|| intanl_tg ../starter/source/fluid/intanl_tg.F
30!|| intgtg ../starter/source/fluid/inte_tg.F
31!|| inthtg ../starter/source/fluid/inte_tg.F
32!|| invert ../starter/source/constraints/general/rbe3/hm_read_rbe3.f
33!||--- uses -----------------------------------------------------
34!|| message_mod ../starter/share/message_module/message_mod.F
35!||====================================================================
36 SUBROUTINE mass_fluid_tg(IFORM, ILVOUT, NNO, NEL, IBUF, ELEM, X,
37 . NORMAL, AF, MFLE, RHO)
38C-----------------------------------------------
39C M o d u l e s
40C-----------------------------------------------
41 USE message_mod
42C-----------------------------------------------
43C I m p l i c i t T y p e s
44C-----------------------------------------------
45#include "implicit_f.inc"
46C-----------------------------------------------
47C C o m m o n B l o c k s
48C-----------------------------------------------
49#include "units_c.inc"
50C-----------------------------------------------
51C D u m m y A r g u m e n t s
52C-----------------------------------------------
53 INTEGER IFORM, ILVOUT, NNO, NEL, IBUF(*), ELEM(3,*)
54 my_real x(3,*), af(*), normal(3,*), mfle(nel,*), rho
55C-----------------------------------------------
56C L o c a l V a r i a b l e s
57C-----------------------------------------------
58 INTEGER IN, JN, KN, N1, N2, N3, NN1, NN2, NN3, NNJ, IEL, JEL, ERR
59 my_real x1, y1, z1, x2, y2, z2, x3, y3, z3, xq, yq, zq, r2,
60 . xp, yp, zp, x12, y12, z12, x13, y13, z13,
61 . nrx, nry, nrz, d2, area, rval, rvlh, rvlg
62 my_real massx, massy, massz, norm
63 my_real, DIMENSION(:,:), ALLOCATABLE :: bbem, cbem, ebem
64 INTEGER :: AERR
65C-------------------------------------------------------------------
66C 1- Compute Bij Cij
67C-------------------------------------------------------------------
68 IF (ilvout>=1) WRITE(istdo,'(A)') ' .. FLUID MASS MATRIX : ASSEMBLY OF INTEGRAL OPERATORS'
69
70 ALLOCATE(bbem(nel,nel), stat = aerr)
71 IF (aerr /= 0) THEN
72 CALL ancmsg(msgid = 1710, anmode=aninfo, msgtype = msgerror)
73 ENDIF
74 ALLOCATE(cbem(nel,nel), stat = aerr)
75 IF (aerr /= 0) THEN
76 CALL ancmsg(msgid = 1710, anmode=aninfo, msgtype = msgerror)
77 ENDIF
78 DO in=1,nel
79 DO jn=1,nel
80 bbem(in,jn)=zero
81 cbem(in,jn)=zero
82 ENDDO
83 ENDDO
84 IF(iform == 1) THEN
85C Collocation BEM - Gauss integration
86 DO iel=1,nel
87C IF (ILVOUT>=2) CALL PROGBAR_C(IEL,NEL)
88 n1=elem(1,iel)
89 n2=elem(2,iel)
90 n3=elem(3,iel)
91 nn1=ibuf(n1)
92 nn2=ibuf(n2)
93 nn3=ibuf(n3)
94 x1=x(1,nn1)
95 x2=x(1,nn2)
96 x3=x(1,nn3)
97 y1=x(2,nn1)
98 y2=x(2,nn2)
99 y3=x(2,nn3)
100 z1=x(3,nn1)
101 z2=x(3,nn2)
102 z3=x(3,nn3)
103C Control point P
104 xp=third*(x1+x2+x3)
105 yp=third*(y1+y2+y3)
106 zp=third*(z1+z2+z3)
107 d2=min((xp-x1)**2+(yp-y1)**2+(zp-z1)**2,
108 . (xp-x2)**2+(yp-y2)**2+(zp-z2)**2,
109 . (xp-x3)**2+(yp-y3)**2+(zp-z3)**2)
110
111 nrx=normal(1,iel)
112 nry=normal(2,iel)
113 nrz=normal(3,iel)
114 DO jel=1,nel
115 IF(iel == jel) THEN
116 bbem(iel,jel)=two*sqrt(pi*af(jel))
117 cbem(iel,jel)=two*pi
118 ELSE
119 n1=elem(1,jel)
120 n2=elem(2,jel)
121 n3=elem(3,jel)
122 nn1=ibuf(n1)
123 nn2=ibuf(n2)
124 nn3=ibuf(n3)
125 x1=x(1,nn1)
126 x2=x(1,nn2)
127 x3=x(1,nn3)
128 y1=x(2,nn1)
129 y2=x(2,nn2)
130 y3=x(2,nn3)
131 z1=x(3,nn1)
132 z2=x(3,nn2)
133 z3=x(3,nn3)
134 area=af(jel)
135C Point source Q
136 xq=third*(x1+x2+x3)
137 yq=third*(y1+y2+y3)
138 zq=third*(z1+z2+z3)
139 CALL inthtg(x1 , y1, z1, x2, y2, z2,
140 . x3, y3, z3, xp, yp, zp,
141 . nrx,nry, nrz, d2, area,
142 . xq, yq, zq, rval )
143 cbem(iel,jel)=rval
144 CALL intgtg(x1 , y1, z1, x2, y2, z2,
145 . x3, y3, z3, xp, yp, zp,
146 . d2, area,
147 . xq, yq, zq, rval )
148 bbem(iel,jel)=rval
149 ENDIF
150 ENDDO
151 ENDDO
152
153 ELSEIF(iform == 2) THEN
154C Analytical integration
155 DO iel=1,nel
156C IF (ILVOUT>=2) CALL PROGBAR_C(IEL,NEL)
157 n1=elem(1,iel)
158 n2=elem(2,iel)
159 n3=elem(3,iel)
160 nn1=ibuf(n1)
161 nn2=ibuf(n2)
162 nn3=ibuf(n3)
163 x1=x(1,nn1)
164 x2=x(1,nn2)
165 x3=x(1,nn3)
166 y1=x(2,nn1)
167 y2=x(2,nn2)
168 y3=x(2,nn3)
169 z1=x(3,nn1)
170 z2=x(3,nn2)
171 z3=x(3,nn3)
172 xp=third*(x1+x2+x3)
173 yp=third*(y1+y2+y3)
174 zp=third*(z1+z2+z3)
175 DO jel=1,nel
176 n1=elem(1,jel)
177 n2=elem(2,jel)
178 n3=elem(3,jel)
179 nn1=ibuf(n1)
180 nn2=ibuf(n2)
181 nn3=ibuf(n3)
182 x1=x(1,nn1)
183 x2=x(1,nn2)
184 x3=x(1,nn3)
185 y1=x(2,nn1)
186 y2=x(2,nn2)
187 y3=x(2,nn3)
188 z1=x(3,nn1)
189 z2=x(3,nn2)
190 z3=x(3,nn3)
191 nrx=normal(1,jel)
192 nry=normal(2,jel)
193 nrz=normal(3,jel)
194 area=af(jel)
195C Matrice C=H
196 CALL intanl_tg(x1, y1, z1, x2, y2, z2,
197 . x3, y3, z3, xp, yp, zp,
198 . nrx, nry, nrz,
199 . area, rvlh, rvlg, iel, jel )
200 IF(iel == jel) THEN
201 cbem(iel,jel)=two*pi
202 ELSE
203 cbem(iel,jel)=rvlh
204 ENDIF
205C Matrice B=G
206 bbem(iel,jel)=rvlg
207 ENDDO
208 ENDDO
209 ENDIF
210
211 IF(ilvout>=3) THEN
212 WRITE (*,'(//,A)') ' BBEM MATRIX'
213 IF(nel < 11) THEN
214 DO in = 1,nel
215 WRITE (*,'(10E13.5)') (bbem(in,jn),jn=1,nel)
216 ENDDO
217 ELSE
218 DO in = 1,10
219 WRITE (*,'(10E13.5)') (bbem(in,jn),jn=1,10)
220 ENDDO
221 WRITE (*,'(//,A)') ' BBEM MATRIX B1J'
222 WRITE (*,'(10E13.5)') (bbem(1,jn),jn=1,nel)
223 ENDIF
224
225 WRITE (*,'(//,A)') ' CBEM MATRIX'
226 IF(nel < 11) THEN
227 DO in = 1,nel
228 WRITE (*,'(10E13.5)') (cbem(in,jn),jn=1,nel)
229 ENDDO
230 ELSE
231 DO in = 1,10
232 WRITE (*,'(10E13.5)') (cbem(in,jn),jn=1,10)
233 ENDDO
234 WRITE (*,'(//,A)') ' CBEM MATRIX C1J'
235 WRITE (*,'(10E13.5)') (cbem(1,jn),jn=1,nel)
236 ENDIF
237 ENDIF
238C---------------------------------------------------
239C 2- Compute fluid mass matrix MFLE
240C---------------------------------------------------
241 IF (ilvout>=1) WRITE(istdo,'(A)') ' .. FLUID MASS MATRIX'
242 ALLOCATE(ebem(nel,nel), stat = aerr)
243 IF (aerr /= 0) THEN
244 CALL ancmsg(msgid = 1710, anmode=aninfo, msgtype = msgerror)
245 ENDIF
246 CALL invert(cbem,ebem,nel,err)
247
248C-----------------------------------------
249C 5 Fluid Mass matrix
250C-----------------------------------------
251 DO in=1,nel
252 DO jn=1,nel
253 mfle(in,jn)=zero
254 cbem(in,jn)=zero
255 ENDDO
256 ENDDO
257 DO in=1,nel
258 DO jn=1,nel
259 DO kn=1,nel
260 cbem(in,jn)=cbem(in,jn)+bbem(in,kn)*ebem(kn,jn)
261 ENDDO
262 ENDDO
263 ENDDO
264 IF(ilvout>=3) THEN
265 WRITE (*,'(//,A)') ' EBEM MATRIX = BBEM*CBEM-1'
266 IF(nel < 11) THEN
267 DO in = 1,nel
268 WRITE (*,'(10E13.5)') (cbem(in,jn),jn=1,nel)
269 ENDDO
270 ELSE
271 DO in = 1,10
272 WRITE (*,'(10E13.5)') (cbem(in,jn),jn=1,10)
273 ENDDO
274 WRITE (*,'(//,A)') ' EBEM MATRIX E1I'
275 WRITE (*,'(10E13.5)') (cbem(1,jn),jn=1,nel)
276 ENDIF
277 ENDIF
278
279 DO in=1,nel
280 DO jn=1,nel
281 mfle(in,jn)=half*rho*af(in)*(cbem(in,jn)+cbem(jn,in))
282 bbem(in,jn)=mfle(in,jn)
283 ENDDO
284 ENDDO
285
286 IF(ilvout>=3) THEN
287 WRITE (*,'(//,A)') ' FLUID MASS MATRIX'
288 IF(nel < 11) THEN
289 DO in = 1,nel
290 WRITE (*,'(10E13.5)') (mfle(in,jn),jn=1,nel)
291 END DO
292 ELSE
293 DO in = 1,10
294 WRITE (*,'(10E13.5)') (mfle(in,jn),jn=1,10)
295 END DO
296 WRITE (*,'(//,A)') ' FLUID MASS MATRIX M1I'
297 WRITE (*,'(10E13.5)') (mfle(1,jn),jn=1,nel)
298 ENDIF
299 ENDIF
300C-----------------------------------------
301C Compute rigid body fluid mass
302C-----------------------------------------
303 massx=zero
304 massy=zero
305 massz=zero
306 DO in=1,nel
307 DO jn=1,nel
308 massx=massx+normal(1,in)*mfle(in,jn)*normal(1,jn)
309 massy=massy+normal(2,in)*mfle(in,jn)*normal(2,jn)
310 massz=massz+normal(3,in)*mfle(in,jn)*normal(3,jn)
311 ENDDO
312 ENDDO
313 WRITE (iout,'(/7X,A,E13.5)') 'DAA : RIGID BODY FLUID MASS XX', massx
314 WRITE (iout,'( 7X,A,E13.5)') 'DAA : RIGID BODY FLUID MASS YY', massy
315 WRITE (iout,'( 7X,A,E13.5)') 'DAA : RIGID BODY FLUID MASS ZZ', massz
316C---------------------------------------------------
317C 3- Compute inverse of fluid mass matrix
318C---------------------------------------------------
319 DO in=1,nel
320 DO jn=1,nel
321 ebem(in,jn)=zero
322 ENDDO
323 ENDDO
324 CALL invert(bbem,ebem,nel,err)
325C
326 IF(ilvout>=3) THEN
327 WRITE (*,'(//,A)') ' INVERSE FLUID MASS MATRIX'
328 IF(nel < 11) THEN
329 DO in = 1,nel
330 WRITE (*,'(10E13.5)') (ebem(in,jn),jn=1,nel)
331 END DO
332 ELSE
333 DO in = 1,10
334 WRITE (*,'(10E13.5)') (ebem(in,jn),jn=1,10)
335 END DO
336 WRITE (*,'(//,A)') ' INVERSE FLUID MASS MATRIX M1I'
337 WRITE (*,'(10E13.5)') (ebem(1,jn),jn=1,nel)
338 ENDIF
339 ENDIF
340
341
342 DO in=1,nel
343 DO jn=1,nel
344 mfle(in,jn)=ebem(in,jn)
345 ENDDO
346 ENDDO
347
348 IF (ALLOCATED(bbem)) DEALLOCATE(bbem)
349 IF (ALLOCATED(ebem)) DEALLOCATE(ebem)
350 IF (ALLOCATED(cbem)) DEALLOCATE(cbem)
351
352 RETURN
353 END
354
355
my_real
#define my_real
Definition cppsort.cpp:32
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
hm_read_bem
subroutine hm_read_bem(igrsurf, iflow, rflow, npc, igrnod, memflow, unitab, x, nom_opt, lgauge, igrv, lsubmodel, iresp)
Definition hm_read_bem.F:274
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
invert
subroutine invert(matrix, inverse, n, errorflag)
Definition hm_read_rbe3.F:1354
hm_read_rbe3
subroutine hm_read_rbe3(irbe3, lrbe3, frbe3, itab, itabm1, igrnod, iskn, lxintd, ikine, iddlevel, nom_opt, itagnd, grnod_uid, unitab, lsubmodel)
Definition hm_read_rbe3.F:49
intanl_tg
subroutine intanl_tg(x1, y1, z1, x2, y2, z2, x3, y3, z3, xp, yp, zp, nrx, nry, nrz, area, rvlh, rvlg, jel, iel)
Definition intanl_tg.F:33
intgtg
subroutine intgtg(x1, y1, z1, x2, y2, z2, x3, y3, z3, xp, yp, zp, d2, jac, xs, ys, zs, rval)
Definition inte_tg.F:145
inthtg
subroutine inthtg(x1, y1, z1, x2, y2, z2, x3, y3, z3, xp, yp, zp, nrx, nry, nrz, d2, jac, xs, ys, zs, rval)
Definition inte_tg.F:33
min
#define min(a, b)
Definition macros.h:20
mass_fluid_tg
subroutine mass_fluid_tg(iform, ilvout, nno, nel, ibuf, elem, x, normal, af, mfle, rho)
Definition mass-fluid_tg.F:38
message_mod
Definition message_mod.F:1249
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
starter
program starter
Definition starter.F:39