OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
pfluid.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!|| pfluid ../engine/source/loads/general/pfluid/pfluid.F
25!||--- called by ------------------------------------------------------
26!|| resol ../engine/source/engine/resol.F
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!|| get_u_numsens ../engine/source/user_interface/usensor.F
30!|| get_u_sens_fpar ../engine/source/user_interface/usensor.F
31!|| get_u_sens_ipar ../engine/source/user_interface/usensor.F
32!|| get_u_sens_value ../engine/source/user_interface/usensor.F
33!|| set_u_sens_value ../engine/source/user_interface/usensor.F
34!||--- uses -----------------------------------------------------
35!|| h3d_mod ../engine/share/modules/h3d_mod.F
36!|| python_funct_mod ../common_source/modules/python_mod.F90
37!|| sensor_mod ../common_source/modules/sensor_mod.F90
38!|| th_surf_mod ../common_source/modules/interfaces/th_surf_mod.F
39!||====================================================================
40 SUBROUTINE pfluid(ILOADP ,RLOAD ,NPC ,TF ,A ,
41 2 V ,X ,XFRAME ,MS ,
42 3 NSENSOR ,SENSOR_TAB ,WFEXC ,WFEXT ,IADC ,
43 4 FSKY ,FSKYV ,LLOADP ,FEXT ,H3D_DATA ,
44 5 TH_SURF ,PYTHON )
45C-----------------------------------------------
46C D e s c r i p t i o n
47C-----------------------------------------------
48C This subroutine is modeling a pressure load on given structural faces (4 or 3 nodes).
49C The pressure load results into a normal force at face centroid. This normal force is then expanded to the nodes composing the face.
50C
51C ILOADP : Integer array related to /LOAD/PBLAST option
52C RLOAD : Real array related to /LOAD/PBLAST option
53C NPC : integer array for /FUNCT options
54C TF : real array for /FUNCT options
55C A : nodal accelerations
56C V : nodal velocities
57C X : nodal coordinates
58C XFRAME : array for /FRAME option
59C AR : rotationnal acceleration
60C VR : rotationnal velocities
61C SENSOR_TAB : data structure for /SENSOR option
62C WEIGHT : -
63C WFEXC : -
64C IADC : (Parith/on only) contains index to be used with FSKY array
65C FSKY : (Parith/on only) Nodal force (directly apply as a nodal acceleration in A array in case of Parith/off)
66C LLOADP : array used to retrieve nodes N1,N2,N3,N4 of the structural face to be loaded
67C FEXT : storage of nodal forces for animation purpose (/ANIM/VECT/FEXT)
68C H3D_DATA : data structure for H3D parameters
69C
70C FORCE STORAGE :
71C Parith/off : Forces are introduced as acceleration directly in A(1:3, 1:NUMNOD) array
72C Parith/on : Forces are saved in FSKY array and will be assembled later using a suitable order
73C
74C SOURCE CODE
75C There are currently 3 parts in this source file
76C -1- "PART-1. PARITH/OFF"
77C -2- "PART-2. PARITH/ON, non vectorial code"
78C -3- "PART-3. PARITH/ON, vectorial code"
79C-----------------------------------------------
80C M o d u l e s
81C-----------------------------------------------
82 USE python_funct_mod
83 USE h3d_mod
84 USE sensor_mod
85 USE th_surf_mod , ONLY : th_surf_
86C-----------------------------------------------
87C I m p l i c i t T y p e s
88C-----------------------------------------------
89#include "implicit_f.inc"
90#include "comlock.inc"
91#include "param_c.inc"
92C-----------------------------------------------
93C C o m m o n B l o c k s
94C-----------------------------------------------
95#include "com01_c.inc"
96#include "com04_c.inc"
97!#include "com06_c.inc"
98#include "com08_c.inc"
99#include "parit_c.inc"
100#include "scr14_c.inc"
101#include "scr16_c.inc"
102#include "tabsiz_c.inc"
103#include "mvsiz_p.inc"
104C-----------------------------------------------
105C E x t e r n a l F u n c t i o n s
106C-----------------------------------------------
107 INTEGER GET_U_NUMSENS,GET_U_SENS_FPAR,GET_U_SENS_IPAR,GET_U_SENS_VALUE,SET_U_SENS_VALUE
108 EXTERNAL GET_U_NUMSENS,GET_U_SENS_FPAR,GET_U_SENS_IPAR,GET_U_SENS_VALUE,SET_U_SENS_VALUE
109C-----------------------------------------------,
110C D u m m y A r g u m e n t s
111C-----------------------------------------------
112 INTEGER, INTENT(IN) :: MS(NUMNOD)
113 INTEGER ,INTENT(IN) :: NSENSOR
114 INTEGER NPC(SNPC),LLOADP(SLLOADP)
115 INTEGER ILOADP(SIZLOADP,NLOADP)
116 INTEGER IADC(DSNCOL)
117 my_real rload(lfacload,nloadp), tf(stf), a(3,numnod), v(3,numnod),
118 . x(3,numnod), xframe(nxframe,numfram+1),wfexc,
119 . fsky(8,sfsky/8), fskyv(sfsky/8,8),fext(3,numnod)
120 TYPE(h3d_database) :: H3D_DATA
121 TYPE (SENSOR_STR_) ,DIMENSION(NSENSOR) ,INTENT(IN) :: SENSOR_TAB
122 TYPE (TH_SURF_) , INTENT(INOUT) :: TH_SURF
123 TYPE(PYTHON_) :: PYTHON
124 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
125C-----------------------------------------------
126C L o c a l V a r i a b l e s
127C-----------------------------------------------
128 INTEGER NL, N1, N2, N3, N4, FUN_HSP, K1, K2, K3, ISENS,K,
129 . IAD,N_OLD,IFRA1,DIR_HSP,I,
130 . FUN_CX,FUN_VEL,DIR_VEL,IFRA2, IANIM,IJK,UP_BOUND,NS,KSURF,NSEGPL
131 my_real NX, NY, NZ, AA,FX, FY, FZ, DYDX, TS,
132 . sixth,wfextt,x_old,xsens,fcx,fcy,
133 . fcx1,fcy1,fcx2,fcy2,vel,vseg,finter,
134 . centroid_depth,pvel,norm,nsign,area,pa,mass,an,vn,pvel0
135 EXTERNAL finter
136 INTEGER :: ISMOOTH
137C-----------------------------------------------
138C S o u r c e L i n e s
139C-----------------------------------------------
140
141 ! init.
142 sixth = one_over_6
143 wfexc = zero
144 wfextt = zero
145 n_old = 0
146 x_old = zero
147 nsegpl = 0
148 ianim = anim_v(5)+outp_v(5) + h3d_data%N_VECT_FINT + anim_v(6)+outp_v(6) + h3d_data%N_VECT_FEXT
149
150C---------------------------------------
151C PART-1. PARITH/OFF
152C---------------------------------------
153
154 IF(iparit == 0) THEN
155
156 DO nl=1,nloadp_f !loop over /LOAD/PFLUID options
157
158 !user parameters (Defined in Starter while reading input file hm_read_pblast.F, and transmitted to Engine with Restart file)
159 fun_hsp = iloadp(7,nl)
160 dir_hsp = iloadp(8,nl)
161 ifra1 = iloadp(9,nl)
162 fcy = rload(1,nl)
163 fcx = rload(2,nl)
164 fun_cx = iloadp(10,nl)
165 fcy1 = rload(3,nl)
166 fcx1 = rload(4,nl)
167 fun_vel = iloadp(11,nl)
168 fcy2 = rload(5,nl)
169 fcx2 = rload(6,nl)
170 dir_vel = max(iloadp(12,nl),1) !To avoid a check bound issue when the velocity options are useless
171 ifra2 = iloadp(13,nl)
172 ifra2 = max(ifra2,1)
173 isens = 0
174
175 !possible sensor
176 xsens = one
177 DO k=1,nsensor
178 IF(iloadp(6,nl) == sensor_tab(k)%SENS_ID) isens=k ! should be moved in Starter to optimize performance
179 ENDDO
180 IF(isens == 0)THEN
181 ts=tt
182 ELSE
183 ts = tt-sensor_tab(isens)%TSTART
184 IF(ts < zero) cycle
185 ENDIF
186
187!$OMP DO SCHEDULE(GUIDED,MVSIZ)
188 DO i = 1,iloadp(1,nl)/4
189 !Structural face is N1,N2,N3,N4
190 n1=lloadp(iloadp(4,nl)+4*(i-1))
191 n2=lloadp(iloadp(4,nl)+4*(i-1)+1)
192 n3=lloadp(iloadp(4,nl)+4*(i-1)+2)
193 n4=lloadp(iloadp(4,nl)+4*(i-1)+3)
194 aa=zero
195 vel=zero
196 pvel=zero
197 !---QUADRAGLE--------------------------------------------
198 IF(n4 /= 0 .AND. n1 /= n2 .AND. n1 /= n3 .AND. n1 /= n4 .AND. n2 /= n3 .AND. n2 /= n4 .AND. n3 /= n4 )THEN !we could optimize performance by moving this check in Starter
199 k1=3*dir_hsp-2
200 k2=3*dir_hsp-1
201 k3=3*dir_hsp
202 ! hydrostatic pressure
203 IF(fun_hsp /= 0)THEN
204 centroid_depth = (xframe(k1,ifra1)*(x(1,n1)+x(1,n2)+x(1,n3)+x(1,n4))/four)
205 centroid_depth = centroid_depth +(xframe(k2,ifra1)*(x(2,n1)+x(2,n2)+x(2,n3)+x(2,n4))/four)
206 centroid_depth = centroid_depth +(xframe(k3,ifra1)*(x(3,n1)+x(3,n2)+x(3,n3)+x(3,n4))/four)
207 ismooth = 0
208 IF (fun_hsp > 0) ismooth = npc(2*nfunct+fun_hsp+1)
209 IF(ismooth < 0) THEN
210 CALL python_call_funct1d(python, -ismooth,(centroid_depth-xframe(9+dir_hsp,ifra1))*fcx , aa)
211 aa = aa*fcy
212 ELSE
213 aa = fcy*finter(fun_hsp,(centroid_depth-xframe(9+dir_hsp,ifra1))*fcx,npc,tf,dydx)
214 ENDIF
215 ENDIF
216 !normal vector
217 nx=(x(2,n3)-x(2,n1))*(x(3,n4)-x(3,n2)) - (x(3,n3)-x(3,n1))*(x(2,n4)-x(2,n2))
218 ny=(x(3,n3)-x(3,n1))*(x(1,n4)-x(1,n2)) - (x(1,n3)-x(1,n1))*(x(3,n4)-x(3,n2))
219 nz=(x(1,n3)-x(1,n1))*(x(2,n4)-x(2,n2)) - (x(2,n3)-x(2,n1))*(x(1,n4)-x(1,n2))
220 norm = sqrt(nx*nx+ny*ny+nz*nz)
221 area = half * norm
222 pa = aa
223 aa = aa * area
224
225 !face velocity
226 k1=3*dir_vel-2
227 k2=3*dir_vel-1
228 k3=3*dir_vel
229 vseg = (xframe(k1,ifra2)*(v(1,n1) + v(1,n2) + v(1,n3) + v(1,n4)) /four)
230 vseg = vseg + (xframe(k2,ifra2)*(v(2,n1) + v(2,n2) + v(2,n3) + v(2,n4)) /four)
231 vseg = vseg + (xframe(k3,ifra2)*(v(3,n1) + v(3,n2) + v(3,n3) + v(3,n4)) /four)
232
233 !VEL,PVEL, and NSIGN
234 IF(fun_vel /= 0)THEN
235 ismooth = 0
236 IF (fun_vel > 0) ismooth = npc(2*nfunct+fun_vel+1)
237 IF(ismooth < 0) THEN
238 CALL python_call_funct1d(python, -ismooth,tt*fcx2, vel)
239 vel = vel*fcy2 - vseg
240 ELSE
241 vel = fcy2*finter(fun_vel,tt*fcx2,npc,tf,dydx) - vseg
242 ENDIF
243 ELSE
244 vel = -vseg
245 ENDIF
246 IF(fun_cx /= 0) THEN
247 pvel = ( (-(nx/norm)*vel*xframe(k1,ifra2)-(ny/norm)*vel*xframe(k2,ifra2)-(nz/norm)*vel*xframe(k3,ifra2))**2 )
248 ismooth = 0
249 IF (fun_cx > 0) ismooth = npc(2*nfunct+fun_cx+1)
250 IF(ismooth < 0) THEN
251 CALL python_call_funct1d(python, -ismooth,tt*fcx1, pvel0)
252 pvel = pvel*pvel0*fcy1/two
253 ELSE
254 pvel=pvel*fcy1*finter(fun_cx,tt*fcx1,npc,tf,dydx)/two
255 ENDIF
256 ENDIF
257 nsign = vel*(nx*xframe(k1,ifra2) + ny*xframe(k2,ifra2) + nz*xframe(k3,ifra2)) ! <Vseg.d,n> where d is flow direction X,Y,Z ! IFRA2=1 is general frame, IFRA2 in [2:NFRAME+1] are for user frames
258 nsign = sign(one,nsign)
259
260 !hydrostatic pressure
261 fx=-aa*(nx/norm)
262 fy=-aa*(ny/norm)
263 fz=-aa*(nz/norm)
264 !drag force
265 fx=fx+pvel*half*nx*nsign
266 fy=fy+pvel*half*ny*nsign
267 fz=fz+pvel*half*nz*nsign
268 !expanding on the 4-node-face
269 fx=fx*fourth
270 fy=fy*fourth
271 fz=fz*fourth
272
273 !External Force Work
274 wfextt=wfextt+dt1*(fx*(v(1,n1)+v(1,n2)+v(1,n3)+v(1,n4))
275 + +fy*(v(2,n1)+v(2,n2)+v(2,n3)+v(2,n4))
276 + +fz*(v(3,n1)+v(3,n2)+v(3,n3)+v(3,n4)))
277
278#include "lockon.inc"
279 !/PARITH/OFF: force is directly added in A array. It will be divided by nodal mass later
280 !-node_1
281 a(1,n1)=a(1,n1)+fx
282 a(2,n1)=a(2,n1)+fy
283 a(3,n1)=a(3,n1)+fz
284 !-node_2
285 a(1,n2)=a(1,n2)+fx
286 a(2,n2)=a(2,n2)+fy
287 a(3,n2)=a(3,n2)+fz
288 !-node_3
289 a(1,n3)=a(1,n3)+fx
290 a(2,n3)=a(2,n3)+fy
291 a(3,n3)=a(3,n3)+fz
292 !-node_4
293 a(1,n4)=a(1,n4)+fx
294 a(2,n4)=a(2,n4)+fy
295 a(3,n4)=a(3,n4)+fz
296#include "lockoff.inc"
297
298 IF(ianim > 0) THEN
299#include "lockon.inc"
300 fext(1,n1) = fext(1,n1)+fx
301 fext(2,n1) = fext(2,n1)+fy
302 fext(3,n1) = fext(3,n1)+fz
303 !
304 fext(1,n2) = fext(1,n2)+fx
305 fext(2,n2) = fext(2,n2)+fy
306 fext(3,n2) = fext(3,n2)+fz
307 !
308 fext(1,n3) = fext(1,n3)+fx
309 fext(2,n3) = fext(2,n3)+fy
310 fext(3,n3) = fext(3,n3)+fz
311 !
312 fext(1,n4) = fext(1,n4)+fx
313 fext(2,n4) = fext(2,n4)+fy
314 fext(3,n4) = fext(3,n4)+fz
315#include "lockoff.inc"
316 ENDIF
317C
318 IF(th_surf%LOADP_FLAG >0 ) THEN
319 !mean velocity (/TH/SURF)
320 vn = v(1,n1)*nx + v(2,n1)*ny + v(3,n1)*nz
321 vn = vn + v(1,n2)*nx + v(2,n2)*ny + v(3,n2)*nz
322 vn = vn + v(1,n3)*nx + v(2,n3)*ny + v(3,n3)*nz
323 vn = vn + v(1,n4)*nx + v(2,n4)*ny + v(3,n4)*nz
324 vn = vn / four
325 !mean acceleration (/TH/SURF)
326 an = fx*nx + fy*ny + fz*nz
327 mass = ( ms(n1)+ms(n2)+ms(n3)+ms(n4) ) / four
328 an = an / mass
329 !TH_SURF%CHANNELS output
330 nsegpl = nsegpl + 1
331#include "lockon.inc"
332 DO ns=th_surf%LOADP_KSEGS(nsegpl) +1,th_surf%LOADP_KSEGS(nsegpl+1)
333 ksurf = th_surf%LOADP_SEGS(ns)
334 th_surf%CHANNELS(3,ksurf)= th_surf%CHANNELS(3,ksurf) + area*vn ! velocity
335 th_surf%CHANNELS(4,ksurf)= th_surf%CHANNELS(4,ksurf) + area*pa ! pressure
336 th_surf%CHANNELS(5,ksurf)= th_surf%CHANNELS(5,ksurf) + area
337 ENDDO
338#include "lockoff.inc"
339 ENDIF
340
341
342 !---TRIANGLE--------------------------------------------
343 ELSE
344 IF(n1 == n2)THEN
345 n2 = n3
346 n3 = n4
347 n4 = 0
348 ELSEIF(n1 == n3)THEN
349 n3 = n4
350 n4 = 0
351 ELSEIF(n1 == n4)THEN
352 n4 = 0
353 ELSEIF(n2 == n3)THEN
354 n3 = n4
355 n4 = 0
356 ELSEIF(n2 == n4)THEN
357 n2 = n3
358 n3 = n4
359 n4 = 0
360 ELSEIF(n3 == n4)THEN
361 n4 = 0
362 ENDIF
363
364 !hydrostatic pressure
365 IF(fun_hsp /= 0)THEN
366 k1=3*dir_hsp-2
367 k2=3*dir_hsp-1
368 k3=3*dir_hsp
369
370 centroid_depth = (xframe(k1,ifra1)*(x(1,n1)+x(1,n2)+x(1,n3))/three)
371 centroid_depth = centroid_depth+(xframe(k2,ifra1)*(x(2,n1)+x(2,n2)+x(2,n3))/three)
372 centroid_depth = centroid_depth+(xframe(k3,ifra1)*(x(3,n1)+x(3,n2)+x(3,n3))/three)
373 ismooth = 0
374 IF (fun_hsp > 0) ismooth = npc(2*nfunct+fun_hsp+1)
375 IF(ismooth < 0) THEN
376 CALL python_call_funct1d(python, -ismooth, (centroid_depth-xframe(9+dir_hsp,ifra1))*fcx ,aa )
377 aa = aa*fcy
378 ELSE
379 aa = fcy*finter(fun_hsp,(centroid_depth-xframe(9+dir_hsp,ifra1))*fcx,npc,tf,dydx)
380 ENDIF
381 ENDIF
382 !normal vector
383 nx = (x(2,n3)-x(2,n1))*(x(3,n3)-x(3,n2)) - (x(3,n3)-x(3,n1))*(x(2,n3)-x(2,n2))
384 ny = (x(3,n3)-x(3,n1))*(x(1,n3)-x(1,n2)) - (x(1,n3)-x(1,n1))*(x(3,n3)-x(3,n2))
385 nz = (x(1,n3)-x(1,n1))*(x(2,n3)-x(2,n2)) - (x(2,n3)-x(2,n1))*(x(1,n3)-x(1,n2))
386 norm = sqrt(nx*nx+ny*ny+nz*nz)
387 area = half * norm
388 pa = aa
389 aa = aa * area
390
391 !face velocity
392 k1=3*dir_vel-2
393 k2=3*dir_vel-1
394 k3=3*dir_vel
395 vseg= (xframe(k1,ifra2)* (v(1,n1) + v(1,n2) + v(1,n3)) /three)
396 vseg=vseg+(xframe(k2,ifra2)* (v(2,n1) + v(2,n2) + v(2,n3)) /three)
397 vseg=vseg+(xframe(k3,ifra2)* (v(3,n1) + v(3,n2) + v(3,n3)) /three)
398
399 ! VEL,PVEL,and NSIGN
400 IF(fun_vel /= 0)THEN
401 ismooth = 0
402 IF (fun_vel > 0) ismooth = npc(2*nfunct+fun_vel+1)
403 IF(ismooth < 0) THEN
404 CALL python_call_funct1d(python, -ismooth,tt*fcx2 , vel )
405 vel = vel*fcy2 - vseg
406 ELSE
407 vel = fcy2*finter(fun_vel,tt*fcx2,npc,tf,dydx) - vseg
408 ENDIF
409 ELSE
410 vel = -vseg
411 ENDIF
412 IF(fun_cx /= 0)THEN
413 pvel = ( (-(nx/norm)*vel*xframe(k1,ifra2)-(ny/norm)*vel*xframe(k2,ifra2)-(nz/norm)*vel*xframe(k3,ifra2))**2 )
414 ismooth = 0
415 IF (fun_cx > 0) ismooth = npc(2*nfunct+fun_cx+1)
416 IF(ismooth < 0) THEN
417 CALL python_call_funct1d(python, -ismooth,tt*fcx1 , pvel0)
418 pvel = pvel*pvel0*fcy1/two
419 ELSE
420 pvel=pvel*fcy1*finter(fun_cx,tt*fcx1,npc,tf,dydx)/two
421 ENDIF
422 ENDIF
423 nsign = vel*(nx*xframe(k1,ifra2) + ny*xframe(k2,ifra2) + nz*xframe(k3,ifra2)) ! <Vseg.d,n> where d is flow direction X,Y,Z ! IFRA2=1 is general frame, IFRA2 in [2:NFRAME+1] are for user frames
424 nsign = sign(one,nsign)
425
426 !hydrostatic pressure
427 fx=-aa*(nx/norm)
428 fy=-aa*(ny/norm)
429 fz=-aa*(nz/norm)
430 !drag force
431 fx=fx+pvel*half*nx*nsign
432 fy=fy+pvel*half*ny*nsign
433 fz=fz+pvel*half*nz*nsign
434 !expanding on the 4-node-face
435 fx=fx*third
436 fy=fy*third
437 fz=fz*third
438
439 !External Force Work
440 wfextt=wfextt+dt1*(fx*(v(1,n1)+v(1,n2)+v(1,n3))
441 + +fy*(v(2,n1)+v(2,n2)+v(2,n3))
442 + +fz*(v(3,n1)+v(3,n2)+v(3,n3)))
443
444
445#include "lockon.inc"
446 !/PARITH/OFF: force is directly added in A array. It will be dividedby nodal mass later
447 !-node_1
448 a(1,n1)=a(1,n1)+fx
449 a(2,n1)=a(2,n1)+fy
450 a(3,n1)=a(3,n1)+fz
451 !-node_2
452 a(1,n2)=a(1,n2)+fx
453 a(2,n2)=a(2,n2)+fy
454 a(3,n2)=a(3,n2)+fz
455 !-node_3
456 a(1,n3)=a(1,n3)+fx
457 a(2,n3)=a(2,n3)+fy
458 a(3,n3)=a(3,n3)+fz
459#include "lockoff.inc"
460
461
462
463 IF(ianim > 0) THEN
464#include "lockon.inc"
465 fext(1,n1) = fext(1,n1)+fx
466 fext(2,n1) = fext(2,n1)+fy
467 fext(3,n1) = fext(3,n1)+fz
468 !
469 fext(1,n2) = fext(1,n2)+fx
470 fext(2,n2) = fext(2,n2)+fy
471 fext(3,n2) = fext(3,n2)+fz
472 !
473 fext(1,n3) = fext(1,n3)+fx
474 fext(2,n3) = fext(2,n3)+fy
475 fext(3,n3) = fext(3,n3)+fz
476#include "lockoff.inc"
477 ENDIF
478
479 IF(th_surf%LOADP_FLAG >0 ) THEN
480 !mean velocity (/TH/SURF)
481 vn = v(1,n1)*nx + v(2,n1)*ny + v(3,n1)*nz
482 vn = vn + v(1,n2)*nx + v(2,n2)*ny + v(3,n2)*nz
483 vn = vn + v(1,n3)*nx + v(2,n3)*ny + v(3,n3)*nz
484 vn = third * vn
485 !mean acceleration (/TH/SURF)
486 an = fx*nx + fy*ny + fz*nz
487 mass = third * ( ms(n1)+ms(n2)+ms(n3) )
488 an = an / mass
489 !TH_SURF%CHANNELS output
490 nsegpl = nsegpl + 1
491#include "lockon.inc"
492 DO ns=th_surf%LOADP_KSEGS(nsegpl) +1,th_surf%LOADP_KSEGS(nsegpl+1)
493 ksurf = th_surf%LOADP_SEGS(ns)
494 th_surf%CHANNELS(3,ksurf)= th_surf%CHANNELS(3,ksurf) + area*vn ! velocity
495 th_surf%CHANNELS(4,ksurf)= th_surf%CHANNELS(4,ksurf) + area*pa ! pressure
496 th_surf%CHANNELS(5,ksurf)= th_surf%CHANNELS(5,ksurf) + area
497 ENDDO
498#include "lockoff.inc"
499 ENDIF
500
501 endif! quadrangle / triangle
502 ENDDO !next I
503!$OMP END DO
504
505 ENDDO !nextNL
506
507!$OMP ATOMIC
508 wfext = wfext + wfextt
509
510
511
512
513 ELSE !IF(IPARIT == 0) THEN
514
515 !otherwise IPARIT /= 0
516
517C---------------------------------------
518C PART-2. PARITH/ON, non vectorial code
519C---------------------------------------
520 IF(ivector == 0) THEN
521
522
523 DO nl=1,nloadp_f !loop over /LOAD/PFLUID options
524
525 !user parameters (Defined in Starter while reading input file hm_read_pblast.F, and transmitted to Engine with Restart file)
526 fun_hsp = iloadp(7,nl)
527 dir_hsp = iloadp(8,nl)
528 ifra1 = iloadp(9,nl)
529 fcy = rload(1,nl)
530 fcx = rload(2,nl)
531 fun_cx = iloadp(10,nl)
532 fcy1 = rload(3,nl)
533 fcx1 = rload(4,nl)
534 fun_vel = iloadp(11,nl)
535 fcy2 = rload(5,nl)
536 fcx2 = rload(6,nl)
537 dir_vel = max(iloadp(12,nl),1) ! To avoid a check bound issue when velocity is useless
538 ifra2 = iloadp(13,nl)
539 isens = 0
540 xsens = one
541
542
543 ! FLUSH fsky array to 0.
544!$OMP DO SCHEDULE(GUIDED,MVSIZ)
545 DO i = 1,iloadp(1,nl)/4
546 !nodes of structural face : N1,N2,N3,N4
547 n1=lloadp(iloadp(4,nl)+4*(i-1))
548 n2=lloadp(iloadp(4,nl)+4*(i-1)+1)
549 n3=lloadp(iloadp(4,nl)+4*(i-1)+2)
550 n4=lloadp(iloadp(4,nl)+4*(i-1)+3)
551 !---QUADRAGLE--------------------------------------------
552 IF(n4 /= 0 .AND. n1 /= n2 .AND. n1 /= n3 .AND. n1 /= n4 .AND.n2 /= n3 .AND. n2 /= n4 .AND. n3 /= n4 )THEN
553 up_bound=4
554 ELSE
555 up_bound=3
556 ENDIF
557 DO ijk=1,up_bound
558 iad = iadc(iloadp(4,nl)+4*(i-1)+(ijk-1))
559 fsky(1:3,iad) = zero
560 ENDDO
561 ENDDO
562!$OMP END DO
563
564 !possible sensor
565 DO k=1,nsensor
566 IF(iloadp(6,nl) == sensor_tab(k)%SENS_ID) isens=k
567 ENDDO
568 IF(isens == 0)THEN
569 ts=tt
570 ELSE
571 ts = tt- sensor_tab(isens)%TSTART
572 IF(ts < zero) cycle
573 ENDIF
574
575!$OMP DO SCHEDULE(GUIDED,MVSIZ)
576 DO i = 1,iloadp(1,nl)/4
577
578 !nodes of structural face : N1,N2,N3,N4
579 n1=lloadp(iloadp(4,nl)+4*(i-1))
580 n2=lloadp(iloadp(4,nl)+4*(i-1)+1)
581 n3=lloadp(iloadp(4,nl)+4*(i-1)+2)
582 n4=lloadp(iloadp(4,nl)+4*(i-1)+3)
583
584 aa = zero
585 vel=zero
586 pvel=zero
587
588 !---QUADRAGLE--------------------------------------------
589 IF(n4 /= 0 .AND. n1 /= n2 .AND. n1 /= n3 .AND. n1 /= n4 .AND.n2 /= n3 .AND. n2 /= n4 .AND. n3 /= n4 )THEN
590
591 !hydrostatic pressure
592 k1=3*dir_hsp-2
593 k2=3*dir_hsp-1
594 k3=3*dir_hsp
595 IF(fun_hsp /= 0) THEN
596 centroid_depth = (xframe(k1,ifra1)*(x(1,n1)+x(1,n2)+x(1,n3)+x(1,n4))/four)
597 centroid_depth = centroid_depth+(xframe(k2,ifra1)*(x(2,n1)+x(2,n2)+x(2,n3)+x(2,n4))/four)
598 centroid_depth = centroid_depth+(xframe(k3,ifra1)*(x(3,n1)+x(3,n2)+x(3,n3)+x(3,n4))/four)
599 ismooth = 0
600 IF (fun_hsp > 0) ismooth = npc(2*nfunct+fun_hsp+1)
601 IF(ismooth < 0) THEN
602 CALL python_call_funct1d(python, -ismooth, (centroid_depth-xframe(9+dir_hsp,ifra1))*fcx , aa )
603 aa = aa*fcy
604 ELSE
605 aa = fcy*finter(fun_hsp,(centroid_depth-xframe(9+dir_hsp,ifra1))*fcx,npc,tf,dydx)
606 ENDIF
607
608 ENDIF
609 !normal vector
610 nx = (x(2,n3)-x(2,n1))*(x(3,n4)-x(3,n2))-(x(3,n3)-x(3,n1))*(x(2,n4)-x(2,n2))
611 ny = (x(3,n3)-x(3,n1))*(x(1,n4)-x(1,n2))-(x(1,n3)-x(1,n1))*(x(3,n4)-x(3,n2))
612 nz = (x(1,n3)-x(1,n1))*(x(2,n4)-x(2,n2))-(x(2,n3)-x(2,n1))*(x(1,n4)-x(1,n2))
613 norm = sqrt(nx*nx+ny*ny+nz*nz)
614 area = half * norm
615 pa = aa
616 aa = aa * area
617
618 !structural face velocity
619 k1=3*dir_vel-2
620 k2=3*dir_vel-1
621 k3=3*dir_vel
622 vseg= (xframe(k1,ifra2)*(v(1,n1) + v(1,n2) + v(1,n3) + v(1,n4)) /four)
623 vseg=vseg+(xframe(k2,ifra2)*(v(2,n1) + v(2,n2) + v(2,n3) + v(2,n4)) /four)
624 vseg=vseg+(xframe(k3,ifra2)*(v(3,n1) + v(3,n2) + v(3,n3) + v(3,n4)) /four)
625
626 !VEL, PVEL, and NSIGN
627 IF(fun_vel /= 0)THEN
628 ismooth = 0
629 IF (fun_vel > 0) ismooth = npc(2*nfunct+fun_vel+1)
630 IF(ismooth < 0) THEN
631 CALL python_call_funct1d(python, -ismooth, tt*fcx2 ,vel )
632 vel = vel*fcy2 - vseg
633 ELSE
634 vel = fcy2*finter(fun_vel,tt*fcx2,npc,tf,dydx) - vseg
635 ENDIF
636 ELSE
637 vel = -vseg
638 ENDIF
639 IF(fun_cx /= 0)THEN
640 pvel = ( (-(nx/norm)*vel*xframe(k1,ifra2) -(ny/norm)*vel*xframe(k2,ifra2)-(nz/norm)*vel*xframe(k3,ifra2))**2 )
641 ismooth = 0
642 IF (fun_cx > 0) ismooth = npc(2*nfunct+fun_cx+1)
643 IF(ismooth < 0) THEN
644 CALL python_call_funct1d(python, -ismooth,tt*fcx1, pvel0)
645 pvel = pvel*pvel0*fcy1/two
646 ELSE
647 pvel=pvel*fcy1*finter(fun_cx,tt*fcx1,npc,tf,dydx)/two
648 ENDIF
649 ENDIF
650 nsign = vel*(nx*xframe(k1,ifra2) + ny*xframe(k2,ifra2) + nz*xframe(k3,ifra2)) ! <Vseg.d,n> where d is flow direction X,Y,Z ! IFRA2=1 is general frame, IFRA2 in [2:NFRAME+1] are for user frames
651 nsign = sign(one,nsign)
652
653 !hydrostatic pressure
654 fx=-aa*(nx/norm)
655 fy=-aa*(ny/norm)
656 fz=-aa*(nz/norm)
657 !drag force
658 fx=fx+pvel*half*nx*nsign
659 fy=fy+pvel*half*ny*nsign
660 fz=fz+pvel*half*nz*nsign
661 !expanding on the 4-node-face
662 fx=fx*fourth
663 fy=fy*fourth
664 fz=fz*fourth
665
666 !NODE_1 : storing nodal force in FSKY array (it will be assembled later)
667 iad = iadc(iloadp(4,nl)+4*(i-1))
668 fsky(1,iad) = fx
669 fsky(2,iad) = fy
670 fsky(3,iad) = fz
671
672 !NODE_2 : storing nodal force in FSKY array (it will be assembled later)
673 iad = iadc(iloadp(4,nl)+4*(i-1)+1)
674 fsky(1,iad) = fx
675 fsky(2,iad) = fy
676 fsky(3,iad) = fz
677
678 !NODE_3 : storing nodal force in FSKY array (it will be assembled later)
679 iad = iadc(iloadp(4,nl)+4*(i-1)+2)
680 fsky(1,iad) = fx
681 fsky(2,iad) = fy
682 fsky(3,iad) = fz
683
684 !NODE_4 : storing nodal force in FSKY array (it will be assembled later)
685 iad = iadc(iloadp(4,nl)+4*(i-1)+3)
686 fsky(1,iad) = fx
687 fsky(2,iad) = fy
688 fsky(3,iad) = fz
689
690 !External force work
691 wfextt=wfextt+dt1*(fx*(v(1,n1)+v(1,n2)+v(1,n3)+v(1,n4))
692 + +fy*(v(2,n1)+v(2,n2)+v(2,n3)+v(2,n4))
693 + +fz*(v(3,n1)+v(3,n2)+v(3,n3)+v(3,n4)))
694
695 IF(ianim > 0) THEN
696#include "lockon.inc"
697 fext(1,n1) = fext(1,n1)+fx
698 fext(2,n1) = fext(2,n1)+fy
699 fext(3,n1) = fext(3,n1)+fz
700 !
701 fext(1,n2) = fext(1,n2)+fx
702 fext(2,n2) = fext(2,n2)+fy
703 fext(3,n2) = fext(3,n2)+fz
704 !
705 fext(1,n3) = fext(1,n3)+fx
706 fext(2,n3) = fext(2,n3)+fy
707 fext(3,n3) = fext(3,n3)+fz
708 !
709 fext(1,n4) = fext(1,n4)+fx
710 fext(2,n4) = fext(2,n4)+fy
711 fext(3,n4) = fext(3,n4)+fz
712#include "lockoff.inc"
713 ENDIF
714C
715 IF(th_surf%LOADP_FLAG >0 ) THEN
716 !mean velocity (/TH/SURF)
717 vn = v(1,n1)*nx + v(2,n1)*ny + v(3,n1)*nz
718 vn = vn + v(1,n2)*nx + v(2,n2)*ny + v(3,n2)*nz
719 vn = vn + v(1,n3)*nx + v(2,n3)*ny + v(3,n3)*nz
720 vn = vn + v(1,n4)*nx + v(2,n4)*ny + v(3,n4)*nz
721 vn = vn / four
722 !mean acceleration (/TH/SURF)
723 an = fx*nx + fy*ny + fz*nz
724 mass = ( ms(n1)+ms(n2)+ms(n3)+ms(n4) ) / four
725 an = an / mass
726 !TH_SURF%CHANNELS output
727 nsegpl = nsegpl + 1
728 DO ns=th_surf%LOADP_KSEGS(nsegpl) +1,th_surf%LOADP_KSEGS(nsegpl+1)
729#include "lockon.inc"
730 ksurf = th_surf%LOADP_SEGS(ns)
731 th_surf%CHANNELS(3,ksurf)= th_surf%CHANNELS(3,ksurf) + area*vn ! velocity
732 th_surf%CHANNELS(4,ksurf)= th_surf%CHANNELS(4,ksurf) + area*pa ! pressure
733 th_surf%CHANNELS(5,ksurf)= th_surf%CHANNELS(5,ksurf) + area
734#include "lockoff.inc"
735 ENDDO
736 ENDIF
737C
738 ELSE
739
740 !---TRIANGLE--------------------------------------------
741 IF(n1 == n2)THEN
742 n2 = n3
743 n3 = n4
744 n4 = 0
745 ELSEIF(n1 == n3)THEN
746 n3 = n4
747 n4 = 0
748 ELSEIF(n1 == n4)THEN
749 n4 = 0
750 ELSEIF(n2 == n3)THEN
751 n3 = n4
752 n4 = 0
753 ELSEIF(n2 == n4)THEN
754 n2 = n3
755 n3 = n4
756 n4 = 0
757 ELSEIF(n3 == n4)THEN
758 n4 = 0
759 ENDIF
760
761 !hydrostatic pressure
762 IF(fun_hsp /= 0)THEN
763 k1=3*dir_hsp-2
764 k2=3*dir_hsp-1
765 k3=3*dir_hsp
766 centroid_depth = (xframe(k1,ifra1)*(x(1,n1)+x(1,n2)+x(1,n3))/three)
767 centroid_depth = centroid_depth + (xframe(k2,ifra1)*(x(2,n1)+x(2,n2)+x(2,n3))/three)
768 centroid_depth = centroid_depth + (xframe(k3,ifra1)*(x(3,n1)+x(3,n2)+x(3,n3))/three)
769 ismooth = 0
770 IF (fun_hsp > 0) ismooth = npc(2*nfunct+fun_hsp+1)
771 IF(ismooth < 0) THEN
772 CALL python_call_funct1d(python, -ismooth,(centroid_depth-xframe(9+dir_hsp,ifra1))*fcx , aa)
773 aa = aa*fcy
774 ELSE
775 aa = fcy*finter(fun_hsp,(centroid_depth-xframe(9+dir_hsp,ifra1))*fcx,npc,tf,dydx)
776 ENDIF
777
778 ENDIF
779
780 !normal vector
781 nx = (x(2,n3)-x(2,n1))*(x(3,n3)-x(3,n2))-(x(3,n3)-x(3,n1))*(x(2,n3)-x(2,n2))
782 ny = (x(3,n3)-x(3,n1))*(x(1,n3)-x(1,n2))-(x(1,n3)-x(1,n1))*(x(3,n3)-x(3,n2))
783 nz = (x(1,n3)-x(1,n1))*(x(2,n3)-x(2,n2))-(x(2,n3)-x(2,n1))*(x(1,n3)-x(1,n2))
784 norm = sqrt(nx*nx+ny*ny+nz*nz)
785 area = half * norm
786 pa = aa
787 aa = aa * area
788
789 !structural face velocity
790 k1=3*dir_vel-2
791 k2=3*dir_vel-1
792 k3=3*dir_vel
793 vseg= (xframe(k1,ifra2)*(v(1,n1) + v(1,n2) + v(1,n3)) /three)
794 vseg=vseg+(xframe(k2,ifra2)*(v(2,n1) + v(2,n2) + v(2,n3)) /three)
795 vseg=vseg+(xframe(k3,ifra2)*(v(3,n1) + v(3,n2) + v(3,n3)) /three)
796
797 !VEL, PVEL, and NSIGN
798 IF(fun_vel /= 0)THEN
799 ismooth = 0
800 IF (fun_vel > 0) ismooth = npc(2*nfunct+fun_vel+1)
801 IF(ismooth < 0) THEN
802 CALL python_call_funct1d(python, -ismooth, tt*fcx2 ,vel )
803 vel = vel*fcy2 - vseg
804 ELSE
805 vel = fcy2*finter(fun_vel,tt*fcx2,npc,tf,dydx) - vseg
806 ENDIF
807 ELSE
808 vel = -vseg
809 ENDIF
810 IF(fun_cx /= 0) THEN
811 pvel = ( (-(nx/norm)*vel*xframe(k1,ifra2)-(ny/norm)*vel*xframe(k2,ifra2)-(nz/norm)*vel*xframe(k3,ifra2))**2 )
812 ismooth = 0
813 IF (fun_cx > 0) ismooth = npc(2*nfunct+fun_cx+1)
814 IF(ismooth < 0) THEN
815 CALL python_call_funct1d(python, -ismooth,tt*fcx1, pvel0)
816 pvel = pvel*pvel0*fcy1/two
817 ELSE
818 pvel=pvel*fcy1*finter(fun_cx,tt*fcx1,npc,tf,dydx)/two
819 ENDIF
820 ENDIF
821 nsign = vel*(nx*xframe(k1,ifra2) + ny*xframe(k2,ifra2) + nz*xframe(k3,ifra2)) ! <Vseg.d,n> where d is flow direction X,Y,Z ! IFRA2=1 is general frame, IFRA2 in [2:NFRAME+1] are for user frames
822 nsign = sign(one,nsign)
823
824 !hydrostatic pressure
825 fx=-aa*(nx/norm)
826 fy=-aa*(ny/norm)
827 fz=-aa*(nz/norm)
828 !drag force
829 fx=fx+pvel*half*nx*nsign
830 fy=fy+pvel*half*ny*nsign
831 fz=fz+pvel*half*nz*nsign
832 !expanding on the 4-node-face
833 fx=fx*third
834 fy=fy*third
835 fz=fz*third
836
837 !NODE_1 : storing nodal force in FSKY array (it will be assembled later)
838 iad = iadc(iloadp(4,nl)+4*(i-1))
839 fsky(1,iad) = fx
840 fsky(2,iad) = fy
841 fsky(3,iad) = fz
842
843 !NODE_2 : storing nodal force in FSKY array (it will be assembled later)
844 iad = iadc(iloadp(4,nl)+4*(i-1)+1)
845 fsky(1,iad) = fx
846 fsky(2,iad) = fy
847 fsky(3,iad) = fz
848
849 !NODE_3 : storing nodal force in FSKY array (it will be assembled later)
850 iad = iadc(iloadp(4,nl)+4*(i-1)+2)
851 fsky(1,iad) = fx
852 fsky(2,iad) = fy
853 fsky(3,iad) = fz
854
855 !external force work
856 wfextt=wfextt+dt1*(fx*(v(1,n1)+v(1,n2)+v(1,n3))
857 + +fy*(v(2,n1)+v(2,n2)+v(2,n3))
858 + +fz*(v(3,n1)+v(3,n2)+v(3,n3)))
859
860 IF(ianim > 0) THEN
861#include "lockon.inc"
862 fext(1,n1) = fext(1,n1)+fx
863 fext(2,n1) = fext(2,n1)+fy
864 fext(3,n1) = fext(3,n1)+fz
865 !
866 fext(1,n2) = fext(1,n2)+fx
867 fext(2,n2) = fext(2,n2)+fy
868 fext(3,n2) = fext(3,n2)+fz
869 !
870 fext(1,n3) = fext(1,n3)+fx
871 fext(2,n3) = fext(2,n3)+fy
872 fext(3,n3) = fext(3,n3)+fz
873#include "lockoff.inc"
874 ENDIF
875C
876 IF(th_surf%LOADP_FLAG >0 ) THEN
877
878 !mean velocity (/TH/SURF)
879 vn = v(1,n1)*nx + v(2,n1)*ny + v(3,n1)*nz
880 vn = vn + v(1,n2)*nx + v(2,n2)*ny + v(3,n2)*nz
881 vn = vn + v(1,n3)*nx + v(2,n3)*ny + v(3,n3)*nz
882 vn = third * vn
883
884 !mean acceleration (/TH/SURF)
885 an = fx*nx + fy*ny + fz*nz
886 mass = third * ( ms(n1)+ms(n2)+ms(n3) )
887 an = an / mass
888
889 nsegpl = nsegpl + 1
890 DO ns=th_surf%LOADP_KSEGS(nsegpl) +1,th_surf%LOADP_KSEGS(nsegpl+1)
891#include "lockon.inc"
892 ksurf = th_surf%LOADP_SEGS(ns)
893 th_surf%CHANNELS(3,ksurf)= th_surf%CHANNELS(3,ksurf) + area*vn ! velocity
894 th_surf%CHANNELS(4,ksurf)= th_surf%CHANNELS(4,ksurf) + area*pa ! pressure
895 th_surf%CHANNELS(5,ksurf)= th_surf%CHANNELS(5,ksurf) + area
896#include "lockoff.inc"
897 ENDDO
898 ENDIF
899C
900 ENDIF
901 enddo!next I
902!$OMP END DO
903 enddo!next NL
904
905!$OMP ATOMIC
906 wfext = wfext + wfextt
907
908 ELSE
909
910C-----------------------------------
911C PART-3. PARITH/ON, vectorial code
912C-----------------------------------
913
914 DO nl=1,nloadp_f !loop over /LOAD/PFLUID options
915
916 !user parameters (Defined in Starter while reading input file hm_read_pblast.F, and transmitted to Engine with Restart file)
917 fun_hsp = iloadp(7,nl)
918 dir_hsp = iloadp(8,nl)
919 ifra1 = iloadp(9,nl)
920 fcy = rload(1,nl)
921 fcx = rload(2,nl)
922 fun_cx = iloadp(10,nl)
923 fcy1 = rload(3,nl)
924 fcx1 = rload(4,nl)
925 fun_vel = iloadp(11,nl)
926 fcy2 = rload(5,nl)
927 fcx2 = rload(6,nl)
928 dir_vel = max(iloadp(12,nl),1)! To avoid a check bound issue when the velocity options is useless
929 ifra2 = iloadp(13,nl)
930 isens = 0
931 xsens = one
932
933 !possible sensor
934 DO k=1,nsensor
935 IF(iloadp(6,nl) == sensor_tab(k)%SENS_ID) isens=k
936 ENDDO
937 IF(isens == 0)THEN
938 ts=tt
939 ELSE
940 ts = tt- sensor_tab(isens)%TSTART
941 IF(ts < zero) cycle
942 ENDIF
943!$OMP DO SCHEDULE(GUIDED,MVSIZ)
944 DO i = 1,iloadp(1,nl)/4
945
946 !nodes of the structural face
947 n1=lloadp(iloadp(4,nl)+4*(i-1))
948 n2=lloadp(iloadp(4,nl)+4*(i-1)+1)
949 n3=lloadp(iloadp(4,nl)+4*(i-1)+2)
950 n4=lloadp(iloadp(4,nl)+4*(i-1)+3)
951
952 aa = zero
953 vel = zero
954 pvel = zero
955
956 IF(n4 /= 0 .AND. n1 /= n2 .AND. n1 /= n3 .AND. n1 /= n4 .AND. n2 /= n3 .AND. n2 /= n4 .AND. n3 /= n4 )THEN
957
958 !hydrostatic pressure
959 k1=3*dir_hsp-2
960 k2=3*dir_hsp-1
961 k3=3*dir_hsp
962 IF(fun_hsp /= 0) THEN
963 centroid_depth = (xframe(k1,ifra1)*(x(1,n1)+x(1,n2)+x(1,n3)+x(1,n4))/four)
964 centroid_depth = centroid_depth+(xframe(k2,ifra1)*(x(2,n1)+x(2,n2)+x(2,n3)+x(2,n4))/four)
965 centroid_depth = centroid_depth+(xframe(k3,ifra1)*(x(3,n1)+x(3,n2)+x(3,n3)+x(3,n4))/four)
966 ismooth = 0
967 IF (fun_hsp > 0) ismooth = npc(2*nfunct+fun_hsp+1)
968 IF(ismooth < 0) THEN
969 CALL python_call_funct1d(python, -ismooth,(centroid_depth-xframe(9+dir_hsp,ifra1))*fcx , aa)
970 aa = aa*fcy
971 ELSE
972 aa = fcy*finter(fun_hsp,(centroid_depth-xframe(9+dir_hsp,ifra1))*fcx,npc,tf,dydx)
973 ENDIF
974 ENDIF
975
976 !normal vector
977 nx = (x(2,n3)-x(2,n1))*(x(3,n4)-x(3,n2)) - (x(3,n3)-x(3,n1))*(x(2,n4)-x(2,n2))
978 ny = (x(3,n3)-x(3,n1))*(x(1,n4)-x(1,n2)) - (x(1,n3)-x(1,n1))*(x(3,n4)-x(3,n2))
979 nz = (x(1,n3)-x(1,n1))*(x(2,n4)-x(2,n2)) - (x(2,n3)-x(2,n1))*(x(1,n4)-x(1,n2))
980 norm = sqrt(nx*nx+ny*ny+nz*nz)
981 area = half * norm
982 pa = aa
983 aa = aa * area
984
985 !structural face velocity
986 k1=3*dir_vel-2
987 k2=3*dir_vel-1
988 k3=3*dir_vel
989 vseg= (xframe(k1,ifra2)*(v(1,n1) + v(1,n2) + v(1,n3) + v(1,n4)) /four)
990 vseg=vseg+(xframe(k2,ifra2)*(v(2,n1) + v(2,n2) + v(2,n3) + v(2,n4)) /four)
991 vseg=vseg+(xframe(k3,ifra2)*(v(3,n1) + v(3,n2) + v(3,n3) + v(3,n4)) /four)
992
993 !VEL,PVEL and NSIGN
994 IF(fun_vel /= 0)THEN
995 ismooth = 0
996 IF (fun_vel > 0) ismooth = npc(2*nfunct+fun_vel+1)
997 IF(ismooth < 0) THEN
998 CALL python_call_funct1d(python, -ismooth, tt*fcx2 ,vel )
999 vel = vel*fcy2 - vseg
1000 ELSE
1001 vel = fcy2*finter(fun_vel,tt*fcx2,npc,tf,dydx) - vseg
1002 ENDIF
1003 ELSE
1004 vel = -vseg
1005 ENDIF
1006 IF(fun_cx /= 0) THEN
1007 pvel = (-nx/norm*vel*xframe(k1,ifra2)- (ny/norm)*vel*xframe(k2,ifra2)-(nz/norm)*vel*xframe(k3,ifra2))**2
1008 ismooth = 0
1009 IF (fun_cx > 0) ismooth = npc(2*nfunct+fun_cx+1)
1010 IF(ismooth < 0) THEN
1011 CALL python_call_funct1d(python, -ismooth,tt*fcx1, pvel0)
1012 pvel = pvel*pvel0*fcy1/two
1013 ELSE
1014 pvel = pvel * fcy1*finter(fun_cx,tt*fcx1,npc,tf,dydx)/two
1015 ENDIF
1016
1017 ENDIF
1018 nsign = vel*(nx*xframe(k1,ifra2) + ny*xframe(k2,ifra2) + nz*xframe(k3,ifra2)) ! <Vseg.d,n> where d is flow direction X,Y,Z ! IFRA2=1 is general frame, IFRA2 in [2:NFRAME+1] are for user frames
1019 nsign = sign(one,nsign)
1020
1021 !hydrostatic pressure
1022 fx=-aa*(nx/norm)
1023 fy=-aa*(ny/norm)
1024 fz=-aa*(nz/norm)
1025 !drag force
1026 fx=fx+pvel*half*nx*nsign
1027 fy=fy+pvel*half*ny*nsign
1028 fz=fz+pvel*half*nz*nsign
1029 !expanding on the 4-node-face
1030 fx=fx*fourth
1031 fy=fy*fourth
1032 fz=fz*fourth
1033
1034 !NODE_1 : storing nodal force in FSKY array (it will be assembled later)
1035 iad = iadc(iloadp(4,nl)+4*(i-1))
1036 fskyv(iad,1) = fx
1037 fskyv(iad,2) = fy
1038 fskyv(iad,3) = fz
1039
1040 !NODE_2 : storing nodal force in FSKY array (it will be assembled later)
1041 iad = iadc(iloadp(4,nl)+4*(i-1)+1)
1042 fskyv(iad,1) = fx
1043 fskyv(iad,2) = fy
1044 fskyv(iad,3) = fz
1045
1046 !NODE_3 : storing nodal force in FSKY array (it will be assembled later)
1047 iad = iadc(iloadp(4,nl)+4*(i-1)+2)
1048 fskyv(iad,1) = fx
1049 fskyv(iad,2) = fy
1050 fskyv(iad,3) = fz
1051
1052 !NODE_4 : storing nodal force in FSKY array (it will be assembled later)
1053 iad = iadc(iloadp(4,nl)+4*(i-1)+3)
1054 fskyv(iad,1) = fx
1055 fskyv(iad,2) = fy
1056 fskyv(iad,3) = fz
1057
1058 !external force work
1059 wfextt=wfextt+dt1*(fx*(v(1,n1)+v(1,n2)+v(1,n3)+v(1,n4))
1060 + +fy*(v(2,n1)+v(2,n2)+v(2,n3)+v(2,n4))
1061 + +fz*(v(3,n1)+v(3,n2)+v(3,n3)+v(3,n4)))
1062
1063 IF(ianim > 0) THEN
1064#include "lockon.inc"
1065 fext(1,n1) = fext(1,n1)+fx
1066 fext(2,n1) = fext(2,n1)+fy
1067 fext(3,n1) = fext(3,n1)+fz
1068 !
1069 fext(1,n2) = fext(1,n2)+fx
1070 fext(2,n2) = fext(2,n2)+fy
1071 fext(3,n2) = fext(3,n2)+fz
1072 !
1073 fext(1,n3) = fext(1,n3)+fx
1074 fext(2,n3) = fext(2,n3)+fy
1075 fext(3,n3) = fext(3,n3)+fz
1076 !
1077 fext(1,n4) = fext(1,n4)+fx
1078 fext(2,n4) = fext(2,n4)+fy
1079 fext(3,n4) = fext(3,n4)+fz
1080#include "lockoff.inc"
1081 ENDIF
1082C
1083 IF(th_surf%LOADP_FLAG >0 ) THEN
1084
1085 !mean velocity (/TH/SURF)
1086 vn = v(1,n1)*nx + v(2,n1)*ny + v(3,n1)*nz
1087 vn = vn + v(1,n2)*nx + v(2,n2)*ny + v(3,n2)*nz
1088 vn = vn + v(1,n3)*nx + v(2,n3)*ny + v(3,n3)*nz
1089 vn = vn + v(1,n4)*nx + v(2,n4)*ny + v(3,n4)*nz
1090 vn = vn / four
1091
1092 !mean acceleration (/TH/SURF)
1093 an = fx*nx + fy*ny + fz*nz
1094 mass = ( ms(n1)+ms(n2)+ms(n3)+ms(n4) ) / four
1095 an = an / mass
1096
1097 nsegpl = nsegpl + 1
1098 DO ns=th_surf%LOADP_KSEGS(nsegpl) +1,th_surf%LOADP_KSEGS(nsegpl+1)
1099#include "lockon.inc"
1100 ksurf = th_surf%LOADP_SEGS(ns)
1101 th_surf%CHANNELS(3,ksurf)= th_surf%CHANNELS(3,ksurf) + area*vn ! velocity
1102 th_surf%CHANNELS(4,ksurf)= th_surf%CHANNELS(4,ksurf) + area*pa ! pressure
1103 th_surf%CHANNELS(5,ksurf)= th_surf%CHANNELS(5,ksurf) + area
1104#include "lockoff.inc"
1105 ENDDO
1106 ENDIF
1107
1108 ELSE
1109
1110 !---TRIANGLE--------------------------------------------
1111 IF(n1 == n2)THEN
1112 n2 = n3
1113 n3 = n4
1114 n4 = 0
1115 ELSEIF(n1 == n3)THEN
1116 n3 = n4
1117 n4 = 0
1118 ELSEIF(n1 == n4)THEN
1119 n4 = 0
1120 ELSEIF(n2 == n3)THEN
1121 n3 = n4
1122 n4 = 0
1123 ELSEIF(n2 == n4)THEN
1124 n2 = n3
1125 n3 = n4
1126 n4 = 0
1127 ELSEIF(n3 == n4)THEN
1128 n4 = 0
1129 ENDIF
1130
1131 !hydrostatic load
1132 k1=3*dir_hsp-2
1133 k2=3*dir_hsp-1
1134 k3=3*dir_hsp
1135 IF(fun_hsp /= 0)THEN
1136 centroid_depth = (xframe(k1,ifra1)*(x(1,n1)+x(1,n2)+x(1,n3))/three)
1137 centroid_depth = centroid_depth+(xframe(k2,ifra1)*(x(2,n1)+x(2,n2)+x(2,n3))/three)
1138 centroid_depth = centroid_depth+(xframe(k3,ifra1)*(x(3,n1)+x(3,n2)+x(3,n3))/three)
1139 ismooth = 0
1140 IF (fun_hsp > 0) ismooth = npc(2*nfunct+fun_hsp+1)
1141 IF(ismooth < 0) THEN
1142 CALL python_call_funct1d(python, -ismooth,(centroid_depth-xframe(9+dir_hsp,ifra1))*fcx , aa)
1143 aa = aa*fcy
1144 ELSE
1145 aa = fcy*finter(fun_hsp,(centroid_depth-xframe(9+dir_hsp,ifra1))*fcx,npc,tf,dydx)
1146 ENDIF
1147
1148 ENDIF
1149
1150 !normal vector
1151 nx = (x(2,n3)-x(2,n1))*(x(3,n3)-x(3,n2)) - (x(3,n3)-x(3,n1))*(x(2,n3)-x(2,n2))
1152 ny = (x(3,n3)-x(3,n1))*(x(1,n3)-x(1,n2)) - (x(1,n3)-x(1,n1))*(x(3,n3)-x(3,n2))
1153 nz = (x(1,n3)-x(1,n1))*(x(2,n3)-x(2,n2)) - (x(2,n3)-x(2,n1))*(x(1,n3)-x(1,n2))
1154 norm = sqrt(nx*nx+ny*ny+nz*nz)
1155 area = half * norm
1156 pa = aa
1157 aa = aa * area
1158
1159 !structural face velocity
1160 k1=3*dir_vel-2
1161 k2=3*dir_vel-1
1162 k3=3*dir_vel
1163 vseg= (xframe(k1,ifra2)* (v(1,n1) + v(1,n2) + v(1,n3)) /three)
1164 vseg=vseg+(xframe(k2,ifra2)* (v(2,n1) + v(2,n2) + v(2,n3)) /three)
1165 vseg=vseg+(xframe(k3,ifra2)* (v(3,n1) + v(3,n2) + v(3,n3)) /three)
1166
1167 !VEL,PVEL, and NSIGN
1168 IF(fun_vel /= 0)THEN
1169 ismooth = 0
1170 IF (fun_vel > 0) ismooth = npc(2*nfunct+fun_vel+1)
1171 IF(ismooth < 0) THEN
1172 CALL python_call_funct1d(python, -ismooth, tt*fcx2 ,vel )
1173 vel = vel*fcy2 - vseg
1174 ELSE
1175 vel = fcy2*finter(fun_vel,tt*fcx2,npc,tf,dydx) - vseg
1176 ENDIF
1177 ELSE
1178 vel = -vseg
1179 ENDIF
1180 IF(fun_cx /= 0)THEN
1181 pvel = ( (-(nx/norm)*vel*xframe(k1,ifra2)-(ny/norm)*vel*xframe(k2,ifra2)-(nz/norm)*vel*xframe(k3,ifra2))**2 )
1182 ismooth = 0
1183 IF (fun_cx > 0) ismooth = npc(2*nfunct+fun_cx+1)
1184 IF(ismooth < 0) THEN
1185 CALL python_call_funct1d(python, -ismooth,tt*fcx1, pvel0)
1186 pvel = pvel*pvel0*fcy1/two
1187 ELSE
1188 pvel = pvel * fcy1* finter(fun_cx,tt*fcx1,npc,tf,dydx)/two
1189 ENDIF
1190
1191 ENDIF
1192 nsign = vel*(nx*xframe(k1,ifra2) + ny*xframe(k2,ifra2) + nz*xframe(k3,ifra2)) ! <Vseg.d,n> where d is flow direction X,Y,Z ! IFRA2=1 is general frame, IFRA2 in [2:NFRAME+1] are for user frames
1193 nsign = sign(one,nsign)
1194
1195 !hydrostatic pressure
1196 fx=-aa*(nx/norm)
1197 fy=-aa*(ny/norm)
1198 fz=-aa*(nz/norm)
1199 !drag force
1200 fx=fx+pvel*half*nx*nsign
1201 fy=fy+pvel*half*ny*nsign
1202 fz=fz+pvel*half*nz*nsign
1203 !expanding on the 4-node-face
1204 fx=fx*third
1205 fy=fy*third
1206 fz=fz*third
1207
1208 !NODE_1 : storing nodal force in FSKY array (it will be assembled later)
1209 iad = iadc(iloadp(4,nl)+4*(i-1))
1210 fskyv(iad,1) = fx
1211 fskyv(iad,2) = fy
1212 fskyv(iad,3) = fz
1213 IF(ianim > 0) THEN
1214#include "lockon.inc"
1215 fext(1,n1) = fext(1,n1)+fx
1216 fext(2,n1) = fext(2,n1)+fy
1217 fext(3,n1) = fext(3,n1)+fz
1218#include "lockoff.inc"
1219 ENDIF
1220
1221 !NODE_2 : storing nodal force in FSKY array (it will be assembled later)
1222 iad = iadc(iloadp(4,nl)+4*(i-1)+1)
1223 fskyv(iad,1) = fx
1224 fskyv(iad,2) = fy
1225 fskyv(iad,3) = fz
1226 IF(ianim > 0) THEN
1227#include "lockon.inc"
1228 fext(1,n2) = fext(1,n2)+fx
1229 fext(2,n2) = fext(2,n2)+fy
1230 fext(3,n2) = fext(3,n2)+fz
1231#include "lockoff.inc"
1232 ENDIF
1233
1234 !NODE_3 : storing nodal force in FSKY array (it will be assembled later)
1235 iad = iadc(iloadp(4,nl)+4*(i-1)+2)
1236 fskyv(iad,1) = fx
1237 fskyv(iad,2) = fy
1238 fskyv(iad,3) = fz
1239 IF(ianim > 0) THEN
1240#include "lockon.inc"
1241 fext(1,n3) = fext(1,n3)+fx
1242 fext(2,n3) = fext(2,n3)+fy
1243 fext(3,n3) = fext(3,n3)+fz
1244#include "lockoff.inc"
1245 ENDIF
1246C
1247 IF(th_surf%LOADP_FLAG >0 ) THEN
1248
1249 !mean velocity (/TH/SURF)
1250 vn = v(1,n1)*nx + v(2,n1)*ny + v(3,n1)*nz
1251 vn = vn + v(1,n2)*nx + v(2,n2)*ny + v(3,n2)*nz
1252 vn = vn + v(1,n3)*nx + v(2,n3)*ny + v(3,n3)*nz
1253 vn = third * vn
1254
1255 !mean acceleration (/TH/SURF)
1256 an = fx*nx + fy*ny + fz*nz
1257 mass = third * ( ms(n1)+ms(n2)+ms(n3)+ms(n4) )
1258 an = an / mass
1259
1260 nsegpl = nsegpl + 1
1261 DO ns=th_surf%LOADP_KSEGS(nsegpl) +1,th_surf%LOADP_KSEGS(nsegpl+1)
1262#include "lockon.inc"
1263 ksurf = th_surf%LOADP_SEGS(ns)
1264 th_surf%CHANNELS(3,ksurf)= th_surf%CHANNELS(3,ksurf) + area*vn ! velocity
1265 th_surf%CHANNELS(4,ksurf)= th_surf%CHANNELS(4,ksurf) + area*pa ! pressure
1266 th_surf%CHANNELS(5,ksurf)= th_surf%CHANNELS(5,ksurf) + area
1267#include "lockoff.inc"
1268 ENDDO
1269 ENDIF
1270
1271 !external force work
1272 wfextt=wfextt+dt1*(fx*(v(1,n1)+v(1,n2)+v(1,n3))
1273 + +fy*(v(2,n1)+v(2,n2)+v(2,n3))
1274 + +fz*(v(3,n1)+v(3,n2)+v(3,n3)))
1275 ENDIF
1276 enddo!next I
1277!$OMP END DO
1278 enddo!next NL
1279C
1280!$OMP ATOMIC
1281 wfext = wfext + wfextt
1282 ENDIF ! IF(IVECTOR == 0)
1283
1284
1285 ENDIF
1286C
1287 RETURN
1288 END
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
max
#define max(a, b)
Definition macros.h:21
h3d_mod
Definition h3d_mod.F:308
th_surf_mod
OPTION /TH/SURF outputs of Pressure and Area needed Tabs.
Definition th_surf_mod.F:60
pfluid
subroutine pfluid(iloadp, rload, npc, tf, a, v, x, xframe, ms, nsensor, sensor_tab, wfexc, wfext, iadc, fsky, fskyv, lloadp, fext, h3d_data, th_surf, python)
Definition pfluid.F:45
h3d_mod::h3d_database
Definition h3d_mod.F:435
th_surf_mod::th_surf_
Definition th_surf_mod.F:75