OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
bemsolv.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!|| bemsolv ../engine/source/fluid/bemsolv.F
25!||--- called by ------------------------------------------------------
26!|| incpflow ../engine/source/fluid/incpflow.F
27!||--- calls -----------------------------------------------------
28!|| arret ../engine/source/system/arret.F
29!|| glbem ../engine/source/fluid/bemsolv.F
30!|| glsinf ../engine/source/fluid/bemsolv.F
31!|| intgtg ../engine/source/fluid/bemsolv.F
32!|| inthtg ../engine/source/fluid/bemsolv.F
33!|| progbar_c ../common_source/sortie/progbar_c.c
34!||====================================================================
35 SUBROUTINE bemsolv(L_ASSEMB, ILVOUT, IFORM, NN , NEL ,
36 . HBEM , GBEM , IBUF , TBEMTG, TELNOR ,
37 . X , Q , PHI , IPIV , PHI_INF)
38C-----------------------------------------------
39C I m p l i c i t T y p e s
40C-----------------------------------------------
41#include "implicit_f.inc"
42C-----------------------------------------------
43C C o m m o n B l o c k s
44C-----------------------------------------------
45#include "units_c.inc"
46C-----------------------------------------------
47C D u m m y A r g u m e n t s
48C-----------------------------------------------
49 INTEGER ILVOUT, IFORM, NN, NEL, IBUF(*), TBEMTG(3,*), IPIV(*)
50 my_real
51 . hbem(nn,*), gbem(nn,*), telnor(3,*), x(3,*), q(*), phi(*),
52 . phi_inf(*)
53 LOGICAL L_ASSEMB
54C-----------------------------------------------
55C L o c a l V a r i a b l e s
56C-----------------------------------------------
57 INTEGER IN, JN, IQ, N1, N2, N3, NN1, NN2, NN3, NNJ, INFO, IEL
58 my_real
59 . x1, y1, z1, x2, y2, z2, x3, y3, z3, x0, y0, z0, d2,
60 . nrx, nry, nrz, area2, rval(3), rsum, s(nn)
61C---------------------------------------------------
62C 1- Assemblage
63C---------------------------------------------------
64#if defined(MYREAL8) && !defined(WITHOUT_LINALG)
65 IF (l_assemb) THEN
66 IF (ilvout>=1) WRITE(istdo,'(A)')
67 . ' ** BEMSOLV : ASSEMBLY OF INTEGRAL OPERATORS'
68 IF (iform==1) THEN
69C Collocation BEM - Elements lineaires en potentiel, constants en flux
70 DO in=1,nn
71 DO jn=1,nn
72 hbem(in,jn)=zero
73 ENDDO
74 DO iq=1,nel
75 gbem(in,iq)=zero
76 ENDDO
77 gbem(in,nel+1)=zero
78 ENDDO
79 DO iel=1,nel
80 IF (ilvout==2) CALL progbar_c(iel,nel)
81 n1=tbemtg(1,iel)
82 n2=tbemtg(2,iel)
83 n3=tbemtg(3,iel)
84 nn1=ibuf(n1)
85 nn2=ibuf(n2)
86 nn3=ibuf(n3)
87 x1=x(1,nn1)
88 x2=x(1,nn2)
89 x3=x(1,nn3)
90 y1=x(2,nn1)
91 y2=x(2,nn2)
92 y3=x(2,nn3)
93 z1=x(3,nn1)
94 z2=x(3,nn2)
95 z3=x(3,nn3)
96 x0=third*(x1+x2+x3)
97 y0=third*(y1+y2+y3)
98 z0=third*(z1+z2+z3)
99 d2=min((x0-x1)**2+(y0-y1)**2+(z0-z1)**2,
100 . (x0-x2)**2+(y0-y2)**2+(z0-z2)**2,
101 . (x0-x3)**2+(y0-y3)**2+(z0-z3)**2)
102 nrx=telnor(1,iel)
103 nry=telnor(2,iel)
104 nrz=telnor(3,iel)
105 area2=sqrt(nrx**2+nry**2+nrz**2)
106 DO jn=1,nn
107 nnj=ibuf(jn)
108 CALL inthtg(nnj, x1 , y1, z1, x2,
109 . y2, z2, x3, y3, z3,
110 . x0, y0, z0, d2, area2,
111 . nrx, nry, nrz, rval, nn1,
112 . nn2, nn3, x )
113 hbem(jn,n1)=hbem(jn,n1)+rval(1)
114 hbem(jn,n2)=hbem(jn,n2)+rval(2)
115 hbem(jn,n3)=hbem(jn,n3)+rval(3)
116 CALL intgtg(nnj, x1 , y1, z1, x2,
117 . y2, z2, x3, y3, z3,
118 . x0, y0, z0, d2, area2,
119 . nrx, nry, nrz, rval, nn1,
120 . nn2, nn3, x )
121 gbem(jn,iel)=gbem(jn,iel)+rval(1)
122 gbem(jn,nel+1)=gbem(jn,nel+1)+rval(1)
123 ENDDO
124 ENDDO
125 DO in=1,nn
126 rsum=zero
127 DO jn=1,nn
128 IF (in/=jn) rsum=rsum+hbem(in,jn)
129 ENDDO
130 hbem(in,in)=-rsum
131 ENDDO
132 ELSEIF (iform==2) THEN
133C Galerkin BEM - Elements lineaires en potentiel, constants en flux
134 DO in=1,nn
135 DO jn=1,nn
136 hbem(in,jn)=zero
137 ENDDO
138 DO iq=1,nel
139 gbem(in,iq)=zero
140 ENDDO
141 gbem(in,nel+1)=zero
142 ENDDO
143 DO iel=1,nel
144 IF (ilvout==2) CALL progbar_c(iel,nel)
145 n1=tbemtg(1,iel)
146 n2=tbemtg(2,iel)
147 n3=tbemtg(3,iel)
148 nn1=ibuf(n1)
149 nn2=ibuf(n2)
150 nn3=ibuf(n3)
151 x1=x(1,nn1)
152 x2=x(1,nn2)
153 x3=x(1,nn3)
154 y1=x(2,nn1)
155 y2=x(2,nn2)
156 y3=x(2,nn3)
157 z1=x(3,nn1)
158 z2=x(3,nn2)
159 z3=x(3,nn3)
160 x0=third*(x1+x2+x3)
161 y0=third*(y1+y2+y3)
162 z0=third*(z1+z2+z3)
163 nrx=telnor(1,iel)
164 nry=telnor(2,iel)
165 nrz=telnor(3,iel)
166 area2=sqrt(nrx**2+nry**2+nrz**2)
167 CALL glbem(x1 , y1 , z1 , x2 , y2 ,
168 . z2 , x3 , y3 , z3 , x0 ,
169 . y0 , z0 , telnor, hbem , gbem,
170 . tbemtg, ibuf, x , n1 , n2 ,
171 . n3 , iel , nn , area2, nel )
172 ENDDO
173 ENDIF
174C
175 DO in=1,nn
176 hbem(in,nn)=-gbem(in,nel+1)
177 ENDDO
178C
179 CALL dgetrf(nn, nn, hbem, nn, ipiv,
180 . info)
181 ENDIF
182C---------------------------------------------------
183C 2- Resolution
184C---------------------------------------------------
185 IF (ilvout>=1) WRITE(istdo,'(A)')
186 .' ** BEMSOLV : LINEAR SYSTEM SOLUTION'
187 CALL dgemv('N', nn, nel+1, one, gbem,
188 . nn, q, 1, zero, s,
189 . 1 )
190C
191 IF (iform==1) THEN
192 DO in=1,nn
193 s(in)=s(in)-phi_inf(in)
194 ENDDO
195 ELSEIF (iform==2) THEN
196 DO iel=1,nel
197 n1=tbemtg(1,iel)
198 n2=tbemtg(2,iel)
199 n3=tbemtg(3,iel)
200 nrx=telnor(1,iel)
201 nry=telnor(2,iel)
202 nrz=telnor(3,iel)
203 area2=sqrt(nrx**2+nry**2+nrz**2)
204 CALL glsinf(n1, n2, n3, area2, s,
205 . phi_inf)
206 ENDDO
207 ENDIF
208C
209 CALL dgetrs('N', nn, 1, hbem, nn, ipiv,
210 . s, nn, info)
211C
212 DO in=1,nn
213 phi(in)=s(in)
214 ENDDO
215#else
216 CALL arret(5)
217#endif
218 RETURN
219 END
220!||====================================================================
221!|| glsinf ../engine/source/fluid/bemsolv.F
222!||--- called by ------------------------------------------------------
223!|| bemsolv ../engine/source/fluid/bemsolv.F
224!||====================================================================
225 SUBROUTINE glsinf(N1 , N2, N3, JAC, S,
226 . PHI_INF)
227C-----------------------------------------------
228C I m p l i c i t T y p e s
229C-----------------------------------------------
230#include "implicit_f.inc"
231C-----------------------------------------------
232C D u m m y A r g u m e n t s
233C-----------------------------------------------
234 INTEGER N1, N2, N3
235 my_real
236 . JAC, S(*), PHI_INF(*)
237C-----------------------------------------------
238C L o c a l V a r i a b l e s
239C-----------------------------------------------
240 INTEGER NPG, IAD, IAD2, IP
241 my_real
242 . PG(50), WPG(25), W, KSIP, ETAP, VAL1, VAL2, VAL3,
243 . phi_infp
244C
245 DATA pg /.33333333,.33333333,
246 . .33333333,.33333333,
247 . .60000000,.20000000,
248 . .20000000,.60000000,
249 . .20000000,.20000000,
250 . .33333333,.33333333,
251 . .79742699,.10128651,
252 . .10128651,.79742699,
253 . .10128651,.10128651,
254 . .05971587,.47014206,
255 . .47014206,.05971587,
256 . .47014206,.47014206,
257 . .06513010,.06513010,
258 . .86973979,.06513010,
259 . .06513010,.86973979,
260 . .31286550,.04869031,
261 . .63844419,.31286550,
262 . .04869031,.63844419,
263 . .63844419,.04869031,
264 . .31286550,.63844419,
265 . .04869031,.31286550,
266 . .26034597,.26034597,
267 . .47930807,.26034597,
268 . .26034597,.47930807,
269 . .33333333,.33333333/
270 DATA wpg /1.00000000,
271 . -.56250000,.52083333,
272 . .52083333,.52083333,
273 . .22500000,.12593918,
274 . .12593918,.12593918,
275 . .13239415,.13239415,
276 . .13239415,
277 . .05334724,.05334724,
278 . .05334724,.07711376,
279 . .07711376,.07711376,
280 . .07711376,.07711376,
281 . .07711376,.17561526,
282 . .17561526,.17561526,
283 . -.14957004/
284C
285 npg=13
286 iad=13
287C
288 iad2=2*(iad-1)+1
289 DO ip=1,npg
290 w=wpg(iad)
291 ksip=pg(iad2)
292 etap=pg(iad2+1)
293 iad=iad+1
294 iad2=iad2+2
295 val1=one-ksip-etap
296 val2=ksip
297 val3=etap
298 phi_infp=phi_inf(n1)*val1+phi_inf(n2)*val2+phi_inf(n3)*val3
299 s(n1)=s(n1)-w*val1*phi_infp*jac
300 s(n2)=s(n2)-w*val2*phi_infp*jac
301 s(n3)=s(n3)-w*val3*phi_infp*jac
302 ENDDO
303C
304 RETURN
305 END
306!||====================================================================
307!|| inthtg ../engine/source/fluid/bemsolv.F
308!||--- called by ------------------------------------------------------
309!|| bemsolv ../engine/source/fluid/bemsolv.F
310!|| bemsolvp ../engine/source/fluid/bemsolvp.F
311!||====================================================================
312 SUBROUTINE inthtg(NNJ, X1 , Y1, Z1, X2,
313 . Y2, Z2, X3, Y3, Z3,
314 . X0, Y0, Z0, D2, JAC,
315 . NRX, NRY, NRZ, RVAL, NN1,
316 . NN2, NN3, X )
317C-----------------------------------------------
318C I m p l i c i t T y p e s
319C-----------------------------------------------
320#include "implicit_f.inc"
321C-----------------------------------------------
322C D u m m y A r g u m e n t s
323C-----------------------------------------------
324 INTEGER NNJ, NN1, NN2, NN3
325 my_real
326 . X1, Y1, Z1, X2, Y2, Z2, X3, Y3, Z3,
327 . X0, Y0, Z0, D2, JAC, NRX, NRY, NRZ, RVAL(*), X(3,*)
328C-----------------------------------------------
329C L o c a l V a r i a b l e s
330C-----------------------------------------------
331 INTEGER NPG, IAD, IAD2, IP, IHG
332 my_real
333 . pg(50), wpg(25), xs, ys, zs, r2,
334 . val1, val2, val3, w, xg, yg, zg, valphi,
335 . eta1, eta2
336C
337 DATA pg /.33333333,.33333333,
338 . .33333333,.33333333,
339 . .60000000,.20000000,
340 . .20000000,.60000000,
341 . .20000000,.20000000,
342 . .33333333,.33333333,
343 . .79742699,.10128651,
344 . .10128651,.79742699,
345 . .10128651,.10128651,
346 . .05971587,.47014206,
347 . .47014206,.05971587,
348 . .47014206,.47014206,
349 . .06513010,.06513010,
350 . .86973979,.06513010,
351 . .06513010,.86973979,
352 . .31286550,.04869031,
353 . .63844419,.31286550,
354 . .04869031,.63844419,
355 . .63844419,.04869031,
356 . .31286550,.63844419,
357 . .04869031,.31286550,
358 . .26034597,.26034597,
359 . .47930807,.26034597,
360 . .26034597,.47930807,
361 . .33333333,.33333333/
362 DATA wpg /1.00000000,
363 . -.56250000,.52083333,
364 . .52083333,.52083333,
365 . .22500000,.12593918,
366 . .12593918,.12593918,
367 . .13239415,.13239415,
368 . .13239415,
369 . .05334724,.05334724,
370 . .05334724,.07711376,
371 . .07711376,.07711376,
372 . .07711376,.07711376,
373 . .07711376,.17561526,
374 . .17561526,.17561526,
375 . -.14957004/
376C
377C DISTANCE A LA SOURCE
378 xs=x(1,nnj)
379 ys=x(2,nnj)
380 zs=x(3,nnj)
381 r2=(x0-xs)**2+(y0-ys)**2+(z0-zs)**2
382C NOMBRE DE POINTS DE GAUSS
383 IF (r2>hundred*d2) THEN
384 npg=1
385 iad=1
386 ELSEIF (r2>twenty5*d2) THEN
387 npg=4
388 iad=2
389 ELSEIF (r2>four*d2) THEN
390 npg=7
391 iad=6
392 ELSE
393 npg=13
394 iad=13
395 ENDIF
396 IF (nn1==nnj) npg=-1
397 IF (nn2==nnj) npg=-1
398 IF (nn3==nnj) npg=-1
399C INTEGRATION
400 rval(1)=zero
401 rval(2)=zero
402 rval(3)=zero
403 IF (npg<0) RETURN
404 iad2=2*(iad-1)+1
405 DO ip=1,npg
406 w=wpg(iad)
407 eta1=pg(iad2)
408 eta2=pg(iad2+1)
409 iad=iad+1
410 iad2=iad2+2
411 val1=one-eta1-eta2
412 val2=eta1
413 val3=eta2
414 xg=val1*x1+val2*x2+val3*x3
415 yg=val1*y1+val2*y2+val3*y3
416 zg=val1*z1+val2*z2+val3*z3
417 r2=(xg-xs)**2+(yg-ys)**2+(zg-zs)**2
418 IF(r2>em20)THEN
419 valphi=-fourth/pi/(r2**three_half)
420 . *(nrx*(xg-xs)+nry*(yg-ys)+nrz*(zg-zs))
421 rval(1)=rval(1)+half*w*val1*valphi
422 rval(2)=rval(2)+half*w*val2*valphi
423 rval(3)=rval(3)+half*w*val3*valphi
424 ENDIF
425C Remarque : la normale (NRX,NRY,NRZ) n'est pas normalisee
426C il faudrait la diviser par JAC et comme on doit multiplier
427C par JAC a la ligne d'en dessous...
428 ENDDO
429C
430 RETURN
431 END
432!||====================================================================
433!|| intgtg ../engine/source/fluid/bemsolv.F
434!||--- called by ------------------------------------------------------
435!|| bemsolv ../engine/source/fluid/bemsolv.F
436!|| bemsolvp ../engine/source/fluid/bemsolvp.F
437!||====================================================================
438 SUBROUTINE intgtg(NNJ, X1 , Y1, Z1, X2,
439 . Y2, Z2, X3, Y3, Z3,
440 . X0, Y0, Z0, D2, JAC,
441 . NRX, NRY, NRZ, RVAL, NN1,
442 . NN2, NN3, X )
443C-----------------------------------------------
444C I m p l i c i t T y p e s
445C-----------------------------------------------
446#include "implicit_f.inc"
447C-----------------------------------------------
448C D u m m y A r g u m e n t s
449C-----------------------------------------------
450 INTEGER NNJ, NN1, NN2, NN3
451 my_real
452 . X1, Y1, Z1, X2, Y2, Z2, X3, Y3, Z3,
453 . X0, Y0, Z0, D2, JAC, NRX, NRY, NRZ, RVAL(*),
454 . X(3,*)
455C-----------------------------------------------
456C L o c a l V a r i a b l e s
457C-----------------------------------------------
458 INTEGER NPG, IAD, IAD2, IP, NNS, IHG
459 my_real
460 . PG(50), WPG(25), XS, YS, ZS, R2,
461 . val, val1, val2, val3, w, xg, yg, zg, valphi,
462 . eta1, eta2, pg4(32), wpg4(16), jjac, ksi1, ksi2
463C
464 DATA pg /.33333333,.33333333,
465 . .33333333,.33333333,
466 . .60000000,.20000000,
467 . .20000000,.60000000,
468 . .20000000,.20000000,
469 . .33333333,.33333333,
470 . .79742699,.10128651,
471 . .10128651,.79742699,
472 . .10128651,.10128651,
473 . .05971587,.47014206,
474 . .47014206,.05971587,
475 . .47014206,.47014206,
476 . .06513010,.06513010,
477 . .86973979,.06513010,
478 . .06513010,.86973979,
479 . .31286550,.04869031,
480 . .63844419,.31286550,
481 . .04869031,.63844419,
482 . .63844419,.04869031,
483 . .31286550,.63844419,
484 . .04869031,.31286550,
485 . .26034597,.26034597,
486 . .47930807,.26034597,
487 . .26034597,.47930807,
488 . .33333333,.33333333/
489 DATA wpg /1.00000000,
490 . -.56250000,.52083333,
491 . .52083333,.52083333,
492 . .22500000,.12593918,
493 . .12593918,.12593918,
494 . .13239415,.13239415,
495 . .13239415,
496 . .05334724,.05334724,
497 . .05334724,.07711376,
498 . .07711376,.07711376,
499 . .07711376,.07711376,
500 . .07711376,.17561526,
501 . .17561526,.17561526,
502 . -.14957004/
503 DATA pg4 /-.861136311594053,-.861136311594053,
504 . -.861136311594053,-.339981034584856,
505 . -.861136311594053, .339981034584856,
506 . -.861136311594053, .861136311594053,
507 . -.339981034584856,-.861136311594053,
508 . -.339981034584856,-.339981034584856,
509 . -.339981034584856, .339981034584856,
510 . -.339981034584856, .861136311594053,
511 . .339981034584856,-.861136311594053,
512 . .339981034584856,-.339981034584856,
513 . .339981034584856, .339981034584856,
514 . .339981034584856, .861136311594053,
515 . .861136311594053,-.861136311594053,
516 . .861136311594053,-.339981034584856,
517 . .861136311594053, .339981034584856,
518 . .861136311594053, .861136311594053/
519C
520 wpg4(1) =.652145154862546 * .652145154862546
521 wpg4(2) =.652145154862546 * .347854845137454
522 wpg4(3) =.652145154862546 * .347854845137454
523 wpg4(4) =.652145154862546 * .652145154862546
524 wpg4(5) =.347854845137454 * .652145154862546
525 wpg4(6) =.347854845137454 * .347854845137454
526 wpg4(7) =.347854845137454 * .347854845137454
527 wpg4(8) =.347854845137454 * .652145154862546
528 wpg4(9) =.347854845137454 * .652145154862546
529 wpg4(10)=.347854845137454 * .347854845137454
530 wpg4(11)=.347854845137454 * .347854845137454
531 wpg4(12)=.347854845137454 * .652145154862546
532 wpg4(13)=.652145154862546 * .652145154862546
533 wpg4(14)=.652145154862546 * .347854845137454
534 wpg4(15)=.652145154862546 * .347854845137454
535 wpg4(16)=.652145154862546 * .652145154862546
536C DISTANCE A LA SOURCE
537 xs=x(1,nnj)
538 ys=x(2,nnj)
539 zs=x(3,nnj)
540 r2=(x0-xs)**2+(y0-ys)**2+(z0-zs)**2
541C NOMBRE DE POINTS DE GAUSS
542 IF (r2>hundred*d2) THEN
543 npg=1
544 iad=1
545 ELSEIF (r2>twenty5*d2) THEN
546 npg=4
547 iad=2
548 ELSEIF (r2>four*d2) THEN
549 npg=7
550 iad=6
551 ELSE
552 npg=13
553 iad=13
554 ENDIF
555C INTEGRATION
556 rval(1)=zero
557 iad2=2*(iad-1)+1
558 DO ip=1,npg
559 w=wpg(iad)
560 eta1=pg(iad2)
561 eta2=pg(iad2+1)
562 iad=iad+1
563 iad2=iad2+2
564 val1=one-eta1-eta2
565 val2=eta1
566 val3=eta2
567 xg=val1*x1+val2*x2+val3*x3
568 yg=val1*y1+val2*y2+val3*y3
569 zg=val1*z1+val2*z2+val3*z3
570 r2=(xg-xs)**2+(yg-ys)**2+(zg-zs)**2
571 IF(r2>em20)THEN
572 valphi=fourth/pi/sqrt(r2)
573 rval(1)=rval(1)+half*w*valphi*jac
574 ENDIF
575 ENDDO
576C
577 RETURN
578 END
579!||====================================================================
580!|| glbem ../engine/source/fluid/bemsolv.F
581!||--- called by ------------------------------------------------------
582!|| bemsolv ../engine/source/fluid/bemsolv.F
583!||--- calls -----------------------------------------------------
584!|| intanl ../engine/source/fluid/bemsolv.F
585!||====================================================================
586 SUBROUTINE glbem(X1 , Y1 , Z1 , X2 , Y2 ,
587 . Z2 , X3 , Y3 , Z3 , X0 ,
588 . Y0 , Z0 , TELNOR, HBEM, GBEM,
589 . TBEMTG, IBUF, X , N1 , N2 ,
590 . N3 , IELS, NN , JAC , NEL )
591C-----------------------------------------------
592C I m p l i c i t T y p e s
593C-----------------------------------------------
594#include "implicit_f.inc"
595C-----------------------------------------------
596C D u m m y A r g u m e n t s
597C-----------------------------------------------
598 INTEGER TBEMTG(3,*), IBUF(*), N1, N2, N3, IELS, NN, NEL
599 my_real
600 . X1, Y1, Z1, X2, Y2, Z2, X3, Y3, Z3, X0, Y0, Z0,
601 . telnor(3,*), x(3,*), jac, hbem(nn,*), gbem(nn,*)
602C-----------------------------------------------
603C L o c a l V a r i a b l e s
604C-----------------------------------------------
605 INTEGER NPG, IAD, JBID, NBID, IDEG, IAD2, IP, IEL, NN1, NN2, NN3
606 my_real
607 . PG(50), WPG(25), W, KSIP, ETAP, VAL1, VAL2, VAL3, XP,
608 . YP, ZP, CP, XX1, YY1, ZZ1, XX2, YY2, ZZ2, XX3, YY3, ZZ3,
609 . XX4, YY4, ZZ4, XX5, YY5, ZZ5, XX6, YY6, ZZ6, XX0, YY0,
610 . ZZ0, NRX, NRY, NRZ, D2, RVL(6), RVLH(3), RVLG
611 .
612C
613 DATA pg /.33333333,.33333333,
614 . .33333333,.33333333,
615 . .60000000,.20000000,
616 . .20000000,.60000000,
617 . .20000000,.20000000,
618 . .33333333,.33333333,
619 . .79742699,.10128651,
620 . .10128651,.79742699,
621 . .10128651,.10128651,
622 . .05971587,.47014206,
623 . .47014206,.05971587,
624 . .47014206,.47014206,
625 . .06513010,.06513010,
626 . .86973979,.06513010,
627 . .06513010,.86973979,
628 . .31286550,.04869031,
629 . .63844419,.31286550,
630 . .04869031,.63844419,
631 . .63844419,.04869031,
632 . .31286550,.63844419,
633 . .04869031,.31286550,
634 . .26034597,.26034597,
635 . .47930807,.26034597,
636 . .26034597,.47930807,
637 . .33333333,.33333333/
638 DATA wpg /1.00000000,
639 . -.56250000,.52083333,
640 . .52083333,.52083333,
641 . .22500000,.12593918,
642 . .12593918,.12593918,
643 . .13239415,.13239415,
644 . .13239415,
645 . .05334724,.05334724,
646 . .05334724,.07711376,
647 . .07711376,.07711376,
648 . .07711376,.07711376,
649 . .07711376,.17561526,
650 . .17561526,.17561526,
651 . -.14957004/
652C
653 npg=13
654 iad=13
655C
656C
657 iad2=2*(iad-1)+1
658 DO ip=1,npg
659 w=wpg(iad)
660 ksip=pg(iad2)
661 etap=pg(iad2+1)
662 iad=iad+1
663 iad2=iad2+2
664 val1=one-ksip-etap
665 val2=ksip
666 val3=etap
667 xp=val1*x1+val2*x2+val3*x3
668 yp=val1*y1+val2*y2+val3*y3
669 zp=val1*z1+val2*z2+val3*z3
670 cp=zero
671 DO iel=1,nel
672 nn1=tbemtg(1,iel)
673 nn2=tbemtg(2,iel)
674 nn3=tbemtg(3,iel)
675 xx1=x(1,ibuf(nn1))
676 yy1=x(2,ibuf(nn1))
677 zz1=x(3,ibuf(nn1))
678 xx2=x(1,ibuf(nn2))
679 yy2=x(2,ibuf(nn2))
680 zz2=x(3,ibuf(nn2))
681 xx3=x(1,ibuf(nn3))
682 yy3=x(2,ibuf(nn3))
683 zz3=x(3,ibuf(nn3))
684 nrx=telnor(1,iel)
685 nry=telnor(2,iel)
686 nrz=telnor(3,iel)
687C Matrice H
688 CALL intanl(xx1 , yy1 , zz1, xx2, yy2,
689 . zz2 , xx3 , yy3, zz3, xp ,
690 . yp , zp , nrx, nry, nrz,
691 . rvlh, rvlg)
692 hbem(n1,nn1)=hbem(n1,nn1)+w*val1*rvlh(1)*jac
693 hbem(n1,nn2)=hbem(n1,nn2)+w*val1*rvlh(2)*jac
694 hbem(n1,nn3)=hbem(n1,nn3)+w*val1*rvlh(3)*jac
695 hbem(n2,nn1)=hbem(n2,nn1)+w*val2*rvlh(1)*jac
696 hbem(n2,nn2)=hbem(n2,nn2)+w*val2*rvlh(2)*jac
697 hbem(n2,nn3)=hbem(n2,nn3)+w*val2*rvlh(3)*jac
698 hbem(n3,nn1)=hbem(n3,nn1)+w*val3*rvlh(1)*jac
699 hbem(n3,nn2)=hbem(n3,nn2)+w*val3*rvlh(2)*jac
700 hbem(n3,nn3)=hbem(n3,nn3)+w*val3*rvlh(3)*jac
701 cp=cp-rvlh(1)-rvlh(2)-rvlh(3)
702C Matrice G
703 gbem(n1,iel)=gbem(n1,iel)+w*val1*rvlg*jac
704 gbem(n2,iel)=gbem(n2,iel)+w*val2*rvlg*jac
705 gbem(n3,iel)=gbem(n3,iel)+w*val3*rvlg*jac
706 gbem(n1,nel+1)=gbem(n1,nel+1)+w*val1*rvlg*jac
707 gbem(n2,nel+1)=gbem(n2,nel+1)+w*val2*rvlg*jac
708 gbem(n3,nel+1)=gbem(n3,nel+1)+w*val3*rvlg*jac
709 ENDDO
710 hbem(n1,n1)=hbem(n1,n1)+w*cp*val1*val1*jac
711 hbem(n1,n2)=hbem(n1,n2)+w*cp*val1*val2*jac
712 hbem(n1,n3)=hbem(n1,n3)+w*cp*val1*val3*jac
713 hbem(n2,n1)=hbem(n2,n1)+w*cp*val2*val1*jac
714 hbem(n2,n2)=hbem(n2,n2)+w*cp*val2*val2*jac
715 hbem(n2,n3)=hbem(n2,n3)+w*cp*val2*val3*jac
716 hbem(n3,n1)=hbem(n3,n1)+w*cp*val3*val1*jac
717 hbem(n3,n2)=hbem(n3,n2)+w*cp*val3*val2*jac
718 hbem(n3,n3)=hbem(n3,n3)+w*cp*val3*val3*jac
719 ENDDO
720C
721 RETURN
722 END
723!||====================================================================
724!|| intanl ../engine/source/fluid/bemsolv.F
725!||--- called by ------------------------------------------------------
726!|| glbem ../engine/source/fluid/bemsolv.F
727!|| glbemp ../engine/source/fluid/bemsolvp.F
728!||====================================================================
729 SUBROUTINE intanl(X1 , Y1 , Z1 , X2 , Y2 ,
730 . Z2 , X3 , Y3 , Z3 , XS ,
731 . YS , ZS , NRX, NRY, NRZ,
732 . RVLH, RVLG)
733C-----------------------------------------------
734C I m p l i c i t T y p e s
735C-----------------------------------------------
736#include "implicit_f.inc"
737C-----------------------------------------------
738C D u m m y A r g u m e n t s
739C-----------------------------------------------
740 my_real
741 . x1, y1, z1, x2, y2, z2, x3, y3, z3, xs, ys, zs,
742 . nrx, nry, nrz, rvlh(*), rvlg
743C-----------------------------------------------
744C L o c a l V a r i a b l e s
745C-----------------------------------------------
746 INTEGER I
747 my_real
748 . vx1, vy1, vz1, vx2, vy2, vz2, s1, s12, s2, nr1, nr2,
749 . x0, y0, z0, ksi(4), eta(4), dksi(3), deta(3), r(4),
750 . xls, yls, zls, xc, yc, s(3), v, fln, arg, i1, i2, i3,
751 . ss11, ss12, ss22, delta, alpha1, beta1, gama1, alpha2,
752 . beta2, gama2, nx, ny, nz, area2, d2, l12, l22, l32,
753 . rx, ry, rz, rr, lm2
754C
755 area2=sqrt(nrx**2+nry**2+nrz**2)
756 nx=nrx/area2
757 ny=nry/area2
758 nz=nrz/area2
759C
760 x0=third*(x1+x2+x3)
761 y0=third*(y1+y2+y3)
762 z0=third*(z1+z2+z3)
763C
764C SIMPLIFICATION SI SOURCE LOIN DE L'ELEMENT
765 d2=(x0-xs)**2+(y0-ys)**2+(z0-zs)**2
766 l12=(x2-x1)**2+(y2-y1)**2+(z2-z1)**2
767 l22=(x3-x2)**2+(y3-y2)**2+(z3-z2)**2
768 l32=(x1-x3)**2+(y1-y3)**2+(z1-z3)**2
769 lm2=max(l12,l22,l32)
770C
771C COORDONNEES LOCALES
772 vx1=x2-x1
773 vy1=y2-y1
774 vz1=z2-z1
775 vx2=x3-x1
776 vy2=y3-y1
777 vz2=z3-z1
778C
779 s1=vx1*vx1+vy1*vy1+vz1*vz1
780 s12=vx1*vx2+vy1*vy2+vz1*vz2
781 nr1=sqrt(s1)
782C
783 vx2=vx2-s12/s1*vx1
784 vy2=vy2-s12/s1*vy1
785 vz2=vz2-s12/s1*vz1
786C
787 s2=vx2*vx2+vy2*vy2+vz2*vz2
788 nr2=sqrt(s2)
789 vx1=vx1/nr1
790 vy1=vy1/nr1
791 vz1=vz1/nr1
792 vx2=vx2/nr2
793 vy2=vy2/nr2
794 vz2=vz2/nr2
795C
796 xls=(xs-x0)*vx1+(ys-y0)*vy1+(zs-z0)*vz1
797 yls=(xs-x0)*vx2+(ys-y0)*vy2+(zs-z0)*vz2
798 zls=(xs-x0)*nx+(ys-y0)*ny+(zs-z0)*nz
799 ksi(1)=(x1-x0)*vx1+(y1-y0)*vy1+(z1-z0)*vz1
800 eta(1)=(x1-x0)*vx2+(y1-y0)*vy2+(z1-z0)*vz2
801 ksi(2)=(x2-x0)*vx1+(y2-y0)*vy1+(z2-z0)*vz1
802 eta(2)=(x2-x0)*vx2+(y2-y0)*vy2+(z2-z0)*vz2
803 ksi(3)=(x3-x0)*vx1+(y3-y0)*vy1+(z3-z0)*vz1
804 eta(3)=(x3-x0)*vx2+(y3-y0)*vy2+(z3-z0)*vz2
805 ksi(4)=ksi(1)
806 eta(4)=eta(1)
807 dksi(1)=ksi(2)-ksi(1)
808 dksi(2)=ksi(3)-ksi(2)
809 dksi(3)=ksi(1)-ksi(3)
810 deta(1)=eta(2)-eta(1)
811 deta(2)=eta(3)-eta(2)
812 deta(3)=eta(1)-eta(3)
813 r(1)=sqrt((xs-x1)**2+(ys-y1)**2+(zs-z1)**2)
814 r(2)=sqrt((xs-x2)**2+(ys-y2)**2+(zs-z2)**2)
815 r(3)=sqrt((xs-x3)**2+(ys-y3)**2+(zs-z3)**2)
816 s(1)=sqrt(l12)
817 s(2)=sqrt(l22)
818 s(3)=sqrt(l32)
819 r(4)=r(1)
820C
821 i1=zero
822 i2=zero
823 i3=zero
824C INTEGRALE DOUBLE COUCHE
825 IF (zls/=zero) THEN
826 DO i=1,3
827 i1=i1
828 . +atan((deta(i)*((xls-ksi(i))**2+zls**2)
829 . -dksi(i)*(xls-ksi(i))*(yls-eta(i)))/(r(i)*zls*dksi(i)))
830 . -atan((deta(i)*((xls-ksi(i+1))**2+zls**2)
831 . -dksi(i)*(xls-ksi(i+1))*(yls-eta(i+1)))/(r(i+1)*zls*dksi(i)))
832 ENDDO
833 ENDIF
834C
835C INTEGRALE SIMPLE COUCHE
836 rvlg=zero
837 DO i=1,3
838 xc=xls-ksi(i)
839 yc=yls-eta(i)
840 v=xc*deta(i)/s(i)-yc*dksi(i)/s(i)
841 arg=(r(i)+r(i+1)-s(i))/(r(i)+r(i+1)+s(i))
842 IF (arg>zero) THEN
843 fln=-log(arg)
844 rvlg=rvlg+v*fln
845 ENDIF
846 ENDDO
847 rvlg=rvlg-zls*i1
848 rvlg=-fourth/pi*rvlg
849C
850C MOMENTS D'ORDRE 1 POUR DOUBLE COUCHE
851 i2=xls*i1
852 i3=yls*i1
853 DO i=1,3
854 xc=xls-ksi(i)
855 yc=yls-eta(i)
856 arg=(r(i)+r(i+1)-s(i))/(r(i)+r(i+1)+s(i))
857 IF (arg>zero) THEN
858 fln=-log(arg)
859 i2=i2+zls*fln*deta(i)/s(i)
860 i3=i3-zls*fln*dksi(i)/s(i)
861 ENDIF
862 ENDDO
863C
864 ss11=
865 .(ksi(2)-ksi(1))*(ksi(2)-ksi(1))+(eta(2)-eta(1))*(eta(2)-eta(1))
866 ss22=
867 .(ksi(3)-ksi(1))*(ksi(3)-ksi(1))+(eta(3)-eta(1))*(eta(3)-eta(1))
868 ss12=
869 .(ksi(2)-ksi(1))*(ksi(3)-ksi(1))+(eta(2)-eta(1))*(eta(3)-eta(1))
870 delta=ss11*ss22-ss12*ss12
871 alpha1=(ss22*(ksi(2)-ksi(1))-ss12*(ksi(3)-ksi(1)))/delta
872 beta1=(ss22*(eta(2)-eta(1))-ss12*(eta(3)-eta(1)))/delta
873 gama1=(ss12*(ksi(1)*(ksi(3)-ksi(1))+eta(1)*(eta(3)-eta(1)))
874 . -ss22*(ksi(1)*(ksi(2)-ksi(1))+eta(1)*(eta(2)-eta(1))))
875 . /delta
876 alpha2=(ss11*(ksi(3)-ksi(1))-ss12*(ksi(2)-ksi(1)))/delta
877 beta2=(ss11*(eta(3)-eta(1))-ss12*(eta(2)-eta(1)))/delta
878 gama2=(ss12*(ksi(1)*(ksi(2)-ksi(1))+eta(1)*(eta(2)-eta(1)))
879 . -ss11*(ksi(1)*(ksi(3)-ksi(1))+eta(1)*(eta(3)-eta(1))))
880 . /delta
881 rvlh(1)=(1.d0-gama1-gama2)*i1-(alpha1+alpha2)*i2
882 . -(beta1+beta2)*i3
883 rvlh(2)=alpha1*i2+beta1*i3+gama1*i1
884 rvlh(3)=alpha2*i2+beta2*i3+gama2*i1
885 rvlh(1)=-rvlh(1)*fourth/pi
886 rvlh(2)=-rvlh(2)*fourth/pi
887 rvlh(3)=-rvlh(3)*fourth/pi
888C
889 RETURN
890 END
intgtg
subroutine intgtg(nnj, x1, y1, z1, x2, y2, z2, x3, y3, z3, x0, y0, z0, d2, jac, nrx, nry, nrz, rval, nn1, nn2, nn3, x)
Definition bemsolv.F:443
inthtg
subroutine inthtg(nnj, x1, y1, z1, x2, y2, z2, x3, y3, z3, x0, y0, z0, d2, jac, nrx, nry, nrz, rval, nn1, nn2, nn3, x)
Definition bemsolv.F:317
bemsolv
subroutine bemsolv(l_assemb, ilvout, iform, nn, nel, hbem, gbem, ibuf, tbemtg, telnor, x, q, phi, ipiv, phi_inf)
Definition bemsolv.F:38
intanl
subroutine intanl(x1, y1, z1, x2, y2, z2, x3, y3, z3, xs, ys, zs, nrx, nry, nrz, rvlh, rvlg)
Definition bemsolv.F:733
glsinf
subroutine glsinf(n1, n2, n3, jac, s, phi_inf)
Definition bemsolv.F:227
glbem
subroutine glbem(x1, y1, z1, x2, y2, z2, x3, y3, z3, x0, y0, z0, telnor, hbem, gbem, tbemtg, ibuf, x, n1, n2, n3, iels, nn, jac, nel)
Definition bemsolv.F:591
my_real
#define my_real
Definition cppsort.cpp:32
alpha2
#define alpha2
Definition eval.h:48
dgetrs
subroutine dgetrs(trans, n, nrhs, a, lda, ipiv, b, ldb, info)
DGETRS
Definition dgetrs.f:121
dgetrf
subroutine dgetrf(m, n, a, lda, ipiv, info)
DGETRF
Definition dgetrf.f:108
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:156
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
progbar_c
void progbar_c(int *icur, int *imax)
Definition progbar_c.c:30
arret
subroutine arret(nn)
Definition arret.F:87