OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mass-fluid_qd.F File Reference
#include "implicit_f.inc"
#include "units_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine mass_fluid_qd (nno, nel, iflow, ibuf, elem, x, normal, af, mfle, cbem, rho, iresp)

Function/Subroutine Documentation

◆ mass_fluid_qd()

subroutine mass_fluid_qd ( integer nno,
integer nel,
integer, dimension(*) iflow,
integer, dimension(*) ibuf,
integer, dimension(5,*) elem,
x,
normal,
af,
mfle,
cbem,
rho,
integer, intent(in) iresp )

Definition at line 39 of file mass-fluid_qd.F.

41C-----------------------------------------------
42C M o d u l e s
43C-----------------------------------------------
44 USE message_mod
45C-----------------------------------------------
46C I m p l i c i t T y p e s
47C-----------------------------------------------
48#include "implicit_f.inc"
49C-----------------------------------------------
50C C o m m o n B l o c k s
51C-----------------------------------------------
52#include "units_c.inc"
53C-----------------------------------------------
54C D u m m y A r g u m e n t s
55C-----------------------------------------------
56 INTEGER NNO, NEL, IFLOW(*), IBUF(*), ELEM(5,*)
57 my_real x(3,*), af(*), normal(3,*), mfle(nel,*), cbem(nel,*), rho
58 INTEGER,INTENT(IN):: IRESP
59C-----------------------------------------------
60C L o c a l V a r i a b l e s
61C-----------------------------------------------
62 INTEGER IFORM, KFORM, ILVOUT, NELMAX
63 INTEGER IN, JN, KN, N1, N2, N3, N4, N5, NN1, NN2, NN3, NN4, NNJ, IEL, JEL, ERR
64 my_real x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, xq, yq, zq, r2,
65 . xp, yp, zp, x13, y13, z13, x24, y24, z24,
66 . nrx, nry, nrz, area, d2, rval, rvlh, rvlg
67 my_real massx, massy, massz, wi(4,2), sum
68 my_real, DIMENSION(:,:), ALLOCATABLE :: bbem, ebem
69 my_real :: alpha, beta
70 INTEGER :: AERR
71 iform = iflow(13)
72 ilvout = iflow(17)
73 kform = iflow(23)
74 nelmax = iflow(28)
75C-------------------------------------------------------------------
76C 1- Compute Bij Cij and Ai
77C-------------------------------------------------------------------
78 IF (ilvout>=1) THEN
79 WRITE(istdo,'(A)') ' .. FLUID MASS MATRIX:
80 . ASSEMBLY OF INTEGRAL OPERATORS'
81 ENDIF
82C Collocation BEM
83 wi(1,1)=fourth
84 wi(2,1)=fourth
85 wi(3,1)=fourth
86 wi(4,1)=fourth
87 wi(1,2)=third
88 wi(2,2)=third
89 wi(3,2)=one_over_6
90 wi(4,2)=one_over_6
91 ALLOCATE(bbem(nel,nel), stat = aerr)
92 IF (aerr /= 0) THEN
93 CALL ancmsg(msgid = 1710, anmode=aninfo, msgtype = msgerror)
94 ENDIF
95 DO in=1,nel
96 DO jn=1,nel
97 bbem(in,jn)=zero
98 cbem(in,jn)=zero
99 ENDDO
100 ENDDO
101
102 IF(iform == 1) THEN
103C Gauss integration
104 DO iel=1,nel
105C IF (ILVOUT>=2) CALL PROGBAR_C(IEL,NEL)
106 n1=elem(1,iel)
107 n2=elem(2,iel)
108 n3=elem(3,iel)
109 n4=elem(4,iel)
110 n5=elem(5,iel)
111 nn1=ibuf(n1)
112 nn2=ibuf(n2)
113 nn3=ibuf(n3)
114 nn4=ibuf(n4)
115 x1=x(1,nn1)
116 x2=x(1,nn2)
117 x3=x(1,nn3)
118 x4=x(1,nn4)
119 y1=x(2,nn1)
120 y2=x(2,nn2)
121 y3=x(2,nn3)
122 y4=x(2,nn4)
123 z1=x(3,nn1)
124 z2=x(3,nn2)
125 z3=x(3,nn3)
126 z4=x(3,nn4)
127C Control point P
128 xp=wi(1,n5)*x1+wi(2,n5)*x2+wi(3,n5)*x3+wi(4,n5)*x4
129 yp=wi(1,n5)*y1+wi(2,n5)*y2+wi(3,n5)*y3+wi(4,n5)*y4
130 zp=wi(1,n5)*z1+wi(2,n5)*z2+wi(3,n5)*z3+wi(4,n5)*z4
131 d2=min((xp-x1)**2+(yp-y1)**2+(zp-z1)**2,
132 . (xp-x2)**2+(yp-y2)**2+(zp-z2)**2,
133 . (xp-x3)**2+(yp-y3)**2+(zp-z3)**2,
134 . (xp-x4)**2+(yp-y4)**2+(zp-z4)**2)
135 nrx=normal(1,iel)
136 nry=normal(2,iel)
137 nrz=normal(3,iel)
138C
139 DO jel=1,nel
140 IF(iel == jel) THEN
141 bbem(iel,jel)=two*sqrt(pi*af(jel))
142 cbem(iel,jel)=two*pi
143 ELSE
144 n1=elem(1,jel)
145 n2=elem(2,jel)
146 n3=elem(3,jel)
147 n4=elem(4,jel)
148 jn=elem(5,jel)
149 nn1=ibuf(n1)
150 nn2=ibuf(n2)
151 nn3=ibuf(n3)
152 nn4=ibuf(n4)
153 x1=x(1,nn1)
154 x2=x(1,nn2)
155 x3=x(1,nn3)
156 x4=x(1,nn4)
157 y1=x(2,nn1)
158 y2=x(2,nn2)
159 y3=x(2,nn3)
160 y4=x(2,nn4)
161 z1=x(3,nn1)
162 z2=x(3,nn2)
163 z3=x(3,nn3)
164 z4=x(3,nn4)
165 area=af(jel)
166C Position de la source Q
167 xq=wi(1,jn)*x1+wi(2,jn)*x2+wi(3,jn)*x3+wi(4,jn)*x4
168 yq=wi(1,jn)*y1+wi(2,jn)*y2+wi(3,jn)*y3+wi(4,jn)*y4
169 zq=wi(1,jn)*z1+wi(2,jn)*z2+wi(3,jn)*z3+wi(4,jn)*z4
170
171 IF(jn == 1) THEN
172 CALL inthqd(x1 , y1, z1, x2, y2, z2,
173 . x3, y3, z3, x4, y4, z4,
174 . xp, yp, zp, xq, yq, zq,
175 . nrx, nry, nrz, d2, area, rval)
176 cbem(iel,jel)=rval
177 CALL intgqd(x1 , y1, z1, x2, y2, z2,
178 . x3, y3, z3, x4, y4, z4,
179 . xp, yp, zp, xq, yq, zq,
180 . d2, area, rval)
181 bbem(iel,jel)=rval
182 ELSEIF(jn == 2) THEN
183 CALL inthtg(x1 , y1, z1, x2, y2, z2,
184 . x3, y3, z3, xp, yp, zp,
185 . nrx,nry, nrz, d2, area,
186 . xq, yq, zq, rval )
187 cbem(iel,jel)=rval
188 CALL intgtg(x1, y1, z1, x2, y2, z2,
189 . x3, y3, z3, xp, yp, zp,
190 . d2, area,
191 . xq, yq, zq, rval )
192 bbem(iel,jel)=rval
193 ENDIF
194 ENDIF
195 ENDDO
196 ENDDO
197
198 ELSEIF(iform == 2) THEN
199C Integration analytique
200 DO iel=1,nel
201C IF (ILVOUT>=2) CALL PROGBAR_C(IEL,NEL)
202 n1=elem(1,iel)
203 n2=elem(2,iel)
204 n3=elem(3,iel)
205 n4=elem(4,iel)
206 in=elem(5,iel)
207 nn1=ibuf(n1)
208 nn2=ibuf(n2)
209 nn3=ibuf(n3)
210 nn4=ibuf(n4)
211 x1=x(1,nn1)
212 x2=x(1,nn2)
213 x3=x(1,nn3)
214 x4=x(1,nn4)
215 y1=x(2,nn1)
216 y2=x(2,nn2)
217 y3=x(2,nn3)
218 y4=x(2,nn4)
219 z1=x(3,nn1)
220 z2=x(3,nn2)
221 z3=x(3,nn3)
222 z4=x(3,nn4)
223C Control point P
224 xp=wi(1,in)*x1+wi(2,in)*x2+wi(3,in)*x3+wi(4,in)*x4
225 yp=wi(1,in)*y1+wi(2,in)*y2+wi(3,in)*y3+wi(4,in)*y4
226 zp=wi(1,in)*z1+wi(2,in)*z2+wi(3,in)*z3+wi(4,in)*z4
227 DO jel=1,nel
228 n1=elem(1,jel)
229 n2=elem(2,jel)
230 n3=elem(3,jel)
231 n4=elem(4,jel)
232 jn=elem(5,jel)
233 nn1=ibuf(n1)
234 nn2=ibuf(n2)
235 nn3=ibuf(n3)
236 nn4=ibuf(n4)
237 x1=x(1,nn1)
238 x2=x(1,nn2)
239 x3=x(1,nn3)
240 x4=x(1,nn4)
241 y1=x(2,nn1)
242 y2=x(2,nn2)
243 y3=x(2,nn3)
244 y4=x(2,nn4)
245 z1=x(3,nn1)
246 z2=x(3,nn2)
247 z3=x(3,nn3)
248 z4=x(3,nn4)
249 nrx=normal(1,jel)
250 nry=normal(2,jel)
251 nrz=normal(3,jel)
252 area=af(jel)
253 IF(jn == 1) THEN
254 CALL intanl_qd(x1, y1, z1, x2, y2, z2, x3, y3, z3,
255 . x4, y4, z4, xp, yp, zp, nrx,nry,nrz,
256 . area, rvlh, rvlg, iel, jel )
257
258 ELSEIF(jn == 2) THEN
259 CALL intanl_tg(x1, y1, z1, x2, y2, z2, x3, y3, z3,
260 . xp, yp, zp, nrx,nry,nrz,
261 . area, rvlh, rvlg, iel, jel )
262 ENDIF
263C Matrice C=H
264 IF(iel == jel) THEN
265 cbem(iel,jel)=two*pi
266 ELSE
267 cbem(iel,jel)=rvlh
268 ENDIF
269C Matrice B=G
270 bbem(iel,jel)=rvlg
271 ENDDO
272 ENDDO
273 ENDIF
274
275 IF(ilvout>=3) THEN
276 WRITE (*,'(//,A)') ' BBEM MATRIX'
277 IF(nel < 11) THEN
278 DO in = 1,nel
279 WRITE (*,'(10E13.5)') (bbem(in,jn),jn=1,nel)
280 ENDDO
281 ELSE
282 DO in = 1,10
283 WRITE (*,'(10E13.5)') (bbem(in,jn),jn=1,10)
284 ENDDO
285 WRITE (*,'(//,A)') ' BBEM MATRIX B1I'
286 WRITE (*,'(10E13.5)') (bbem(1,jn),jn=1,nel)
287 ENDIF
288
289 WRITE (*,'(//,A)') ' CBEM MATRIX'
290 IF(nel < 11) THEN
291 DO in = 1,nel
292 WRITE (*,'(10E13.5)') (cbem(in,jn),jn=1,nel)
293 ENDDO
294 ELSE
295 DO in = 1,10
296 WRITE (*,'(10E13.5)') (cbem(in,jn),jn=1,10)
297 ENDDO
298 WRITE (*,'(//,A)') ' CBEM MATRIX C1I'
299 WRITE (*,'(10E13.5)') (cbem(1,jn),jn=1,nel)
300 ENDIF
301 ENDIF
302C------------------------------------------------------------------------
303C 2- Compute matrix MFLE = CBEMt * BBEM + BBEMt * CBEM (Case DAA KFORM=1)
304C------------------------------------------------------------------------
305 IF(kform == 1 .AND. nel >= nelmax) THEN
306c MFLE = CBEMt*BBEM ! DGEMM C := alpha*op( A )*op( B ) +beta*C
307C C = MFLE ; beta = 0
308C op = T A = CBEM
309C op = N , B = BBEM
310c MFLE = MFLE + BBEMt*CBEM
311C C = MFLE ; beta = 1
312C op = N A = BBEMt
313C op = T , B = CBEM
314
315C CMAT(1:NEL,1:NEL) = TRANSPOSE(BBEM(1:NEL,1:NEL))
316C MFLE(1:NEL,1:NEL) = MATMUL(CMAT(1:NEL,1:NEL),CBEM(1:NEL,1:NEL))
317C CMAT(1:NEL,1:NEL) = TRANSPOSE(CBEM(1:NEL,1:NEL))
318C MFLE(1:NEL,1:NEL) = MFLE(1:NEL,1:NEL) + MATMUL(CMAT(1:NEL,1:NEL),BBEM(1:NEL,1:NEL))
319C CMAT(1:NEL,1:NEL) = CBEM(1:NEL,1:NEL)
320
321 IF (iresp == 0)THEN
322 ! Double Precision version : use either Lapack / MKP or ARMPL
323
324 alpha = 1.0d0
325 beta = 0.0d0
326 CALL dgemm('T','N',nel,nel,nel,alpha,cbem,nel,bbem,nel,beta,mfle,nel)
327 alpha = 1.0d0
328 beta = 1.0d0
329 CALL dgemm('N','T',nel,nel,nel,alpha,bbem,nel,cbem,nel,beta,mfle,nel)
330 ELSE
331 ! Single Precision version / Bad performance
332 DO jn=1,nel
333 DO in=1,nel
334 sum=zero
335 DO kn=1,nel
336 sum=sum+bbem(kn,in)*cbem(kn,jn)+cbem(kn,in)*bbem(kn,jn)
337 ENDDO
338 mfle(in,jn)=sum
339 ENDDO
340 ENDDO
341 ENDIF
342 RETURN
343 ENDIF
344C---------------------------------------------------
345C 2- Compute fluid mass matrix MFLE
346C---------------------------------------------------
347 IF (ilvout>=1) WRITE(istdo,'(A)') ' .. FLUID MASS MATRIX'
348 ALLOCATE(ebem(nel,nel), stat = aerr)
349 IF (aerr /= 0) THEN
350 CALL ancmsg(msgid = 1710, anmode=aninfo, msgtype = msgerror)
351 ENDIF
352 CALL invert(cbem,ebem,nel,err)
353
354 IF(ilvout>=3) THEN
355 WRITE (*,'(//,A)') ' CBEM-1 MATRIX'
356 IF(nel < 11) THEN
357 DO in = 1,nel
358 WRITE (*,'(10E13.5)') (ebem(in,jn),jn=1,nel)
359 END DO
360 ELSE
361 DO in = 1,10
362 WRITE (*,'(10E13.5)') (ebem(in,jn),jn=1,10)
363 END DO
364 WRITE (*,'(//,A)') ' CBEM-1 MATRIX C1I'
365 WRITE (*,'(10E13.5)') (ebem(1,jn),jn=1,nel)
366 ENDIF
367 ENDIF
368
369 DO in=1,nel
370 DO jn=1,nel
371 mfle(in,jn)=zero
372 cbem(in,jn)=zero
373 ENDDO
374 ENDDO
375 cbem(1:nel, 1:nel) = matmul(bbem(1:nel, 1:nel), ebem(1:nel, 1:nel))
376c$$$ DO IN=1,NEL
377c$$$ DO JN=1,NEL
378c$$$ DO KN=1,NEL
379c$$$ CBEM(IN,JN)=CBEM(IN,JN)+BBEM(IN,KN)*EBEM(KN,JN)
380c$$$ ENDDO
381c$$$ ENDDO
382c$$$ ENDDO
383
384 IF(ilvout>=3) THEN
385 WRITE (*,'(//,A)') ' EBEM MATRIX = BBEM*CBEM-1'
386 IF(nel < 11) THEN
387 DO in = 1,nel
388 WRITE (*,'(10E13.5)') (cbem(in,jn),jn=1,nel)
389 ENDDO
390 ELSE
391 DO in = 1,10
392 WRITE (*,'(10E13.5)') (cbem(in,jn),jn=1,10)
393 ENDDO
394 WRITE (*,'(//,A)') ' EBEM MATRIX E1I'
395 WRITE (*,'(10E13.5)') (cbem(1,jn),jn=1,nel)
396 ENDIF
397 ENDIF
398
399 DO in=1,nel
400 DO jn=1,nel
401 mfle(in,jn)=half*rho*af(in)*(cbem(in,jn)+cbem(jn,in))
402 bbem(in,jn)=mfle(in,jn)
403 ENDDO
404 ENDDO
405C
406 IF(ilvout>=3) THEN
407 WRITE (*,'(//,A)') ' FLUID MASS MATRIX'
408 IF(nel < 11) THEN
409 DO in = 1,nel
410 WRITE (*,'(10E13.5)') (mfle(in,jn),jn=1,nel)
411 END DO
412 ELSE
413 DO in = 1,10
414 WRITE (*,'(10E13.5)') (mfle(in,jn),jn=1,10)
415 END DO
416 WRITE (*,'(//,A)') ' FLUID MASS MATRIX M1I'
417 WRITE (*,'(10E13.5)') (mfle(1,jn),jn=1,nel)
418 ENDIF
419 ENDIF
420C---------------------------------
421C COMPUTE RIGD BODY FLUID MASS
422C---------------------------------
423 massx=zero
424 massy=zero
425 massz=zero
426 DO in=1,nel
427 DO jn=1,nel
428 massx=massx+normal(1,in)*mfle(in,jn)*normal(1,jn)
429 massy=massy+normal(2,in)*mfle(in,jn)*normal(2,jn)
430 massz=massz+normal(3,in)*mfle(in,jn)*normal(3,jn)
431 ENDDO
432 ENDDO
433 WRITE (iout,'(/7X,A,E13.5)') 'DAA : RIGID BODY FLUID MASS XX', massx
434 WRITE (iout,'( 7X,A,E13.5)') 'DAA : RIGID BODY FLUID MASS YY', massy
435 WRITE (iout,'( 7X,A,E13.5)') 'DAA : RIGID BODY FLUID MASS ZZ', massz
436
437C---------------------------------------------------
438C 3- Compute inverse of fluid mass matrix
439C---------------------------------------------------
440 IF(kform==1) THEN
441 DO in=1,nel
442 DO jn=1,nel
443 ebem(in,jn)=zero
444 ENDDO
445 ENDDO
446 CALL invert(bbem,ebem,nel,err)
447
448C
449 IF(ilvout>=3) THEN
450 WRITE (*,'(//,A)') ' INVERSE FLUID MASS MATRIX'
451 IF(nel < 11) THEN
452 DO in = 1,nel
453 WRITE (*,'(10E13.5)') (ebem(in,jn),jn=1,nel)
454 END DO
455 ELSE
456 DO in = 1,10
457 WRITE (*,'(10E13.5)') (ebem(in,jn),jn=1,10)
458 END DO
459 WRITE (*,'(//,A)') ' INVERSE FLUID MASS MATRIX MII'
460 WRITE (*,'(10E13.5)') (ebem(jn,jn),jn=1,nel)
461 ENDIF
462 ENDIF
463
464 DO in=1,nel
465 DO jn=1,nel
466 mfle(in,jn)=ebem(in,jn)
467 ENDDO
468 ENDDO
469 ENDIF
470 IF (ALLOCATED(bbem)) DEALLOCATE(bbem)
471 IF (ALLOCATED(ebem))DEALLOCATE(ebem)
472
473 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:187
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
intanl_qd
subroutine intanl_qd(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, xp, yp, zp, nrx, nry, nrz, area, rvlh, rvlg, jel, iel)
Definition intanl_qd.F:31
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
intgqd
subroutine intgqd(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, xp, yp, zp, xs, ys, zs, d2, jac, rval)
Definition inte_qd.F:124
inthqd
subroutine inthqd(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, xp, yp, zp, xs, ys, zs, nrx, nry, nrz, d2, jac, rval)
Definition inte_qd.F:32
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
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