OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
incpflow.F File Reference
#include "implicit_f.inc"
#include "com01_c.inc"
#include "com04_c.inc"
#include "com08_c.inc"
#include "task_c.inc"
#include "flowcom.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine incpflow (nno, nel, ninout, nni, iflow, ibuf, elem, iinout, ibufi, ibufr, ibufc, ibufl, rflow, phi, pres, u, rinout, x, v, a, npc, tf, nsensor, sensor_tab, cnp, itagel, ibufelr, ibufelc, ipiv, hbem, gbem, ibufil, cnpi, ipint, python, wfext)
subroutine trgrad (x1, x2, x3, y1, y2, y3, z1, z2, z3, phi1, phi2, phi3, grx, gry, grz)

Function/Subroutine Documentation

◆ incpflow()

subroutine incpflow ( integer nno,
integer nel,
integer ninout,
integer nni,
integer, dimension(*) iflow,
integer, dimension(*) ibuf,
integer, dimension(3,*) elem,
integer, dimension(niioflow,*) iinout,
integer, dimension(*) ibufi,
integer, dimension(*) ibufr,
integer, dimension(*) ibufc,
integer, dimension(*) ibufl,
rflow,
phi,
pres,
u,
rinout,
x,
v,
a,
integer, dimension(*) npc,
tf,
integer, intent(in) nsensor,
type (sensor_str_), dimension(nsensor) sensor_tab,
integer, dimension(*) cnp,
integer, dimension(*) itagel,
integer, dimension(*) ibufelr,
integer, dimension(*) ibufelc,
integer, dimension(*) ipiv,
hbem,
gbem,
integer, dimension(*) ibufil,
integer, dimension(*) cnpi,
integer ipint,
type(python_), intent(inout) python,
double precision, intent(inout) wfext )

Definition at line 38 of file incpflow.F.

