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