37
38
39
40
41
42
43
44
45
46
47
48
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 "com08_c.inc"
68
69
70
71 INTEGER :: IXS(NIXS,*), NV46, IPM(NPROPMI,*),NEL
72 my_real :: pm(npropm,*), flux(6,*), flu1(*), x(3,*),w(3,*)
73 TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECT
74
75
76
77 INTEGER MAT(MVSIZ), NC1(MVSIZ), NC2(MVSIZ), NC3(MVSIZ), NC4(MVSIZ), NC5(MVSIZ), NC6(MVSIZ)
78
79
80
81 . x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz), x5(mvsiz), x6(mvsiz), x7(mvsiz), x8(mvsiz),
82 . y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz), y5(mvsiz), y6(mvsiz), y7(mvsiz), y8(mvsiz), z1(mvsiz),
83 . z2(mvsiz), z3(mvsiz), z4(mvsiz), z5(mvsiz), z6(mvsiz), z7(mvsiz), z8(mvsiz), n1x(mvsiz),
84 . n2x(mvsiz), n3x(mvsiz), n4x(mvsiz), n5x(mvsiz), n6x(mvsiz), n1y(mvsiz), n2y(mvsiz), n3y(mvsiz),
85 . n4y(mvsiz), n5y(mvsiz), n6y(mvsiz), n1z(mvsiz), n2z(mvsiz), n3z(mvsiz), n4z(mvsiz), n5z(mvsiz),
86 . n6z(mvsiz), flux1(mvsiz), flux2(mvsiz), flux3(mvsiz), flux4(mvsiz), flux5(mvsiz),
87 . flux6(mvsiz), vx1(mvsiz), vx2(mvsiz), vx3(mvsiz), vx4(mvsiz), vx5(mvsiz), vx6(mvsiz),
88 . vy1(mvsiz), vy2(mvsiz), vy3(mvsiz), vy4(mvsiz), vy5(mvsiz), vy6(mvsiz), vz1(mvsiz), vz2(mvsiz),
89 . vz3(mvsiz), vz4(mvsiz), vz5(mvsiz), vz6(mvsiz),
90 . reduc,upwl(6,mvsiz),valvois(6,nv46,mvsiz), valel(6,mvsiz),vec(3,6),cf(3,mvsiz),
91 . z(3),zadj(3),zcf(3),zzadj(3),zcf_,zzadj_,cos, ps, denom, denom1,denom2,lambda,xsi,wface(3,6,mvsiz),sr1,sr2,
92 . surf(6), term2, n(3,6,mvsiz)
93 LOGICAL debug_outp
94 INTEGER :: idbf, idbl, NC(8,MVSIZ), IAD2, LGTH
95
96 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/
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
118 imat = ixs(1,1+nft)
119 ialefvm_flg = ipm(251,imat)
120 IF(ialefvm_flg <= 1)RETURN
121
122
123
124
125 DO i=1,nel
126 ii = i + nft
127 mat(i) = ixs(1,ii)
128
129 nc1(i)=ixs(2,ii)
130 nc2(i)=ixs(3,ii)
131 nc3(i)=ixs(4,ii)
132 nc4(i)=ixs(5,ii)
133 nc5(i)=ixs(6,ii)
134 nc6(i)=ixs(7,ii)
135 nc7(i)=ixs(8,ii)
136 nc8(i)=ixs(9,ii)
137
138
139 x1(i)=x(1,nc1(i))
140 y1(i)=x(2,nc1(i))
141 z1(i)=x(3,nc1(i))
142
143 x2(i)=x(1,nc2(i))
144 y2(i)=x(2,nc2(i))
145 z2(i)=x(3,nc2(i))
146
147 x3(i)=x(1,nc3(i))
148 y3(i)=x(2,nc3(i))
149 z3(i)=x(3,nc3(i))
150
151 x4(i)=x(1,nc4(i))
152 y4(i)=x(2,nc4(i))
153 z4(i)=x(3,nc4(i))
154
155 x5(i)=x(1,nc5(i))
156 y5(i)=x(2,nc5(i))
157 z5(i)=x(3,nc5(i))
158
159 x6(i)=x(1,nc6(i))
160 y6(i)=x(2,nc6(i))
161 z6(i)=x(3,nc6(i))
162
163 x7(i)=x(1,nc7(i))
164 y7(i)=x(2,nc7(i))
165 z7(i)=x(3,nc7(i))
166
167 x8(i)=x(1,nc8(i))
168 y8(i)=x(2,nc8(i))
169 z8(i)=x(3,nc8(i))
170 ENDDO
171
172
173
174
175 DO i=1,nel
176 wface(1,1,i) = fourth*(w(1,nc1(i))+w(1,nc2(i))+w(1,nc3(i))+w(1,nc4(i)))
177 wface(2,1,i) = fourth*(w(2,nc1(i))+w(2,nc2(i))+w(2,nc3(i))+w(2,nc4(i)))
178 wface(3,1,i) = fourth*(w(3,nc1(i))+w(3,nc2(i))+w(3,nc3(i))+w(3,nc4(i)))
179
180 wface(1,2,i) = fourth*(w(1,nc3(i))+w(1,nc4(i))+w(1,nc7(i))+w(1,nc8(i)))
181 wface(2,2,i) = fourth*(w(2,nc3(i))+w(2,nc4(i))+w(2,nc7(i))+w(2,nc8(i)))
182 wface(3,2,i) = fourth*(w(3,nc3(i))+w(3,nc4(i))+w(3,nc7(i))+w(3,nc8(i)))
183
184
185 wface(1,3,i) = fourth*(w(1,nc5(i))+w(1,nc6(i))+w(1,nc7(i))+w(1,nc8(i)))
186 wface(2,3,i) = fourth*(w(2,nc5(i))+w(2,nc6(i))+w(2,nc7(i))+w(2,nc8(i)))
187 wface(3,3,i) = fourth*(w(3,nc5(i))+w(3,nc6(i))+w(3,nc7(i))+w(3,nc8(i)))
188
189 wface(1,4,i) = fourth*(w(1,nc1(i))+w(1,nc2(i))+w(1,nc5(i))+w(1,nc6(i)))
190 wface(2,4,i) = fourth*(w(2,nc1(i))+w(2,nc2(i))+w(2,nc5(i))+w(2,nc6(i
191 wface(3,4,i) = fourth*(w(3,nc1(i))+w(3,nc2(i))+w(3,nc5(i))+w(3,nc6(i)))
192
193 wface(1,5,i) = fourth*(w(1,nc2(i))+w(1,nc3(i))+w(1,nc6(i))+w(1,nc7(i)))
194 wface(2,5,i) = fourth*(w(2,nc2(i))+w(2,nc3(i))+w(2,nc6(i))+w(2,nc7(i)))
195 wface(3,5,i) = fourth*(w(3,nc2(i))+w(3,nc3(i))+w(3,nc6(i))+w(3,nc7(i)))
196
197 wface(1,6,i) = fourth*(w(1,nc1(i))+w(1,nc4(i))+w(1,nc5(i))+w(1,nc8(i)))
198 wface(2,6,i) = fourth*(w(2,nc1(i))+w(2,nc4(i))+w(2,nc5(i))+w(2,nc8(i)))
199 wface(3,6,i) = fourth*(w(3,nc1(i))+w(3,nc4(i))+w(3,nc5(i))+w(3,nc8(i)))
200 ENDDO
201
202
203
204
205
206 DO i=1,nel
207 ii = i + nft
214 iad2 = ale_connect%ee_connect%iad_connect(ii)
215 lgth = ale_connect%ee_connect%iad_connect(ii+1) - iad2
216 DO j=1,lgth
217 idv = ale_connect%ee_connect%connected(iad2 + j - 1)
218 IF(idv <= 0) THEN
225 ELSE
232 ENDIF
233 enddo
234 enddo
235
236
237
238
239
240
241
242 DO i=1,nel
243
244 n1x(i)=(y3(i)-y1(i))*(z2(i)-z4(i)) - (z3(i)-z1(i))*(y2(i)-y4(i))
245 n1y(i)=(z3(i)-z1(i))*(x2(i)-x4(i)) - (x3(i)-x1(i))*(z2(i)-z4(i))
246 n1z(i)=(x3(i)-x1(i))*(y2(i)-y4(i)) - (y3(i)-y1(i))*(x2(i)-x4(i))
247
248 n2x(i)=(y7(i)-y4(i))*(z3(i)-z8(i)) - (z7(i)-z4(i
249 n2y(i)=(z7(i)-z4(i))*(x3(i)-x8
250 n2z(i)=(x7(i)-x4(i))*(y3(i)-y8(i)) -
251
252 n3x(i)=(y6(i)-y8(i))*(z7(i)-z5(i))
253 n3y(i)=(z6(i)-z8(i))*(x7(i)-x5
254 n3z(i)=(x6(i)-x8(i))*(y7(i)-y5(i)) - (y6(i)-y8(i))*(x7(i)-x5(i))
255
256 n4x(i)=(y2(i)-y5(i))*(z6(i)-z1(i)) - (z2(i)-z5(i))*(y6(i)-y1(i))
257 n4y(i)=(z2(i)-z5(i))*(x6(i)-x1(i)) - (x2(i)-x5(i))*(z6(i)-z1(i))
258 n4z(i)=(x2(i)-x5(i))*(y6(i)-y1(i)) - (y2(i)-y5(i))*(x6(i)-x1(i))
259
260 n5x(i)=(y7(i)-y2(i))*(z6(i)-z3(i)) - (z7(i)-z2(i))*(y6(i)-y3
261 n5y(i)=(z7(i)-z2(i))*(x6(i)-x3(i)) - (x7(i)-x2(i))*(z6(i)-z3(i))
262 n5z(i)=(x7(i)-x2(i))*(y6(i)-y3(i)) - (y7(i)-y2(i))*(x6(i)-x3(i))
263
264 n6x(i)=(y8(i)-y1(i))*(z4(i)-z5(i)) - (z8(i)-z1(i))*(y4(i)-y5(i))
265 n6y(i)=(z8(i)-z1(i))*(x4(i)-x5(i)) - (x8(i)-x1(i))*(z4(i)-z5(i))
266 n6z(i)=(x8(i)-x1(i))*(y4(i)-y5(i)) - (y8(i)-y1(i))*(x4(i)-x5(i))
267
268 n(1,1,i) = n1x(i)
269 n(2,1,i) = n1y(i)
270 n(3,1,i) = n1z(i)
271
272 n(1,2,i) = n2x(i)
273 n(2,2,i) = n2y(i)
274 n(3,2,i) = n2z(i)
275 !
276 n(1,3,i) = n3x(i)
277 n(2,3,i) = n3y(i)
278 n(3,3,i) = n3z(i)
279
280 n(1,4,i) = n4x(i)
281 n(2,4,i) = n4y(i)
282 n(3,4,i) = n4z(i)
283
284 n(1,5,i) = n5x(i)
285 n(2,5,i) = n5y(i)
286 n(3,5,i) = n5z(i)
287
288 n(1,6,i) = n6x(i)
289 n(2,6,i) = n6y(i)
290 n(3,6,i) = n6z(i)
291 END DO
292
293
294
295
296
297
298
300
302
303
304 CASE(1)
305
306 DO i=1,nel
307 ii = i + nft
308 DO j=1,nv46
309 vec(1,j) = half * (valel(1,i)/valel(4,i) + valvois(1,j,i)/valvois(4,j,i))
310 vec(2,j) = half * (valel(2,i)/valel(4,i) + valvois(2,j,i)/valvois(
311 vec(3,j) = half * (valel(3,i)/valel(4,i) + valvois(3,j,i)/valvois(4,j,i
312 vec(1,j) = vec(1,j) - wface(1,j,i)
313 vec(2,j) = vec(2,j) - wface(2,j,i)
314 vec(3,j) = vec(3,j) - wface(3,j,i)
315 enddo
316 DO j=1,nv46
320 ENDDO
321
322 vx1(i) = vec(1,1)
323 vy1(i) = vec(2,1)
324 vz1(i) = vec(3,1)
325
326 vx2(i) = vec(1,2)
327 vy2(i) = vec(2,2)
328 vz2(i) = vec(3,2)
329
330 vx3(i) = vec(1,3)
331 vy3(i) = vec(2,3)
332 vz3(i) = vec(3,3)
333
334 vx4(i) = vec(1,4)
335 vy4(i) = vec(2,4)
336 vz4(i) = vec(3,4)
337
338 vx5(i) = vec(1,5)
339 vy5(i) = vec(2,5)
340 vz5(i) = vec(3,5)
341
342 vx6(i) = vec(1,6)
343 vy6(i) = vec(2,6)
344 vz6(i) = vec(3,6)
345 enddo
346
347 CASE(2)
348
349 DO i=1,nel
350 ii = i + nft
351 DO j=1,nv46
352 vec(1,j) = (valel(1,i) + valvois(1,j,i)) / (valel(4,i)+valvois(4,j,i))
353 vec(2,j) = (valel(2,i) + valvois(2,j,i)) / (valel(4,i)+valvois(4,j,i))
354 vec(3,j) = (valel(3,i) + valvois(3,j,i)) / (valel(4,i)+valvois(4,j,i))
355 vec(1,j) = vec(1,j) - wface(1,j,i)
356 vec(2,j) = vec(2,j) - wface(2,j,i)
357 vec(3,j) = vec(3,j) - wface(3,j,i)
358 enddo
359 DO j=1,nv46
363 ENDDO
364
365 vx1(i) = vec(1,1)
366 vy1(i) = vec(2,1)
367 vz1(i) = vec(3,1)
368
369 vx2(i) = vec(1,2)
370 vy2(i) = vec(2,2)
371 vz2(i) = vec(3,2)
372
373 vx3(i) = vec(1,3)
374 vy3(i) = vec(2,3)
375 vz3(i) = vec(3,3)
376
377 vx4(i) = vec(1,4)
378 vy4(i) = vec(2,4)
379 vz4(i) = vec(3,4)
380
381 vx5(i) = vec(1,5)
382 vy5(i) = vec(2,5)
383 vz5(i) = vec(3,5)
384
385 vx6(i) = vec(1,6)
386 vy6(i) = vec(2,6)
387 vz6(i) = vec(3,6)
388 enddo
389
390 CASE(3)
391
392 DO i=1,nel
393 ii = i + nft
394 sr1 = sqrt(valel(4,i))
395 DO j=1,nv46
396 sr2 = sqrt(valvois(4,j,i))
397 vec(1,j) = (valel(1,i)/sr1 + valvois(1,j,i)/sr2) / (sr1+sr2)
398 vec(2,j) = (valel(2,i)/sr1 + valvois(2,j,i)/sr2) / (sr1+sr2)
399 vec(3,j) = (valel(3,i)/sr1 + valvois(3,j,i)/sr2) / (sr1+sr2)
400 vec(1,j) = vec(1,j) - wface(1,j,i)
401 vec(2,j) = vec(2,j) - wface(2,j,i)
402 vec(3,j) = vec(3,j) - wface(3,j,i)
403 enddo
404 DO j=1,nv46
408 ENDDO
409
410 vx1(i) = vec(1,1)
411 vy1(i) = vec(2,1)
412 vz1(i) = vec(3,1)
413
414 vx2(i) = vec(1,2)
415 vy2(i) = vec(2,2)
416 vz2(i) = vec(3,2)
417
418 vx3(i) = vec(1,3)
419 vy3(i) = vec(2,3)
420 vz3(i) = vec(3,3)
421
422 vx4(i) = vec(1,4)
423 vy4(i) = vec(2,4)
424 vz4(i) = vec(3,4)
425
426 vx5(i) = vec(1,5)
427 vy5(i) = vec(2,5)
428 vz5(i) = vec(3,5)
429
430 vx6(i) = vec(1,6)
431 vy6(i) = vec(2,6)
432 vz6(i) = vec(3,6)
433 enddo
434
435 CASE(4)
436
437 DO i=1,nel
438 ii = i + nft
439 DO j=1,nv46
440 IF(dt1==zero)THEN
441 vec(1,j) = (valel(1,i) + valvois(1,j,i) )
442 . / (valel(4,i) + valvois(4,j,i) )
443 vec(2,j) = (valel(2,i) + valvois(2,j,i) )
444 . / (valel(4,i) + valvois(4,j,i) )
445 vec(3,j) = (valel(3,i) + valvois(3,j,i) )
446 . / (valel(4,i) + valvois(4,j,i) )
447 ELSE
448 vec(1,j) = (valel(1,i)*valel(5,i) + valvois(1,j,i)*valvois(5,j,i))
449 . / (valel(4,i)*valel(5,i) + valvois(4,j,i)*valvois(5,j,i))
450 vec(2,j) = (valel(2,i)*valel(5,i) + valvois(2,j,i)*valvois(5,j,i))
451 . / (valel(4,i)*valel(5,i) + valvois(4,j,i)*valvois(5,j,i))
452 vec(3,j) = (valel(3,i)*valel(5,i) + valvois(3,j,i)*valvois(5,j,i))
453 . / (valel(4,i)*valel(5,i) + valvois(4,j,i)*valvois(5,j,i))
454 ENDIF
455 vec(1,j) = vec(1,j) - wface(1,j,i)
456 vec(2,j) = vec(2,j) - wface(2,j,i)
457 vec(3,j) = vec(3,j) - wface(3,j,i)
458 enddo
459 DO j=1,nv46
463 ENDDO
464
465 vx1(i) = vec(1,1)
466 vy1(i) = vec(2,1)
467 vz1(i) = vec(3,1)
468
469 vx2(i) = vec(1,2)
470 vy2(i) = vec(2,2)
471 vz2(i) = vec(3,2)
472
473 vx3(i) = vec(1,3)
474 vy3(i) = vec(2,3)
475 vz3(i) = vec(3,3)
476
477 vx4(i) = vec(1,4)
478 vy4(i) = vec(2,4)
479 vz4(i) = vec(3,4)
480
481 vx5(i) = vec(1,5)
482 vy5(i) = vec(2,5)
483 vz5(i) = vec(3,5)
484
485 vx6(i) = vec(1,6)
486 vy6(i) = vec(2,6)
487 vz6(i) = vec(3,6)
488 enddo
489
490 CASE(5)
491
492 DO i=1,nel
493 ii = i + nft
494 IF(dt1==zero)THEN
495 DO j=1,nv46
496 vec(1,j) = (valel(1,i) + valvois
497 . / (valel(4,i) + valvois(4,j,i))
498 vec(2,j) = (valel(2,i) + valvois(2,j,i))
499 . / (valel(4,i) + valvois(4,j,i))
500 vec(3,j) = (valel(3,i) + valvois(3,j,i))
501 . / (valel(4,i) + valvois(4,j,i))
502 vec(1,j) = vec(1,j) - wface(1,j,i)
503 vec(2,j) = vec(2,j) - wface(2,j,i)
504 vec(3,j) = vec(3,j) - wface(3,j
505 enddo
506 ELSE
507 surf(1) = one/sqrt(n1x(i)*n1x(i)+n1y(i)*n1y(i)+n1z(i)*n1z(i))
508 surf(2) = one/sqrt(n2x(i)*n2x(i)+n2y(i)*n2y(i)+n2z(i)*n2z(i))
509 surf(3) = one/sqrt(n3x(i)*n3x(i)+n3y(i)*n3y
510 surf(4) = one/sqrt(n4x(i)*n4x(i)+n4y(i)*n4y(i)+n4z(i)*n4z(i))
511 surf(5) = one/sqrt(n5x(i)*n5x(i)+n5y(i)*n5y(i)+n5z(i)*n5z(i))
512 surf(6) = one/sqrt(n6x(i)*n6x(i)+n6y(i)*n6y(i)+n6z(i)*n6z(i))
513 DO j=1,nv46
514 term2 = surf(j) * ( valel(6,i)-valvois(6,j,i) ) / (valel(4,i)*valel(5,i)+valvois(4,j,i)*valvois(5,j,i))
515
516 vec(1,j) = (valel(1,i)*valel(5,i) + valvois(1,j,i)*valvois(5,j,i))
517 . / (valel(4,i)*valel(5,i)+valvois(4,j,i)*valvois(5,j,i)) + term2*n(1,j,i)
518 vec(2,j) = (valel(2,i)*valel(5,i) + valvois(2,j,i)*valvois(5,j,i))
519 . / (valel(4,i)*valel(5,i)+valvois(4,j,i)*valvois(5,j,i)) + term2*n(2,j,i)
520 vec(3,j) = (valel(3,i)*valel(5,i) + valvois(3,j,i)*valvois(5,j,i))
521 . / (valel(4,i)*valel(5,i)+valvois(4,j,i)*valvois(5,j,i)) + term2*n(3,j,i)
522 vec(1,j) = vec(1,j) - wface(1,j,i)
523 vec(2,j) = vec(2,j) - wface(2,j,i)
524 vec(3,j) = vec(3,j) - wface(3,j,i)
525 enddo
526 ENDIF
527 DO j=1,nv46
531 ENDDO
532 !---face-1
533 vx1(i) = vec(1,1)
534 vy1(i) = vec(2,1)
535 vz1(i) = vec(3,1)
536
537 vx2(i) = vec(1,2)
538 vy2(i) = vec(2,2)
539 vz2(i) = vec(3,2)
540
541 vx3(i) = vec(1,3)
542 vy3(i) = vec(2,3)
543 vz3(i) = vec(3,3)
544
545 vx4(i) = vec(1,4)
546 vy4(i) = vec(2,4)
547 vz4(i) = vec(3,4)
548
549 vx5(i) = vec(1,5)
550 vy5(i) = vec(2,5)
551 vz5(i) = vec(3,5)
552
553 vx6(i) = vec(1,6)
554 vy6(i) = vec(2,6)
555 vz6(i) = vec(3,6)
556 enddo
557
558 CASE(6)
559
560
561 DO i=1,nel
562 ii = i + nft
563 nc(1,i) = ixs(2,ii)
564 nc(2,i) = ixs(3,ii)
565 nc(3,i) = ixs(4,ii)
566 nc(4,i) = ixs(5,ii)
567 nc(5,i) = ixs(6,ii)
568 nc(6,i) = ixs(7,ii)
569 nc(7,i) = ixs(8,ii)
570 nc(8,i) = ixs(9,ii)
571 z(1) = one_over_8*sum(x(1,nc(1:8,i)))
572 z(2) = one_over_8*sum(x(2,nc(1:8,i)))
573 z(3) = one_over_8*sum(x(3,nc(1:8,i)))
574 iad2 = ale_connect%ee_connect%iad_connect(ii)
575 lgth = ale_connect%ee_connect%iad_connect(ii + 1) - iad2
576 DO j=1,lgth
577 idv = ale_connect%ee_connect%connected(iad2 + j - 1)
578 IF(idv<=0)THEN
579 vec(1:3,j) = zero
580 cycle
581 ENDIF
582 ix1=ixs(icf(1,j)+1,ii)
583 ix2=ixs(icf(2,j)+1,ii)
584 ix3=ixs(icf(3,j)+1,ii)
585 ix4=ixs(icf(4,j)+1,ii)
586 cf(1,i) = fourth*(x(1,ix1)+x(1,ix2)+x(1,ix3)+x(1,ix4))
587 cf(2,i) = fourth*(x(2,ix1)+x(2,ix2)+x(2,ix3)+x(2,ix4))
588 cf(3,i) = fourth*(x(3,ix1)+x(3,ix2)+x(3,ix3)+x(3,ix4))
589 nc(1,i)=ixs(2,idv)
590 nc(2,i)=ixs(3,idv)
591 nc(3,i)=ixs(4,idv)
592 nc(4,i)=ixs(5,idv)
593 nc(5,i)=ixs(6,idv)
594 nc(6,i)=ixs(7,idv)
595 nc(7,i)=ixs(8,idv)
596 nc(8,i)=ixs(9,idv)
597 zadj(1) = one_over_8*sum(x(1,nc(1:8,i)))
598 zadj(2) = one_over_8*sum(x(2,nc(1:8,i)))
599 zadj(3) = one_over_8*sum(x(3,nc(1:8,i)))
600 zzadj(1) = zadj(1)-z(1)
601 zzadj(2) = zadj(2)-z(2)
602 zzadj(3) = zadj(3)-z(3)
603 zcf(1) = cf(1,i) - z(1)
604 zcf(2) = cf(2,i) - z(2)
605 zcf(3) = cf(3,i) - z(3)
606 ps = zcf(1)*zzadj(1) + zcf(2)*zzadj(2) + zcf(3)*zzadj(3)
607 zcf_ = zcf(1)**2 + zcf(2)**2 + zcf(3)**2
608 zzadj_ = zzadj(1)**2 + zzadj(2)**2 + zzadj(3)**2
609 denom1 = sqrt(zcf_)
610 denom2 = sqrt(zzadj_)
611 denom = denom1*denom2
612 cos = ps/denom
613 xsi = sqrt(zcf_) * cos
614 lambda = xsi / denom2
615
616 vec(1,j) = valel(1,i) + lambda*(valvois(1,j,i)/valvois(4,j,i)-valel(1,i)/valel(4,i)) - wface(1,j,i)
617 vec(2,j) = valel(2,i) + lambda*(valvois(2,j,i)/valvois(4,j,i)-valel
618 vec(3,j) = valel(3,i) + lambda*(valvois(3,j,i)/valvois(4,j,i)-valel(3,i)/valel(4,i)) - wface(3,j,i)
619 enddo
620 DO j=1,nv46
624 ENDDO
625
626 vx1(i) = vec(1,1)
627 vy1(i) = vec(2,1)
628 vz1(i) = vec(3,1)
629
630 vx2(i) = vec(1,2)
631 vy2(i) = vec(2,2)
632 vz2(i) = vec(3,2)
633
634 vx3(i) = vec(1,3)
635 vy3(i) = vec(2,3)
636 vz3(i) = vec(3,3)
637
638 vx4(i) = vec(1,4)
639 vy4(i) = vec(2,4)
640 vz4(i) = vec(3,4)
641
642 vx5(i) = vec(1,5)
643 vy5(i) = vec(2,5)
644 vz5(i) = vec(3,5)
645
646 vx6(i) = vec(1,6)
647 vy6(i) = vec(2,6)
648 vz6(i) = vec(3,6)
649 enddo
650
651 END SELECT
652
653
654
655
656
657 DO i=1,nel
658 flux1(i)=half*(vx1(i)*n1x(i)+vy1(i)*n1y(i)+vz1(i)*n1z(i))
659 flux2(i)=half*(vx2(i)*n2x(i)+vy2(i)*n2y(i)+vz2(i)*n2z(i))
660 flux3(i)=half*(vx3(i)*n3x(i)+vy3(i)*n3y(i)+vz3(i)*n3z(i))
661 flux4(i)=half*(vx4(i)*n4x(i)+vy4(i)*n4y(i)+vz4(i)*n4z(i))
662 flux5(i)=half*(vx5(i)*n5x(i)+vy5(i)*n5y(i)+vz5(i)*n5z(i))
663 flux6(i)=half*(vx6(i)*n6x(i)+vy6(i)*n6y(i)+vz6(i)*n6z(i))
664 ENDDO
665
666
668 debug_outp = .false.
670 do i=lft,llt
671 ii = nft + i
673 debug_outp = .true.
674 idbf = i
675 idbl = i
676 EXIT
677 endif
678 enddo
680 debug_outp=.true.
681 idbf = lft
682 idbl = llt
683 endif
684
685 if(debug_outp)then
686
687 print *, " |----alefvm_aflux3.F-----|"
688 print *, " | THREAD INFORMATION |"
689 print *, " |------------------------|"
690 print *, " NCYCLE =", ncycle
691 do i=idbf,idbl
692 ii = nft + i
693 print *, " brique=", ixs(11,nft+i)
694 write (*,fmt='(A,6E26.14)') " Flux(1:6) =", flux1(i),flux2(i),flux3(i),flux4(i),flux5(i),flux6
695 write(*,fmt='(A24,1A20)') " ",
696 . "#--------- cell-----------"
697 write (*,fmt=
'(A,1E26.14)')
" V_cell-X =",
alefvm_buffer%FCELL(1,ii)
698 write (*,fmt=
'(A,1E26.14)')
" V_cell-Y =",
alefvm_buffer%FCELL(2,ii)
699 write (*,fmt=
'(A,1E26.14)')
" V_cell-Z =",
alefvm_buffer%FCELL(3,ii)
700 write(*,fmt='(A24,6A26)') " ",
701 . "#--------- adj_1 ---------","#--------- adj_2 ---------",
702 . "#--------- adj_3 ---------","#--------- adj_4 ---------",
703 . "#--------- adj_5 ---------","#--------- adj_6 --------#"
704 write (*,fmt=
'(A,6E26.14)')
" V_faces-X =",
alefvm_buffer%F_FACE(1,1:6,ii)
705 write (*,fmt=
'(A,6E26.14)')
" V_faces-Y =",
alefvm_buffer%F_FACE(2,1:6,ii)
706 write (*,fmt=
'(A,6E26.14)')
" V_faces-Z =",
alefvm_buffer%F_FACE(3,1:6,ii)
707 print *, " "
708 enddo
709
710 endif
711 endif
712
713
714
715
716
717
718
719 IF(nint(pm(19,mat(1)))==51)THEN
720 DO i=1,nel
721 flux(1,i)=flux1(i)
722 flux(2,i)=flux2(i)
723 flux(3,i)=flux3(i)
724 flux(4,i)=flux4(i)
725 flux(5,i)=flux5(i)
726 flux(6,i)=flux6(i)
727 ENDDO
728 RETURN
729 ENDIF
730
731
732
733
734
735 DO j=1,6
736 DO i=1,nel
737 upwl(j,i)=pm(16,mat(i))
738 ENDDO
739 ENDDO
740
741 DO i=1,nel
742 reduc=pm(92,mat(i))
743 iad2 = ale_connect%ee_connect%iad_connect(i + nft)
744
745 ii=ale_connect%ee_connect%connected(iad2 + 1 - 1)
746 IF(ii==0)THEN
747 flux1(i)=flux1(i)*reduc
748 ENDIF
749
750 ii=ale_connect%ee_connect%connected(iad2 + 2 - 1)
751 IF(ii==0)THEN
752 flux2(i)=flux2(i)*reduc
753 ENDIF
754
755 ii=ale_connect%ee_connect%connected(iad2 + 3 - 1)
756 IF(ii==0)THEN
757 flux3(i)=flux3(i)*reduc
758 ENDIF
759
760 ii=ale_connect%ee_connect%connected(iad2 + 4 - 1)
761 IF(ii==0)THEN
762 flux4(i)=flux4(i)*reduc
763 ENDIF
764
765 ii=ale_connect%ee_connect%connected(iad2 + 5 - 1)
766 IF(ii==0)THEN
767 flux5(i)=flux5(i)*reduc
768 ENDIF
769
770 ii=ale_connect%ee_connect%connected(iad2 + 6 - 1)
771 IF(ii==0)THEN
772 flux6(i)=flux6(i)*reduc
773 ENDIF
774 ENDDO
775
776 DO i=1,nel
777 flux(1,i)=flux1(i)-upwl(1,i)*abs(flux1(i))
778 flux(2,i)=flux2(i)-upwl(2,i)*abs(flux2(i))
779 flux(3,i)=flux3(i)-upwl(3,i)*abs(flux3(i))
780 flux(4,i)=flux4(i)-upwl(4,i)*abs(flux4(i))
781 flux(5,i)=flux5(i)-upwl(5,i)*abs(flux5(i))
782 flux(6,i)=flux6(i)-upwl(6,i)*abs(flux6(i))
783
784 flu1(i) =flux1(i)+upwl(1,i)*abs(flux1(i))
785 . +flux2(i)+upwl(2,i)*abs(flux2(i))
786 . +flux3(i)+upwl(3,i)*abs(flux3(i))
787 . +flux4(i)+upwl(4,i)*abs(flux4(i))
788 . +flux5(i)+upwl(5,i)*abs(flux5(i))
789 . +flux6(i)+upwl(6,i)*abs(flux6(i))
790 ENDDO
791
792
793
794
795
796 RETURN
type(alefvm_buffer_), target alefvm_buffer
type(alefvm_param_), target alefvm_param