46C-----------------------------------------------
47C M o d u l e s
48C-----------------------------------------------
49 USE sensor_mod
50 USE python_funct_mod
51C-----------------------------------------------
52C I m p l i c i t T y p e s
53C-----------------------------------------------
54#include "implicit_f.inc"
55C-----------------------------------------------
56C C o m m o n B l o c k s
57C-----------------------------------------------
58#include "com01_c.inc"
59#include "com04_c.inc"
60!#include "com06_c.inc"
61#include "com08_c.inc"
62#include "task_c.inc"
63#include "flowcom.inc"
64C-----------------------------------------------
65C D u m m y A r g u m e n t s
66C-----------------------------------------------
67 INTEGER ,INTENT(IN) :: NSENSOR
68 INTEGER NNO, NEL, NINOUT, NNI, IFLOW(*), IBUF(*), ELEM(3,*),
69 . IINOUT(NIIOFLOW,*), IBUFI(*), IBUFR(*), IBUFC(*),
70 . IBUFL(*), NPC(*),CNP(*),
71 . ITAGEL(*), IBUFELR(*), IBUFELC(*), IPIV(*), IBUFIL(*),
72 . CNPI(*), IPINT, ISMOOTH
73 my_real
74 . rflow(*), phi(*), pres(*), u(3,*), rinout(nrioflow,*),
75 . x(3,*), v(3,*), a(3,*), tf(*), hbem(*),
76 . gbem(*)
77 TYPE (SENSOR_STR_) ,DIMENSION(NSENSOR) :: SENSOR_TAB
78 TYPE(PYTHON_) , INTENT(INOUT) :: PYTHON
79 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
80C-----------------------------------------------
81C L o c a l V a r i a b l e s
82C-----------------------------------------------
83 INTEGER ILVOUT, I, II, NNO_L, III, LENBUF, PP, NNO_L2, J,
84 . IBUFL2(NNO), JJ, N1, N2, N3, NN1, NN2, NN3, NNRP, NNCP,
85 . CR_LOC_GLOB(NNO), CR_GLOB_LOC(NNO), CC_LOC_GLOB(NNO),
86 . CC_GLOB_LOC(NNO), NREIP, NCEIP, NEP, NR1, NR2, NR3, NC1,
87 . NC2, NC3, CREI_LOC_GLOB(NEL), CCEI_LOC_GLOB(NEL),
88 . CE_LOC_GLOB(NEL), CE_GLOB_LOC(NEL), NPROW, NPCOL, NBLOC,
89 . IFORM, IFREE, ISENS, IFUNC, IIO, NERP, IC, NNI_L2,
90 . IBUFIL2(NNI), ITEST, ITAGNI(NNI), NNI_L, IFT, ILT, IR, NL,
91 . IFPA, IFPRES, NELF, LELF(NEL), K, IAD, IAD2, IFVINI
92 my_real dtsub, tl, xl(3,nno+nni), vl(3,nno+nni), x1, y1, z1, x2,
93 . y2, z2, x3, y3, z3, x12, y12, z12, x13, y13, z13, x23,
94 . y23, z23, nrx, nry, nrz, norm(3,nel),q(nel), ql(nel),
95 . area, vx, vy, vz, qi(ninout), vel, tsg, tstart, dydx,
96 . sfree, isq, vfree, elarea(nel), scalt, nodarea(nno),
97 . phi1, phi2, phi3, ux, uy, uz, arean1, arean2, arean3,
98 . xmin, xmax, ymin, ymax, zmin, zmax, xx, yy, zz, st,
99 . vx1, vy1, vz1, vx2, vy2, vz2, vx3, vy3, vz3, nx1, ny1,
100 . nz1, nx2, ny2, nz2, nx3, ny3, nz3, rr1, rr2, rr3, ss1,
101 . ss2, ss3, sarea, nx, ny, nz, xc, yc ,zc, ss, phim, qm,
102 . r, fac, fac2, dpsdx, dpsdy, dpsdz, dqsdx, dqsdy, dqsdz,
103 . grx, gry, grz, tole, xm, ym, zm, phis, qs,
104 . sfpa, sfpres, rho, pa, pio, phi_old(nno+nni), uu,
105 . coef, area2, w, ksi, eta, val1, val2, val3, x0, y0, z0,
106 . r2, d2, vgx_old, vgy_old, vgz_old, vgx, vgy, vgz, agx,
107 . agy, agz, sfvini, vv, phip(nno+nni), phi_inf(nno+nni)
108 INTEGER NPG
109 my_real pg(50), wpg(25)
110 DATA pg /.33333333,.33333333,
111 . .33333333,.33333333,
112 . .60000000,.20000000,
113 . .20000000,.60000000,
114 . .20000000,.20000000,
115 . .33333333,.33333333,
116 . .79742699,.10128651,
117 . .10128651,.79742699,
118 . .10128651,.10128651,
119 . .05971587,.47014206,
120 . .47014206,.05971587,
121 . .47014206,.47014206,
122 . .06513010,.06513010,
123 . .86973979,.06513010,
124 . .06513010,.86973979,
125 . .31286550,.04869031,
126 . .63844419,.31286550,
127 . .04869031,.63844419,
128 . .63844419,.04869031,
129 . .31286550,.63844419,
130 . .04869031,.31286550,
131 . .26034597,.26034597,
132 . .47930807,.26034597,
133 . .26034597,.47930807,
134 . .33333333,.33333333/
135 DATA wpg /1.00000000,
136 . -.56250000,.52083333,
137 . .52083333,.52083333,
138 . .22500000,.12593918,
139 . .12593918,.12593918,
140 . .13239415,.13239415,
141 . .13239415,
142 . .05334724,.05334724,
143 . .05334724,.07711376,
144 . .07711376,.07711376,
145 . .07711376,.07711376,
146 . .07711376,.17561526,
147 . .17561526,.17561526,
148 . -.14957004/
149C
150 my_real, ALLOCATABLE :: sbuf(:), rbuf(:)
151 LOGICAL L_ASSEMB
152 INTEGER :: FID ! function id
153
154 my_real finter,finter_smooth
155 EXTERNAL finter,finter_smooth
156C-----------------------------------------------
157 vv = zero
158 pa = zero
159 ifree = -huge(ifree)
160 iform=iflow(13)
161 ilvout=iflow(17)
162C Subcycling
163 l_assemb=.true.
164 dtsub=rflow(3)
165 IF (dtsub>zero) THEN
166 l_assemb=.false.
167 tl=rflow(4)
168 IF (tt>tl) THEN
169 l_assemb=.true.
170 tl=max(tl+dtsub,tt)
171 rflow(4)=tl
172 ENDIF
173 ENDIF
174C---------------------------------------------------
175C 0- Preliminaires
176C---------------------------------------------------
177C coordinates and local velocities SPMD
178 IF (nspmd == 1) THEN
179 DO i=1,nno
180 ii=ibuf(i)
181 xl(1,i)=x(1,ii)
182 xl(2,i)=x(2,ii)
183 xl(3,i)=x(3,ii)
184 vl(1,i)=v(1,ii)
185 vl(2,i)=v(2,ii)
186 vl(3,i)=v(3,ii)
187 ENDDO
188 DO i=1,nni
189 ii=ibufi(i)
190 xl(1,nno+i)=x(1,ii)
191 xl(2,nno+i)=x(2,ii)
192 xl(3,nno+i)=x(3,ii)
193 vl(1,nno+i)=v(1,ii)
194 vl(2,nno+i)=v(2,ii)
195 vl(3,nno+i)=v(3,ii)
196 ENDDO
197 ELSE
198 ALLOCATE(sbuf(6*(nno+nni)), rbuf(6*(nno+nni)))
199
200 lenbuf=6*(nno+nni)
201 DO i=1,lenbuf
202 sbuf(i)=zero
203 rbuf(i)=zero
204 ENDDO
205 nno_l=iflow(16)
206
207 DO i=1,nno_l
208 ii=ibufl(i)
209 iii=ibuf(ii)
210 sbuf(6*(ii-1)+1)=x(1,iii)/cnp(ii)
211 sbuf(6*(ii-1)+2)=x(2,iii)/cnp(ii)
212 sbuf(6*(ii-1)+3)=x(3,iii)/cnp(ii)
213 sbuf(6*(ii-1)+4)=v(1,iii)/cnp(ii)
214 sbuf(6*(ii-1)+5)=v(2,iii)/cnp(ii)
215 sbuf(6*(ii-1)+6)=v(3,iii)/cnp(ii)
216 ENDDO
217 nni_l=iflow(22)
218 DO i=1,nni_l
219 ii=ibufil(i)
220 iii=ibufi(ii)
221 sbuf(6*nno+6*(ii-1)+1)=x(1,iii)/cnpi(ii)
222 sbuf(6*nno+6*(ii-1)+2)=x(2,iii)/cnpi(ii)
223 sbuf(6*nno+6*(ii-1)+3)=x(3,iii)/cnpi(ii)
224 sbuf(6*nno+6*(ii-1)+4)=v(1,iii)/cnpi(ii)
225 sbuf(6*nno+6*(ii-1)+5)=v(2,iii)/cnpi(ii)
226 sbuf(6*nno+6*(ii-1)+6)=v(3,iii)/cnpi(ii)
227 ENDDO
228 CALL spmd_fl_sum(sbuf, lenbuf, rbuf)
229 DO i=1,nno+nni
230 xl(1,i)=rbuf(6*(i-1)+1)
231 xl(2,i)=rbuf(6*(i-1)+2)
232 xl(3,i)=rbuf(6*(i-1)+3)
233 vl(1,i)=rbuf(6*(i-1)+4)
234 vl(2,i)=rbuf(6*(i-1)+5)
235 vl(3,i)=rbuf(6*(i-1)+6)
236 ENDDO
237C
238 DEALLOCATE(sbuf, rbuf)
239 ENDIF
240C outgoing unitary normal vectors
241 DO i=1,nno
242 nodarea(i)=zero
243 ENDDO
244C
245 DO i=1,nel
246 n1=elem(1,i)
247 n2=elem(2,i)
248 n3=elem(3,i)
249 x1=xl(1,n1)
250 x2=xl(1,n2)
251 x3=xl(1,n3)
252 y1=xl(2,n1)
253 y2=xl(2,n2)
254 y3=xl(2,n3)
255 z1=xl(3,n1)
256 z2=xl(3,n2)
257 z3=xl(3,n3)
258 x12=x2-x1
259 y12=y2-y1
260 z12=z2-z1
261 x13=x3-x1
262 y13=y3-y1
263 z13=z3-z1
264 x23=x3-x2
265 y23=y3-y2
266 z23=z3-z2
267 nrx=y12*z13-z12*y13
268 nry=z12*x13-x12*z13
269 nrz=x12*y13-y12*x13
270 norm(1,i)=nrx
271 norm(2,i)=nry
272 norm(3,i)=nrz
273 area=sqrt(nrx**2+nry**2+nrz**2)
274 elarea(i)=half*area
275 nodarea(n1)=nodarea(n1)+one_over_6*area
276 nodarea(n2)=nodarea(n2)+one_over_6*area
277 nodarea(n3)=nodarea(n3)+one_over_6*area
278 ENDDO
279
280 DO i=1,nno+nni
281 phi_old(i)=phi(i)
282 ENDDO
283C velocity in far field
284 ifvini=iflow(24)
285 vgx_old=rflow(12)*rflow(9)
286 vgy_old=rflow(12)*rflow(10)
287 vgz_old=rflow(12)*rflow(11)
288 vgx=zero
289 vgy=zero
290 vgz=zero
291 agx=zero
292 agy=zero
293 agz=zero
294 IF (ifvini>0) THEN
295 tstart=zero
296 sfvini=rflow(7)
297 scalt=rflow(8)
298 ismooth = 0
299 IF (ifvini > 0) ismooth = npc(2*nfunct+ifvini+1)
300 IF (tt>tstart.AND.dt1>zero) THEN
301 tsg=(tt-tstart)*scalt
302 IF (ismooth == 0) THEN
303 vv=sfvini*finter(ifvini, tsg, npc, tf, dydx)
304 ELSE IF(ismooth > 0) THEN
305 vv=sfvini*finter_smooth(ifvini, tsg, npc, tf, dydx)
306 ELSE
307 fid = -ismooth ! the id the python function is saved in the position of ISMOOTH in the NPC array
308 CALL python_call_funct1d(python, fid,tsg,vv)
309 vv = vv * sfvini
310 ENDIF
311 ENDIF
312 rflow(12)=vv
313 vgx=rflow(9)*vv
314 vgy=rflow(10)*vv
315 vgz=rflow(11)*vv
316 IF (dt1>zero) THEN
317 agx=(vgx-vgx_old)/dt1
318 agy=(vgy-vgy_old)/dt1
319 agz=(vgz-vgz_old)/dt1
320 ENDIF
321C potentiel at far field
322 DO i=1,nno+nni
323 xx=xl(1,i)
324 yy=xl(2,i)
325 zz=xl(3,i)
326 phi_inf(i)=vgx*(xx-xl(1,nno)) + vgy*(yy-xl(2,nno)) + vgz*(zz-xl(3,nno))
327 ENDDO
328 ELSE
329 DO i=1,nno+nni
330 phi_inf(i)=zero
331 ENDDO
332 ENDIF
333C---------------------------------------------------
334C 1- Surface flux
335C---------------------------------------------------
336 DO i=1,nel
337 n1=elem(1,i)
338 n2=elem(2,i)
339 n3=elem(3,i)
340 nrx=norm(1,i)
341 nry=norm(2,i)
342 nrz=norm(3,i)
343 area=elarea(i)
344 IF(area>em20)THEN
345 area2=two*area
346 nrx=nrx/area2
347 nry=nry/area2
348 nrz=nrz/area2
349 vx=third*(vl(1,n1)+vl(1,n2)+vl(1,n3))
350 vy=third*(vl(2,n1)+vl(2,n2)+vl(2,n3))
351 vz=third*(vl(3,n1)+vl(3,n2)+vl(3,n3))
352 q(i)=vx*nrx+vy*nry+vz*nrz
353 ELSE
354 q(i)=zero
355 ENDIF
356
357 ENDDO
358
359 DO i=1,ninout
360 qi(i)=zero
361 IF (iinout(3,i)==0) THEN
362 ifree=i
363 ELSE
364 isens=iinout(2,i)
365 ifunc=iinout(3,i)
366 vel=rinout(1,i)
367 scalt=rinout(3,i)
368 ismooth = 0
369 IF (ifunc > 0) ismooth = npc(2*nfunct+ifunc+1)
370 IF (isens==0) THEN
371 tstart=zero
372 ELSE
373 tstart=sensor_tab(isens)%TSTART
374 ENDIF
375 IF (tt>=tstart.AND.dt1>zero) THEN
376 tsg=(tt-tstart)*scalt
377 IF (ismooth == 0) THEN
378 qi(i)=vel*finter(ifunc, tsg, npc, tf, dydx)
379 ELSE IF(ismooth > 0 ) THEN
380 qi(i)=vel*finter_smooth(ifunc, tsg, npc, tf, dydx)
381 ELSE
382 fid = -ismooth ! the id the python function is saved in the position of ISMOOTH in the NPC array
383 CALL python_call_funct1d(python, fid,tsg,qi(i))
384 qi(i) = qi(i) * vel
385 ENDIF
386 ENDIF
387 ENDIF
388 ENDDO
389
390 IF (ninout>zero) THEN
391 sfree=zero
392 isq=zero
393 DO i=1,nel
394 iio=itagel(i)
395 IF (iio/=0) THEN
396 IF (iio==ifree) THEN
397 sfree=sfree+elarea(i)
398 ELSE
399 q(i)=q(i)+qi(iio)
400 ENDIF
401 ENDIF
402 isq=isq+q(i)*elarea(i)
403 ENDDO
404
405 IF(sfree/=zero)THEN
406 vfree=-isq/sfree
407 ELSE
408 vfree=zero
409 ENDIF
410 qi(ifree)=vfree
411 DO i=1,nel
412 iio=itagel(i)
413 IF (iio==ifree) q(i)=q(i)+vfree
414 ENDDO
415 ENDIF
416
417 IF (nspmd > 1) THEN
418 DO i=1,nel
419 ql(i)=zero
420 ENDDO
421 nerp=0
422 ic=ibufelc(1)
423 IF (ic>0) THEN
424 DO i=1,nel
425 IF (ibufelr(i)>0) THEN
426 nerp=nerp+1
427 ql(nerp)=q(i)
428 ENDIF
429 ENDDO
430 ENDIF
431 ENDIF
432C---------------------------------------------------
433C 2- BEM
434C---------------------------------------------------
435 IF (nspmd == 1) THEN
436 IF (dt1>zero)
437 . CALL bemsolv(l_assemb, ilvout, iform, nno, nel,
438 . hbem, gbem, ibuf, elem, norm,
439 . x, q, phi, ipiv, phi_inf)
440 ELSE
441C MApping Tables
442 nnrp=0
443 nncp=0
444 DO i=1,nno
445 cr_loc_glob(i)=0
446 cr_glob_loc(i)=0
447 cc_loc_glob(i)=0
448 cc_glob_loc(i)=0
449 ENDDO
450 DO i=1,nel
451 crei_loc_glob(i)=0
452 ccei_loc_glob(i)=0
453 ce_loc_glob(i)=0
454 ce_glob_loc(i)=0
455 ENDDO
456
457 DO i=1,nno
458 IF (ibufr(i)>0) THEN
459 nnrp=nnrp+1
460 cr_loc_glob(nnrp)=i
461 cr_glob_loc(i)=nnrp
462 ENDIF
463 IF (ibufc(i)>0) THEN
464 nncp=nncp+1
465 cc_loc_glob(nncp)=i
466 cc_glob_loc(i)=nncp
467 ENDIF
468 ENDDO
469 nreip=0
470 nceip=0
471 nep=0
472 DO i=1,nel
473 n1=elem(1,i)
474 n2=elem(2,i)
475 n3=elem(3,i)
476 nr1=ibufr(n1)
477 nr2=ibufr(n2)
478 nr3=ibufr(n3)
479 nc1=ibufc(n1)
480 nc2=ibufc(n2)
481 nc3=ibufc(n3)
482 IF (nr1/=0.OR.nr2/=0.OR.nr3/=0) THEN
483 nreip=nreip+1
484 crei_loc_glob(nreip)=i
485 ENDIF
486 IF (nc1/=0.OR.nc2/=0.OR.nc3/=0) THEN
487 nceip=nceip+1
488 ccei_loc_glob(nceip)=i
489 ENDIF
490 IF (ibufelc(i)>0) THEN
491 nep=nep+1
492 ce_loc_glob(nep)=i
493 ce_glob_loc(i)=nep
494 ENDIF
495 ENDDO
496
497 nprow=iflow(18)
498 npcol=iflow(19)
499 nbloc=iflow(12)
500 IF (dt1>zero)
501 . CALL bemsolvp(
502 . elem, norm, xl, nnrp, cr_loc_glob,
503 . cr_glob_loc, nncp, cc_loc_glob, cc_glob_loc, nceip,
504 . ccei_loc_glob, nep, ce_loc_glob, ce_glob_loc, nreip,
505 . crei_loc_glob, nno, nel, hbem, gbem,
506 . ql, phi, iform, l_assemb, nprow,
507 . npcol, nbloc, ipiv, nerp, ilvout,
508 . phi_inf )
509 ENDIF
510C---------------------------------------------------
511C 3- Velocities
512C---------------------------------------------------
513 IF (nspmd == 1) THEN
514 phi(nno)=zero
515 ! on surface
516 DO i=1,nno
517 u(1,i)=zero
518 u(2,i)=zero
519 u(3,i)=zero
520 ENDDO
521
522 DO i=1,nel
523 n1=elem(1,i)
524 n2=elem(2,i)
525 n3=elem(3,i)
526 x1=xl(1,n1)
527 y1=xl(2,n1)
528 z1=xl(3,n1)
529 x2=xl(1,n2)
530 y2=xl(2,n2)
531 z2=xl(3,n2)
532 x3=xl(1,n3)
533 y3=xl(2,n3)
534 z3=xl(3,n3)
535 phi1=phi(n1)
536 phi2=phi(n2)
537 phi3=phi(n3)
538 CALL trgrad(x1, x2, x3, y1, y2,
539 . y3, z1, z2, z3, phi1,
540 . phi2, phi3, grx, gry, grz )
541 nrx=norm(1,i)
542 nry=norm(2,i)
543 nrz=norm(3,i)
544 area=elarea(i)
545 area2=two*area
546 ux=grx+q(i)*nrx/area2
547 uy=gry+q(i)*nry/area2
548 uz=grz+q(i)*nrz/area2
549 arean1=nodarea(n1)
550 arean2=nodarea(n2)
551 arean3=nodarea(n3)
552 u(1,n1)=u(1,n1)+ux*third*area/arean1
553 u(2,n1)=u(2,n1)+uy*third*area/arean1
554 u(3,n1)=u(3,n1)+uz*third*area/arean1
555 u(1,n2)=u(1,n2)+ux*third*area/arean2
556 u(2,n2)=u(2,n2)+uy*third*area/arean2
557 u(3,n2)=u(3,n2)+uz*third*area/arean2
558 u(1,n3)=u(1,n3)+ux*third*area/arean3
559 u(2,n3)=u(2,n3)+uy*third*area/arean3
560 u(3,n3)=u(3,n3)+uz*third*area/arean3
561 ENDDO
562C velocities at auxiliary nodes
563 IF (nni>0.AND.ipint==1) THEN
564 itest=iflow(21)
565 IF (itest>0) THEN
566 tole=rflow(6)
567 xmin=ep30
568 xmax=-ep30
569 ymin=ep30
570 ymax=-ep30
571 zmin=ep30
572 zmax=-ep30
573 DO i=1,nno
574 xx=xl(1,i)
575 yy=xl(2,i)
576 zz=xl(3,i)
577 xmin=min(xmin,xx)
578 xmax=max(xmax,xx)
579 ymin=min(ymin,yy)
580 ymax=max(ymax,yy)
581 zmin=min(zmin,zz)
582 zmax=max(zmax,zz)
583 ENDDO
584
585 DO i=1,nni
586 itagni(i)=0
587 xx=xl(1,nno+i)
588 yy=xl(2,nno+i)
589 zz=xl(3,nno+i)
590 IF (xx>xmax.OR.xx<xmin.OR.
591 . yy>ymax.OR.yy<ymin.OR.
592 . zz>zmax.OR.zz<zmin) cycle
593 st=zero
594 DO j=1,nel
595 n1=elem(1,j)
596 n2=elem(2,j)
597 n3=elem(3,j)
598 x1=xl(1,n1)
599 y1=xl(2,n1)
600 z1=xl(3,n1)
601 x2=xl(1,n2)
602 y2=xl(2,n2)
603 z2=xl(3,n2)
604 x3=xl(1,n3)
605 y3=xl(2,n3)
606 z3=xl(3,n3)
607 vx1=x1-xx
608 vy1=y1-yy
609 vz1=z1-zz
610 vx2=x2-xx
611 vy2=y2-yy
612 vz2=z2-zz
613 vx3=x3-xx
614 vy3=y3-yy
615 vz3=z3-zz
616
617 nx1=vy1*vz2-vz1*vy2
618 ny1=vz1*vx2-vx1*vz2
619 nz1=vx1*vy2-vy1*vx2
620 nx2=vy2*vz3-vz2*vy3
621 ny2=vz2*vx3-vx2*vz3
622 nz2=vx2*vy3-vy2*vx3
623 nx3=vy3*vz1-vz3*vy1
624 ny3=vz3*vx1-vx3*vz1
625 nz3=vx3*vy1-vy3*vx1
626 rr1=sqrt(nx1**2+ny1**2+nz1**2)
627 rr2=sqrt(nx2**2+ny2**2+nz2**2)
628 rr3=sqrt(nx3**2+ny3**2+nz3**2)
629 ss1=-(nx1*nx2+ny1*ny2+nz1*nz2)/rr1/rr2
630 ss2=-(nx2*nx3+ny2*ny3+nz2*nz3)/rr2/rr3
631 ss3=-(nx3*nx1+ny3*ny1+nz3*nz1)/rr3/rr1
632 IF (ss1>one) ss1=one
633 IF (ss1<-one) ss1=-one
634 IF (ss2>one) ss2=one
635 IF (ss2<-one) ss2=-one
636 IF (ss3>one) ss3=one
637 IF (ss3<-one) ss3=-one
638 sarea=acos(ss1)+acos(ss2)+acos(ss3)-pi
639
640 vx1=x2-x1
641 vy1=y2-y1
642 vz1=z2-z1
643 vx2=x3-x1
644 vy2=y3-y1
645 vz2=z3-z1
646 nx=vy1*vz2-vz1*vy2
647 ny=vz1*vx2-vx1*vz2
648 nz=vx1*vy2-vy1*vx2
649 xc=third*(x1+x2+x3)
650 yc=third*(y1+y2+y3)
651 zc=third*(z1+z2+z3)
652 ss=nx*(xc-xx)+ny*(yc-yy)+nz*(zc-zz)
653 IF (ss<zero) sarea=-sarea
654 st=st+sarea
655 ENDDO
656 IF (itest==1.AND.abs(st)>tole) itagni(i)=1
657 IF (itest==2.AND.abs(st)<=tole) itagni(i)=1
658 ENDDO
659 ELSE
660 DO i=1,nni
661 itagni(i)=1
662 ENDDO
663 ENDIF
664
665 DO i=1,nni
666 phi(nno+i)=zero
667 u(1,nno+i)=zero
668 u(2,nno+i)=zero
669 u(3,nno+i)=zero
670 IF (itagni(i)==0) cycle
671 xx=xl(1,nno+i)
672 yy=xl(2,nno+i)
673 zz=xl(3,nno+i)
674 DO j=1,nel
675 n1=elem(1,j)
676 n2=elem(2,j)
677 n3=elem(3,j)
678 x1=xl(1,n1)
679 y1=xl(2,n1)
680 z1=xl(3,n1)
681 x2=xl(1,n2)
682 y2=xl(2,n2)
683 z2=xl(3,n2)
684 x3=xl(1,n3)
685 y3=xl(2,n3)
686 z3=xl(3,n3)
687 x0=third*(x1+x2+x3)
688 y0=third*(y1+y2+y3)
689 z0=third*(z1+z2+z3)
690 r2=(x0-xx)**2+(y0-yy)**2+(z0-zz)**2
691 d2=min((x0-x1)**2+(y0-y1)**2+(z0-z1)**2,
692 . (x0-x2)**2+(y0-y2)**2+(z0-z2)**2,
693 . (x0-x3)**2+(y0-y3)**2+(z0-z3)**2)
694C number of Gauss point on the triangle
695 IF (r2>hundred*d2) THEN
696 npg=1
697 iad=1
698 ELSEIF (r2>twenty5*d2) THEN
699 npg=4
700 iad=2
701 ELSEIF (r2>four*d2) THEN
702 npg=7
703 iad=6
704 ELSE
705 npg=13
706 iad=13
707 ENDIF
708
709 qm=q(j)
710 area=elarea(j)
711 area2=two*area
712 nrx=norm(1,j)/area2
713 nry=norm(2,j)/area2
714 nrz=norm(3,j)/area2
715 iad2=2*(iad-1)+1
716 DO k=1,npg
717 w=wpg(iad)
718 ksi=pg(iad2)
719 eta=pg(iad2+1)
720 iad=iad+1
721 iad2=iad2+2
722 val1=one-ksi-eta
723 val2=ksi
724 val3=eta
725 xm=x1*val1+x2*val2+x3*val3
726 ym=y1*val1+y2*val2+y3*val3
727 zm=z1*val1+z2*val2+z3*val3
728 phim=phi(n1)*val1+phi(n2)*val2+phi(n3)*val3
729 r=sqrt((xx-xm)**2+(yy-ym)**2+(zz-zm)**2)
730 fac=fourth/pi/r**3
731 fac2=three/r**2* ((xx-xm)*nrx+(yy-ym)*nry+(zz-zm)*nrz)
732
733 phis=fourth/pi/r
734 qs=fac*((xx-xm)*nrx+(yy-ym)*nry+(zz-zm)*nrz)
735 dpsdx=-fac*(xx-xm)
736 dpsdy=-fac*(yy-ym)
737 dpsdz=-fac*(zz-zm)
738 dqsdx=fac*(nrx-(xx-xm)*fac2)
739 dqsdy=fac*(nry-(yy-ym)*fac2)
740 dqsdz=fac*(nrz-(zz-zm)*fac2)
741
742 phi(nno+i)=phi(nno+i)+area*w*(qm*phis-phim*qs)
743 u(1,nno+i)=u(1,nno+i)+area*w*(qm*dpsdx-phim*dqsdx)
744 u(2,nno+i)=u(2,nno+i)+area*w*(qm*dpsdy-phim*dqsdy)
745 u(3,nno+i)=u(3,nno+i)+area*w*(qm*dpsdz-phim*dqsdz)
746 ENDDO
747 ENDDO
748 phi(nno+i)=phi(nno+i)+phi_inf(nno+i)
749 u(1,nno+i)=u(1,nno+i)+vgx
750 u(2,nno+i)=u(2,nno+i)+vgy
751 u(3,nno+i)=u(3,nno+i)+vgz
752 ENDDO
753 ENDIF
754 ELSE
755C exchange of POtentials
756 ALLOCATE(sbuf(nno), rbuf(nno))
757 DO i=1,nno
758 sbuf(i)=zero
759 rbuf(i)=zero
760 ENDDO
761 ic=ibufc(1)
762 IF (ic>0) THEN
763 ii=0
764 DO i=1,nno
765 ir=ibufr(i)
766 IF (ir>0) THEN
767 ii=ii+1
768 sbuf(i)=phi(ii)
769 ENDIF
770 ENDDO
771 ENDIF
772 lenbuf=nno
773 CALL spmd_fl_sum(sbuf, lenbuf, rbuf)
774 DO i=1,nno
775 phi(i)=rbuf(i)
776 ENDDO
777 DEALLOCATE(sbuf, rbuf)
778 phi(nno)=zero
779 !velocity on surface
780 DO i=1,nno
781 u(1,i)=zero
782 u(2,i)=zero
783 u(3,i)=zero
784 ENDDO
785
786 DO i=1,nel
787 n1=elem(1,i)
788 n2=elem(2,i)
789 n3=elem(3,i)
790 x1=xl(1,n1)
791 y1=xl(2,n1)
792 z1=xl(3,n1)
793 x2=xl(1,n2)
794 y2=xl(2,n2)
795 z2=xl(3,n2)
796 x3=xl(1,n3)
797 y3=xl(2,n3)
798 z3=xl(3,n3)
799 phi1=phi(n1)
800 phi2=phi(n2)
801 phi3=phi(n3)
802 CALL trgrad(x1, x2, x3, y1, y2,
803 . y3, z1, z2, z3, phi1,
804 . phi2, phi3, grx, gry, grz )
805 nrx=norm(1,i)
806 nry=norm(2,i)
807 nrz=norm(3,i)
808 area=elarea(i)
809 IF(area>em20)THEN
810 area2=two*area
811 ux=grx+q(i)*nrx/area2
812 uy=gry+q(i)*nry/area2
813 uz=grz+q(i)*nrz/area2
814 arean1=nodarea(n1)
815 arean2=nodarea(n2)
816 arean3=nodarea(n3)
817 u(1,n1)=u(1,n1)+ux*third*area/arean1
818 u(2,n1)=u(2,n1)+uy*third*area/arean1
819 u(3,n1)=u(3,n1)+uz*third*area/arean1
820 u(1,n2)=u(1,n2)+ux*third*area/arean2
821 u(2,n2)=u(2,n2)+uy*third*area/arean2
822 u(3,n2)=u(3,n2)+uz*third*area/arean2
823 u(1,n3)=u(1,n3)+ux*third*area/arean3
824 u(2,n3)=u(2,n3)+uy*third*area/arean3
825 u(3,n3)=u(3,n3)+uz*third*area/arean3
826 ENDIF
827
828 ENDDO
829C auxiliary nodes are dispatched on nodes
830 IF (nni>0.AND.ipint==1) THEN
831 nl=nni/nspmd+1
832 ift=ispmd*nl+1
833 ilt=min(ift+nl-1,nni)
834 itest=iflow(21)
835 IF (itest>0) THEN
836 tole=rflow(6)
837 xmin=ep30
838 xmax=-ep30
839 ymin=ep30
840 ymax=-ep30
841 zmin=ep30
842 zmax=-ep30
843 DO i=1,nno
844 xx=xl(1,i)
845 yy=xl(2,i)
846 zz=xl(3,i)
847 xmin=min(xmin,xx)
848 xmax=max(xmax,xx)
849 ymin=min(ymin,yy)
850 ymax=max(ymax,yy)
851 zmin=min(zmin,zz)
852 zmax=max(zmax,zz)
853 ENDDO
854
855 DO i=ift,ilt
856 itagni(i)=0
857 xx=xl(1,nno+i)
858 yy=xl(2,nno+i)
859 zz=xl(3,nno+i)
860 IF (xx>xmax.OR.xx<xmin.OR.
861 . yy>ymax.OR.yy<ymin.OR.
862 . zz>zmax.OR.zz<zmin) cycle
863 st=zero
864 DO j=1,nel
865 n1=elem(1,j)
866 n2=elem(2,j)
867 n3=elem(3,j)
868 x1=xl(1,n1)
869 y1=xl(2,n1)
870 z1=xl(3,n1)
871 x2=xl(1,n2)
872 y2=xl(2,n2)
873 z2=xl(3,n2)
874 x3=xl(1,n3)
875 y3=xl(2,n3)
876 z3=xl(3,n3)
877 vx1=x1-xx
878 vy1=y1-yy
879 vz1=z1-zz
880 vx2=x2-xx
881 vy2=y2-yy
882 vz2=z2-zz
883 vx3=x3-xx
884 vy3=y3-yy
885 vz3=z3-zz
886
887 nx1=vy1*vz2-vz1*vy2
888 ny1=vz1*vx2-vx1*vz2
889 nz1=vx1*vy2-vy1*vx2
890 nx2=vy2*vz3-vz2*vy3
891 ny2=vz2*vx3-vx2*vz3
892 nz2=vx2*vy3-vy2*vx3
893 nx3=vy3*vz1-vz3*vy1
894 ny3=vz3*vx1-vx3*vz1
895 nz3=vx3*vy1-vy3*vx1
896 rr1=sqrt(nx1**2+ny1**2+nz1**2)
897 rr2=sqrt(nx2**2+ny2**2+nz2**2)
898 rr3=sqrt(nx3**2+ny3**2+nz3**2)
899 ss1=-(nx1*nx2+ny1*ny2+nz1*nz2)/rr1/rr2
900 ss2=-(nx2*nx3+ny2*ny3+nz2*nz3)/rr2/rr3
901 ss3=-(nx3*nx1+ny3*ny1+nz3*nz1)/rr3/rr1
902 IF (ss1>one) ss1=one
903 IF (ss1<-one) ss1=-one
904 IF (ss2>one) ss2=one
905 IF (ss2<-one) ss2=-one
906 IF (ss3>one) ss3=one
907 IF (ss3<-one) ss3=-one
908 sarea=acos(ss1)+acos(ss2)+acos(ss3)-pi
909
910 vx1=x2-x1
911 vy1=y2-y1
912 vz1=z2-z1
913 vx2=x3-x1
914 vy2=y3-y1
915 vz2=z3-z1
916 nx=vy1*vz2-vz1*vy2
917 ny=vz1*vx2-vx1*vz2
918 nz=vx1*vy2-vy1*vx2
919 xc=third*(x1+x2+x3)
920 yc=third*(y1+y2+y3)
921 zc=third*(z1+z2+z3)
922 ss=nx*(xc-xx)+ny*(yc-yy)+nz*(zc-zz)
923 IF (ss<zero) sarea=-sarea
924 st=st+sarea
925 ENDDO
926 IF (itest==1.AND.abs(st)>tole) itagni(i)=1
927 IF (itest==2.AND.abs(st)<=tole) itagni(i)=1
928 ENDDO
929 ELSE
930 DO i=ift,ilt
931 itagni(i)=1
932 ENDDO
933 ENDIF
934
935 DO i=ift,ilt
936 phi(nno+i)=zero
937 u(1,nno+i)=zero
938 u(2,nno+i)=zero
939 u(3,nno+i)=zero
940 IF (itagni(i)==0) cycle
941 xx=xl(1,nno+i)
942 yy=xl(2,nno+i)
943 zz=xl(3,nno+i)
944 DO j=1,nel
945 n1=elem(1,j)
946 n2=elem(2,j)
947 n3=elem(3,j)
948 x1=xl(1,n1)
949 y1=xl(2,n1)
950 z1=xl(3,n1)
951 x2=xl(1,n2)
952 y2=xl(2,n2)
953 z2=xl(3,n2)
954 x3=xl(1,n3)
955 y3=xl(2,n3)
956 z3=xl(3,n3)
957 x0=third*(x1+x2+x3)
958 y0=third*(y1+y2+y3)
959 z0=third*(z1+z2+z3)
960 r2=(x0-xx)**2+(y0-yy)**2+(z0-zz)**2
961 d2=min((x0-x1)**2+(y0-y1)**2+(z0-z1)**2,
962 . (x0-x2)**2+(y0-y2)**2+(z0-z2)**2,
963 . (x0-x3)**2+(y0-y3)**2+(z0-z3)**2)
964C number of gauss point on the triangle
965 IF (r2>hundred*d2) THEN
966 npg=1
967 iad=1
968 ELSEIF (r2>twenty5*d2) THEN
969 npg=4
970 iad=2
971 ELSEIF (r2>four*d2) THEN
972 npg=7
973 iad=6
974 ELSE
975 npg=13
976 iad=13
977 ENDIF
978
979 qm=q(j)
980 area=elarea(j)
981 IF(area>em20)THEN
982 area2=two*area
983 nrx=norm(1,j)/area2
984 nry=norm(2,j)/area2
985 nrz=norm(3,j)/area2
986 ELSE
987 nrx=zero
988 nry=zero
989 nrz=zero
990 ENDIF
991 iad2=2*(iad-1)+1
992 DO k=1,npg
993 w=wpg(iad)
994 ksi=pg(iad2)
995 eta=pg(iad2+1)
996 iad=iad+1
997 iad2=iad2+2
998 val1=one-ksi-eta
999 val2=ksi
1000 val3=eta
1001 xm=x1*val1+x2*val2+x3*val3
1002 ym=y1*val1+y2*val2+y3*val3
1003 zm=z1*val1+z2*val2+z3*val3
1004 phim=phi(n1)*val1+phi(n2)*val2+phi(n3)*val3
1005 r=sqrt((xx-xm)**2+(yy-ym)**2+(zz-zm)**2)
1006 IF(r>em20)THEN
1007 fac=fourth/pi/r**3
1008 fac2=three/r**2* ((xx-xm)*nrx+(yy-ym)*nry+(zz-zm)*nrz)
1009 phis=fourth/pi/r
1010 ELSE
1011 fac = zero
1012 fac2 = zero
1013 phis = zero
1014 ENDIF
1015 qs=fac*((xx-xm)*nrx+(yy-ym)*nry+(zz-zm)*nrz)
1016 dpsdx=-fac*(xx-xm)
1017 dpsdy=-fac*(yy-ym)
1018 dpsdz=-fac*(zz-zm)
1019 dqsdx=fac*(nrx-(xx-xm)*fac2)
1020 dqsdy=fac*(nry-(yy-ym)*fac2)
1021 dqsdz=fac*(nrz-(zz-zm)*fac2)
1022
1023 phi(nno+i)=phi(nno+i)+area*w*(qm*phis-phim*qs)
1024 u(1,nno+i)=u(1,nno+i)+area*w*(qm*dpsdx-phim*dqsdx)
1025 u(2,nno+i)=u(2,nno+i)+area*w*(qm*dpsdy-phim*dqsdy)
1026 u(3,nno+i)=u(3,nno+i)+area*w*(qm*dpsdz-phim*dqsdz)
1027 ENDDO
1028 ENDDO
1029 phi(nno+i)=phi(nno+i)+phi_inf(nno+i)
1030 u(1,nno+i)=u(1,nno+i)+vgx
1031 u(2,nno+i)=u(2,nno+i)+vgy
1032 u(3,nno+i)=u(3,nno+i)+vgz
1033 ENDDO
1034C Exchange of potentials and velocities on internal nodes
1035 ALLOCATE(sbuf(4*nni), rbuf(4*nni))
1036 lenbuf=4*nni
1037 DO i=1,lenbuf
1038 sbuf(i)=zero
1039 rbuf(i)=zero
1040 ENDDO
1041 DO i=ift,ilt
1042 sbuf(4*(i-1)+1)=phi(nno+i)
1043 sbuf(4*(i-1)+2)=u(1,nno+i)
1044 sbuf(4*(i-1)+3)=u(2,nno+i)
1045 sbuf(4*(i-1)+4)=u(3,nno+i)
1046 ENDDO
1047 CALL spmd_fl_sum(sbuf, lenbuf, rbuf)
1048 DO i=1,nni
1049 phi(nno+i)=rbuf(4*(i-1)+1)
1050 u(1,nno+i)=rbuf(4*(i-1)+2)
1051 u(2,nno+i)=rbuf(4*(i-1)+3)
1052 u(3,nno+i)=rbuf(4*(i-1)+4)
1053 ENDDO
1054 ENDIF
1055 ENDIF
1056C---------------------------------------------------
1057C 4- Pressure
1058C---------------------------------------------------
1059 rho=rflow(5)
1060 ifpa=iflow(23)
1061 IF (ifpa/=0) THEN
1062C imposed stagnation pressure
1063 sfpa=rflow(1)
1064 scalt=rflow(2)
1065 tstart=zero
1066 ismooth = 0
1067 IF (ifpa > 0) ismooth = npc(2*nfunct+ifpa+1)
1068 IF (tt>=tstart.AND.dt1>zero) THEN
1069 tsg=(tt-tstart)*scalt
1070 IF (ismooth == 0) THEN
1071 pa=sfpa*finter(ifpa, tsg, npc, tf, dydx)
1072 ELSE IF(ismooth > 0) THEN
1073 pa=sfpa*finter_smooth(ifpa, tsg, npc, tf, dydx)
1074 ELSE
1075 fid = -ismooth ! the id the python function is saved in the position of ISMOOTH in the NPC array
1076 CALL python_call_funct1d(python, fid,tsg,pa)
1077 pa = pa * sfpa
1078 ENDIF ! IF (ISMOOTH == 0)
1079 ENDIF
1080 ELSE
1081C pressure imposed at inlet/outlet
1082 DO i=1,ninout
1083 IF (iinout(4,i)/=0) THEN
1084 ifpres=iinout(4,i)
1085 sfpres=rinout(2,i)
1086 scalt=rinout(3,i)
1087 tstart=0
1088 pio=zero
1089 ismooth = 0
1090 IF (ifpres > 0) ismooth = npc(2*nfunct+ifpres+1)
1091 IF (tt>tstart.AND.dt1>zero) THEN
1092 tsg=(tt-tstart)*scalt
1093 IF (ismooth == 0) THEN
1094 pio=sfpres*finter(ifpres, tsg, npc, tf, dydx)
1095 ELSEIF(ismooth > 0) THEN
1096 pio=sfpres*finter_smooth(ifpres, tsg, npc, tf, dydx)
1097 ELSE
1098 fid = -ismooth ! the id the python function is saved in the position of ISMOOTH in the NPC array
1099 CALL python_call_funct1d(python, fid,tsg,pio)
1100 pio = pio * sfpres
1101 ENDIF ! IF (ISMOOTH == 0)
1102 ENDIF
1103 pa=pio+half*rho*qi(i)**2
1104 ENDIF
1105 ENDDO
1106 ENDIF
1107
1108 IF (dt1>zero) THEN
1109 IF (ncycle>1) THEN
1110 DO i=1,nno+nni
1111 xx=xl(1,i)
1112 yy=xl(2,i)
1113 zz=xl(3,i)
1114 phip(i)=(phi(i)-phi_old(i))/dt1
1115 . +agx*(xx-xl(1,nno))+agy*(yy-xl(2,nno))
1116 . +agz*(zz-xl(3,nno))
1117 ENDDO
1118 ELSE
1119 DO i=1,nno+nni
1120 phip(i)=zero
1121 ENDDO
1122 ENDIF
1123
1124 DO i=1,nno+nni
1125 uu=u(1,i)**2+u(2,i)**2+u(3,i)**2
1126 pres(i)=pa-rho*phip(i)-half*rho*uu
1127 ENDDO
1128 ELSE
1129 DO i=1,nno+nni
1130 pres(i)=zero
1131 ENDDO
1132 ENDIF
1133C---------------------------------------------------
1134C 5- Force on the boundary
1135C---------------------------------------------------
1136 IF (nspmd == 1) THEN
1137 DO i=1,nel
1138 n1=elem(1,i)
1139 n2=elem(2,i)
1140 n3=elem(3,i)
1141 coef=one_over_6*third*(pres(n1)+pres(n2)+pres(n3))
1142 nrx=norm(1,i)
1143 nry=norm(2,i)
1144 nrz=norm(3,i)
1145 nn1=ibuf(n1)
1146 nn2=ibuf(n2)
1147 nn3=ibuf(n3)
1148 a(1,nn1)=a(1,nn1)+coef*nrx
1149 a(2,nn1)=a(2,nn1)+coef*nry
1150 a(3,nn1)=a(3,nn1)+coef*nrz
1151 a(1,nn2)=a(1,nn2)+coef*nrx
1152 a(2,nn2)=a(2,nn2)+coef*nry
1153 a(3,nn2)=a(3,nn2)+coef*nrz
1154 a(1,nn3)=a(1,nn3)+coef*nrx
1155 a(2,nn3)=a(2,nn3)+coef*nry
1156 a(3,nn3)=a(3,nn3)+coef*nrz
1157 wfext=wfext+coef*nrx*v(1,nn1)*dt1
1158 . +coef*nry*v(2,nn1)*dt1
1159 . +coef*nrz*v(3,nn1)*dt1
1160 wfext=wfext+coef*nrx*v(1,nn2)*dt1
1161 . +coef*nry*v(2,nn2)*dt1
1162 . +coef*nrz*v(3,nn2)*dt1
1163 wfext=wfext+coef*nrx*v(1,nn3)*dt1
1164 . +coef*nry*v(2,nn3)*dt1
1165 . +coef*nrz*v(3,nn3)*dt1
1166 ENDDO
1167 ELSE
1168 nelf=0
1169 DO i=1,nel
1170 n1=elem(1,i)
1171 n2=elem(2,i)
1172 n3=elem(3,i)
1173 IF (ibuf(n1)/=0.OR.ibuf(n2)/=0.OR.ibuf(n3)/=0) THEN
1174 nelf=nelf+1
1175 lelf(nelf)=i
1176 ENDIF
1177 ENDDO
1178 DO i=1,nelf
1179 ii=lelf(i)
1180 n1=elem(1,ii)
1181 n2=elem(2,ii)
1182 n3=elem(3,ii)
1183 coef=one_over_6*third*(pres(n1)+pres(n2)+pres(n3))
1184 nrx=norm(1,ii)
1185 nry=norm(2,ii)
1186 nrz=norm(3,ii)
1187 IF (ibuf(n1)/=0) THEN
1188 nn1=ibuf(n1)
1189 a(1,nn1)=a(1,nn1)+coef*nrx
1190 a(2,nn1)=a(2,nn1)+coef*nry
1191 a(3,nn1)=a(3,nn1)+coef*nrz
1192 wfext=wfext+coef*nrx*v(1,nn1)*dt1
1193 . +coef*nry*v(2,nn1)*dt1
1194 . +coef*nrz*v(3,nn1)*dt1
1195 ENDIF
1196 IF (ibuf(n2)/=0) THEN
1197 nn2=ibuf(n2)
1198 a(1,nn2)=a(1,nn2)+coef*nrx
1199 a(2,nn2)=a(2,nn2)+coef*nry
1200 a(3,nn2)=a(3,nn2)+coef*nrz
1201 wfext=wfext+coef*nrx*v(1,nn2)*dt1
1202 . +coef*nry*v(2,nn2)*dt1
1203 . +coef*nrz*v(3,nn2)*dt1
1204 ENDIF
1205 IF (ibuf(n3)/=0) THEN
1206 nn3=ibuf(n3)
1207 a(1,nn3)=a(1,nn3)+coef*nrx
1208 a(2,nn3)=a(2,nn3)+coef*nry
1209 a(3,nn3)=a(3,nn3)+coef*nrz
1210 wfext=wfext+coef*nrx*v(1,nn3)*dt1
1211 . +coef*nry*v(2,nn3)*dt1
1212 . +coef*nrz*v(3,nn3)*dt1
1213 ENDIF
1214 ENDDO
1215 ENDIF
1216
1217 RETURN
bemsolv
subroutine bemsolv(l_assemb, ilvout, iform, nn, nel, hbem, gbem, ibuf, tbemtg, telnor, x, q, phi, ipiv, phi_inf)
Definition bemsolv.F:38
bemsolvp
subroutine bemsolvp(elem, norm, x, nnrp, cr_loc_glob, cr_glob_loc, nncp, cc_loc_glob, cc_glob_loc, nceip, ccei_loc_glob, nep, ce_loc_glob, ce_glob_loc, nreip, crei_loc_glob, nno, nel, hbem, gbem, q, phi, iform, l_assemb, nprow, npcol, nbloc, ipiv, nerp, ilvout, phi_inf)
Definition bemsolvp.F:43
my_real
#define my_real
Definition cppsort.cpp:32
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
trgrad
subroutine trgrad(x1, x2, x3, y1, y2, y3, z1, z2, z3, phi1, phi2, phi3, grx, gry, grz)
Definition incpflow.F:1227
ymax
subroutine ymax(idn, fac, npc, pld, stiffmin, stiffmax, stiffini, stiffavg)
Definition law100_upd.F:272
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
spmd_fl_sum
subroutine spmd_fl_sum(lsum, len, lsumt)
Definition spmd_fl_sum.F:37
nl
character *2 function nl()
Definition message.F:2354

