OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
volpresp.F File Reference
#include "implicit_f.inc"
#include "param_c.inc"
#include "com04_c.inc"
#include "com08_c.inc"
#include "scr14_c.inc"
#include "scr16_c.inc"
#include "parit_c.inc"

Go to the source code of this file.

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ volprep()

subroutine volprep ( integer, dimension(*) ivolu,
rvolu,
integer njet,
integer, dimension(nibjet,*) ibagjet,
rbagjet,
integer, intent(in) nsensor,
type (sensor_str_), dimension(nsensor), intent(in) sensor_tab,
x,
v,
a,
normal,
integer, dimension(*) npc,
tf,
integer nn,
integer, dimension(nn,4) surf_nodes,
integer, dimension(4,*) iadmv,
fsky,
fskyv,
fext,
type(h3d_database) h3d_data,
integer, dimension(nn) surf_eltyp,
integer, dimension(nn) surf_elem,
integer, intent(in) node_number,
integer, intent(in) total_contribution_number,
integer, dimension(nn,4), intent(in) contribution_index,
integer, dimension(node_number+1), intent(in) contribution_number,
integer, dimension(node_number), intent(in) node_id,
intent(inout) contribution,
double precision, intent(inout) wfext,
type(python_), intent(inout) python )
Parameters
[in]nsensornumber of sensor
[in]node_numbernumber of node of the airbag
[in]total_contribution_numbertotal number of contribution
[in]contribution_indexindex to the array CONTRIBUTION
[in]contribution_numbernumber of contribution per node
[in]node_idnode id

Definition at line 33 of file volpresp.F.

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
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