OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
load_pressure.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!|| load_pressure ../engine/source/loads/general/load_pressure/load_pressure.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!|| pblast_mod ../common_source/modules/loads/pblast_mod.F90
37!|| python_funct_mod ../common_source/modules/python_mod.F90
38!|| sensor_mod ../common_source/modules/sensor_mod.F90
39!|| th_surf_mod ../common_source/modules/interfaces/th_surf_mod.F
40!||====================================================================
41 SUBROUTINE load_pressure (ILOADP ,LOADP ,LLOADP ,NPC ,TF ,
42 2 A ,V ,X ,SKEW ,SENSOR_TAB,
43 3 IADC ,FSKY ,FEXT ,TAGNCONT ,NSENSOR ,
44 4 LOADP_HYD_INTER , H3D_DATA , PYTHON,
45 5 NPRESLOAD ,LOADP_TAGDEL,TH_SURF,PBLAST ,WFEXT)
46C-----------------------------------------------
47C M o d u l e s
48C-----------------------------------------------
49 USE python_funct_mod
50 USE h3d_mod
51 USE sensor_mod
52 USE th_surf_mod , ONLY : th_surf_
53 USE pblast_mod
54C-----------------------------------------------
55C I m p l i c i t T y p e s
56C-----------------------------------------------
57#include "implicit_f.inc"
58#include "comlock.inc"
59#include "param_c.inc"
60C-----------------------------------------------
61C C o m m o n B l o c k s
62C-----------------------------------------------
63#include "com04_c.inc"
64!#include "com06_c.inc"
65#include "com08_c.inc"
66#include "impl1_c.inc"
67#include "scr14_c.inc"
68#include "scr16_c.inc"
69#include "parit_c.inc"
70#include "tabsiz_c.inc"
71C-----------------------------------------------
72C E x t e r n a l F u n c t i o n s
73C-----------------------------------------------
74 INTEGER,EXTERNAL :: GET_U_NUMSENS,GET_U_SENS_FPAR,GET_U_SENS_IPAR,GET_U_SENS_VALUE,SET_U_SENS_VALUE
75C-----------------------------------------------,
76C D u m m y A r g u m e n t s
77C-----------------------------------------------
78 TYPE(python_), INTENT(inout) :: PYTHON
79 INTEGER, INTENT(IN) :: NPRESLOAD,NSENSOR
80 INTEGER NPC(*),LLOADP(SLLOADP)
81 INTEGER ILOADP(SIZLOADP,*),TAGNCONT(NLOADP_HYD_INTER,NUMNOD),LOADP_HYD_INTER(NLOADP_HYD),IADC(*)
82 my_real loadp(sloadp,*), tf(*), a(3,*), v(3,*), x(3,*), skew(lskew,*),fsky(8,lsky), fext(3,*)
83 TYPE(h3d_database) :: H3D_DATA
84 INTEGER, INTENT(IN) :: LOADP_TAGDEL(NPRESLOAD)
85 TYPE (sensor_str_) ,DIMENSION(NSENSOR) ,INTENT(IN) :: sensor_tab
86 TYPE (TH_SURF_) , INTENT(INOUT) :: TH_SURF
87 TYPE(pblast_), INTENT(IN) :: PBLAST
88 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
89C-----------------------------------------------
90C L o c a l V a r i a b l e s
91C-----------------------------------------------
92 INTEGER N1, ISK, N2, N3, N4, ISENS,
93 . IAD ,NP ,IFUNC ,NPRES ,NINTERP ,IDIR ,INORM ,NPL ,IANIM,IFLOAD,
94 . SEGCONT ,N ,IADN , IJK ,UP_BOUND,TAGN1,TAGN2,TAGN3,TAGN4,NUMPRESLOAD,
95 . IDEL, NSEGPL, NS, KSURF,NIDXLOAD
96 INTEGER :: IS_TABULATED
97 my_real NX, NY, NZ, AA, FX, FY, FZ, DYDX, TS, WFEXTT,F1, FCX,FCY,NORM,AREA
98 my_real,EXTERNAL :: FINTER
99C---------------------------------------------------------------------
100C D e s c r i p t i o n
101C----------------------------------------------------------------------
102C This subroutine is applying pressure load to
103C a defined surface on dynamic zones which is not in contact
104C---------------------------------------------------------------------
105 wfextt = zero
106 fx = zero
107 fy = zero
108 fz = zero
109 ianim = anim_v(5)+outp_v(5)+h3d_data%N_VECT_FINT+
110 . anim_v(6)+outp_v(6)+h3d_data%N_VECT_FEXT
111C
112 numpresload = 0
113 nsegpl = th_surf%NSEGLOADPF+th_surf%NSEGLOADPB
114 nidxload = nloadp_f+pblast%NLOADP_B
115 IF(iparit==0) THEN
116C code SPMD Parith/OFF or SMP
117 DO np= 1,nloadp_hyd
118
119 npres = iloadp(1,nidxload+np)
120 ifunc = iloadp(3,nidxload+np)
121 iad = iloadp(4,nidxload+np)
122 ninterp = iloadp(5,nidxload+np)
123 idir = iloadp(6,nidxload+np)
124 isens = iloadp(7,nidxload+np)
125 isk = iloadp(8,nidxload+np)
126 inorm = iloadp(9,nidxload+np)
127 ifload = iloadp(10,nidxload+np)
128 fcy = loadp(1,nidxload+np)
129 fcx = loadp(2,nidxload+np)
130C
131 IF(isens==0)THEN
132 ts=tt
133 ELSE
134 ts = tt-sensor_tab(isens)%TSTART
135 IF(ts<zero)cycle
136 ENDIF
137
138
139 is_tabulated = npc(2*nfunct+ifunc+1)
140 IF(is_tabulated >= 0) THEN
141 f1 = finter(ifunc,ts*fcx,npc,tf,dydx)
142 ELSE
143 is_tabulated = -is_tabulated
144 CALL python_call_funct1d(python, is_tabulated,ts*fcx, f1)
145 ENDIF
146
147 aa = fcy*f1
148
149 DO n=1, npres/4
150C
151 n1 = lloadp(iad+4*(n-1))
152 n2 = lloadp(iad+4*(n-1)+1)
153 n3 = lloadp(iad+4*(n-1)+2)
154 n4 = lloadp(iad+4*(n-1)+3)
155 numpresload = numpresload + 1
156 nsegpl = nsegpl + 1
157
158 idel = loadp_tagdel(numpresload)
159
160 IF(idel > 0 ) cycle ! SEGMENT DELETED
161
162C----------------
163C Check if segment is in contact
164C----------------
165 segcont = 0
166
167 tagn1 = 0
168 tagn2 = 0
169 tagn3 = 0
170 tagn4 = 0
171
172 IF(ninterp > 0 ) THEN
173 npl = loadp_hyd_inter(np)
174 IF(n4/=0) THEN
175 segcont = tagncont(npl,n1) + tagncont(npl,n2) +
176 . tagncont(npl,n3)+tagncont(npl,n4)
177 IF(segcont >= 2 .AND.ifload==1) THEN
178 segcont = 1
179 ELSEIF(segcont <= 1.AND.ifload==2) THEN
180 segcont = 1
181 ELSE
182 segcont = 0
183 ENDIF
184 ELSE
185 segcont = tagncont(npl,n1) + tagncont(npl,n2) +
186 . tagncont(npl,n3)
187 IF(segcont >= 2 .AND.ifload==1) THEN
188 segcont = 1
189 ELSEIF(segcont <= 1.AND.ifload==2) THEN
190 segcont = 1
191 ELSE
192 segcont = 0
193 ENDIF
194 ENDIF
195 ENDIF
196
197 IF(segcont == 0) THEN
198
199 IF(n4/=0)THEN
200
201 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))
202 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))
203 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))
204 IF(inorm == 1) THEN
205 fx = aa*nx*one_over_8
206 fy = aa*ny*one_over_8
207 fz = aa*nz*one_over_8
208 ELSEIF(inorm==2) THEN
209 norm = sqrt(nx*nx+ny*ny+nz*nz)
210 IF(isk == 0) THEN
211 IF(idir == 1 ) THEN
212 fx = aa*one_over_8*norm
213 fy = zero
214 fz = zero
215 ELSEIF(idir==2) THEN
216 fy = aa*one_over_8*norm
217 fx = zero
218 fz = zero
219 ELSEIF(idir==3) THEN
220 fz = aa*one_over_8*norm
221 fx = zero
222 fy = zero
223 ENDIF
224 ELSE
225 IF(idir == 1 ) THEN
226 fx = aa*one_over_8*norm*skew(1,isk)
227 fy = zero
228 fz = zero
229 ELSEIF(idir==2) THEN
230 fy = aa*one_over_8*norm*skew(2,isk)
231 fx = zero
232 fz = zero
233 ELSEIF(idir==3) THEN
234 fz = aa*one_over_8*norm*skew(3,isk)
235 fx = zero
236 fy = zero
237 ENDIF
238 ENDIF
239 ELSEIF(inorm==3) THEN
240
241 IF(isk == 0) THEN
242 IF(idir == 1 ) THEN
243 fx = aa*one_over_8*nx
244 fy = zero
245 fz = zero
246 ELSEIF(idir==2) THEN
247 fy = aa*one_over_8*ny
248 fx = zero
249 fz = zero
250 ELSEIF(idir==3) THEN
251 fz = aa*one_over_8*nz
252 fx = zero
253 fy = zero
254 ENDIF
255 ELSE
256 IF(idir == 1 ) THEN
257 fx = aa*one_over_8*skew(1,isk)*nx
258 fy = zero
259 fz = zero
260 ELSEIF(idir==2) THEN
261 fy = aa*one_over_8*skew(2,isk)*ny
262 fx = zero
263 fz = zero
264 ELSEIF(idir==3) THEN
265 fz = aa*one_over_8*skew(3,isk)*nz
266 fx = zero
267 fy = zero
268 ENDIF
269 ENDIF
270 ENDIF
271C
272 IF(tagn1 == 0) THEN
273 a(1,n1)=a(1,n1)+fx
274 a(2,n1)=a(2,n1)+fy
275 a(3,n1)=a(3,n1)+fz
276 IF(ianim >0 .AND.impl_s==0) THEN
277 fext(1,n1) = fext(1,n1)+fx
278 fext(2,n1) = fext(2,n1)+fy
279 fext(3,n1) = fext(3,n1)+fz
280 ENDIF
281 ENDIF
282C
283 IF(tagn2 == 0) THEN
284 a(1,n2)=a(1,n2)+fx
285 a(2,n2)=a(2,n2)+fy
286 a(3,n2)=a(3,n2)+fz
287 IF(ianim >0 .AND.impl_s==0) THEN
288 fext(1,n2) = fext(1,n2)+fx
289 fext(2,n2) = fext(2,n2)+fy
290 fext(3,n2) = fext(3,n2)+fz
291 ENDIF
292 ENDIF
293C
294 IF(tagn3 == 0) THEN
295 a(1,n3)=a(1,n3)+fx
296 a(2,n3)=a(2,n3)+fy
297 a(3,n3)=a(3,n3)+fz
298 IF(ianim >0 .AND.impl_s==0) THEN
299 fext(1,n3) = fext(1,n3)+fx
300 fext(2,n3) = fext(2,n3)+fy
301 fext(3,n3) = fext(3,n3)+fz
302 ENDIF
303 ENDIF
304C
305 IF(tagn4 == 0) THEN
306 a(1,n4)=a(1,n4)+fx
307 a(2,n4)=a(2,n4)+fy
308 a(3,n4)=a(3,n4)+fz
309 IF(ianim >0 .AND.impl_s==0) THEN
310 fext(1,n4) = fext(1,n4)+fx
311 fext(2,n4) = fext(2,n4)+fy
312 fext(3,n4) = fext(3,n4)+fz
313 ENDIF
314 ENDIF
315C
316 IF(th_surf%LOADP_FLAG >0 ) THEN
317 area = half*sqrt(nx*nx+ny*ny+nz*nz)
318 DO ns=th_surf%LOADP_KSEGS(nsegpl) +1,th_surf%LOADP_KSEGS(nsegpl+1)
319 ksurf = th_surf%LOADP_SEGS(ns)
320 th_surf%channels(4,ksurf)= th_surf%channels(4,ksurf) + area*aa
321 th_surf%channels(5,ksurf)= th_surf%channels(5,ksurf) + area
322 ENDDO
323 ENDIF
324C
325 wfextt=wfextt+dt1*(fx*(v(1,n1)+v(1,n2)+v(1,n3)+v(1,n4))
326 + +fy*(v(2,n1)+v(2,n2)+v(2,n3)+v(2,n4))
327 + +fz*(v(3,n1)+v(3,n2)+v(3,n3)+v(3,n4)))
328C
329 ELSE
330 ! true triangles.
331 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))
332 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))
333 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))
334
335 IF(inorm == 1) THEN
336 fx = aa*nx*one_over_6
337 fy = aa*ny*one_over_6
338 fz = aa*nz*one_over_6
339 ELSEIF(inorm==2) THEN
340 norm = sqrt(nx*nx+ny*ny+nz*nz)
341 IF(isk == 0) THEN
342 IF(idir == 1 ) THEN
343 fx = aa*one_over_6*norm
344 fy = zero
345 fz = zero
346 ELSEIF(idir==2) THEN
347 fy = aa*one_over_6*norm
348 fx = zero
349 fz = zero
350 ELSEIF(idir==3) THEN
351 fz = aa*one_over_6*norm
352 fx = zero
353 fy = zero
354 ENDIF
355 ELSE
356 IF(idir == 1 ) THEN
357 fx = aa*one_over_6*norm*skew(1,isk)
358 fy = zero
359 fz = zero
360 ELSEIF(idir==2) THEN
361 fy = aa*one_over_6*norm*skew(2,isk)
362 fx = zero
363 fz = zero
364 ELSEIF(idir==3) THEN
365 fz = aa*one_over_6*norm*skew(3,isk)
366 fx = zero
367 fy = zero
368 ENDIF
369 ENDIF
370 ELSEIF(inorm==3) THEN
371
372 IF(isk == 0) THEN
373 IF(idir == 1 ) THEN
374 IF(nx /= zero) fx = aa*one_over_6*nx
375 fy = zero
376 fz = zero
377 ELSEIF(idir==2) THEN
378 IF(ny /= zero) fy = aa*one_over_6*ny
379 fx = zero
380 fz = zero
381 ELSEIF(idir==3) THEN
382 IF(nz /= zero) fz = aa*one_over_6*nz
383 fx = zero
384 fy = zero
385 ENDIF
386 ELSE
387 IF(idir == 1 ) THEN
388 fx = aa*one_over_6*skew(1,isk)*nx
389 fy = zero
390 fz = zero
391 ELSEIF(idir==2) THEN
392 fy = aa*one_over_6*skew(2,isk)*ny
393 fx = zero
394 fz = zero
395 ELSEIF(idir==3) THEN
396 fz = aa*one_over_6*skew(3,isk)*nz
397 fx = zero
398 fy = zero
399 ENDIF
400 ENDIF
401 ENDIF
402C
403 IF(tagn1 == 0) THEN
404 a(1,n1)=a(1,n1)+fx
405 a(2,n1)=a(2,n1)+fy
406 a(3,n1)=a(3,n1)+fz
407 IF(ianim >0 .AND.impl_s==0) THEN
408 fext(1,n1) = fext(1,n1)+fx
409 fext(2,n1) = fext(2,n1)+fy
410 fext(3,n1) = fext(3,n1)+fz
411 ENDIF
412
413 ENDIF
414C
415 IF(tagn2 == 0) THEN
416 a(1,n2)=a(1,n2)+fx
417 a(2,n2)=a(2,n2)+fy
418 a(3,n2)=a(3,n2)+fz
419 IF(ianim >0 .AND.impl_s==0) THEN
420 fext(1,n2) = fext(1,n2)+fx
421 fext(2,n2) = fext(2,n2)+fy
422 fext(3,n2) = fext(3,n2)+fz
423 ENDIF
424 ENDIF
425C
426 IF(tagn3 == 0) THEN
427 a(1,n3)=a(1,n3)+fx
428 a(2,n3)=a(2,n3)+fy
429 a(3,n3)=a(3,n3)+fz
430 IF(ianim >0 .AND.impl_s==0) THEN
431 fext(1,n3) = fext(1,n3)+fx
432 fext(2,n3) = fext(2,n3)+fy
433 fext(3,n3) = fext(3,n3)+fz
434 ENDIF
435 ENDIF
436C
437 IF(th_surf%LOADP_FLAG >0 ) THEN
438 area = half*sqrt(nx*nx+ny*ny+nz*nz)
439 DO ns=th_surf%LOADP_KSEGS(nsegpl) +1,th_surf%LOADP_KSEGS(nsegpl+1)
440 ksurf = th_surf%LOADP_SEGS(ns)
441 th_surf%channels(4,ksurf)= th_surf%channels(4,ksurf) + area*aa
442 th_surf%channels(5,ksurf)= th_surf%channels(5,ksurf) + area
443 ENDDO
444 ENDIF
445C
446 wfextt=wfextt+dt1*(fx*(v(1,n1)+v(1,n2)+v(1,n3))
447 + + fy*(v(2,n1)+v(2,n2)+v(2,n3))
448 + + fz*(v(3,n1)+v(3,n2)+v(3,n3)))
449 ENDIF
450C
451 ENDIF
452C
453 ENDDO
454C
455 ENDDO
456C
457!$OMP ATOMIC
458 wfext = wfext + wfextt
459C
460 ELSE
461C-------------------------
462C code SPMD Parith/ON
463C-------------------------
464
465 DO np= 1,nloadp_hyd
466
467 npres = iloadp(1,nidxload+np)
468 ifunc = iloadp(3,nidxload+np)
469 iad = iloadp(4,nidxload+np)
470 ninterp = iloadp(5,nidxload+np)
471 idir = iloadp(6,nidxload+np)
472 isens = iloadp(7,nidxload+np)
473 isk = iloadp(8,nidxload+np)
474 inorm = iloadp(9,nidxload+np)
475 ifload = iloadp(10,nidxload+np)
476 fcy = loadp(1,nidxload+np)
477 fcx = loadp(2,nidxload+np)
478C
479 IF(isens==0)THEN
480 ts=tt
481 ELSE
482 ts = tt-sensor_tab(isens)%TSTART
483 IF(ts<zero)cycle
484 ENDIF
485
486 ! -------------
487 ! flush fsky array to 0.
488 DO n = 1,npres/4
489 n1=lloadp(iad+4*(n-1))
490 n2=lloadp(iad+4*(n-1)+1)
491 n3=lloadp(iad+4*(n-1)+2)
492 n4=lloadp(iad+4*(n-1)+3)
493
494 IF(n4/=0 .AND. n1/=n2 .AND. n1/=n3 .AND. n1/=n4 .AND.
495 . n2/=n3 .AND. n2/=n4 .AND. n3/=n4 )THEN
496 up_bound=4
497 ELSE
498 up_bound=3
499 ENDIF
500 DO ijk=1,up_bound
501 iadn = iadc(iad + 4*(n-1)+(ijk-1))
502 fsky(1:3,iadn) = zero
503 ENDDO
504 ENDDO
505 ! -------------
506
507 is_tabulated = npc(2*nfunct+ifunc+1)
508 IF(is_tabulated >= 0) THEN
509 f1 = finter(ifunc,ts*fcx,npc,tf,dydx)
510 ELSE
511 is_tabulated = -is_tabulated
512 CALL python_call_funct1d(python, is_tabulated,ts*fcx, f1)
513 ENDIF
514
515 aa = fcy*f1
516
517 DO n=1, npres/4
518 n1 = lloadp(iad+4*(n-1))
519 n2 = lloadp(iad+4*(n-1)+1)
520 n3 = lloadp(iad+4*(n-1)+2)
521 n4 = lloadp(iad+4*(n-1)+3)
522
523 numpresload = numpresload + 1
524 nsegpl = nsegpl + 1
525
526 idel = loadp_tagdel(numpresload)
527
528 IF(idel > 0 ) cycle ! SEGMENT DELETED
529
530C----------------
531C Check if segment is in contact
532C----------------
533 segcont = 0
534
535 tagn1 = 0
536 tagn2 = 0
537 tagn3 = 0
538 tagn4 = 0
539
540 IF(ninterp > 0 ) THEN
541 npl = loadp_hyd_inter(np)
542 IF(n4/=0) THEN
543 segcont = tagncont(npl,n1) + tagncont(npl,n2) +
544 . tagncont(npl,n3)+tagncont(npl,n4)
545 IF(segcont >= 2 .AND.ifload==1) THEN
546 segcont = 1
547 ELSEIF(segcont <= 1.AND.ifload==2) THEN
548 segcont = 1
549 ELSE
550 segcont = 0
551 ENDIF
552 ELSE
553 segcont = tagncont(npl,n1) + tagncont(npl,n2) +
554 . tagncont(npl,n3)
555 IF(segcont >= 2 .AND.ifload==1) THEN
556 segcont = 1
557 ELSEIF(segcont <= 1.AND.ifload==2) THEN
558 segcont = 1
559 ELSE
560 segcont = 0
561 ENDIF
562 ENDIF
563 ENDIF
564
565 IF(segcont == 0) THEN
566
567 IF(n4/=0)THEN
568
569 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))
570 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))
571 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))
572 IF(inorm == 1) THEN
573 fx = aa*nx*one_over_8
574 fy = aa*ny*one_over_8
575 fz = aa*nz*one_over_8
576 ELSEIF(inorm==2) THEN
577 norm = sqrt(nx*nx+ny*ny+nz*nz)
578 IF(isk == 0) THEN
579 IF(idir == 1 ) THEN
580 fx = aa*one_over_8*norm
581 fy = zero
582 fz = zero
583 ELSEIF(idir==2) THEN
584 fy = aa*one_over_8*norm
585 fx = zero
586 fz = zero
587 ELSEIF(idir==3) THEN
588 fz = aa*one_over_8*norm
589 fx = zero
590 fy = zero
591 ENDIF
592 ELSE
593 IF(idir == 1 ) THEN
594 fx = aa*one_over_8*norm*skew(1,isk)
595 fy = zero
596 fz = zero
597 ELSEIF(idir==2) THEN
598 fy = aa*one_over_8*norm*skew(2,isk)
599 fx = zero
600 fz = zero
601 ELSEIF(idir==3) THEN
602 fz = aa*one_over_8*norm*skew(3,isk)
603 fx = zero
604 fy = zero
605 ENDIF
606 ENDIF
607 ELSEIF(inorm==3) THEN
608
609 IF(isk == 0) THEN
610 IF(idir == 1 ) THEN
611 fx = aa*one_over_8*nx
612 fy = zero
613 fz = zero
614 ELSEIF(idir==2) THEN
615 fy = aa*one_over_8*ny
616 fx = zero
617 fz = zero
618 ELSEIF(idir==3) THEN
619 fz = aa*one_over_8*nz
620 fx = zero
621 fy = zero
622 ENDIF
623 ELSE
624 IF(idir == 1 ) THEN
625 fx = aa*one_over_8*skew(1,isk)*nx
626 fy = zero
627 fz = zero
628 ELSEIF(idir==2) THEN
629 fy = aa*one_over_8*skew(2,isk)*ny
630 fx = zero
631 fz = zero
632 ELSEIF(idir==3) THEN
633 fz = aa*one_over_8*skew(3,isk)*nz
634 fx = zero
635 fy = zero
636 ENDIF
637 ENDIF
638 ENDIF
639C
640 IF(tagn1 == 0) THEN
641 iadn = iadc(iad+4*(n-1))
642 fsky(1,iadn) = fx
643 fsky(2,iadn) = fy
644 fsky(3,iadn) = fz
645 IF(ianim >0 .AND.impl_s==0) THEN
646 fext(1,n1) = fext(1,n1) + fx
647 fext(2,n1) = fext(2,n1) + fy
648 fext(3,n1) = fext(3,n1) + fz
649 ENDIF
650 ENDIF
651C
652 IF(tagn2 == 0) THEN
653 iadn = iadc(iad+4*(n-1)+1)
654 fsky(1,iadn) = fx
655 fsky(2,iadn) = fy
656 fsky(3,iadn) = fz
657 IF(ianim >0 .AND.impl_s==0) THEN
658 fext(1,n2) = fext(1,n2) + fx
659 fext(2,n2) = fext(2,n2) + fy
660 fext(3,n2) = fext(3,n2) + fz
661 ENDIF
662 ENDIF
663C
664 IF(tagn3 == 0) THEN
665 iadn = iadc(iad+4*(n-1)+2)
666 fsky(1,iadn) = fx
667 fsky(2,iadn) = fy
668 fsky(3,iadn) = fz
669 IF(ianim >0 .AND.impl_s==0) THEN
670 fext(1,n3) = fext(1,n3) + fx
671 fext(2,n3) = fext(2,n3) + fy
672 fext(3,n3) = fext(3,n3) + fz
673 ENDIF
674 ENDIF
675C
676 IF(tagn4 == 0) THEN
677 iadn = iadc(iad+4*(n-1)+3)
678 fsky(1,iadn) = fx
679 fsky(2,iadn) = fy
680 fsky(3,iadn) = fz
681 IF(ianim >0 .AND.impl_s==0) THEN
682 fext(1,n4) = fext(1,n4) + fx
683 fext(2,n4) = fext(2,n4) + fy
684 fext(3,n4) = fext(3,n4) + fz
685 ENDIF
686 ENDIF
687C
688 IF(th_surf%LOADP_FLAG >0 ) THEN
689 area = half*sqrt(nx*nx+ny*ny+nz*nz)
690 DO ns=th_surf%LOADP_KSEGS(nsegpl) +1,th_surf%LOADP_KSEGS(nsegpl+1)
691 ksurf = th_surf%LOADP_SEGS(ns)
692 th_surf%channels(4,ksurf)= th_surf%channels(4,ksurf) + area*aa
693 th_surf%channels(5,ksurf)= th_surf%channels(5,ksurf) + area
694 ENDDO
695 ENDIF
696 wfextt=wfextt+dt1*(fx*(v(1,n1)+v(1,n2)+v(1,n3)+v(1,n4))
697 + +fy*(v(2,n1)+v(2,n2)+v(2,n3)+v(2,n4))
698 + +fz*(v(3,n1)+v(3,n2)+v(3,n3)+v(3,n4)))
699C
700 ELSE
701
702 ! true triangles.
703 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))
704 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))
705 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))
706
707 IF(inorm == 1) THEN
708 fx = aa*nx*one_over_6
709 fy = aa*ny*one_over_6
710 fz = aa*nz*one_over_6
711 ELSEIF(inorm==2) THEN
712 norm = sqrt(nx*nx+ny*ny+nz*nz)
713 IF(isk == 0) THEN
714 IF(idir == 1 ) THEN
715 fx = aa*one_over_6*norm
716 fy = zero
717 fz = zero
718 ELSEIF(idir==2) THEN
719 fy = aa*one_over_6*norm
720 fx = zero
721 fz = zero
722 ELSEIF(idir==3) THEN
723 fz = aa*one_over_6*norm
724 fx = zero
725 fy = zero
726 ENDIF
727 ELSE
728 IF(idir == 1 ) THEN
729 fx = aa*one_over_6*norm*skew(1,isk)
730 fy = zero
731 fz = zero
732 ELSEIF(idir==2) THEN
733 fy = aa*one_over_6*norm*skew(2,isk)
734 fx = zero
735 fz = zero
736 ELSEIF(idir==3) THEN
737 fz = aa*one_over_6*norm*skew(3,isk)
738 fx = zero
739 fy = zero
740 ENDIF
741 ENDIF
742 ELSEIF(inorm==3) THEN
743
744 IF(isk == 0) THEN
745 IF(idir == 1 ) THEN
746 IF(nx /= zero) fx = aa*one_over_6*nx
747 fy = zero
748 fz = zero
749 ELSEIF(idir==2) THEN
750 IF(ny /= zero) fy = aa*one_over_6*ny
751 fx = zero
752 fz = zero
753 ELSEIF(idir==3) THEN
754 IF(nz /= zero) fz = aa*one_over_6*nz
755 fx = zero
756 fy = zero
757 ENDIF
758 ELSE
759 IF(idir == 1 ) THEN
760 fx = aa*one_over_6*skew(1,isk)*nx
761 fy = zero
762 fz = zero
763 ELSEIF(idir==2) THEN
764 fy = aa*one_over_6*skew(2,isk)*ny
765 fx = zero
766 fz = zero
767 ELSEIF(idir==3) THEN
768 fz = aa*one_over_6*skew(3,isk)*nz
769 fx = zero
770 fy = zero
771 ENDIF
772 ENDIF
773 ENDIF
774C
775 IF(tagn1 == 0) THEN
776 iadn = iadc(iad+4*(n-1))
777 fsky(1,iadn) = fx
778 fsky(2,iadn) = fy
779 fsky(3,iadn) = fz
780 IF(ianim >0 .AND.impl_s==0) THEN
781 fext(1,n1) = fext(1,n1) + fx
782 fext(2,n1) = fext(2,n1) + fy
783 fext(3,n1) = fext(3,n1) + fz
784 ENDIF
785 ENDIF
786C
787 IF(tagn2 == 0) THEN
788 iadn = iadc(iad+4*(n-1)+1)
789 fsky(1,iadn) = fx
790 fsky(2,iadn) = fy
791 fsky(3,iadn) = fz
792 IF(ianim >0 .AND.impl_s==0) THEN
793 fext(1,n2) = fext(1,n2) + fx
794 fext(2,n2) = fext(2,n2) + fy
795 fext(3,n2) = fext(3,n2) + fz
796 ENDIF
797 ENDIF
798C
799 IF(tagn3 == 0) THEN
800 iadn = iadc(iad+4*(n-1)+2)
801 fsky(1,iadn) = fx
802 fsky(2,iadn) = fy
803 fsky(3,iadn) = fz
804 IF(ianim >0 .AND.impl_s==0) THEN
805 fext(1,n3) = fext(1,n3) + fx
806 fext(2,n3) = fext(2,n3) + fy
807 fext(3,n3) = fext(3,n3) + fz
808 ENDIF
809 ENDIF
810
811C
812 IF(th_surf%LOADP_FLAG >0 ) THEN
813 area = half*sqrt(nx*nx+ny*ny+nz*nz)
814 DO ns=th_surf%LOADP_KSEGS(nsegpl) +1,th_surf%LOADP_KSEGS(nsegpl+1)
815 ksurf = th_surf%LOADP_SEGS(ns)
816 th_surf%channels(4,ksurf)= th_surf%channels(4,ksurf) + area*aa
817 th_surf%channels(5,ksurf)= th_surf%channels(5,ksurf) + area
818 ENDDO
819 ENDIF
820
821 wfextt=wfextt+dt1*(fx*(v(1,n1)+v(1,n2)+v(1,n3))
822 + + fy*(v(2,n1)+v(2,n2)+v(2,n3))
823 + + fz*(v(3,n1)+v(3,n2)+v(3,n3)))
824
825 ENDIF
826C
827 ENDIF
828C
829 ENDDO
830C
831 ENDDO
832C
833!$omp atomic
834 wfext = wfext + wfextt
835C
836 ENDIF
837 RETURN
838 END
my_real
#define my_real
Definition cppsort.cpp:32
load_pressure
subroutine load_pressure(iloadp, loadp, lloadp, npc, tf, a, v, x, skew, sensor_tab, iadc, fsky, fext, tagncont, nsensor, loadp_hyd_inter, h3d_data, python, npresload, loadp_tagdel, th_surf, pblast, wfext)
Definition load_pressure.F:46
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
h3d_mod::h3d_database
Definition h3d_mod.F:435
th_surf_mod::th_surf_
Definition th_surf_mod.F:75