◆ trgrad()

subroutine trgrad ( x1,
x2,
x3,
y1,
y2,
y3,
z1,
z2,
z3,
phi1,
phi2,
phi3,
grx,
gry,
grz )

Definition at line 1224 of file incpflow.F.

1227C-----------------------------------------------
1228C I m p l i c i t T y p e s
1229C-----------------------------------------------
1230#include "implicit_f.inc"
1231C-----------------------------------------------
1232C D u m m y A r g u m e n t s
1233C-----------------------------------------------
1234 my_real
1235 . x1, x2, x3, y1, y2, y3, z1, z2, z3,
1236 . phi1, phi2, phi3, grx, gry, grz
1237C-----------------------------------------------
1238C L o c a l V a r i a b l e s
1239C-----------------------------------------------
1240 my_real
1241 . vx1, vy1, vz1, vx2, vy2, vz2, ex1, ex2, ey2,
1242 . xl1, yl1, xl2, yl2, xl3, yl3, jac, jj11, jj12, jj21, jj22,
1243 . grlx, grly
1244
1245 vx1=x2-x1
1246 vy1=y2-y1
1247 vz1=z2-z1
1248 ex1=sqrt(vx1**2+vy1**2+vz1**2)
1249 IF(ex1>em20)THEN
1250 vx1=vx1/ex1
1251 vy1=vy1/ex1
1252 vz1=vz1/ex1
1253 ELSE
1254 vx1=zero
1255 vy1=zero
1256 vz1=zero
1257 ENDIF
1258
1259 vx2=x3-x1
1260 vy2=y3-y1
1261 vz2=z3-z1
1262 ex2=vx1*vx2+vy1*vy2+vz1*vz2
1263 vx2=vx2-ex2*vx1
1264 vy2=vy2-ex2*vy1
1265 vz2=vz2-ex2*vz1
1266 ey2=sqrt(vx2**2+vy2**2+vz2**2)
1267 IF(ey2>em20)THEN
1268 vx2=vx2/ey2
1269 vy2=vy2/ey2
1270 vz2=vz2/ey2
1271 ELSE
1272 vx2=zero
1273 vy2=zero
1274 vz2=zero
1275 ENDIF
1276
1277 xl1=zero
1278 yl1=zero
1279 xl2=ex1
1280 yl2=zero
1281 xl3=ex2
1282 yl3=ey2
1283 jac=(xl2-xl1)*(yl3-yl1)-(yl2-yl1)*(xl3-xl1)
1284 IF(abs(jac)>em20)THEN
1285 jj11=(yl3-yl1)/jac
1286 jj12=-(xl3-xl1)/jac
1287 jj21=-(yl2-yl1)/jac
1288 jj22=(xl2-xl1)/jac
1289 ELSE
1290 jj11=zero
1291 jj12=zero
1292 jj21=zero
1293 jj22=zero
1294 ENDIF
1295
1296 grlx=(phi2-phi1)*jj11+(phi3-phi1)*jj21
1297 grly=(phi2-phi1)*jj12+(phi3-phi1)*jj22
1298
1299 grx=grlx*vx1+grly*vx2
1300 gry=grlx*vy1+grly*vy2
1301 grz=grlx*vz1+grly*vz2
1302
1303 RETURN