36
37
38
39
40
41
42
43
44
45
46
47
52 use element_mod , only : nixs
53
54
55
56#include "implicit_f.inc"
57
58
59
60#include "mvsiz_p.inc"
61
62
63
64#include "vect01_c.inc"
65#include "param_c.inc"
66#include "com01_c.inc"
67#include "com04_c.inc"
68#include "com08_c.inc"
69
70
71
72 INTEGER :: IXS(NIXS,NUMELS), NV46, IPM(NPROPMI,*),NEL
73 my_real :: pm(npropm,nummat), flux(6,*), flu1(*), veul(lveul,*),x(3,numnod)
74 TYPE(t_ale_connectivity)INTENT(IN) :: ALE_CONNECT
75
76
77
78 INTEGER MAT(MVSIZ), I, II,,IMAT,IALEFVM_FLG,IDV,ICF(4,6),IX1,IX2,IX3,IX4
80 . n2x(mvsiz), n3x(mvsiz), n4x(mvsiz), n5x(mvsiz), n6x(mvsiz), n1y(mvsiz), n2y(mvsiz), n3y(mvsiz),
81 . n4y(mvsiz), n5y(mvsiz), n6y(mvsiz), n1z(mvsiz),
82 . n6z(mvsiz), flux1(mvsiz), flux2
83 . flux6(mvsiz), vx1(mvsiz), vx2(mvsiz), vx3(mvsiz), vx4(mvsiz), vx5(mvsiz), vx6(mvsiz),
84 . vy1(mvsiz), vy2(mvsiz), vy3(mvsiz), vy4(mvsiz), vy5(mvsiz), vy6(mvsiz), vz1(mvsiz), vz2
85 . vz3(mvsiz), vz4(mvsiz), vz5(mvsiz), vz6(mvsiz),
86 . reduc,upwl(6,mvsiz),valvois(6,nv46,mvsiz), valel(6,mvsiz), vec(3,6),cf(3,mvsiz),
87 . z(3),zadj(3),zcf(3),zzadj(3),zzadj_,ps,lambda,sr1,sr2, srf, surf(6),term2
88 LOGICAL debug_outp
89 INTEGER :: idbf,idbl,NC(8,MVSIZ),IAD2,LGTH
90
91 DATA icf/1,4,3,2,3,4,8,7,5,6,7,8,1,2,6,5,2,3,7,6,1,5,8,4/
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
113 imat = ixs(1,1+nft)
114 ialefvm_flg = ipm(251,imat)
115 IF(ialefvm_flg <= 1)RETURN
116
117
118
119
120 DO i=1,nel
121 ii = i + nft
122 mat(i) = ixs(1,ii)
123 ENDDO
124
125
126
127
128
129 DO i=1,nel
130 ii = i + nft
131 iad2 = ale_connect%ee_connect%iad_connect(ii)
132 lgth = ale_connect%ee_connect%iad_connect(ii+1) - iad2
139 DO j=1,lgth
140 idv = ale_connect%ee_connect%connected(iad2 + j - 1)
141 IF(idv <= 0) THEN
142
149 ELSE
150
157 ENDIF
158
159 enddo
160
161
162 enddo
163
164
165
166
167
168
169
170 DO i=1,nel
171 ii=i+nft
172 n1x(i)=veul(14,ii)
173 n2x(i)=veul(15,ii)
174 n3x(i)=veul(16,ii)
175 n4x(i)=veul(17,ii)
176 n5x(i)=veul(18,ii)
177 n6x(i)=veul(19,ii)
178
179 n1y(i)=veul(20,ii)
180 n2y(i)=veul(21,ii)
181 n3y(i)=veul(22,ii)
182 n4y(i)=veul(23,ii)
183 n5y(i)=veul(24,ii)
184 n6y(i)=veul(25,ii)
185
186 n1z(i)=veul(26,ii)
187 n2z(i)=veul(27,ii)
188 n3z(i)=veul(28,ii)
189 n4z(i)=veul(29,ii)
190 n5z(i)=veul(30,ii)
191 n6z(i)=veul(31,ii)
192 ENDDO
193
194
195
196
197
198
200
201
203
204 CASE(1)
205
206 DO i=1,nel
207 ii = i + nft
208 DO j=1,nv46
209 vec(1,j) = half * (valel(1,i)/valel(4,i) + valvois(1,j,i)/valvois(4,j,i))
210 vec(2,j) = half * (valel(2,i)/valel(4,i) + valvois(2,j,i)/valvois(4,j,i))
211 vec(3,j) = half * (valel(3,i)/valel(4,i) + valvois(3,j,i)/valvois(4,j,i))
212 enddo
213 DO j=1,nv46
217 ENDDO
218
219 vx1(i) = vec(1,1)
220 vy1(i) = vec(2,1)
221 vz1(i) = vec(3,1)
222
223 vx2(i) = vec(1,2)
224 vy2(i) = vec(2,2)
225 vz2(i) = vec(3,2)
226
227 vx3(i) = vec(1,3)
228 vy3(i) = vec(2,3)
229 vz3(i) = vec(3,3)
230
231 vx4(i) = vec(1,4)
232 vy4(i) = vec(2,4)
233 vz4(i) = vec(3,4)
234
235 vx5(i) = vec(1,5)
236 vy5(i) = vec(2,5)
237 vz5(i) = vec(3,5)
238
239 vx6(i) = vec(1,6)
240 vy6(i) = vec(2,6)
241 vz6(i) = vec(3,6)
242 enddo
243
244 CASE(2)
245
246 DO i=1,nel
247 ii = i + nft
248 DO j=1,nv46
249 vec(1,j) = (valel(1,i) + valvois(1,j,i)) / (valel(4,i)+valvois(4,j,i))
250 vec(2,j) = (valel(2,i) + valvois(2,j,i)) / (valel(4,i)+valvois(4,j,i))
251 vec(3,j) = (valel(3,i) + valvois(3,j,i)) / (valel(4,i)+valvois(4,j,i))
252 enddo
253 DO j=1,nv46
257 ENDDO
258
259 vx1(i) = vec(1,1)
260 vy1(i) = vec(2,1)
261 vz1(i) = vec(3,1)
262
263 vx2(i) = vec(1,2)
264 vy2(i) = vec(2,2)
265 vz2(i) = vec(3,2)
266
267 vx3(i) = vec(1,3)
268 vy3(i) = vec(2,3)
269 vz3(i) = vec(3,3)
270
271 vx4(i) = vec(1,4)
272 vy4(i) = vec(2,4)
273 vz4(i) = vec(3,4)
274
275 vx5(i) = vec(1,5)
276 vy5(i) = vec(2,5)
277 vz5(i) = vec(3,5)
278
279 vx6(i) = vec(1,6)
280 vy6(i) = vec(2,6)
281 vz6(i) = vec(3,6)
282 enddo
283
284 CASE(3)
285
286 DO i=1,nel
287 ii = i + nft
288 sr1 = sqrt(valel(4,i))
289 DO j=1,nv46
290 sr2 = sqrt(valvois(4,j,i))
291 vec(1,j) = (valel(1,i)/sr1 + valvois(1,j,i)/sr2) / (sr1+sr2)
292 vec(2,j) = (valel(2,i)/sr1 + valvois(2,j,i)/sr2) / (sr1+sr2)
293 vec(3,j) = (valel(3,i)/sr1 + valvois(3,j,i)/sr2) / (sr1+sr2)
294 enddo
295 DO j=1,nv46
299 ENDDO
300
301 vx1(i) = vec(1,1)
302 vy1(i) = vec(2,1)
303 vz1(i) = vec(3,1)
304
305 vx2(i) = vec(1,2)
306 vy2(i) = vec(2,2)
307 vz2(i) = vec(3,2)
308
309 vx3(i) = vec(1,3)
310 vy3(i) = vec(2,3)
311 vz3(i) = vec(3,3)
312
313 vx4(i) = vec(1,4)
314 vy4(i) = vec(2,4)
315 vz4(i) = vec(3,4)
316
317 vx5(i) = vec(1,5)
318 vy5(i) = vec(2,5)
319 vz5(i) = vec(3,5)
320
321 vx6(i) = vec(1,6)
322 vy6(i) = vec(2,6)
323 vz6(i) = vec(3,6)
324 enddo
325
326 CASE(4)
327
328 DO i=1,nel
329 ii = i + nft
330 DO j=1,nv46
331 IF(dt1==zero)THEN
332 vec(1,j) = (valel(1,i) + valvois(1,j,i) )
333 . / (valel(4,i) + valvois(4,j,i) )
334 vec(2,j) = (valel(2,i) + valvois(2,j,i) )
335 . / (valel(4,i) + valvois(4,j,i) )
336 vec(3,j) = (valel(3,i) + valvois(3,j,i) )
337 . / (valel(4,i) + valvois(4,j,i) )
338 ELSE
339 vec(1,j) = (valel(1,i)*valel(5,i) + valvois(1,j,i)*valvois(5,j,i))
340 . / (valel(4,i)*valel(5,i) + valvois(4,j,i)*valvois(5,j,i))
341 vec(2,j) = (valel(2,i)*valel(5,i) + valvois(2,j,i)*valvois(5,j,i))
342 . / (valel(4,i)*valel(5,i) + valvois(4,j,i)*valvois
343 vec(3,j) = (valel(3,i)*valel(5,i) + valvois(3,j,i)*valvois(5,j,i))
344 . / (valel(4,i)*valel(5,i) + valvois(4,j,i)*valvois(5,j,i))
345 ENDIF
346 enddo
347 DO j=1,nv46
351 ENDDO
352
353 vx1(i) = vec(1,1)
354 vy1(i) = vec(2,1)
355 vz1(i) = vec(3,1)
356
357 vx2(i) = vec(1,2)
358 vy2(i) = vec(2,2)
359 vz2(i) = vec(3,2)
360
361 vx3(i) = vec(1,3)
362 vy3(i) = vec(2,3)
363 vz3(i) = vec(3,3)
364
365 vx4(i) = vec(1,4)
366 vy4(i) = vec(2,4)
367 vz4(i) = vec(3,4)
368
369 vx5(i) = vec(1,5)
370 vy5(i) = vec(2,5)
371 vz5(i) = vec(3,5)
372
373 vx6(i) = vec(1,6)
374 vy6(i) = vec(2,6)
375 vz6(i) = vec(3,6)
376 enddo
377
378 CASE(5)
379
380 DO i=1,nel
381 ii = i + nft
382 IF(dt1==zero)THEN
383 DO j=1,nv46
384 vec(1,j) = (valel(1,i) + valvois(1,j,i))
385 . / (valel(4,i) + valvois(4,j,i))
386 vec(2,j) = (valel(2,i) + valvois(2,j,i))
387 . / (valel(4,i) + valvois(4,j,i))
388 vec(3,j) = (valel(3,i) + valvois(3,j,i))
389 . / (valel(4,i) + valvois(4,j,i))
390 enddo
391
392 ELSE
393 surf(1) = one/sqrt(n1x(i)*n1x(i)+n1y(i)*n1y(i)+n1z(i)*n1z(i))
394 surf(2) = one/sqrt(n2x(i)*n2x(i)+n2y(i)*n2y(i)+n2z(i)*n2z(i))
395 surf(3) = one/sqrt(n3x(i)*n3x(i)+n3y(i)*n3y(i)+n3z(i)*n3z(i))
396 surf(4) = one/sqrt(n4x(i)*n4x(i)+n4y(i)*n4y(i)+n4z(i)*n4z(i))
397 surf(5) = one/sqrt(n5x(i)*n5x(i)+n5y(i)*n5y(i)+n5z(i)*n5z(i))
398 surf(6) = one/sqrt(n6x(i)*n6x(i)+n6y(i)*n6y(i)+n6z(i)*n6z(i))
399 DO j=1,nv46
400
401 term2 = +surf(j) * ( valel(6,i)-valvois(6,j,i) ) / (valel(4,i
402 vec(1,j) = (valel(1,i)*valel(5,i) + valvois(1,j,i)*valvois(5,j,i))
403 . / (valel(4,i)*valel(5,i)+valvois(4,j,i)*valvois(5,j,i)) + term2*veul(14-1+j,ii)
404 vec(2,j) = (valel(2,i)*valel(5,i) + valvois(2,j,i)*valvois(5,j,i))
405 . / (valel(4,i)*valel(5,i)+valvois(4,j,i)*valvois(5,j,i)) + term2*veul(20-1+j,ii)
406 vec(3,j) = (valel(3,i)*valel(5,i) + valvois(3,j,i)*valvois(5,j,i))
407 . / (valel(4,i)*valel(5,i)+valvois(4,j,i)*valvois(5,j
408 enddo
409 ENDIF
410 DO j=1,nv46
414 ENDDO
415
416 vx1(i) = vec(1,1)
417 vy1(i) = vec(2,1)
418 vz1(i) = vec(3,1)
419
420 vx2(i) = vec(1,2)
421 vy2(i) = vec
422 vz2(i) = vec(3,2)
423
424 vx3(i) = vec(1,3)
425 vy3(i) = vec(2,3)
426 vz3(i) = vec(3,3)
427
428 vx4(i) = vec(1,4)
429 vy4(i) = vec(2,4)
430 vz4(i) = vec(3,4)
431
432 vx5(i) = vec(1,5)
433 vy5(i) = vec(2,5)
434 vz5(i) = vec(3,5)
435
436 vx6(i) = vec(1,6)
437 vy6(i) = vec(2,6)
438 vz6(i) = vec(3,6)
439 enddo
440
441 CASE(6)
442
443
444 DO i=1,nel
445 ii = i + nft
446 nc(1,i) = ixs(2,ii)
447 nc(2,i) = ixs(3,ii)
448 nc(3,i) = ixs(4,ii)
449 nc(4,i) = ixs(5,ii)
450 nc(5,i) = ixs(6,ii)
451 nc(6,i) = ixs(7,ii)
452 nc(7,i) = ixs(8,ii)
453 nc(8,i) = ixs(9,ii)
454 z(1) = one_over_8*sum(x(1,nc(1:8,i)))
455 z(2) = one_over_8*sum(x(2,nc(1:8,i)))
456 z(3) = one_over_8*sum(x(3,nc(1:8,i)))
457 iad2 = ale_connect%ee_connect%iad_connect(ii)
458 lgth = ale_connect%ee_connect%iad_connect(ii+1) - iad2
459 DO j=1,lgth
460 idv = ale_connect%ee_connect%connected(iad2 + j - 1)
461 IF(idv<=0)THEN
462 vec(1:3,j) = zero
463 cycle
464 ENDIF
465 ix1=ixs(icf(1,j)+1,ii)
466 ix2=ixs(icf(2,j)+1,ii)
467 ix3=ixs(icf(3,j)+1,ii)
468 ix4=ixs(icf(4,j)+1,ii)
469 cf(1,i) = fourth*(x(1,ix1)+x(1,ix2)+x(1,ix3)+x(1,ix4))
470 cf(2,i) = fourth*(x(2,ix1)+x(2,ix2)+x(2,ix3)+x(2,ix4))
471 cf(3,i) = fourth*(x(3,ix1)+x(3,ix2)+x
472 nc(1,i)=ixs(2,idv)
473 nc(2,i)=ixs(3,idv)
474 nc(3,i)=ixs(4,idv)
475 nc(4,i)=ixs(5,idv)
476 nc(5,i)=ixs(6,idv)
477 nc(6,i)=ixs(7,idv)
478 nc(7,i)=ixs(8,idv)
479 nc(8,i)=ixs(9,idv)
480 zadj(1) = one_over_8*sum(x(1,nc(1:8,i)))
481 zadj(2) = one_over_8*sum(x(2,nc(1:8,i)))
482 zadj(3) = one_over_8
483 zzadj(1) = zadj(1)-z(1)
484 zzadj(2) = zadj(2)-z(2)
485 zzadj(3) = zadj(3)-z(3)
486 zcf(1) = cf(1,i) - z(1)
487 zcf(2) = cf(2,i) - z(2)
488 zcf(3) = cf(3,i) - z(3)
489 ps = zcf(1)*zzadj(1) + zcf(2)*zzadj(2) + zcf(3)*zzadj(3)
490 zzadj_ = zzadj(1)**2 + zzadj(2)**2 + zzadj(3)**2
491 lambda = ps /
max(em20,zzadj_)
492
493
494
495
496 lambda =
min(
max(zero,lambda) , one)
497 lambda = sin(half*3.14159265358979d00*lambda)
498 lambda = lambda * lambda
499
500 sr1 = valel(4,i)
501 sr2 = valvois(4,j,i)
502 srf = sr1 + lambda*(sr2-sr1)
503
504 vec(1,j) = valel(1,i) + lambda*(valvois(1,j,i)-valel(1,i))
505 vec(2,j) = valel(2,i) + lambda*(valvois(2,j,i)-valel(2,i))
506 vec(3,j) = valel(3,i) + lambda*(valvois(3,j,i)-valel(3,i))
507
508 vec(1,j) = vec(1,j) / srf
509 vec(2,j) = vec(2,j) / srf
510 vec(3,j) = vec(3,j) / srf
511 enddo
512 DO j=1,nv46
516 ENDDO
517
518 vx1(i) = vec(1,1)
519 vy1(i) = vec(2,1)
520 vz1(i) = vec(3,1)
521
522 vx2(i) = vec(1,2)
523 vy2(i) = vec(2,2)
524 vz2(i) = vec(3,2)
525
526 vx3(i) = vec(1,3)
527 vy3(i) = vec(2,3)
528 vz3(i) = vec(3,3)
529
530 vx4(i) = vec(1,4)
531 vy4(i) = vec(2,4)
532 vz4(i) = vec(3,4)
533
534 vx5(i) = vec(1,5)
535 vy5(i) = vec(2,5)
536 vz5(i) = vec(3,5)
537
538 vx6(i) = vec(1,6)
539 vy6(i) = vec(2,6)
540 vz6(i) = vec(3,6)
541 enddo
542
543 END SELECT
544
545
546
547
548
549 DO i=1,nel
550 flux1(i)=half*(vx1(i)*n1x(i)+vy1(i)*n1y(i)+vz1(i)*n1z(i))
551 flux2(i)=half*(vx2(i)*n2x(i)+vy2(i)*n2y(i)+vz2(i)*n2z(i))
552 flux3(i)=half*(vx3(i)*n3x(i)+vy3(i)*n3y(i)+vz3(i)*n3z(i))
553 flux4(i)=half*(vx4(i)*n4x(i)+vy4(i)*n4y(i)+vz4(i)*n4z(i))
554 flux5(i)=half*(vx5(i)*n5x(i)+vy5(i)*n5y(i)+vz5(i)*n5z(i))
555 flux6(i)=half*(vx6(i)*n6x(i)+vy6(i)*n6y(i)+vz6(i)*n6z(i))
556 ENDDO
557
558
559
560
562 debug_outp = .false.
564 do i=lft,llt
565 ii = nft + i
567 debug_outp = .true.
568 idbf = i
569 idbl = i
570 EXIT
571 endif
572 enddo
574 debug_outp=.true.
575 idbf = lft
576 idbl = llt
577 endif
578
579 if(debug_outp)then
580
581 print *, " |----alefvm_eflux3.F-----|"
582 print *, " | THREAD INFORMATION |"
583 print *, " |------------------------|"
584 print *, " NCYCLE =", ncycle
585 do i=idbf,idbl
586 ii = nft + i
587 print *, " brique=", ixs(11,nft+i)
588 write (*,fmt='(A,6E26.14)') " Flux(1:6) =", flux1(i),flux2(i),flux3(i),flux4(i),flux5(i),flux6(i)
589 write(*,fmt='(A24,1A20)') " ",
590 . "#--------- cell-----------"
591 write (*,fmt=
'(A,1E26.14)')
" V_cell-X =",
alefvm_buffer%FCELL(1,ii)
592 write (*,fmt=
'(A,1E26.14)')
" V_cell-Y =",
alefvm_buffer%FCELL(2,ii)
593 write (*,fmt=
'(A,1E26.14)')
" V_cell-Z =",
alefvm_buffer%FCELL(3,ii)
594 write(*,fmt='(A24,6A26)') " ",
595 . "#--------- adj_1 ---------","#--------- adj_2 ---------",
596 . "#--------- adj_3 ---------","#--------- adj_4 ---------",
597 . "#--------- adj_5 ---------","#--------- adj_6 --------#"
598 write (*,fmt=
'(A,6E26.14)')
" V_faces-X =",
alefvm_buffer%F_FACE(1,1:6,ii)
599 write (*,fmt=
'(A,6E26.14)')
" V_faces-Y =",
alefvm_buffer%F_FACE(2,1:6,ii)
600 write (*,fmt=
'(A,6E26.14)')
" V_faces-Z =",
alefvm_buffer%F_FACE(3,1:
601 print *, " "
602 enddo
603
604 endif
605 endif
606
607
608
609
610
611
612
613 IF(nint(pm(19,mat(1)))==51)THEN
614 DO i=1,nel
615 flux(1,i)=flux1(i)
616 flux(2,i)=flux2(i)
617 flux(3,i)=flux3(i)
618 flux(4,i)=flux4(i)
619 flux(5,i)=flux5(i)
620 flux(6,i)=flux6(i)
621 ENDDO
622 RETURN
623 ENDIF
624
625
626
627
628
629 DO j=1,6
630 DO i=1,nel
631 upwl(j,i)=pm(16,mat(i))
632 ENDDO
633 ENDDO
634
635 DO i=1,nel
636 iad2 = ale_connect%ee_connect%iad_connect(i + nft)
637 reduc=pm(92,mat(i))
638
639 ii=ale_connect%ee_connect%connected(iad2 + 1 - 1)
640 IF(ii==0)THEN
641 flux1(i)=flux1(i)*reduc
642 ENDIF
643
644 ii=ale_connect%ee_connect%connected(iad2 + 2 - 1)
645 IF(ii==0)THEN
646 flux2(i)=flux2(i)*reduc
647 ENDIF
648
649 ii=ale_connect%ee_connect%connected(iad2 + 3 - 1)
650 IF(ii==0)THEN
651 flux3(i)=flux3(i)*reduc
652 ENDIF
653
654 ii=ale_connect%ee_connect%connected(iad2 + 4 - 1)
655 IF(ii==0)THEN
656 flux4(i)=flux4(i)*reduc
657 ENDIF
658
659 ii=ale_connect%ee_connect%connected(iad2 + 5 - 1)
660 IF(ii==0)THEN
661 flux5(i)=flux5(i)*reduc
662 ENDIF
663
664 ii=ale_connect%ee_connect%connected(iad2 + 6 - 1)
665 IF(ii==0)THEN
666 flux6(i)=flux6(i)*reduc
667 ENDIF
668 ENDDO
669
670 DO i=1,nel
671 flux(1,i)=flux1(i)-upwl(1,i)*abs(flux1(i))
672 flux(2,i)=flux2(i)-upwl(2,i)*abs(flux2(i))
673 flux(3,i)=flux3(i)-upwl(3,i)*abs(flux3(i))
674 flux(4,i)=flux4(i)-upwl(4,i)*abs(flux4(i))
675 flux(5,i)=flux5(i)-upwl(5,i)*abs(flux5(i))
676 flux(6,i)=flux6(i)-upwl(6,i)*abs(flux6(i))
677
678 flu1(i) =flux1(i)+upwl(1,i)*abs(flux1(i))
679 . +flux2(i)+upwl(2,i)*abs(flux2(i))
680 . +flux3(i)+upwl(3,i)*abs(flux3(i))
681 . +flux4(i)+upwl(4,i)*abs(flux4(i))
682 . +flux5(i)+upwl(5,i)*abs(flux5(i))
683 . +flux6(i)+upwl(6,i)*abs(flux6(i))
684 ENDDO
685
686
687
688
689
690 RETURN
type(alefvm_buffer_), target alefvm_buffer
type(alefvm_param_), target alefvm_param