OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
volpresp.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!|| volprep ../engine/source/airbag/volpresp.F
25!||--- called by ------------------------------------------------------
26!|| monvol0 ../engine/source/airbag/monvol0.F
27!||--- uses -----------------------------------------------------
28!|| finter_mixed_mod ../engine/source/tools/finter_mixed.F90
29!|| h3d_mod ../engine/share/modules/h3d_mod.F
30!|| python_funct_mod ../common_source/modules/python_mod.F90
31!|| sensor_mod ../common_source/modules/sensor_mod.F90
32!||====================================================================
33 SUBROUTINE volprep(IVOLU ,RVOLU ,NJET ,IBAGJET ,RBAGJET ,
34 2 NSENSOR,SENSOR_TAB,X ,V ,A ,
35 3 NORMAL ,NPC ,TF ,NN ,SURF_NODES,
36 4 IADMV ,FSKY ,FSKYV ,FEXT ,H3D_DATA ,
37 5 SURF_ELTYP,SURF_ELEM,NODE_NUMBER,TOTAL_CONTRIBUTION_NUMBER,
38 6 CONTRIBUTION_INDEX,CONTRIBUTION_NUMBER,
39 7 NODE_ID,CONTRIBUTION,WFEXT,PYTHON)
40C-----------------------------------------------
41C M o d u l e s
42C-----------------------------------------------
43 USE python_funct_mod
44 USE finter_mixed_mod
45 USE h3d_mod
46 USE sensor_mod
47C-----------------------------------------------
48C AIRBAGS INPUT FORMAT 4.4
49C-----------------------------------------------
50C I m p l i c i t T y p e s
51C-----------------------------------------------
52#include "implicit_f.inc"
53C-----------------------------------------------
54C o m m o n B l o c k s
55C-----------------------------------------------
56#include "param_c.inc"
57#include "com04_c.inc"
58#include "com08_c.inc"
59#include "scr14_c.inc"
60#include "scr16_c.inc"
61#include "parit_c.inc"
62C-----------------------------------------------
63C D u m m y A r g u m e n t s
64C-----------------------------------------------
65 INTEGER ,INTENT(IN) :: NSENSOR !< number of sensor
66 INTEGER ,INTENT(IN) :: NODE_NUMBER !< number of node of the airbag
67 INTEGER ,INTENT(IN) :: TOTAL_CONTRIBUTION_NUMBER !< total number of contribution
68 INTEGER IVOLU(*),NJET,IBAGJET(NIBJET,*),
69 . NPC(*),NN,IADMV(4,*),SURF_NODES(NN,4),SURF_ELTYP(NN),
70 . SURF_ELEM(NN)
71 INTEGER, DIMENSION(NN,4), INTENT(in) :: CONTRIBUTION_INDEX !< index to the array CONTRIBUTION
72 INTEGER, DIMENSION(NODE_NUMBER+1), INTENT(in) :: CONTRIBUTION_NUMBER !< number of contribution per node
73 INTEGER, DIMENSION(NODE_NUMBER), INTENT(in) :: NODE_ID !< node id
74 my_real, DIMENSION(TOTAL_CONTRIBUTION_NUMBER,3), INTENT(inout) :: contribution !< contribution array, the force are saved here
75 my_real rvolu(*),rbagjet(nrbjet,*),x(3,*), v(3,*), a(3,*), normal(3,nn),tf(*),fsky(8,lsky),fskyv(lsky,8),fext(3,*)
76 TYPE(h3d_database) :: H3D_DATA
77 TYPE (SENSOR_STR_) ,DIMENSION(NSENSOR) ,INTENT(IN) :: SENSOR_TAB
78 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
79 TYPE(python_), intent(inout) :: PYTHON
80C-----------------------------------------------
81C L o c a l V a r i a b l e s
82C-----------------------------------------------
83 INTEGER I,J,II,NJ1,NJ2,NJ3,IPT,IPA,IPZ,N1,N2,N3,N4,IJ,ISENS,K,ITYP
84 INTEGER :: ADDRESS,NOD,LOCAL_CONTRIBUTION_NUMBER
85 my_real :: pres,fx,fy,fz,p,pext,dp,q,pj,pj0,xx,yy,zz,
86 . theta,nx1,nx2,ny1,ny2,nz1,nz2,an1,an,
87 . xm,ym,zm,nx3,ny3,nz3,
88 . ps1,ps2,ps3,an2,an3,det,mu,facr,
89 . api,dist,tstart,ts,tmpo,fpt,fpa,fpz,scalt,scala,scald
90 my_real, DIMENSION(3) :: tmp
91C-----------------------------------------------
92 api=hundred80/pi
93C-----------------------------------------------
94!$OMP PARALLEL PRIVATE(N1,N2,N3,N4,P,II,FX,FY,FZ,ADDRESS,NOD,LOCAL_CONTRIBUTION_NUMBER,TMP)
95 pres =rvolu(12)
96 q =rvolu(23)
97 pext =rvolu(3)
98 dp =q+pres-pext
99 scalt =rvolu(26)
100 scala =rvolu(29)
101 scald =rvolu(30)
102!$OMP DO SCHEDULE(guided)
103 DO i=1,total_contribution_number
104 contribution(i,1) = zero
105 contribution(i,2) = zero
106 contribution(i,3) = zero
107 ENDDO
108!$OMP END DO
109
110!$OMP DO SCHEDULE(guided)
111 DO i=1,nn
112C
113 n1 = surf_nodes(i,1)
114 n2 = surf_nodes(i,2)
115 n3 = surf_nodes(i,3)
116 IF(surf_eltyp(i) /= 7) n4 = surf_nodes(i,4)
117C
118 IF(surf_eltyp(i) == 3)THEN
119 p=dp*fourth
120 ii=surf_elem(i)
121C
122 fx=p*normal(1,i)
123 fy=p*normal(2,i)
124 fz=p*normal(3,i)
125C
126 fsky(1,iadmv(1,i))=fx
127 fsky(2,iadmv(1,i))=fy
128 fsky(3,iadmv(1,i))=fz
129C
130 fsky(1,iadmv(2,i))=fx
131 fsky(2,iadmv(2,i))=fy
132 fsky(3,iadmv(2,i))=fz
133C
134 fsky(1,iadmv(3,i))=fx
135 fsky(2,iadmv(3,i))=fy
136 fsky(3,iadmv(3,i))=fz
137C
138 fsky(1,iadmv(4,i))=fx
139 fsky(2,iadmv(4,i))=fy
140 fsky(3,iadmv(4,i))=fz
141 IF(anim_v(5)+outp_v(5)+h3d_data%N_VECT_FINT+
142 . anim_v(6)+outp_v(6)+h3d_data%N_VECT_FEXT >0) THEN
143 address = contribution_index(i,1)
144 contribution(address,1) = fx
145 contribution(address,2) = fy
146 contribution(address,3) = fz
147 address = contribution_index(i,2)
148 contribution(address,1) = fx
149 contribution(address,2) = fy
150 contribution(address,3) = fz
151 address = contribution_index(i,3)
152 contribution(address,1) = fx
153 contribution(address,2) = fy
154 contribution(address,3) = fz
155 address = contribution_index(i,4)
156 contribution(address,1) = fx
157 contribution(address,2) = fy
158 contribution(address,3) = fz
159 ENDIF
160 ELSEIF(surf_eltyp(i)==7)THEN
161 p=dp*third
162 ii=surf_elem(i) + numelc
163C
164 fx=p*normal(1,i)
165 fy=p*normal(2,i)
166 fz=p*normal(3,i)
167C
168 fsky(1,iadmv(1,i))=fx
169 fsky(2,iadmv(1,i))=fy
170 fsky(3,iadmv(1,i))=fz
171C
172 fsky(1,iadmv(2,i))=fx
173 fsky(2,iadmv(2,i))=fy
174 fsky(3,iadmv(2,i))=fz
175C
176 fsky(1,iadmv(3,i))=fx
177 fsky(2,iadmv(3,i))=fy
178 fsky(3,iadmv(3,i))=fz
179 IF(anim_v(5)+outp_v(5)+h3d_data%N_VECT_FINT+
180 . anim_v(6)+outp_v(6)+h3d_data%N_VECT_FEXT >0) THEN
181 address = contribution_index(i,1)
182 contribution(address,1) = fx
183 contribution(address,2) = fy
184 contribution(address,3) = fz
185 address = contribution_index(i,2)
186 contribution(address,1) = fx
187 contribution(address,2) = fy
188 contribution(address,3) = fz
189 address = contribution_index(i,3)
190 contribution(address,1) = fx
191 contribution(address,2) = fy
192 contribution(address,3) = fz
193 ENDIF
194 ELSE
195 p=dp*fourth
196 ii=i+numelc+numeltg
197C
198 fx=p*normal(1,i)
199 fy=p*normal(2,i)
200 fz=p*normal(3,i)
201C
202 fsky(1,iadmv(1,i))=fx
203 fsky(2,iadmv(1,i))=fy
204 fsky(3,iadmv(1,i))=fz
205C
206 fsky(1,iadmv(2,i))=fx
207 fsky(2,iadmv(2,i))=fy
208 fsky(3,iadmv(2,i))=fz
209C
210 fsky(1,iadmv(3,i))=fx
211 fsky(2,iadmv(3,i))=fy
212 fsky(3,iadmv(3,i))=fz
213C
214 fsky(1,iadmv(4,i))=fx
215 fsky(2,iadmv(4,i))=fy
216 fsky(3,iadmv(4,i))=fz
217 IF(anim_v(5)+outp_v(5)+h3d_data%N_VECT_FINT+
218 . anim_v(6)+outp_v(6)+h3d_data%N_VECT_FEXT >0) THEN
219 address = contribution_index(i,1)
220 contribution(address,1) = fx
221 contribution(address,2) = fy
222 contribution(address,3) = fz
223 address = contribution_index(i,2)
224 contribution(address,1) = fx
225 contribution(address,2) = fy
226 contribution(address,3) = fz
227 address = contribution_index(i,3)
228 contribution(address,1) = fx
229 contribution(address,2) = fy
230 contribution(address,3) = fz
231 address = contribution_index(i,4)
232 contribution(address,1) = fx
233 contribution(address,2) = fy
234 contribution(address,3) = fz
235 ENDIF
236 ENDIF
237 ENDDO
238
239!$OMP END DO
240
241 IF(anim_v(5)+outp_v(5)+h3d_data%N_VECT_FINT+
242 . anim_v(6)+outp_v(6)+h3d_data%N_VECT_FEXT >0) THEN
243!$OMP DO SCHEDULE(guided)
244 DO i=1,node_number
245 nod = node_id(i)
246 address = contribution_number(i)
247 local_contribution_number = contribution_number(i+1) - contribution_number(i)
248 tmp(1) = zero
249 tmp(2) = zero
250 tmp(3) = zero
251 DO j=1,local_contribution_number
252 tmp(1) = tmp(1) + contribution(address+j,1)
253 tmp(2) = tmp(2) + contribution(address+j,2)
254 tmp(3) = tmp(3) + contribution(address+j,3)
255 ENDDO
256 fext(1,nod) = fext(1,nod)+ tmp(1)
257 fext(2,nod) = fext(2,nod)+ tmp(2)
258 fext(3,nod) = fext(3,nod)+ tmp(3)
259 ENDDO
260!$OMP END DO
261 ENDIF
262
263!$OMP END PARALLEL
264C
265 ityp=ivolu(2)
266 IF(ityp/=4.AND.ityp/=5.AND.ityp/=7.AND.ityp/=9)THEN
267 RETURN
268 ENDIF
269C-------------------------------------------
270C INJECTEURS
271C-------------------------------------------
272 DO ij =1,njet
273 nj1 = ibagjet( 5,ij)
274 nj2 = ibagjet( 6,ij)
275 nj3 = ibagjet( 7,ij)
276 ipt = ibagjet( 8,ij)
277 ipa = ibagjet( 9,ij)
278 ipz = ibagjet(10,ij)
279 fpt = rbagjet(12,ij)
280 fpa = rbagjet(13,ij)
281 fpz = rbagjet(14,ij)
282
283 isens=ibagjet(4,ij)
284 IF(isens==0)THEN
285 tstart=zero
286 ELSE
287 tstart=sensor_tab(isens)%TSTART
288 ENDIF
289
290 IF(nj1/=0.AND.tt>=tstart)THEN
291 ts=tt-tstart
292 pj0=fpt*finter_mixed(python,nfunct,ipt,ts*scalt,npc,tf)
293!$OMP DO SCHEDULE(guided)
294 DO i=1,total_contribution_number
295 contribution(i,1) = zero
296 contribution(i,2) = zero
297 contribution(i,3) = zero
298 ENDDO
299!$OMP END DO
300
301!$OMP PARALLEL PRIVATE(N1,N2,N3,N4,PJ,II,XX,YY,ZZ,XM,YM,ZM,NX1,NY1,NZ1,NX2,NY2,NZ2,NX3,NY3,NZ3,
302!$OMP+ PS1,PS2,PS3,AN2,AN3,DET,MU,FACR,AN1,AN,TMPO,THETA,DIST,FX,FY,FZ,ADDRESS,
303!$OMP+ NOD,LOCAL_CONTRIBUTION_NUMBER,TMP)
304
305!$OMP DO SCHEDULE(guided) REDUCTION(+:WFEXT)
306 DO i=1,nn
307
308 IF(surf_eltyp(i)==3)THEN
309 pj=fourth*pj0
310
311 n1 = surf_nodes(i,1)
312 n2 = surf_nodes(i,2)
313 n3 = surf_nodes(i,3)
314 n4 = surf_nodes(i,4)
315
316 ii=surf_elem(i)
317
318 xx = fourth*(x(1,n1)+x(1,n2)+x(1,n3)+x(1,n4))
319 yy = fourth*(x(2,n1)+x(2,n2)+x(2,n3)+x(2,n4))
320 zz = fourth*(x(3,n1)+x(3,n2)+x(3,n3)+x(3,n4))
321
322 xm=half*(x(1,nj1)+x(1,nj3))
323 ym=half*(x(2,nj1)+x(2,nj3))
324 zm=half*(x(3,nj1)+x(3,nj3))
325
326 nx1 = xx-xm
327 ny1 = yy-ym
328 nz1 = zz-zm
329
330 ! decomposition de (M,P) sur (M,N2) et (M,N3)
331 nx2 = x(1,nj2)-xm
332 ny2 = x(2,nj2)-ym
333 nz2 = x(3,nj2)-zm
334
335 nx3 = x(1,nj3)-xm
336 ny3 = x(2,nj3)-ym
337 nz3 = x(3,nj3)-zm
338
339 ps1 = nx1*nx2+ny1*ny2+nz1*nz2
340 ps2 = nx2*nx3+ny2*ny3+nz2*nz3
341 ps3 = nx1*nx3+ny1*ny3+nz1*nz3
342 an2 = nx2*nx2+ny2*ny2+nz2*nz2
343 an3 = nx3*nx3+ny3*ny3+nz3*nz3
344 det = ps2*ps2-an2*an3
345
346 mu = (ps2*ps1-ps3*an2)/sign(max(em30,abs(det)),det)
347
348 facr =min(one,max(-one,mu))
349 nx1=nx1-facr*nx3
350 ny1=ny1-facr*ny3
351 nz1=nz1-facr*nz3
352
353 an1 = (nx1**2+ny1**2+nz1**2)
354 an = max(em30,sqrt((normal(1,i)**2+normal(2,i)**2+normal(3,i)**2)*an1))
355 pj=pj*max(zero,(nx1*normal(1,i)+ny1*normal(2,i)+nz1*normal(3,i))/an)
356 an = max(em30,sqrt((nx2**2+ny2**2+nz2**2)*an1))
357 tmpo = (nx1*nx2+ny1*ny2+nz1*nz2)/an
358 tmpo = sign(min(one,abs(tmpo)),tmpo)
359 theta=api*acos(tmpo)
360 pj=pj*fpa*finter_mixed(python,nfunct,ipa,theta*scala,npc,tf)
361 dist = sqrt(an1)
362 IF(ipz/=0)pj=pj*fpz*finter_mixed(python,nfunct,ipz,dist*scald,npc,tf)
363
364 fx=pj*normal(1,i)
365 fy=pj*normal(2,i)
366 fz=pj*normal(3,i)
367
368 fsky(1,iadmv(1,i))=fsky(1,iadmv(1,i))+fx
369 fsky(2,iadmv(1,i))=fsky(2,iadmv(1,i))+fy
370 fsky(3,iadmv(1,i))=fsky(3,iadmv(1,i))+fz
371
372 fsky(1,iadmv(2,i))=fsky(1,iadmv(2,i))+fx
373 fsky(2,iadmv(2,i))=fsky(2,iadmv(2,i))+fy
374 fsky(3,iadmv(2,i))=fsky(3,iadmv(2,i))+fz
375
376 fsky(1,iadmv(3,i))=fsky(1,iadmv(3,i))+fx
377 fsky(2,iadmv(3,i))=fsky(2,iadmv(3,i))+fy
378 fsky(3,iadmv(3,i))=fsky(3,iadmv(3,i))+fz
379
380 fsky(1,iadmv(4,i))=fsky(1,iadmv(4,i))+fx
381 fsky(2,iadmv(4,i))=fsky(2,iadmv(4,i))+fy
382 fsky(3,iadmv(4,i))=fsky(3,iadmv(4,i))+fz
383
384 IF(anim_v(5)+outp_v(5)+h3d_data%N_VECT_FINT+
385 . anim_v(6)+outp_v(6)+h3d_data%N_VECT_FEXT >0) THEN
386 address = contribution_index(i,1)
387 contribution(address,1) = fx
388 contribution(address,2) = fy
389 contribution(address,3) = fz
390 address = contribution_index(i,2)
391 contribution(address,1) = fx
392 contribution(address,2) = fy
393 contribution(address,3) = fz
394 address = contribution_index(i,3)
395 contribution(address,1) = fx
396 contribution(address,2) = fy
397 contribution(address,3) = fz
398 address = contribution_index(i,4)
399 contribution(address,1) = fx
400 contribution(address,2) = fy
401 contribution(address,3) = fz
402 ENDIF
403
404 wfext=wfext+dt1*(fx*(v(1,n1)+v(1,n2)+v(1,n3)+v(1,n4))
405 + +fy*(v(2,n1)+v(2,n2)+v(2,n3)+v(2,n4))
406 + +fz*(v(3,n1)+v(3,n2)+v(3,n3)+v(3,n4)))
407 ELSEIF(surf_eltyp(i)==7)THEN
408 pj=pj0*third
409
410 n1 = surf_nodes(i,1)
411 n2 = surf_nodes(i,2)
412 n3 = surf_nodes(i,3)
413
414 ii=surf_elem(i) + numelc
415
416 xx = (x(1,n1)+x(1,n2)+x(1,n3))*third
417 yy = (x(2,n1)+x(2,n2)+x(2,n3))*third
418 zz = (x(3,n1)+x(3,n2)+x(3,n3))*third
419
420 xm=half*(x(1,nj1)+x(1,nj3))
421 ym=half*(x(2,nj1)+x(2,nj3))
422 zm=half*(x(3,nj1)+x(3,nj3))
423
424 nx1 = xx-xm
425 ny1 = yy-ym
426 nz1 = zz-zm
427
428 ! decomposition de (M,P) sur (M,N2) et (M,N3)
429 nx2 = x(1,nj2)-xm
430 ny2 = x(2,nj2)-ym
431 nz2 = x(3,nj2)-zm
432
433 nx3 = x(1,nj3)-xm
434 ny3 = x(2,nj3)-ym
435 nz3 = x(3,nj3)-zm
436
437 ps1 = nx1*nx2+ny1*ny2+nz1*nz2
438 ps2 = nx2*nx3+ny2*ny3+nz2*nz3
439 ps3 = nx1*nx3+ny1*ny3+nz1*nz3
440 an2 = nx2*nx2+ny2*ny2+nz2*nz2
441 an3 = nx3*nx3+ny3*ny3+nz3*nz3
442 det = ps2*ps2-an2*an3
443
444 mu = (ps2*ps1-ps3*an2)/sign(max(em30,abs(det)),det)
445
446 facr =min(one,max(-one,mu))
447 nx1=nx1-facr*nx3
448 ny1=ny1-facr*ny3
449 nz1=nz1-facr*nz3
450
451 an1 = (nx1**2+ny1**2+nz1**2)
452 an = max(em30,sqrt((normal(1,i)**2+normal(2,i)**2+normal(3,i)**2)*an1))
453 pj=pj*max(zero,(nx1*normal(1,i)+ny1*normal(2,i)+nz1*normal(3,i))/an)
454 an = max(em30,sqrt((nx2**2+ny2**2+nz2**2)*an1))
455 tmpo = (nx1*nx2+ny1*ny2+nz1*nz2)/an
456 tmpo = sign(min(one,abs(tmpo)),tmpo)
457 theta=api*acos(tmpo)
458 pj=pj*fpa*finter_mixed(python,nfunct,ipa,theta*scala,npc,tf)
459 dist = sqrt(an1)
460 IF(ipz/=0)pj=pj*fpz*finter_mixed(python,nfunct,ipz,dist*scald,npc,tf)
461
462 fx=pj*normal(1,i)
463 fy=pj*normal(2,i)
464 fz=pj*normal(3,i)
465
466 fsky(1,iadmv(1,i))=fsky(1,iadmv(1,i))+fx
467 fsky(2,iadmv(1,i))=fsky(2,iadmv(1,i))+fy
468 fsky(3,iadmv(1,i))=fsky(3,iadmv(1,i))+fz
469
470 fsky(1,iadmv(2,i))=fsky(1,iadmv(2,i))+fx
471 fsky(2,iadmv(2,i))=fsky(2,iadmv(2,i))+fy
472 fsky(3,iadmv(2,i))=fsky(3,iadmv(2,i))+fz
473
474 fsky(1,iadmv(3,i))=fsky(1,iadmv(3,i))+fx
475 fsky(2,iadmv(3,i))=fsky(2,iadmv(3,i))+fy
476 fsky(3,iadmv(3,i))=fsky(3,iadmv(3,i))+fz
477
478 IF(anim_v(5)+outp_v(5)+h3d_data%N_VECT_FINT+
479 . anim_v(6)+outp_v(6)+h3d_data%N_VECT_FEXT >0) THEN
480 address = contribution_index(i,1)
481 contribution(address,1) = fx
482 contribution(address,2) = fy
483 contribution(address,3) = fz
484 address = contribution_index(i,2)
485 contribution(address,1) = fx
486 contribution(address,2) = fy
487 contribution(address,3) = fz
488 address = contribution_index(i,3)
489 contribution(address,1) = fx
490 contribution(address,2) = fy
491 contribution(address,3) = fz
492 ENDIF
493
494 wfext=wfext+dt1*(fx*(v(1,n1)+v(1,n2)+v(1,n3))
495 + +fy*(v(2,n1)+v(2,n2)+v(2,n3))
496 + +fz*(v(3,n1)+v(3,n2)+v(3,n3)))
497 ELSE
498 pj=fourth*pj0
499
500 n1 = surf_nodes(i,1)
501 n2 = surf_nodes(i,2)
502 n3 = surf_nodes(i,3)
503 n4 = surf_nodes(i,4)
504
505 ii=i+numelc+numeltg
506
507 xx = fourth*(x(1,n1)+x(1,n2)+x(1,n3)+x(1,n4))
508 yy = fourth*(x(2,n1)+x(2,n2)+x(2,n3)+x(2,n4))
509 zz = fourth*(x(3,n1)+x(3,n2)+x(3,n3)+x(3,n4))
510
511 xm=half*(x(1,nj1)+x(1,nj3))
512 ym=half*(x(2,nj1)+x(2,nj3))
513 zm=half*(x(3,nj1)+x(3,nj3))
514
515 nx1 = xx-xm
516 ny1 = yy-ym
517 nz1 = zz-zm
518
519 ! decomposition de (M,P) sur (M,N2) et (M,N3)
520 nx2 = x(1,nj2)-xm
521 ny2 = x(2,nj2)-ym
522 nz2 = x(3,nj2)-zm
523
524 nx3 = x(1,nj3)-xm
525 ny3 = x(2,nj3)-ym
526 nz3 = x(3,nj3)-zm
527
528 ps1 = nx1*nx2+ny1*ny2+nz1*nz2
529 ps2 = nx2*nx3+ny2*ny3+nz2*nz3
530 ps3 = nx1*nx3+ny1*ny3+nz1*nz3
531 an2 = nx2*nx2+ny2*ny2+nz2*nz2
532 an3 = nx3*nx3+ny3*ny3+nz3*nz3
533 det = ps2*ps2-an2*an3
534
535 mu = (ps2*ps1-ps3*an2)/sign(max(em30,abs(det)),det)
536
537 facr =min(one,max(-one,mu))
538 nx1=nx1-facr*nx3
539 ny1=ny1-facr*ny3
540 nz1=nz1-facr*nz3
541
542 an1 = (nx1**2+ny1**2+nz1**2)
543 an = max(em30,sqrt((normal(1,i)**2+normal(2,i)**2+normal(3,i)**2)*an1))
544 pj=pj*max(zero,(nx1*normal(1,i)+ny1*normal(2,i)+nz1*normal(3,i))/an)
545 an = max(em30,sqrt((nx2**2+ny2**2+nz2**2)*an1))
546 tmpo = (nx1*nx2+ny1*ny2+nz1*nz2)/an
547 tmpo = sign(min(one,abs(tmpo)),tmpo)
548 theta=api*acos(tmpo)
549 pj=pj*fpa*finter_mixed(python,nfunct,ipa,theta*scala,npc,tf)
550 dist = sqrt(an1)
551 IF(ipz/=0)pj=pj*fpz*finter_mixed(python,nfunct,ipz,dist*scald,npc,tf)
552
553 fx=pj*normal(1,i)
554 fy=pj*normal(2,i)
555 fz=pj*normal(3,i)
556
557 fsky(1,iadmv(1,i))=fsky(1,iadmv(1,i))+fx
558 fsky(2,iadmv(1,i))=fsky(2,iadmv(1,i))+fy
559 fsky(3,iadmv(1,i))=fsky(3,iadmv(1,i))+fz
560
561 fsky(1,iadmv(2,i))=fsky(1,iadmv(2,i))+fx
562 fsky(2,iadmv(2,i))=fsky(2,iadmv(2,i))+fy
563 fsky(3,iadmv(2,i))=fsky(3,iadmv(2,i))+fz
564
565 fsky(1,iadmv(3,i))=fsky(1,iadmv(3,i))+fx
566 fsky(2,iadmv(3,i))=fsky(2,iadmv(3,i))+fy
567 fsky(3,iadmv(3,i))=fsky(3,iadmv(3,i))+fz
568
569 fsky(1,iadmv(4,i))=fsky(1,iadmv(4,i))+fx
570 fsky(2,iadmv(4,i))=fsky(2,iadmv(4,i))+fy
571 fsky(3,iadmv(4,i))=fsky(3,iadmv(4,i))+fz
572
573 IF(anim_v(5)+outp_v(5)+h3d_data%N_VECT_FINT+
574 . anim_v(6)+outp_v(6)+h3d_data%N_VECT_FEXT >0) THEN
575 address = contribution_index(i,1)
576 contribution(address,1) = fx
577 contribution(address,2) = fy
578 contribution(address,3) = fz
579 address = contribution_index(i,2)
580 contribution(address,1) = fx
581 contribution(address,2) = fy
582 contribution(address,3) = fz
583 address = contribution_index(i,3)
584 contribution(address,1) = fx
585 contribution(address,2) = fy
586 contribution(address,3) = fz
587 address = contribution_index(i,4)
588 contribution(address,1) = fx
589 contribution(address,2) = fy
590 contribution(address,3) = fz
591 ENDIF
592
593 wfext=wfext+dt1*(fx*(v(1,n1)+v(1,n2)+v(1,n3)+v(1,n4))
594 + +fy*(v(2,n1)+v(2,n2)+v(2,n3)+v(2,n4))
595 + +fz*(v(3,n1)+v(3,n2)+v(3,n3)+v(3,n4)))
596 ENDIF
597 ENDDO
598!$OMP END DO
599
600 IF(anim_v(5)+outp_v(5)+h3d_data%N_VECT_FINT+
601 . anim_v(6)+outp_v(6)+h3d_data%N_VECT_FEXT >0) THEN
602!$OMP DO SCHEDULE(guided)
603 DO i=1,node_number
604 nod = node_id(i)
605 address = contribution_number(i)
606 local_contribution_number = contribution_number(i+1) - contribution_number(i)
607 tmp(1) = zero
608 tmp(2) = zero
609 tmp(3) = zero
610 DO j=1,local_contribution_number
611 tmp(1) = tmp(1) + contribution(address+j,1)
612 tmp(2) = tmp(2) + contribution(address+j,2)
613 tmp(3) = tmp(3) + contribution(address+j,3)
614 ENDDO
615 fext(1,nod) = fext(1,nod)+ tmp(1)
616 fext(2,nod) = fext(2,nod)+ tmp(2)
617 fext(3,nod) = fext(3,nod)+ tmp(3)
618 ENDDO
619!$omp END DO
620 ENDIF
621
622!$OMP END PARALLEL
623 ENDIF
624 ENDDO
625
626
627 RETURN
628 END
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
h3d_mod
Definition h3d_mod.F:308
h3d_mod::h3d_database
Definition h3d_mod.F:435
volprep
subroutine volprep(ivolu, rvolu, njet, ibagjet, rbagjet, nsensor, sensor_tab, x, v, a, normal, npc, tf, nn, surf_nodes, iadmv, fsky, fskyv, fext, h3d_data, surf_eltyp, surf_elem, node_number, total_contribution_number, contribution_index, contribution_number, node_id, contribution, wfext, python)
Definition volpresp.F:40