SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
fem_sparsekit.f90
Go to the documentation of this file.
1 MODULE fem_s_m
2 
3 CONTAINS
4 
5 
6  SUBROUTINE qs_00_m(mesh, alpha, ia, ja, a0)
7  !=================================================
8 
9  ! alpha < w, _ > ===> A0 incremental accumulation
10 
11  USE gauss_points
12 
13  IMPLICIT NONE
14 
15  REAL(KIND=8), INTENT(IN) :: alpha
16  INTEGER, DIMENSION(:), INTENT(IN) :: ia
17  INTEGER, DIMENSION(:), INTENT(IN) :: ja
18  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
19 
20  INTEGER :: m, l, ni, nj, i, j, p
21  REAL(KIND=8) :: al, x
22 
23  TYPE(mesh_type), TARGET :: mesh
24  INTEGER, DIMENSION(:,:), POINTER :: jj
25  INTEGER, POINTER :: me
26 
27  CALL gauss(mesh)
28  jj => mesh%jj
29  me => mesh%me
30 
31  DO m = 1, me
32  DO l = 1, l_g
33 
34  al = alpha * rj(l,m)
35 
36  DO nj = 1, n_w; j = jj(nj, m)
37  DO ni = 1, n_w; i = jj(ni, m)
38 
39  ! IF (j >= i) THEN
40  x = ww(nj,l) * al * ww(ni,l)
41  DO p = ia(i), ia(i+1) - 1
42  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
43  ENDIF
44  ENDDO
45  ! ENDIF
46 
47  ENDDO
48  ENDDO
49 
50  ENDDO
51  ENDDO
52 
53  END SUBROUTINE qs_00_m
54 
55  SUBROUTINE qs_100_m(mesh, vv, ia, ja, a0)
56  !=================================================
57 
58  ! < D w, vv _ > ===> A0 incremental accumulation
59 
60  USE gauss_points
61 
62  IMPLICIT NONE
63 
64  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: vv
65  INTEGER, DIMENSION(:), INTENT(IN) :: ia
66  INTEGER, DIMENSION(:), INTENT(IN) :: ja
67  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
68 
69  TYPE(mesh_type), TARGET :: mesh
70  INTEGER, DIMENSION(:,:), POINTER :: jj
71  INTEGER, POINTER :: me
72 
73  INTEGER :: m, l, ni, nj, i, j, p, k
74  REAL(KIND=8) :: x
75  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: vl
76 
77  CALL gauss(mesh)
78  jj => mesh%jj
79  me => mesh%me
80 
81  DO m = 1, me
82  DO l = 1, l_g
83 
84  vl = 0
85  DO ni = 1, n_w
86  vl = vl + vv(:,jj(ni,m)) * ww(ni,l)
87  END DO
88 
89 
90  vl = vl * rj(l,m)
91 
92  DO nj = 1, n_w; j = jj(nj, m)
93  DO ni = 1, n_w; i = jj(ni, m)
94 
95  x = 0
96  DO k = 1, k_d
97  x = x + vl(k) * dw(k, ni, l, m)
98  END DO
99  x = x * ww(nj,l)
100 
101  DO p = ia(i), ia(i+1) - 1
102  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
103  ENDIF
104  ENDDO
105 
106  ENDDO
107  ENDDO
108 
109  ENDDO
110  ENDDO
111 
112  END SUBROUTINE qs_100_m
113 
114  !------------------------------------------------------------------------------
115  SUBROUTINE qs_000_m(mesh, ff, ia, ja, a0)
116  !=================================================
117 
118  ! alpha < w, f _ > ===> A0 incremental accumulation
119 
120  USE gauss_points
121 
122  IMPLICIT NONE
123 
124  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
125  INTEGER, DIMENSION(:), INTENT(IN) :: ia
126  INTEGER, DIMENSION(:), INTENT(IN) :: ja
127  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
128 
129  INTEGER :: m, l, ni, nj, n, i, j, p
130  REAL(KIND=8) :: al, x , ffl
131 
132  TYPE(mesh_type), TARGET :: mesh
133  INTEGER, DIMENSION(:,:), POINTER :: jj
134  INTEGER, POINTER :: me
135 
136  CALL gauss(mesh)
137  jj => mesh%jj
138  me => mesh%me
139 
140  DO m = 1, me
141 
142  DO l = 1, l_g
143 
144  ffl =0.d0
145  DO n = 1, n_w
146  ffl = ffl + ff(jj(n,m)) * ww(n,l)
147  ENDDO
148 
149  al = ffl * rj(l,m)
150 
151  DO nj = 1, n_w; j = jj(nj, m)
152  DO ni = 1, n_w; i = jj(ni, m)
153 
154  ! IF (j >= i) THEN
155  x = ww(nj,l) * al * ww(ni,l)
156  DO p = ia(i), ia(i+1) - 1
157  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
158  ENDIF
159  ENDDO
160  ! ENDIF
161 
162  ENDDO
163  ENDDO
164 
165  ENDDO
166  ENDDO
167 
168  END SUBROUTINE qs_000_m
169 
170 
171  SUBROUTINE qs_11_m (mesh, alpha, ia, ja, a0)
172  !=================================================
173 
174  ! alpha << (Dw), (D_) >> ===> A0 incremental accumulation
175 
176  USE gauss_points
177 
178  IMPLICIT NONE
179 
180  REAL(KIND=8), INTENT(IN) :: alpha
181  INTEGER, DIMENSION(:), INTENT(IN) :: ia
182  INTEGER, DIMENSION(:), INTENT(IN) :: ja
183  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
184 
185  INTEGER :: m, l, ni, nj, i, j, p
186  REAL(KIND=8) :: al, x
187 
188  TYPE(mesh_type), TARGET :: mesh
189  INTEGER, DIMENSION(:,:), POINTER :: jj
190  INTEGER, POINTER :: me
191 
192  CALL gauss(mesh)
193  jj => mesh%jj
194  me => mesh%me
195 
196  DO m = 1, me
197 
198  DO l = 1, l_g
199 
200  al = alpha * rj(l,m)
201 
202  DO nj = 1, n_w; j = jj(nj, m)
203  DO ni = 1, n_w; i = jj(ni, m)
204 
205  ! IF (j >= i) THEN
206  x = al * sum(dw(:,nj,l,m) * dw(:,ni,l,m))
207  DO p = ia(i), ia(i+1) - 1
208  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
209  ENDIF
210  ENDDO
211  ! ENDIF
212 
213  ENDDO
214  ENDDO
215 
216  ENDDO
217  ENDDO
218 
219  END SUBROUTINE qs_11_m
220 
221  SUBROUTINE qs_11_bb_p1_m (mesh, alpha, bcd)
222  !=================================================
223 
224  ! alpha << (Dw_h^H), (D(_)_h^H) >> ===> bcd incremental accumulation
225 
226  USE gauss_points
227 
228  IMPLICIT NONE
229 
230  REAL(KIND=8), INTENT(IN) :: alpha
231  REAL(KIND=8), DIMENSION(:,:), INTENT(INOUT) :: bcd
232 
233  INTEGER :: m, m_mother, l, ni, nj, i, j
234  REAL(KIND=8) :: al, x, mest
235 
236  TYPE(mesh_type), TARGET :: mesh
237  INTEGER, DIMENSION(:,:), POINTER :: jj
238  INTEGER, POINTER :: me
239 
240  CALL gauss(mesh)
241  jj => mesh%jj
242  me => mesh%me
243 
244  DO m = 1, me
245 
246  m_mother = (m-1)/n_w + 1
247 
248  mest = 0
249  DO l = 1, l_g
250  mest = mest + rj(l,m)
251  END DO
252  mest = alpha*(n_w*mest)**(1.d0/k_d)
253 
254  DO l = 1, l_g
255 
256  al = mest * rj(l,m)
257 
258  nj = n_w; j = jj(nj, m)
259  ni = n_w; i = jj(ni, m)
260 
261  x = al * sum(dw(:,nj,l,m) * dw(:,ni,l,m))
262  bcd(m_mother,n_w+1) = bcd(m_mother,n_w+1) + x
263 
264  ENDDO
265  ENDDO
266 
267  END SUBROUTINE qs_11_bb_p1_m
268 
269 
270  !=====================================================
271  SUBROUTINE qs_v_grad_v_grad_w (mesh, gg, ia, ja, a0)
272  !=====================================================
273 
274  ! << (g.D)w, (g.D)_ >> ===> A0 incremental accumulation
275 
276  USE gauss_points
277 
278  IMPLICIT NONE
279  !--------------------------------------------------------------------------
280  ! Declaration des variables globales
281  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
282  INTEGER, DIMENSION(:), INTENT(IN) :: ia
283  INTEGER, DIMENSION(:), INTENT(IN) :: ja
284  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
285 
286  TYPE(mesh_type), TARGET :: mesh
287  INTEGER, DIMENSION(:,:), POINTER :: jj
288  INTEGER, POINTER :: me
289 
290  !--------------------------------------------------------------------------
291  ! Declaration des variables locales
292  INTEGER :: l, k, ni, nj, p, m, i, j
293  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: gl
294  REAL(KIND=8) :: x
295  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: y
296 
297  CALL gauss(mesh)
298  jj => mesh%jj
299  me => mesh%me
300 
301  boucle_mm : DO m = 1, me
302 
303  boucle_l : DO l = 1, l_g
304 
305  gl = 0
306  boucle_k : DO k = 1, k_d
307  boucle_ni : DO ni =1 ,n_w
308  gl(k) = gl(k) + gg(k, jj(ni,m)) * ww(ni,l)
309  END DO boucle_ni
310  ENDDO boucle_k
311 
312  y = 0
313  boucle_ni_2 : DO ni = 1, n_w
314  boucle_k_2 : DO k = 1, k_d
315  y(ni) = y(ni) + gl(k) * dw(k,ni,l,m)
316  END DO boucle_k_2
317  END DO boucle_ni_2
318 
319  boucle_ni_3 : DO ni = 1, n_w
320  i = jj(ni, m)
321  boucle_nj : DO nj = 1, n_w
322  j = jj(nj, m)
323  x = rj(l,m) * y(nj) * y(ni)
324  boucle_p : DO p = ia(i), ia(i+1) - 1
325  IF (ja(p) == j) THEN
326  a0(p) = a0(p) + x
327  EXIT
328  ENDIF
329  ENDDO boucle_p
330  ENDDO boucle_nj
331  ENDDO boucle_ni_3
332 
333  ENDDO boucle_l
334 
335  ENDDO boucle_mm
336 
337  END SUBROUTINE qs_v_grad_v_grad_w
338 
339 
340  !=====================================================
341  SUBROUTINE qs_ls_mass_adv (mesh, alpha, gg, ia, ja, a0)
342  !=====================================================
343 
344  ! <<alpha w + (g.D)w, alpha_ +(g.D)_ >> ===> A0 incremental accumulation
345 
346  USE gauss_points
347 
348  IMPLICIT NONE
349  !--------------------------------------------------------------------------
350  ! Declaration des variables globales
351  REAL(KIND=8), INTENT(IN) :: alpha
352  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
353  INTEGER, DIMENSION(:), INTENT(IN) :: ia
354  INTEGER, DIMENSION(:), INTENT(IN) :: ja
355  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
356 
357  TYPE(mesh_type), TARGET :: mesh
358  INTEGER, DIMENSION(:,:), POINTER :: jj
359  INTEGER, POINTER :: me
360 
361  !--------------------------------------------------------------------------
362  ! Declaration des variables locales
363  INTEGER :: l, k, ni, nj, p, m, i, j
364  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: gl
365  REAL(KIND=8) :: x
366  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: y
367 
368  !--------------------------------------------------------------------------
369 
370  CALL gauss(mesh)
371  jj => mesh%jj
372  me => mesh%me
373 
374  boucle_mm : DO m = 1, me
375 
376  boucle_l : DO l = 1, l_g
377 
378  gl = 0
379  boucle_k : DO k = 1, k_d
380  boucle_ni : DO ni =1 ,n_w
381  gl(k) = gl(k) + gg(k, jj(ni,m)) * ww(ni,l)
382  END DO boucle_ni
383  ENDDO boucle_k
384 
385  y = alpha*ww(:,l)
386  boucle_ni_2 : DO ni = 1, n_w
387  boucle_k_2 : DO k = 1, k_d
388  y(ni) = y(ni) + gl(k) * dw(k,ni,l,m)
389  END DO boucle_k_2
390  END DO boucle_ni_2
391 
392  boucle_ni_3 : DO ni = 1, n_w
393  i = jj(ni, m)
394  boucle_nj : DO nj = 1, n_w
395  j = jj(nj, m)
396  x = rj(l,m) * y(nj) * y(ni)
397  boucle_p : DO p = ia(i), ia(i+1) - 1
398  IF (ja(p) == j) THEN
399  a0(p) = a0(p) + x
400  EXIT
401  ENDIF
402  ENDDO boucle_p
403  ENDDO boucle_nj
404  ENDDO boucle_ni_3
405 
406  ENDDO boucle_l
407 
408  ENDDO boucle_mm
409 
410  END SUBROUTINE qs_ls_mass_adv
411 
412 
413  !=====================================================
414  SUBROUTINE qs_gals_mass_adv_m (mesh, param, alpha, gg, ia, ja, a0)
415  !=====================================================
416 
417  ! <<w + h(alpha w + (g.D)w), alpha_ +(g.D)_ >> ===> A0 incremental accumulation
418 
419  USE gauss_points
420 
421  IMPLICIT NONE
422  !--------------------------------------------------------------------------
423  ! Declaration des variables globales
424  REAL(KIND=8), INTENT(IN) :: alpha, param
425  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
426  INTEGER, DIMENSION(:), INTENT(IN) :: ia
427  INTEGER, DIMENSION(:), INTENT(IN) :: ja
428  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
429 
430  TYPE(mesh_type), TARGET :: mesh
431  INTEGER, DIMENSION(:,:), POINTER :: jj
432  INTEGER, POINTER :: me
433 
434  !--------------------------------------------------------------------------
435  ! Declaration des variables locales
436  INTEGER :: l, k, ni, nj, p, m, i, j
437  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: gl
438  REAL(KIND=8) :: x, mest, hloc
439  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: y
440 
441  !--------------------------------------------------------------------------
442 
443  CALL gauss(mesh)
444  jj => mesh%jj
445  me => mesh%me
446 
447  DO m = 1, me
448 
449  mest = 0
450  DO l = 1, l_g
451  mest = mest + rj(l,m)
452  END DO
453  hloc = param*mest**(1.d0/k_d)
454 
455  DO l = 1, l_g
456 
457  gl = 0
458  DO k = 1, k_d
459  DO ni =1 ,n_w
460  gl(k) = gl(k) + gg(k, jj(ni,m)) * ww(ni,l)
461  END DO
462  ENDDO
463 
464  y = alpha*ww(:,l)
465  DO ni = 1, n_w
466  DO k = 1, k_d
467  y(ni) = y(ni) + gl(k) * dw(k,ni,l,m)
468  END DO
469  END DO
470 
471  DO ni = 1, n_w
472  i = jj(ni, m)
473  DO nj = 1, n_w
474  j = jj(nj, m)
475  x = rj(l,m) * y(nj) * (ww(ni,l) + hloc*y(ni))
476  DO p = ia(i), ia(i+1) - 1
477  IF (ja(p) == j) THEN
478  a0(p) = a0(p) + x
479  EXIT
480  ENDIF
481  ENDDO
482  ENDDO
483  ENDDO
484 
485  ENDDO
486 
487  ENDDO
488 
489  END SUBROUTINE qs_gals_mass_adv_m
490 
491  !=====================================================
492  SUBROUTINE qs_h_v_grad_v_grad_w (mesh, gg, ia, ja, a0)
493  !=====================================================
494 
495  ! << h (g.D)w, (g.D)_ >> ===> A0 incremental accumulation
496 
497  USE gauss_points
498 
499  IMPLICIT NONE
500  !--------------------------------------------------------------------------
501  ! Declaration des variables globales
502  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
503  INTEGER, DIMENSION(:), INTENT(IN) :: ia
504  INTEGER, DIMENSION(:), INTENT(IN) :: ja
505  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
506 
507  TYPE(mesh_type), TARGET :: mesh
508  INTEGER, DIMENSION(:,:), POINTER :: jj
509  INTEGER, POINTER :: me
510 
511  !--------------------------------------------------------------------------
512  ! Declaration des variables locales
513  INTEGER :: l, k, ni, nj, p, m, i, j
514  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: gl
515  REAL(KIND=8) :: x, mest, hloc
516  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: y
517 
518  !--------------------------------------------------------------------------
519 
520  CALL gauss(mesh)
521  jj => mesh%jj
522  me => mesh%me
523 
524  boucle_mm : DO m = 1, me
525 
526  mest = 0
527  DO l = 1, l_g
528  mest = mest + rj(l,m)
529  END DO
530  hloc = mest**(1.d0/k_d)
531 
532  boucle_l : DO l = 1, l_g
533 
534  gl = 0
535  boucle_k : DO k = 1, k_d
536  boucle_ni : DO ni =1 ,n_w
537  gl(k) = gl(k) + gg(k, jj(ni,m)) * ww(ni,l)
538  END DO boucle_ni
539  ENDDO boucle_k
540 
541  y = 0
542  boucle_ni_2 : DO ni = 1, n_w
543  boucle_k_2 : DO k = 1, k_d
544  y(ni) = y(ni) + gl(k) * dw(k,ni,l,m)
545  END DO boucle_k_2
546  END DO boucle_ni_2
547 
548  boucle_ni_3 : DO ni = 1, n_w
549  i = jj(ni, m)
550  boucle_nj : DO nj = 1, n_w
551  j = jj(nj, m)
552  x = rj(l,m) * y(nj) * y(ni)
553  boucle_p : DO p = ia(i), ia(i+1) - 1
554  IF (ja(p) == j) THEN
555  a0(p) = a0(p) + x*hloc
556  EXIT
557  ENDIF
558  ENDDO boucle_p
559  ENDDO boucle_nj
560  ENDDO boucle_ni_3
561 
562  ENDDO boucle_l
563 
564  ENDDO boucle_mm
565 
566  END SUBROUTINE qs_h_v_grad_v_grad_w
567 
568 
569  SUBROUTINE qs_dif_mass_adv_m(mesh, visco, alpha, gg, ia, ja, a0)
570  !==========================================================
571 
572  ! << visco (Dw), D_) >>
573  ! + alpha * < w, _ >
574  ! + < w, (g.D)_ > + < w, (D.g) _ > / 2
575  ! ===> a0
576 
577  USE gauss_points
578 
579  REAL(KIND=8), INTENT(IN) :: visco, alpha
580  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
581  INTEGER, DIMENSION(:), INTENT(IN) :: ia
582  INTEGER, DIMENSION(:), INTENT(IN) :: ja
583  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
584 
585  TYPE(mesh_type), TARGET :: mesh
586  INTEGER, DIMENSION(:,:), POINTER :: jj
587  INTEGER, POINTER :: me
588 
589  INTEGER :: k, l, m, ni, nj, i, j, p
590  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: gl
591  REAL(KIND=8) :: dg, xij, masslm, viscolm
592  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
593  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: aij
594  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: y
595 
596  CALL gauss(mesh)
597  jj => mesh%jj
598  me => mesh%me
599 
600  DO l = 1, l_g
601  DO ni = 1, n_w
602  DO nj = 1, n_w
603  wwprod(ni,nj,l) = ww(ni,l)*ww(nj,l)
604  END DO
605  END DO
606  END DO
607 
608  DO m = 1, me
609  aij = 0
610  DO l = 1, l_g
611 
612  dg = 0.
613  gl = 0.
614  DO k = 1, k_d
615  DO ni =1 ,n_w
616  gl(k) = gl(k) + gg(k, jj(ni,m)) * ww(ni,l)
617  dg = dg + gg(k, jj(ni,m)) * dw(k,ni,l,m)
618  END DO
619  ENDDO
620  viscolm = visco*rj(l,m)
621  masslm = (alpha+0.5*dg)*rj(l,m)
622 
623  y = 0.
624  DO ni = 1, n_w;
625  DO k = 1, k_d
626  y(ni) = y(ni) + gl(k) * dw(k,ni,l,m)
627  END DO
628  END DO
629  y = rj(l,m)*y
630 
631  DO ni = 1, n_w; i = jj(ni, m)
632  DO nj = 1, n_w; j = jj(nj, m)
633 
634  xij = 0.
635  DO k = 1, k_d
636  xij = xij + dw(k,nj,l,m) * dw(k,ni,l,m)
637  END DO
638 
639  aij(ni,nj) = aij(ni,nj) + viscolm*xij + masslm*wwprod(ni,nj,l) + &
640  y(nj)*ww(ni,l)
641 
642  ENDDO
643  ENDDO
644 
645  ENDDO
646 
647  DO ni = 1, n_w; i = jj(ni, m)
648  DO nj = 1, n_w; j = jj(nj, m)
649  DO p = ia(i), ia(i+1) - 1
650  IF (ja(p) == j) then; a0(p) = a0(p) + aij(ni,nj); exit;
651  ENDIF
652  ENDDO
653  END DO
654  END DO
655 
656  ENDDO
657 
658  END SUBROUTINE qs_dif_mass_adv_m
659 
660  SUBROUTINE qs_diff_mass_adv_m(mesh, visco, alpha, gg, ia, ja, a0)
661  !==========================================================
662 
663  ! << visco (Dw), D_) >>
664  ! + alpha * < w, _ >
665  ! + < w, (g.D)_ >
666  ! ===> a0
667 
668  USE gauss_points
669 
670  IMPLICIT NONE
671 
672  REAL(KIND=8), INTENT(IN) :: visco, alpha
673  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
674  INTEGER, DIMENSION(:), INTENT(IN) :: ia
675  INTEGER, DIMENSION(:), INTENT(IN) :: ja
676  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
677 
678  TYPE(mesh_type), TARGET :: mesh
679  INTEGER, DIMENSION(:,:), POINTER :: jj
680  INTEGER, POINTER :: me
681 
682  INTEGER :: k, l, m, ni, nj, i, j, p
683  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: gl
684  REAL(KIND=8) :: xij, masslm, viscolm
685  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
686  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: aij
687  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: y
688 
689  CALL gauss(mesh)
690 
691  jj => mesh%jj
692  me => mesh%me
693 
694  DO l = 1, l_g
695  DO ni = 1, n_w
696  DO nj = 1, n_w
697  wwprod(ni,nj,l) = ww(ni,l)*ww(nj,l)
698  END DO
699  END DO
700  END DO
701 
702  DO m = 1, me
703  aij = 0
704  DO l = 1, l_g
705 
706  gl = 0.
707  DO k = 1, k_d
708  DO ni =1 ,n_w
709  gl(k) = gl(k) + gg(k, jj(ni,m)) * ww(ni,l)
710  END DO
711  ENDDO
712  viscolm = visco*rj(l,m)
713  masslm = alpha*rj(l,m)
714 
715  y = 0.
716  DO ni = 1, n_w;
717  DO k = 1, k_d
718  y(ni) = y(ni) + gl(k) * dw(k,ni,l,m)
719  END DO
720  END DO
721  y = rj(l,m)*y
722 
723  DO ni = 1, n_w; i = jj(ni, m)
724  DO nj = 1, n_w; j = jj(nj, m)
725 
726  xij = 0.
727  DO k = 1, k_d
728  xij = xij + dw(k,nj,l,m) * dw(k,ni,l,m)
729  END DO
730 
731  aij(ni,nj) = aij(ni,nj) + viscolm*xij + masslm*wwprod(ni,nj,l) + &
732  y(nj)*ww(ni,l)
733 
734  ENDDO
735  ENDDO
736 
737  ENDDO
738 
739  DO ni = 1, n_w; i = jj(ni, m)
740  DO nj = 1, n_w; j = jj(nj, m)
741  DO p = ia(i), ia(i+1) - 1
742  IF (ja(p) == j) then; a0(p) = a0(p) + aij(ni,nj); exit;
743  ENDIF
744  ENDDO
745  END DO
746  END DO
747 
748  ENDDO
749 
750  END SUBROUTINE qs_diff_mass_adv_m
751 
752  SUBROUTINE qs_dif_mass_adv_skew_m(mesh, visco, alpha, gg, ia, ja, a0)
753  !==========================================================
754 
755  ! << visco (Dw), D_) >>
756  ! + alpha * < w, _ >
757  ! + < w, (g.D)_ >/2 - < _, (g.D) w > / 2
758  ! ===> a0
759 
760  USE gauss_points
761 
762  IMPLICIT NONE
763  REAL(KIND=8), INTENT(IN) :: visco, alpha
764  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
765  INTEGER, DIMENSION(:), INTENT(IN) :: ia
766  INTEGER, DIMENSION(:), INTENT(IN) :: ja
767  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
768 
769  TYPE(mesh_type), TARGET :: mesh
770  INTEGER, DIMENSION(:,:), POINTER :: jj
771  INTEGER, POINTER :: me
772 
773  INTEGER :: k, l, m, ni, nj, i, j, p
774  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: gl
775  REAL(KIND=8) :: x, xij, masslm, viscolm
776  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
777  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: y
778 
779  CALL gauss(mesh)
780  jj => mesh%jj
781  me => mesh%me
782 
783  DO l = 1, l_g
784  DO ni = 1, n_w
785  DO nj = 1, n_w
786  wwprod(ni,nj,l) = ww(ni,l)*ww(nj,l)
787  END DO
788  END DO
789  END DO
790 
791  DO m = 1, me
792  DO l = 1, l_g
793 
794  viscolm = visco*rj(l,m)
795  masslm = alpha*rj(l,m)
796 
797  gl = 0.
798  DO k = 1, k_d
799  DO ni =1 ,n_w
800  gl(k) = gl(k) + gg(k, jj(ni,m)) * ww(ni,l)
801  END DO
802  ENDDO
803 
804  y = 0.
805  DO ni = 1, n_w;
806  DO k = 1, k_d
807  y(ni) = y(ni) + gl(k) * dw(k,ni,l,m)
808  END DO
809  END DO
810  y = 0.5*rj(l,m)*y
811 
812  DO ni = 1, n_w; i = jj(ni, m)
813  DO nj = 1, n_w; j = jj(nj, m)
814 
815  xij = 0.
816  DO k = 1, k_d
817  xij = xij + dw(k,nj,l,m) * dw(k,ni,l,m)
818  END DO
819 
820  x = viscolm*xij + masslm*wwprod(ni,nj,l) + &
821  y(nj)*ww(ni,l) - ww(nj,l)*y(ni)
822 
823  DO p = ia(i), ia(i+1) - 1
824  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
825  ENDIF
826  ENDDO
827 
828  ENDDO
829  ENDDO
830 
831  ENDDO
832  ENDDO
833 
834  END SUBROUTINE qs_dif_mass_adv_skew_m
835 
836  SUBROUTINE qs_dif_mass_m(mesh, visco, alpha, ia, ja, a0)
837  !==========================================================
838 
839  ! << visco (Dw), D_) >>
840  ! + alpha * < w, _ >
841  ! ===> a0
842 
843  USE gauss_points
844 
845  IMPLICIT NONE
846 
847  REAL(KIND=8), INTENT(IN) :: visco, alpha
848  INTEGER, DIMENSION(:), INTENT(IN) :: ia
849  INTEGER, DIMENSION(:), INTENT(IN) :: ja
850  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
851 
852  TYPE(mesh_type), TARGET :: mesh
853  INTEGER, DIMENSION(:,:), POINTER :: jj
854  INTEGER, POINTER :: me
855 
856  INTEGER :: k, l, m, ni, nj, i, j, p
857  REAL(KIND=8) :: x, xij, masslm, viscolm
858  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
859 
860  CALL gauss(mesh)
861  jj => mesh%jj
862  me => mesh%me
863 
864  DO l = 1, l_g
865  DO ni = 1, n_w
866  DO nj = 1, n_w
867  wwprod(ni,nj,l) = ww(ni,l)*ww(nj,l)
868  END DO
869  END DO
870  END DO
871 
872  DO m = 1, me
873  DO l = 1, l_g
874 
875  viscolm = visco*rj(l,m)
876  masslm = alpha*rj(l,m)
877 
878  DO ni = 1, n_w; i = jj(ni, m)
879  DO nj = 1, n_w; j = jj(nj, m)
880 
881  xij = 0.
882  DO k = 1, k_d
883  xij = xij + dw(k,nj,l,m) * dw(k,ni,l,m)
884  END DO
885 
886  x = viscolm*xij + masslm*wwprod(ni,nj,l)
887 
888  DO p = ia(i), ia(i+1) - 1
889  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
890  ENDIF
891  ENDDO
892 
893  ENDDO
894  ENDDO
895 
896  ENDDO
897  ENDDO
898 
899  END SUBROUTINE qs_dif_mass_m
900 
901  SUBROUTINE qs_adv_m_bis(mesh, gg, ia, ja, a0)
902  !==========================================================
903 
904  ! + < w, (g.D)_ > + < w, (D.g) _ > / 2
905  ! ===> a0
906 
907  USE gauss_points
908 
909  IMPLICIT NONE
910  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
911  INTEGER, DIMENSION(:), INTENT(IN) :: ia
912  INTEGER, DIMENSION(:), INTENT(IN) :: ja
913  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
914 
915  TYPE(mesh_type), TARGET :: mesh
916  INTEGER, DIMENSION(:,:), POINTER :: jj
917  INTEGER, POINTER :: me
918 
919  INTEGER :: k, l, m, ni, nj, i, j, p
920  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: gl
921  REAL(KIND=8) :: dg, x, masslm
922  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
923  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: y
924 
925  CALL gauss(mesh)
926  jj => mesh%jj
927  me => mesh%me
928 
929  DO l = 1, l_g
930  DO ni = 1, n_w
931  DO nj = 1, n_w
932  wwprod(ni,nj,l) = ww(ni,l)*ww(nj,l)
933  END DO
934  END DO
935  END DO
936 
937  DO m = 1, me
938  DO l = 1, l_g
939 
940  dg = 0.
941  gl = 0.
942  DO k = 1, k_d
943  DO ni =1 ,n_w
944  gl(k) = gl(k) + gg(k, jj(ni,m)) * ww(ni,l)
945  dg = dg + gg(k, jj(ni,m)) * dw(k,ni,l,m)
946  END DO
947  ENDDO
948  masslm = 0.5*dg*rj(l,m)
949 
950  y = 0.
951  DO ni = 1, n_w;
952  DO k = 1, k_d
953  y(ni) = y(ni) + gl(k) * dw(k,ni,l,m)
954  END DO
955  END DO
956  y = rj(l,m)*y
957 
958  DO ni = 1, n_w; i = jj(ni, m)
959  DO nj = 1, n_w; j = jj(nj, m)
960 
961  x = masslm*wwprod(ni,nj,l) + y(nj)*ww(ni,l)
962 
963  DO p = ia(i), ia(i+1) - 1
964  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
965  ENDIF
966  ENDDO
967 
968 
969  ENDDO
970  ENDDO
971 
972  ENDDO
973  ENDDO
974  END SUBROUTINE qs_adv_m_bis
975 
976  SUBROUTINE qs_adv_m(mesh, gg, ia, ja, a0)
977  !==========================================================
978 
979  ! + < w, (g.D)_ > + < w, (D.g) _ > / 2
980  ! ===> a0
981 
982  USE gauss_points
983 
984  IMPLICIT NONE
985  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
986  INTEGER, DIMENSION(:), INTENT(IN) :: ia
987  INTEGER, DIMENSION(:), INTENT(IN) :: ja
988  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
989 
990  TYPE(mesh_type), TARGET :: mesh
991  INTEGER, DIMENSION(:,:), POINTER :: jj
992  INTEGER, POINTER :: me
993 
994  INTEGER :: k, l, m, ni, nj, i, j, p
995  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: gl
996  REAL(KIND=8) :: dg, x, masslm, rjlm
997  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
998  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
999  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: y
1000  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
1001 
1002  CALL gauss(mesh)
1003  jj => mesh%jj
1004  me => mesh%me
1005 
1006  DO l = 1, l_g
1007  DO ni = 1, n_w
1008  DO nj = 1, n_w
1009  wwprod(ni,nj,l) = 0.5d0*ww(ni,l)*ww(nj,l)
1010  END DO
1011  END DO
1012  END DO
1013 
1014  DO m = 1, me
1015  DO ni = 1, n_w
1016  j_loc(ni) = jj(ni,m)
1017  END DO
1018 
1019  DO l = 1, l_g
1020 
1021  rjlm = rj(l,m)
1022  dg = 0.
1023  DO k = 1, k_d
1024  gl(k) = 0.
1025  DO ni =1 ,n_w
1026  dw_loc(k,ni) = dw(k,ni,l,m)
1027  gl(k) = gl(k) + gg(k, j_loc(ni)) * ww(ni,l)
1028  dg = dg + gg(k, j_loc(ni)) * dw_loc(k,ni)
1029  END DO
1030  gl(k) = gl(k)*rjlm
1031  ENDDO
1032  masslm = dg*rjlm
1033 
1034  DO ni = 1, n_w
1035  y(ni) = 0.
1036  DO k = 1, k_d
1037  y(ni) = y(ni) + gl(k) * dw_loc(k,ni)
1038  END DO
1039  END DO
1040 
1041  DO ni = 1, n_w; i = j_loc(ni)
1042  DO nj = 1, n_w; j = j_loc(nj)
1043 
1044  x = masslm*wwprod(ni,nj,l) + y(nj)*ww(ni,l)
1045 
1046  DO p = ia(i), ia(i+1) - 1
1047  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
1048  ENDIF
1049  ENDDO
1050 
1051 
1052  ENDDO
1053  ENDDO
1054 
1055  ENDDO
1056  ENDDO
1057  END SUBROUTINE qs_adv_m
1058 
1059  SUBROUTINE qs_001_m(mesh, gg, ia, ja, a0)
1060  !==========================================================
1061 
1062  ! < w, (g.D)_ > ===> a0
1063 
1064  USE gauss_points
1065 
1066  IMPLICIT NONE
1067  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
1068  INTEGER, DIMENSION(:), INTENT(IN) :: ia
1069  INTEGER, DIMENSION(:), INTENT(IN) :: ja
1070  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
1071 
1072  TYPE(mesh_type), TARGET :: mesh
1073  INTEGER, DIMENSION(:,:), POINTER :: jj
1074  INTEGER, POINTER :: me
1075 
1076  INTEGER :: k, l, m, ni, nj, i, j, p
1077  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: gl
1078  REAL(KIND=8) :: x
1079  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: y
1080 
1081  CALL gauss(mesh)
1082  jj => mesh%jj
1083  me => mesh%me
1084 
1085  DO m = 1, me
1086 
1087  DO l = 1, l_g
1088 
1089  gl = 0.
1090  DO k = 1, k_d
1091  DO ni =1 ,n_w
1092  gl(k) = gl(k) + gg(k, jj(ni,m)) * ww(ni,l)
1093  END DO
1094  ENDDO
1095 
1096  y = 0.
1097  DO ni = 1, n_w;
1098  DO k = 1, k_d
1099  y(ni) = y(ni) + gl(k) * dw(k,ni,l,m)
1100  END DO
1101  END DO
1102  y = rj(l,m)*y
1103 
1104  DO ni = 1, n_w; i = jj(ni, m)
1105  DO nj = 1, n_w; j = jj(nj, m)
1106 
1107  x = y(nj)*ww(ni,l)
1108 
1109  DO p = ia(i), ia(i+1) - 1
1110  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
1111  ENDIF
1112  ENDDO
1113 
1114 
1115  ENDDO
1116  ENDDO
1117 
1118  ENDDO
1119  ENDDO
1120  END SUBROUTINE qs_001_m
1121 
1122  SUBROUTINE qs_h_100_m(mesh, alpha, gg, ia, ja, a0)
1123  !==========================================================
1124 
1125  ! < H_loc alpha _, (g.D)w > ===> a0
1126 
1127  USE gauss_points
1128 
1129  IMPLICIT NONE
1130  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
1131  REAL(KIND=8), INTENT(IN) :: alpha
1132  INTEGER, DIMENSION(:), INTENT(IN) :: ia
1133  INTEGER, DIMENSION(:), INTENT(IN) :: ja
1134  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
1135 
1136  TYPE(mesh_type), TARGET :: mesh
1137  INTEGER, DIMENSION(:,:), POINTER :: jj
1138  INTEGER, POINTER :: me
1139 
1140  INTEGER :: k, l, m, ni, nj, i, j, p
1141  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: gl
1142  REAL(KIND=8) :: x, mest, hloc
1143  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: y
1144 
1145  CALL gauss(mesh)
1146  jj => mesh%jj
1147  me => mesh%me
1148 
1149  DO m = 1, me
1150 
1151  mest = 0
1152  DO l = 1, l_g
1153  mest = mest + rj(l,m)
1154  END DO
1155  hloc = alpha*mest**(1.d0/k_d)
1156 
1157  DO l = 1, l_g
1158 
1159  gl = 0
1160  DO k = 1, k_d
1161  DO ni =1 ,n_w
1162  gl(k) = gl(k) + gg(k, jj(ni,m)) * ww(ni,l)
1163  END DO
1164  ENDDO
1165 
1166  y = 0
1167  DO ni = 1, n_w;
1168  DO k = 1, k_d
1169  y(ni) = y(ni) + gl(k) * dw(k,ni,l,m)
1170  END DO
1171  END DO
1172  y = rj(l,m)*y
1173 
1174  DO ni = 1, n_w; i = jj(ni, m)
1175  DO nj = 1, n_w; j = jj(nj, m)
1176 
1177  x = y(ni)*ww(nj,l)
1178 
1179  DO p = ia(i), ia(i+1) - 1
1180  IF (ja(p) == j) then; a0(p) = a0(p) + hloc*x; exit;
1181  ENDIF
1182  ENDDO
1183 
1184 
1185  ENDDO
1186  ENDDO
1187 
1188  ENDDO
1189  ENDDO
1190  END SUBROUTINE qs_h_100_m
1191 
1192  SUBROUTINE qs_1_sgs_1_m (mesh, ff, ia, ja, a0)
1193  !=======================================
1194 
1195  ! << (Dw), h* f (D_) >> ===> A0 incremental accumulation
1196 
1197  USE gauss_points
1198 
1199  IMPLICIT NONE
1200  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
1201  INTEGER, DIMENSION(:), INTENT(IN) :: ia
1202  INTEGER, DIMENSION(:), INTENT(IN) :: ja
1203  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
1204 
1205  INTEGER :: m, l, ni, nj, i, j, k, p
1206  REAL(KIND=8) :: fl, x, xij, h, exp
1207 
1208  TYPE(mesh_type), TARGET :: mesh
1209  INTEGER, DIMENSION(:,:), POINTER :: jj
1210  INTEGER, POINTER :: me
1211 
1212  CALL gauss(mesh)
1213  jj => mesh%jj
1214  me => mesh%me
1215 
1216  m = SIZE(jj,1)
1217  SELECT CASE(m)
1218  CASE(3); exp = 1.d0/2
1219  CASE(6); exp = 1.d0/2
1220  CASE(4); exp = 1.d0/3
1221  CASE(10); exp = 1.d0/3
1222  END SELECT
1223 
1224  DO m = 1, me
1225 
1226  h = 0
1227  DO l = 1, l_g
1228  h = h + rj(l,m)
1229  END DO
1230  h = h**exp
1231 
1232  DO l = 1, l_g
1233 
1234  fl = 0.
1235  DO ni =1 ,n_w
1236  fl = fl + ff(jj(ni,m)) * ww(ni,l)
1237  END DO
1238  fl = h*fl*rj(l,m)
1239 
1240  DO ni = 1, n_w; i = jj(ni, m)
1241  DO nj = 1, n_w; j = jj(nj, m)
1242 
1243  xij = 0.
1244  DO k = 1, k_d
1245  xij = xij + dw(k,nj,l,m) * dw(k,ni,l,m)
1246  END DO
1247 
1248  x = fl*xij
1249  DO p = ia(i), ia(i+1) - 1
1250  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
1251  ENDIF
1252  ENDDO
1253 
1254  ENDDO
1255  ENDDO
1256 
1257  ENDDO
1258  ENDDO
1259 
1260  END SUBROUTINE qs_1_sgs_1_m
1261 
1262  SUBROUTINE qs_101_m (mesh, ff, ia, ja, a0)
1263  !=======================================
1264 
1265  ! << (Dw), f (D_) >> ===> A0 incremental accumulation
1266 
1267  USE gauss_points
1268 
1269  IMPLICIT NONE
1270  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
1271  INTEGER, DIMENSION(:), INTENT(IN) :: ia
1272  INTEGER, DIMENSION(:), INTENT(IN) :: ja
1273  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
1274 
1275  INTEGER :: m, l, ni, nj, i, j, k, p
1276  REAL(KIND=8) :: fl, x, xij
1277 
1278  TYPE(mesh_type), TARGET :: mesh
1279  INTEGER, DIMENSION(:,:), POINTER :: jj
1280  INTEGER, POINTER :: me
1281 
1282  CALL gauss(mesh)
1283  jj => mesh%jj
1284  me => mesh%me
1285 
1286  DO m = 1, me
1287 
1288  DO l = 1, l_g
1289 
1290  fl = 0.
1291  DO ni =1 ,n_w
1292  fl = fl + ff(jj(ni,m)) * ww(ni,l)
1293  END DO
1294  fl = fl*rj(l,m)
1295 
1296  DO ni = 1, n_w; i = jj(ni, m)
1297  DO nj = 1, n_w; j = jj(nj, m)
1298 
1299  xij = 0.
1300  DO k = 1, k_d
1301  xij = xij + dw(k,nj,l,m) * dw(k,ni,l,m)
1302  END DO
1303 
1304  x = fl*xij
1305  DO p = ia(i), ia(i+1) - 1
1306  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
1307  ENDIF
1308  ENDDO
1309 
1310  ENDDO
1311  ENDDO
1312 
1313  ENDDO
1314  ENDDO
1315 
1316  END SUBROUTINE qs_101_m
1317 
1318  SUBROUTINE qs_v_plap_m (mesh, p, vstar, ff, ia, ja, a0)
1319  !=======================================
1320 
1321  ! << (Dw), |vstar.D(ff)|^(p-2) (D_) >> ===> A0 incremental accumulation
1322 
1323  USE gauss_points
1324 
1325  IMPLICIT NONE
1326  REAL(KIND=8), INTENT(IN) :: p
1327  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: vstar
1328  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
1329  INTEGER, DIMENSION(:), INTENT(IN) :: ia
1330  INTEGER, DIMENSION(:), INTENT(IN) :: ja
1331  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
1332 
1333  TYPE(mesh_type), TARGET :: mesh
1334  INTEGER, DIMENSION(:,:), POINTER :: jj
1335  INTEGER, POINTER :: me
1336 
1337  REAL(KIND=8) :: x, xrjlm, xij
1338  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: vl, gl
1339  INTEGER :: ni, nj, i, j, k, l, m, q
1340 
1341  CALL gauss(mesh)
1342  jj => mesh%jj
1343  me => mesh%me
1344 
1345  DO m = 1, me
1346 
1347  DO l = 1, l_g
1348 
1349  vl = 0.
1350  gl = 0.
1351  DO k= 1, k_d
1352  DO ni =1 ,n_w
1353  vl(k) = vl(k) + vstar(k,jj(ni,m)) * ww(ni,l)
1354  gl(k) = gl(k) + ff(jj(ni,m)) * dw(k,ni,l,m)
1355  END DO
1356  END DO
1357 
1358  x = 0.
1359  DO k= 1, k_d
1360  x = x + vl(k)*gl(k)
1361  END DO
1362  IF (abs(x) .LT. 1.d-12) THEN
1363  xrjlm = 0.
1364  ELSE
1365  xrjlm = (abs(x)**(p-2.d0))* rj(l,m)
1366  ENDIF
1367 
1368  DO ni = 1, n_w; i = jj(ni, m)
1369  DO nj = 1, n_w; j = jj(nj, m)
1370 
1371  xij = 0.
1372  DO k = 1, k_d
1373  xij = xij + dw(k,nj,l,m) * dw(k,ni,l,m)
1374  END DO
1375 
1376  x = xrjlm * xij
1377  DO q = ia(i), ia(i+1) - 1
1378  IF (ja(q) == j) then; a0(q) = a0(q) + x; exit;
1379  ENDIF
1380  ENDDO
1381 
1382  ENDDO
1383  ENDDO
1384 
1385  ENDDO
1386  ENDDO
1387 
1388  END SUBROUTINE qs_v_plap_m
1389 
1390  SUBROUTINE qs_plap_m (mesh, p, ff, ia, ja, a0)
1391  !=======================================
1392 
1393  ! << (Dw), |D(ff)|^(p-2) (D_) >> ===> A0 incremental accumulation
1394 
1395  USE gauss_points
1396 
1397  IMPLICIT NONE
1398  REAL(KIND=8), INTENT(IN) :: p
1399  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
1400  INTEGER, DIMENSION(:), INTENT(IN) :: ia
1401  INTEGER, DIMENSION(:), INTENT(IN) :: ja
1402  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
1403 
1404  TYPE(mesh_type), TARGET :: mesh
1405  INTEGER, DIMENSION(:,:), POINTER :: jj
1406  INTEGER, POINTER :: me
1407 
1408  REAL(KIND=8) :: x, xrjlm, xij
1409  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: gl
1410  INTEGER :: ni, nj, i, j, k, l, m, q
1411 
1412  CALL gauss(mesh)
1413  jj => mesh%jj
1414  me => mesh%me
1415 
1416  DO m = 1, me
1417 
1418  DO l = 1, l_g
1419 
1420  gl = 0.
1421  DO k= 1, k_d
1422  DO ni =1 ,n_w
1423  gl(k) = gl(k) + ff(jj(ni,m)) * dw(k,ni,l,m)
1424  END DO
1425  END DO
1426 
1427  x = 0.
1428  DO k= 1, k_d
1429  x = x + gl(k)*gl(k)
1430  END DO
1431  IF (abs(x) .LT. 1.d-20) THEN
1432  xrjlm = 0.
1433  ELSE
1434  xrjlm = (sqrt(x)**(p-2.d0))* rj(l,m)
1435  ENDIF
1436 
1437  DO ni = 1, n_w; i = jj(ni, m)
1438  DO nj = 1, n_w; j = jj(nj, m)
1439 
1440  xij = 0.
1441  DO k = 1, k_d
1442  xij = xij + dw(k,nj,l,m) * dw(k,ni,l,m)
1443  END DO
1444 
1445  x = xrjlm * xij
1446  DO q = ia(i), ia(i+1) - 1
1447  IF (ja(q) == j) then; a0(q) = a0(q) + x; exit;
1448  ENDIF
1449  ENDDO
1450 
1451  ENDDO
1452  ENDDO
1453 
1454  ENDDO
1455  ENDDO
1456 
1457  END SUBROUTINE qs_plap_m
1458 
1459  SUBROUTINE qs_adv_stab_m(mesh, gg, stab, istab, ia, ja, a0)
1460  !==========================================================
1461 
1462  ! + < w, (g.D)_ > + < w, (D.g) _ > / 2
1463  !if (istab .eq. 1) + c*h^k << (Dw), |D(gg)|^p (D_) >>
1464  !if (istab .eq. 2) + c*h^k << (gg.Dw), |gg.D(gg)|^p (gg.D_) >>
1465  ! ===> a0
1466  !
1467 
1468  USE gauss_points
1469 
1470  IMPLICIT NONE
1471  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
1472  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: stab
1473  INTEGER, INTENT(IN) :: istab
1474  INTEGER, DIMENSION(:), INTENT(IN) :: ia
1475  INTEGER, DIMENSION(:), INTENT(IN) :: ja
1476  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
1477 
1478  TYPE(mesh_type), TARGET :: mesh
1479  INTEGER, DIMENSION(:,:), POINTER :: jj
1480  INTEGER, POINTER :: me
1481 
1482  INTEGER :: k, kp, l, m, ni, nj, i, j, p
1483  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: vl
1484  REAL(KIND=8) :: dg, x, xij, masslm, rjlm
1485  REAL(KIND=8) :: h,hmu,xrjlm, exponent
1486  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
1487  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
1488  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%k_d) :: gradv
1489  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: vdw
1490  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
1491 
1492  CALL gauss(mesh)
1493  jj => mesh%jj
1494  me => mesh%me
1495 
1496  exponent = 1./k_d
1497 
1498  DO l = 1, l_g
1499  DO ni = 1, n_w
1500  DO nj = 1, n_w
1501  wwprod(ni,nj,l) = 0.5*ww(ni,l)*ww(nj,l)
1502  END DO
1503  END DO
1504  END DO
1505 
1506  DO m = 1, me
1507 
1508  DO ni = 1, n_w
1509  j_loc(ni) = jj(ni,m)
1510  END DO
1511 
1512  h = 0.
1513  DO l = 1, l_g
1514  h = h + rj(l,m)
1515  END DO
1516  ! h = h**exponent
1517  hmu = stab(1)*h**(stab(2)*exponent)
1518 
1519  DO l = 1, l_g
1520 
1521  rjlm = rj(l,m)
1522  dg = 0.
1523  gradv = 0.
1524 
1525  dw_loc = dw(:,:,l,m)
1526 
1527  DO k = 1, k_d
1528  vl(k) = 0.
1529  DO ni =1 ,n_w
1530  vl(k) = vl(k) + gg(k, j_loc(ni)) * ww(ni,l)
1531  DO kp = 1, k_d
1532  gradv(k,kp) = gradv(k,kp) + gg(k, j_loc(ni)) * dw_loc(kp,ni)
1533  END DO
1534  END DO
1535  dg = dg + gradv(k,k)
1536  ENDDO
1537 
1538  masslm = dg*rjlm
1539 
1540  x = 0.
1541  IF (istab .EQ. 1) THEN
1542  DO k= 1, k_d
1543  DO kp = 1, k_d
1544  x = x + (gradv(k,kp) + gradv(kp,k))**2
1545  END DO
1546  END DO
1547  ELSE IF (istab .EQ. 2) THEN
1548  DO k= 1, k_d
1549  DO kp = 1, k_d
1550  x = x +(vl(kp)*gradv(k,kp))**2
1551  END DO
1552  END DO
1553  ELSE
1554  WRITE(*,*) 'qs_adv_stab_M: istab > 2'
1555  stop
1556  END IF
1557 
1558  IF (abs(x) .LT. 1.d-12) THEN
1559  xrjlm = 0.
1560  ELSE
1561  xrjlm = hmu*(sqrt(x)**(stab(3)))* rjlm
1562  ENDIF
1563 
1564  DO ni = 1, n_w
1565  vdw(ni) = 0.
1566  DO k = 1, k_d
1567  vdw(ni) = vdw(ni) + vl(k)*dw_loc(k,ni)
1568  END DO
1569  END DO
1570 
1571  DO ni = 1, n_w; i = j_loc(ni)
1572  DO nj = 1, n_w; j = j_loc(nj)
1573 
1574 
1575  xij = 0.
1576  IF (istab .EQ. 1) THEN
1577  DO k = 1, k_d
1578  xij = xij + dw_loc(k,nj) * dw_loc(k,ni)
1579  END DO
1580  ELSE IF (istab .EQ. 2) THEN
1581  xij = vdw(ni)*vdw(nj)
1582  END IF
1583 
1584  x = masslm*wwprod(ni,nj,l) + vdw(nj)*ww(ni,l)*rjlm + xrjlm * xij
1585  ! x = x + xrjlm * xij
1586 
1587  DO p = ia(i), ia(i+1) - 1
1588  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
1589  ENDIF
1590  ENDDO
1591 
1592 
1593  ENDDO
1594  ENDDO
1595 
1596  ENDDO
1597  ENDDO
1598 
1599  END SUBROUTINE qs_adv_stab_m
1600 
1601  SUBROUTINE bs_101_m (mesh, ff, ia, ja, a0)
1602  !==============================================
1603 
1604  ! - < J(w,_), f > ===> a0 incremental accumulation
1605  ! < (Dw) x k, ff (D_) > ===> a0 incremental accumulation
1606 
1607  USE gauss_points
1608 
1609  IMPLICIT NONE
1610  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
1611  INTEGER, DIMENSION(:), INTENT(IN) :: ia, ja
1612  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
1613 
1614  REAL(KIND=8) :: f, x, s
1615  INTEGER :: m, l, n, n1, i, j, p
1616 
1617  TYPE(mesh_type), TARGET :: mesh
1618  INTEGER, DIMENSION(:,:), POINTER :: jj
1619  INTEGER, POINTER :: me
1620 
1621  CALL gauss(mesh)
1622  jj => mesh%jj
1623  me => mesh%me
1624 
1625  DO m = 1, me
1626 
1627  DO l = 1, l_g
1628 
1629  f = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m)
1630 
1631  DO n = 1, n_w; i = jj(n, m)
1632  DO n1 = 1, n_w; j = jj(n1, m)
1633  s = dw(1,n,l,m) * dw(2,n1,l,m) - dw(2,n,l,m) * dw(1,n1,l,m)
1634  x = -s * f
1635 
1636  DO p = ia(i), ia(i+1) - 1
1637  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
1638  ENDIF
1639  END DO
1640 
1641  ENDDO
1642  ENDDO
1643 
1644  ENDDO
1645  ENDDO
1646 
1647 
1648  END SUBROUTINE bs_101_m
1649 
1650 
1651  !==============================================
1652 
1653  SUBROUTINE qs_00_s_m (mesh, fs, ia, ja, a0)
1654  !==========================================
1655 
1656  ! < ws, fs _ >_s ===> A0 incremental accumulation of boundary terms
1657 
1658  USE gauss_points
1659 
1660  IMPLICIT NONE
1661  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: fs
1662  INTEGER, DIMENSION(:), INTENT(IN) :: ia
1663  INTEGER, DIMENSION(:), INTENT(IN) :: ja
1664  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
1665 
1666  INTEGER :: ms, ls, ns, ns1, i, j, q
1667  REAL(KIND=8) :: fls, x
1668 
1669  TYPE(mesh_type), TARGET :: mesh
1670  INTEGER, DIMENSION(:,:), POINTER :: js, is
1671  INTEGER, POINTER :: mes
1672 
1673  CALL gauss(mesh)
1674  js => mesh%jjs
1675  is => mesh%iis
1676  mes => mesh%mes
1677 
1678  DO ms = 1, mes
1679 
1680  DO ls = 1, l_gs
1681  fls = 0
1682  DO ns = 1, n_ws
1683  fls = fls + fs(is(ns,ms)) * wws(ns,ls) * rjs(ls,ms)
1684  END DO
1685  fls = fls * rjs(ls,ms)
1686 
1687  DO ns = 1, n_ws; i = js(ns, ms)
1688  DO ns1 = 1, n_ws; j = js(ns1, ms)
1689  x = wws(ns,ls) * fls * wws(ns1,ls)
1690 
1691  DO q = ia(i), ia(i+1) - 1
1692  IF (ja(q) == j) then; a0(q) = a0(q) + x; exit;
1693  ENDIF
1694  ENDDO
1695 
1696  ENDDO
1697  ENDDO
1698  ENDDO
1699  ENDDO
1700 
1701  END SUBROUTINE qs_00_s_m
1702 
1703  SUBROUTINE cv_11_m (mesh, alpha, ia, ja, a0)
1704  !=======================================
1705 
1706  ! alpha << (D x w), (D x _) >> ===> a0 incremental acc. in 3D
1707 
1708  USE gauss_points
1709 
1710  IMPLICIT NONE
1711  REAL(KIND=8), INTENT(IN) :: alpha
1712  INTEGER, DIMENSION(:), INTENT(IN) :: ia
1713  INTEGER, DIMENSION(:), INTENT(IN) :: ja
1714  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
1715 
1716  INTEGER :: m, l, n, q, n0, k, k1, k2, h, h1, h2, i, j, i_b, j_b, &
1717  bloc_size
1718  REAL(KIND=8) :: fla, x
1719 
1720  TYPE(mesh_type), TARGET :: mesh
1721  INTEGER, DIMENSION(:,:), POINTER :: jj
1722  INTEGER, POINTER :: me
1723 
1724  CALL gauss(mesh)
1725  jj => mesh%jj
1726  me => mesh%me
1727 
1728  bloc_size = (SIZE(ia) - 1)/3
1729 
1730  DO m = 1, me
1731  DO l = 1, l_g
1732 
1733  fla = alpha * rj(l,m)
1734 
1735  DO n = 1, n_w; i_b = jj(n, m)
1736  DO k = 1, k_d; k1 = modulo(k,k_d) + 1; k2 = modulo(k+1,k_d) + 1
1737  i = i_b + (k-1)*bloc_size
1738 
1739  DO n0 = 1, n_w; j_b = jj(n0,m)
1740  DO h = 1, k_d; h1 = modulo(h,k_d) + 1; h2 = modulo(h+1,k_d) + 1
1741  j = j_b + (h-1)*bloc_size
1742 
1743  IF (h.EQ.k) THEN
1744  x = dw(k1,n,l,m)*dw(h1,n0,l,m) + &
1745  dw(k2,n,l,m)*dw(h2,n0,l,m)
1746  ELSE
1747  x = - dw(h,n,l,m)*dw(k,n0,l,m)
1748  END IF
1749 
1750  x = x * fla
1751 
1752  DO q = ia(i), ia(i+1) - 1
1753  IF (ja(q) == j) then; a0(q) = a0(q) + x; exit;
1754  ENDIF
1755  ENDDO
1756 
1757  ENDDO
1758  ENDDO
1759 
1760  ENDDO
1761  ENDDO
1762 
1763  ENDDO
1764  ENDDO
1765 
1766  END SUBROUTINE cv_11_m
1767 
1768  SUBROUTINE cc_101_m (mesh, gg, ia, ja, a0)
1769  !=======================================
1770 
1771  ! << (D x w), gg x (D x _) >> ===> a0 incremental acc. in 3D
1772 
1773  USE gauss_points
1774 
1775  IMPLICIT NONE
1776  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
1777  INTEGER, DIMENSION(:), INTENT(IN) :: ia
1778  INTEGER, DIMENSION(:), INTENT(IN) :: ja
1779  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
1780 
1781  INTEGER :: m, l, n, q, n0, k, k1, k2, h, h1, h2, i, j, i_b, j_b, &
1782  bloc_size, s , o
1783  REAL(KIND=8) :: x
1784  REAL(KIND=8), DIMENSION(3) :: gl
1785 
1786  TYPE(mesh_type), TARGET :: mesh
1787  INTEGER, DIMENSION(:,:), POINTER :: jj
1788  INTEGER, POINTER :: me
1789 
1790  CALL gauss(mesh)
1791  jj => mesh%jj
1792  me => mesh%me
1793 
1794  bloc_size = (SIZE(ia) - 1)/3
1795 
1796  DO m = 1, me
1797  DO l = 1, l_g
1798 
1799  DO k = 1, k_d
1800  gl(k) = sum(gg(k, jj(:,m)) * ww(:,l)) * rj(l,m)
1801  ENDDO
1802 
1803  DO n = 1, n_w; i_b = jj(n, m)
1804  DO k = 1, k_d; k1 = modulo(k,k_d) + 1; k2 = modulo(k+1,k_d) + 1
1805  i = i_b + (k-1)*bloc_size
1806 
1807  DO n0 = 1, n_w; j_b = jj(n0,m)
1808  DO h = 1, k_d; h1 = modulo(h,k_d) + 1; h2 = modulo(h+1,k_d) + 1
1809  j = j_b + (h-1)*bloc_size
1810 
1811  IF (h == k) THEN
1812  x = gl(h) * ( dw(k2,n,l,m) * dw(h1,n0,l,m) &
1813  - dw(k1,n,l,m) * dw(h2,n0,l,m) )
1814  ELSE
1815  o = 6 - k - h
1816  s = 0; IF (k > h) s = 1
1817  x = (-1)**(k+h+s) * &
1818  ( dw(o,n,l,m) * (gl(h1)*dw(h1,n0,l,m) + gl(h2)*dw(h2,n0,l,m)) &
1819  + dw(h,n,l,m) * gl(h) * dw(o,n0,l,m) )
1820  END IF
1821 
1822  x = x * rj(l,m)
1823 
1824  DO q = ia(i), ia(i+1) - 1
1825  IF (ja(q) == j) then; a0(q) = a0(q) + x; exit;
1826  ENDIF
1827  ENDDO
1828 
1829  ENDDO
1830  ENDDO
1831 
1832  ENDDO
1833  ENDDO
1834 
1835  ENDDO
1836  ENDDO
1837 
1838  END SUBROUTINE cc_101_m
1839 
1840 
1841  SUBROUTINE cv_00_s_m (mesh, fs, ia, ja, a0)
1842  !==========================================
1843 
1844  ! < ws x ns, fs (_ x ns) >_s ===> A0 incremental accumulation of boundary terms
1845 
1846  USE gauss_points
1847 
1848  IMPLICIT NONE
1849  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: fs
1850  INTEGER, DIMENSION(:), INTENT(IN) :: ia
1851  INTEGER, DIMENSION(:), INTENT(IN) :: ja
1852  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
1853 
1854  INTEGER :: ms, ls, ns, ns1, i, j, q, i_b, j_b, bloc_size, h, k
1855  REAL(KIND=8) :: fls, x, wwprod
1856  REAL(KIND=8), DIMENSION(3,3) :: tab
1857 
1858  TYPE(mesh_type), TARGET :: mesh
1859  INTEGER, DIMENSION(:,:), POINTER :: js, is
1860  INTEGER, POINTER :: mes
1861 
1862  CALL gauss(mesh)
1863  js => mesh%jjs
1864  is => mesh%iis
1865  mes => mesh%mes
1866 
1867  bloc_size = (SIZE(ia) - 1)/3
1868 
1869  DO ms = 1, mes
1870  DO ls = 1, l_gs
1871 
1872  fls = sum(fs(is(:,ms)) * wws(:,ls)) * rjs(ls,ms)
1873 
1874  tab(1,1) = rnorms(2,ls,ms)**2 + rnorms(3,ls,ms)**2
1875  tab(2,2) = rnorms(1,ls,ms)**2 + rnorms(3,ls,ms)**2
1876  tab(3,3) = rnorms(1,ls,ms)**2 + rnorms(2,ls,ms)**2
1877 
1878  tab(1,2) = -rnorms(2,ls,ms) * rnorms(1,ls,ms)
1879  tab(1,3) = -rnorms(3,ls,ms) * rnorms(1,ls,ms)
1880  tab(2,1) = -rnorms(1,ls,ms) * rnorms(2,ls,ms)
1881  tab(2,3) = -rnorms(3,ls,ms) * rnorms(2,ls,ms)
1882  tab(3,1) = -rnorms(1,ls,ms) * rnorms(3,ls,ms)
1883  tab(3,2) = -rnorms(2,ls,ms) * rnorms(3,ls,ms)
1884 
1885 
1886  DO ns = 1, n_ws; i_b = js(ns, ms)
1887  DO ns1 = 1, n_ws; j_b = js(ns1, ms)
1888 
1889  wwprod = wws(ns, ls) * wws(ns1, ls) * fls
1890 
1891  DO k = 1, k_d
1892  DO h = 1, k_d
1893 
1894  i = i_b + (k-1)*bloc_size; j = j_b + (h-1)*bloc_size
1895 
1896  x = wwprod * tab(k,h)
1897 
1898  DO q = ia(i), ia(i+1) - 1
1899  IF (ja(q) == j) then; a0(q) = a0(q) + x; exit;
1900  ENDIF
1901  ENDDO
1902 
1903  ENDDO
1904  ENDDO
1905  ENDDO
1906  ENDDO
1907  ENDDO
1908  ENDDO
1909 
1910  END SUBROUTINE cv_00_s_m
1911 
1912  SUBROUTINE cc_00_s_m (mesh, alpha, ia, ja, a0)
1913  !==========================================
1914 
1915  ! < ws x ns, (_ x ns) >_s ===> A0 incremental accumulation of boundary terms
1916 
1917  USE gauss_points
1918 
1919  IMPLICIT NONE
1920  REAL(KIND=8), INTENT(IN) :: alpha
1921  INTEGER, DIMENSION(:), INTENT(IN) :: ia
1922  INTEGER, DIMENSION(:), INTENT(IN) :: ja
1923  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
1924 
1925  INTEGER :: ms, ls, ns, ns1, i, j, q, i_b, j_b, bloc_size, h, k
1926  REAL(KIND=8) :: x, wwprod
1927  REAL(KIND=8), DIMENSION(3,3) :: tab
1928 
1929  TYPE(mesh_type), TARGET :: mesh
1930  INTEGER, DIMENSION(:,:), POINTER :: js
1931  INTEGER, POINTER :: mes
1932 
1933  CALL gauss(mesh)
1934  js => mesh%jjs
1935  mes => mesh%mes
1936 
1937  bloc_size = (SIZE(ia) - 1)/3
1938 
1939  DO ms = 1, mes
1940  DO ls = 1, l_gs
1941 
1942  tab(1,1) = rnorms(2,ls,ms)**2 + rnorms(3,ls,ms)**2
1943  tab(2,2) = rnorms(1,ls,ms)**2 + rnorms(3,ls,ms)**2
1944  tab(3,3) = rnorms(1,ls,ms)**2 + rnorms(2,ls,ms)**2
1945 
1946  tab(1,2) = -rnorms(2,ls,ms) * rnorms(1,ls,ms)
1947  tab(1,3) = -rnorms(3,ls,ms) * rnorms(1,ls,ms)
1948  tab(2,1) = -rnorms(1,ls,ms) * rnorms(2,ls,ms)
1949  tab(2,3) = -rnorms(3,ls,ms) * rnorms(2,ls,ms)
1950  tab(3,1) = -rnorms(1,ls,ms) * rnorms(3,ls,ms)
1951  tab(3,2) = -rnorms(2,ls,ms) * rnorms(3,ls,ms)
1952 
1953  tab = tab * rjs(ls,ms) * alpha
1954 
1955  DO ns = 1, n_ws; i_b = js(ns, ms)
1956  DO ns1 = 1, n_ws; j_b = js(ns1, ms)
1957 
1958  wwprod = wws(ns, ls) * wws(ns1, ls)
1959 
1960  DO k = 1, k_d
1961  DO h = 1, k_d
1962 
1963  i = i_b + (k-1)*bloc_size; j = j_b + (h-1)*bloc_size
1964 
1965  x = wwprod * tab(k,h)
1966 
1967  DO q = ia(i), ia(i+1) - 1
1968  IF (ja(q) == j) then; a0(q) = a0(q) + x; exit;
1969  ENDIF
1970  ENDDO
1971 
1972  ENDDO
1973  ENDDO
1974  ENDDO
1975  ENDDO
1976  ENDDO
1977  ENDDO
1978 
1979  END SUBROUTINE cc_00_s_m
1980 
1981  SUBROUTINE qs_dif_massvar_m (mesh, visco, ff, ia, ja, a0)
1982  !===========================================================
1983 
1984  ! << visco (Dw), D_) >>
1985  ! + < w, ff _ >
1986  ! ===> a0 ! incremental accumulation
1987 
1988  USE gauss_points
1989 
1990  IMPLICIT NONE
1991  REAL(KIND=8), INTENT(IN) :: visco
1992  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
1993  INTEGER, DIMENSION(:), INTENT(IN) :: ia, ja
1994  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
1995 
1996  TYPE(mesh_type), TARGET :: mesh
1997  INTEGER, DIMENSION(:,:), POINTER :: jj
1998  INTEGER, POINTER :: me
1999 
2000  INTEGER :: k, l, m, ni, nj, i, j, p
2001  REAL(KIND=8) :: x, xij, masslm, viscolm
2002  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
2003 
2004  CALL gauss(mesh)
2005  jj => mesh%jj
2006  me => mesh%me
2007 
2008  DO l = 1, l_g
2009  DO ni = 1, n_w
2010  DO nj = 1, n_w
2011  wwprod(ni,nj,l) = ww(ni,l)*ww(nj,l)
2012  END DO
2013  END DO
2014  END DO
2015 
2016  DO m = 1, me
2017  DO l = 1, l_g
2018 
2019  viscolm = visco*rj(l,m)
2020 
2021  masslm = 0
2022  DO ni = 1, n_w
2023  masslm = masslm + ff(jj(ni,m))*ww(ni,l)
2024  END DO
2025  masslm = masslm * rj(l,m)
2026 
2027  DO ni = 1, n_w; i = jj(ni, m)
2028  DO nj = 1, n_w; j = jj(nj, m)
2029 
2030  xij = 0.
2031  DO k = 1, k_d
2032  xij = xij + dw(k,nj,l,m) * dw(k,ni,l,m)
2033  END DO
2034 
2035  x = viscolm*xij + masslm*wwprod(ni,nj,l)
2036 
2037  DO p = ia(i), ia(i+1) - 1
2038  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
2039  ENDIF
2040  ENDDO
2041 
2042  ENDDO
2043  ENDDO
2044 
2045  ENDDO
2046  ENDDO
2047 
2048  END SUBROUTINE qs_dif_massvar_m
2049 
2050  SUBROUTINE qs_massvar_adv_m (mesh, ff, gg, ia, ja, a0, a_skew)
2051  !===================================================================
2052 
2053  ! + < w, ff _ >
2054  ! + < w, (g.D)_ > ! + skew * < w, (D.g) _ >
2055  ! ===> a0 ! incremental accumulation
2056 
2057  USE gauss_points
2058 
2059  IMPLICIT NONE
2060  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
2061  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
2062  INTEGER, DIMENSION(:), INTENT(IN) :: ia, ja
2063  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
2064  REAL(KIND=8), OPTIONAL, INTENT(IN) :: a_skew
2065 
2066  TYPE(mesh_type), TARGET :: mesh
2067  INTEGER, DIMENSION(:,:), POINTER :: jj
2068  INTEGER, POINTER :: me
2069 
2070  INTEGER :: k, l, m, ni, nj, i, j, p
2071  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: gl
2072  REAL(KIND=8) :: fl, dg, x, masslm, skew=1
2073  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
2074  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: y
2075  ! REAL(KIND=8), DIMENSION(k_d,n_w) :: dw_loc
2076 
2077  CALL gauss(mesh)
2078  jj => mesh%jj
2079  me => mesh%me
2080 
2081  DO l = 1, l_g
2082  DO ni = 1, n_w
2083  DO nj = 1, n_w
2084  wwprod(ni,nj,l) = ww(ni,l)*ww(nj,l)
2085  END DO
2086  END DO
2087  END DO
2088 
2089  IF (present(a_skew)) THEN
2090  skew = a_skew
2091  ELSE
2092  skew = 0.5d0
2093  END IF
2094 
2095  DO m = 1, me
2096  DO l = 1, l_g
2097 
2098  ! dw_loc = dw(:,:,l,m)
2099 
2100  dg = 0.
2101  gl = 0.
2102  fl = 0
2103  DO ni =1 ,n_w
2104  fl = fl + ff(jj(ni,m)) * ww(ni,l)
2105  DO k = 1, k_d
2106  gl(k) = gl(k) + gg(k, jj(ni,m)) * ww(ni,l)
2107  dg = dg + gg(k, jj(ni,m)) * dw(k,ni,l,m) ! divergence
2108  END DO
2109  ENDDO
2110  masslm = (fl + skew*dg)*rj(l,m) ! skew form
2111  ! masslm = fl * rj(l,m)
2112 
2113  y = 0.
2114  DO ni = 1, n_w;
2115  DO k = 1, k_d
2116  y(ni) = y(ni) + gl(k) * dw(k,ni,l,m)
2117  END DO
2118  END DO
2119  y = rj(l,m)*y
2120 
2121 
2122  DO ni = 1, n_w; i = jj(ni, m)
2123  DO nj = 1, n_w; j = jj(nj, m)
2124 
2125  x = masslm*wwprod(ni,nj,l) + y(nj)*ww(ni,l)
2126 
2127  DO p = ia(i), ia(i+1) - 1
2128  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
2129  ENDIF
2130  ENDDO
2131 
2132  ENDDO
2133  ENDDO
2134 
2135  ENDDO
2136  ENDDO
2137 
2138  END SUBROUTINE qs_massvar_adv_m
2139 
2140  SUBROUTINE qs_varden_adv_m (mesh, mass, ff, gg, ia, ja, a0, a_skew)
2141  !===================================================================
2142 
2143  ! + < w, ff _ >
2144  ! + < w, ( ff*g.D)_ > ! + skew * < w, ff*(D.g) _ >
2145  ! ===> a0 ! incremental accumulation
2146 
2147  USE gauss_points
2148 
2149  IMPLICIT NONE
2150  REAL(KIND=8), INTENT(IN) :: mass
2151  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
2152  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
2153  INTEGER, DIMENSION(:), INTENT(IN) :: ia, ja
2154  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
2155  REAL(KIND=8), OPTIONAL, INTENT(IN) :: a_skew
2156 
2157  TYPE(mesh_type), TARGET :: mesh
2158  INTEGER, DIMENSION(:,:), POINTER :: jj
2159  INTEGER, POINTER :: me
2160 
2161  INTEGER :: k, l, m, ni, nj, i, j, p
2162  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: gl
2163  REAL(KIND=8) :: fl, dg, x, masslm, skew=1
2164  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
2165  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: y
2166 
2167 
2168  CALL gauss(mesh)
2169  jj => mesh%jj
2170  me => mesh%me
2171 
2172  DO l = 1, l_g
2173  DO ni = 1, n_w
2174  DO nj = 1, n_w
2175  wwprod(ni,nj,l) = ww(ni,l)*ww(nj,l)
2176  END DO
2177  END DO
2178  END DO
2179 
2180  IF (present(a_skew)) THEN
2181  skew = a_skew
2182  ELSE
2183  skew = 0.25d0
2184  END IF
2185 
2186  DO m = 1, me
2187  DO l = 1, l_g
2188 
2189  dg = 0.
2190  gl = 0.
2191  fl = 0
2192  DO ni =1 ,n_w
2193  fl = fl + ff(jj(ni,m)) * ww(ni,l)
2194  DO k = 1, k_d
2195  gl(k) = gl(k) + gg(k, jj(ni,m)) * ww(ni,l)
2196  dg = dg + gg(k, jj(ni,m)) * dw(k,ni,l,m) ! divergence
2197  END DO
2198  ENDDO
2199  masslm = (mass*fl + fl*skew*dg)*rj(l,m) ! skew form
2200 
2201  y = 0.
2202  DO ni = 1, n_w;
2203  DO k = 1, k_d
2204  y(ni) = y(ni) + gl(k) * dw(k,ni,l,m)
2205  END DO
2206  END DO
2207  y = rj(l,m)*fl*y
2208 
2209 
2210  DO ni = 1, n_w; i = jj(ni, m)
2211  DO nj = 1, n_w; j = jj(nj, m)
2212 
2213  x = masslm*wwprod(ni,nj,l) + y(nj)*ww(ni,l)
2214 
2215  DO p = ia(i), ia(i+1) - 1
2216  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
2217  ENDIF
2218  ENDDO
2219 
2220  ENDDO
2221  ENDDO
2222 
2223  ENDDO
2224  ENDDO
2225 
2226  END SUBROUTINE qs_varden_adv_m
2227 
2228  SUBROUTINE qs_mass_adv_m (mesh, alpha, gg, ia, ja, a0, a_skew)
2229  !===================================================================
2230 
2231  ! + alpha * < w, _ >
2232  ! + < w, (g.D)_ > + < w, (D.g) _ > / 2
2233  ! ===> a0 ! incremental accumulation
2234 
2235  USE gauss_points
2236 
2237  IMPLICIT NONE
2238  REAL(KIND=8), INTENT(IN) :: alpha
2239  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
2240  INTEGER, DIMENSION(:), INTENT(IN) :: ia, ja
2241  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
2242  REAL(KIND=8), OPTIONAL, INTENT(IN) :: a_skew
2243 
2244  TYPE(mesh_type), TARGET :: mesh
2245  INTEGER, DIMENSION(:,:), POINTER :: jj
2246  INTEGER, POINTER :: me
2247 
2248  INTEGER :: k, l, m, ni, nj, i, j, p
2249  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: gl
2250  REAL(KIND=8) :: dg, x, masslm, skew
2251  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
2252  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: y
2253 
2254  CALL gauss(mesh)
2255  jj => mesh%jj
2256  me => mesh%me
2257 
2258  DO l = 1, l_g
2259  DO ni = 1, n_w
2260  DO nj = 1, n_w
2261  wwprod(ni,nj,l) = ww(ni,l)*ww(nj,l)
2262  END DO
2263  END DO
2264  END DO
2265 
2266  IF (present(a_skew)) THEN
2267  skew = a_skew
2268  ELSE
2269  skew = 0.5d0
2270  END IF
2271 
2272 
2273  DO m = 1, me
2274  DO l = 1, l_g
2275 
2276 
2277  dg = 0.
2278  gl = 0.
2279  DO k = 1, k_d
2280  DO ni =1 ,n_w
2281  gl(k) = gl(k) + gg(k, jj(ni,m)) * ww(ni,l)
2282  dg = dg + gg(k, jj(ni,m)) * dw(k,ni,l,m)
2283  END DO
2284  ENDDO
2285  masslm = (alpha+skew*dg)*rj(l,m)
2286 
2287  y = 0.
2288  DO ni = 1, n_w;
2289  DO k = 1, k_d
2290  y(ni) = y(ni) + gl(k) * dw(k,ni,l,m)
2291  END DO
2292  END DO
2293  y = rj(l,m)*y
2294 
2295  DO ni = 1, n_w; i = jj(ni, m)
2296  DO nj = 1, n_w; j = jj(nj, m)
2297 
2298  x = masslm*wwprod(ni,nj,l) + &
2299  y(nj)*ww(ni,l)
2300 
2301  DO p = ia(i), ia(i+1) - 1
2302  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
2303  ENDIF
2304  ENDDO
2305 
2306  ENDDO
2307  ENDDO
2308 
2309  ENDDO
2310  ENDDO
2311 
2312  END SUBROUTINE qs_mass_adv_m
2313 
2314 
2315  SUBROUTINE qs_mass_div_m (mesh, alpha, gg, ia, ja, a0)
2316  !===================================================================
2317 
2318  ! + alpha * < w, _ >
2319  ! + < w, (g.D)_ > + < w, (D.g) _ >
2320  ! ===> a0 ! incremental accumulation
2321 
2322  USE gauss_points
2323 
2324  IMPLICIT NONE
2325  REAL(KIND=8), INTENT(IN) :: alpha
2326  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
2327  INTEGER, DIMENSION(:), INTENT(IN) :: ia, ja
2328  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
2329 
2330  TYPE(mesh_type), TARGET :: mesh
2331  INTEGER, DIMENSION(:,:), POINTER :: jj
2332  INTEGER, POINTER :: me
2333 
2334  INTEGER :: k, l, m, ni, nj, i, j, p
2335  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: gl
2336  REAL(KIND=8) :: dg, x, masslm
2337  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
2338  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: y
2339  ! REAL(KIND=8), DIMENSION(k_d,n_w) :: dw_loc
2340 
2341  CALL gauss(mesh)
2342  jj => mesh%jj
2343  me => mesh%me
2344 
2345  DO l = 1, l_g
2346  DO ni = 1, n_w
2347  DO nj = 1, n_w
2348  wwprod(ni,nj,l) = ww(ni,l)*ww(nj,l)
2349  END DO
2350  END DO
2351  END DO
2352 
2353  DO m = 1, me
2354  DO l = 1, l_g
2355 
2356  ! dw_loc = dw(:,:,l,m)
2357 
2358  dg = 0.
2359  gl = 0.
2360  DO k = 1, k_d
2361  DO ni =1 ,n_w
2362  gl(k) = gl(k) + gg(k, jj(ni,m)) * ww(ni,l)
2363  dg = dg + gg(k, jj(ni,m)) * dw(k,ni,l,m)
2364  END DO
2365  ENDDO
2366  masslm = (alpha+dg)*rj(l,m)
2367 
2368  y = 0.
2369  DO ni = 1, n_w;
2370  DO k = 1, k_d
2371  y(ni) = y(ni) + gl(k) * dw(k,ni,l,m)
2372  END DO
2373  END DO
2374  y = rj(l,m)*y
2375 
2376  DO ni = 1, n_w; i = jj(ni, m)
2377  DO nj = 1, n_w; j = jj(nj, m)
2378 
2379  x = masslm*wwprod(ni,nj,l) + &
2380  y(nj)*ww(ni,l)
2381 
2382  DO p = ia(i), ia(i+1) - 1
2383  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
2384  ENDIF
2385  ENDDO
2386 
2387  ENDDO
2388  ENDDO
2389 
2390  ENDDO
2391  ENDDO
2392 
2393  END SUBROUTINE qs_mass_div_m
2394 
2395  SUBROUTINE elast_m (mesh, alpha, ia, ja, a0)
2396  !=======================================
2397 
2398  ! << (Dw), (D_) >>
2399  ! << (Dw), (D_)^t >>
2400  ! alpha << (D.w), (D._) >>
2401  ! ===> a0 incremental acc. in k_d Dimensions
2402 
2403  USE gauss_points
2404 
2405  IMPLICIT NONE
2406  REAL(KIND=8), INTENT(IN) :: alpha
2407  INTEGER, DIMENSION(:), INTENT(IN) :: ia
2408  INTEGER, DIMENSION(:), INTENT(IN) :: ja
2409  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
2410 
2411  INTEGER :: m, l, p, ni, nj, k, k1, h, i, j, i_b, j_b, &
2412  bloc_size
2413  REAL(KIND=8) :: x
2414 
2415  TYPE(mesh_type), TARGET :: mesh
2416  INTEGER, DIMENSION(:,:), POINTER :: jj
2417  INTEGER, POINTER :: me
2418 
2419  CALL gauss(mesh)
2420  jj => mesh%jj
2421  me => mesh%me
2422 
2423  bloc_size = (SIZE(ia) - 1)/k_d
2424 
2425  DO m = 1, me
2426  DO l = 1, l_g
2427 
2428  DO ni = 1, n_w; i_b = jj(ni, m)
2429  DO k = 1, k_d;
2430  i = i_b + (k-1)*bloc_size
2431 
2432  DO nj = 1, n_w; j_b = jj(nj,m)
2433  DO h = 1, k_d;
2434  j = j_b + (h-1)*bloc_size
2435 
2436  x = dw(h,ni,l,m)*dw(k,nj,l,m) &
2437  + alpha*dw(k,ni,l,m)*dw(h,nj,l,m)
2438  IF (h.EQ.k) THEN
2439  DO k1 = 1, k_d
2440  x = x + dw(k1,ni,l,m)*dw(k1,nj,l,m)
2441  END DO
2442  END IF
2443 
2444  x = x * rj(l,m)
2445 
2446  DO p = ia(i), ia(i+1) - 1
2447  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
2448  ENDIF
2449  ENDDO
2450 
2451  ENDDO
2452  ENDDO
2453 
2454  ENDDO
2455  ENDDO
2456 
2457  ENDDO
2458  ENDDO
2459 
2460  END SUBROUTINE elast_m
2461 
2462  SUBROUTINE curl_div_2d_m (mesh, alpha, stabh, expdiv, ia, ja, a0)
2463  !=======================================
2464  ! << w , _ >>
2465  ! + << (Dxw), (Dx_) >>
2466  ! + alpha << (D.w), (D._) >>
2467  ! ===> a0 incremental acc. in k_d Dimensions
2468 
2469  USE gauss_points
2470 
2471  IMPLICIT NONE
2472  REAL(KIND=8), INTENT(IN) :: stabh, expdiv
2473  INTEGER, DIMENSION(:), INTENT(IN) :: ia
2474  INTEGER, DIMENSION(:), INTENT(IN) :: ja
2475  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
2476  REAL(KIND=8), INTENT(IN) :: alpha
2477  INTEGER :: m, l, p, ni, nj, k, k1, h, i, j, i_b, j_b, &
2478  bloc_size
2479  REAL(KIND=8) :: x, alphah
2480 
2481  TYPE(mesh_type), TARGET :: mesh
2482  INTEGER, DIMENSION(:,:), POINTER :: jj
2483  INTEGER, POINTER :: me
2484 
2485  CALL gauss(mesh)
2486  jj => mesh%jj
2487  me => mesh%me
2488 
2489  bloc_size = (SIZE(ia) - 1)/2
2490 
2491  DO m = 1, me
2492  alphah = stabh*sum(rj(:,m))**(expdiv/2)
2493  DO l = 1, l_g
2494 
2495  DO ni = 1, n_w; i_b = jj(ni, m)
2496  DO k = 1, 2
2497  i = i_b + (k-1)*bloc_size
2498  k1 = modulo(k,2) + 1
2499 
2500  DO nj = 1, n_w; j_b = jj(nj,m)
2501  DO h = 1, 2
2502  j = j_b + (h-1)*bloc_size
2503 
2504  x = alphah*dw(k,ni,l,m)*dw(h,nj,l,m)
2505  IF (h==k) THEN
2506  x = x + alpha*ww(ni,l)*ww(nj,l) + dw(k1,ni,l,m)*dw(k1,nj,l,m)
2507  ELSE
2508  x = x - dw(h,ni,l,m)*dw(k,nj,l,m)
2509  END IF
2510 
2511  x = x * rj(l,m)
2512 
2513  DO p = ia(i), ia(i+1) - 1
2514  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
2515  ENDIF
2516  ENDDO
2517 
2518  ENDDO
2519  ENDDO
2520 
2521  ENDDO
2522  ENDDO
2523 
2524  ENDDO
2525  ENDDO
2526 
2527  END SUBROUTINE curl_div_2d_m
2528 
2529  SUBROUTINE curl_grad_2d_m (mesh, alpha, stab, exp_sta, ia, ja, a0)
2530  !=======================================
2531  ! alpha(phij, phii) + (Curl phij, Curl phii) + (Grad psij, phii) + sta*h^(2*exp_sta)*(div phij, Div phii)
2532  ! - (Grad phij, psii) + h^(2(1-exp_sta))*(Grad psij, Grad psii)
2533  ! ===> a0 incremental acc. in k_d Dimensions
2534 
2535  USE gauss_points
2536 
2537  IMPLICIT NONE
2538  REAL(KIND=8), INTENT(IN) :: stab, exp_sta
2539  INTEGER, DIMENSION(:), INTENT(IN) :: ia
2540  INTEGER, DIMENSION(:), INTENT(IN) :: ja
2541  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
2542  REAL(KIND=8), INTENT(IN) :: alpha
2543  INTEGER :: m, l, p, ni, nj, k, k1, h, i, j, i_b, j_b, &
2544  bloc_size, type_fe
2545  REAL(KIND=8) :: x, alphah, betah, hloc
2546 
2547  TYPE(mesh_type), TARGET :: mesh
2548  INTEGER, DIMENSION(:,:), POINTER :: jj
2549  INTEGER, POINTER :: me
2550 
2551  CALL gauss(mesh)
2552  jj => mesh%jj
2553  me => mesh%me
2554 
2555  bloc_size = (SIZE(ia) - 1)/3
2556  IF (mesh%gauss%n_w == 3) THEN
2557  type_fe = 1
2558  ELSE IF (mesh%gauss%n_w == 6) THEN
2559  type_fe = 2
2560  ELSE
2561  WRITE(*,*) ' BUG, bad FE'
2562  END IF
2563 
2564  DO m = 1, me
2565  hloc = sqrt(sum(rj(:,m)))/type_fe
2566  alphah = stab*hloc**(2*exp_sta)
2567  betah = hloc**(2*(1-exp_sta))*(1.d0/stab)
2568  DO l = 1, l_g
2569 
2570  DO ni = 1, n_w; i_b = jj(ni, m)
2571  DO k = 1, 2
2572  i = i_b + (k-1)*bloc_size
2573  k1 = modulo(k,2) + 1
2574 
2575  DO nj = 1, n_w; j_b = jj(nj,m)
2576  DO h = 1, 2
2577  j = j_b + (h-1)*bloc_size
2578 
2579  x = alphah*dw(k,ni,l,m)*dw(h,nj,l,m)
2580  IF (h==k) THEN
2581  x = x + alpha*ww(ni,l)*ww(nj,l) + dw(k1,ni,l,m)*dw(k1,nj,l,m)
2582  ELSE
2583  x = x - dw(h,ni,l,m)*dw(k,nj,l,m)
2584  END IF
2585  x = x * rj(l,m)
2586  DO p = ia(i), ia(i+1) - 1
2587  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
2588  ENDIF
2589  ENDDO
2590 
2591  END DO
2592  ENDDO
2593 
2594  ENDDO
2595  ENDDO
2596 
2597  DO ni = 1, n_w; i_b = jj(ni, m)
2598  DO k = 1, 2
2599  i = i_b + (k-1)*bloc_size
2600  DO nj = 1, n_w; j_b = jj(nj,m)
2601  h = 3
2602  j = j_b + (h-1)*bloc_size
2603  x = dw(k,nj,l,m)*ww(ni,l)* rj(l,m)
2604  DO p = ia(i), ia(i+1) - 1
2605  IF (ja(p) == j) then; a0(p) = a0(p) + x;
2606  EXIT
2607  ENDIF
2608  ENDDO
2609 
2610  ENDDO
2611  ENDDO
2612  ENDDO
2613 
2614  DO ni = 1, n_w; i_b = jj(ni, m)
2615  k = 3
2616  i = i_b + (k-1)*bloc_size
2617  DO nj = 1, n_w; j_b = jj(nj,m)
2618  DO h = 1, 3
2619  j = j_b + (h-1)*bloc_size
2620  IF (h<3) THEN
2621  x = -ww(nj,l) * dw(h,ni,l,m)* rj(l,m)
2622  ELSE
2623  x = betah*sum(dw(:,ni,l,m) * dw(:,nj,l,m))* rj(l,m)
2624  END IF
2625  DO p = ia(i), ia(i+1) - 1
2626  IF (ja(p) == j) then; a0(p) = a0(p) + x;
2627  EXIT
2628  ENDIF
2629  ENDDO
2630  END DO
2631  END DO
2632 
2633  ENDDO
2634  ENDDO
2635  END DO
2636  END SUBROUTINE curl_grad_2d_m
2637 
2638  SUBROUTINE curl_surf_2d_m (mesh, stab_b, exp_b, consist, adj, ia, ja, a0)
2639  !=======================================
2640  ! alpha*h << wxn , _xn >>
2641  ! - << (Dxw)xn, _ >>
2642  ! ===> a0 incremental acc. in k_d Dimensions
2643 
2644  USE gauss_points
2645 
2646  IMPLICIT NONE
2647  INTEGER, DIMENSION(:), INTENT(IN) :: ia
2648  INTEGER, DIMENSION(:), INTENT(IN) :: ja
2649  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
2650  REAL(KIND=8), INTENT(IN) :: stab_b, consist, adj, exp_b
2651  INTEGER :: m, ms, ls, p, ni, nj, k, k1, h, h1, i, j, i_b, j_b, &
2652  bloc_size
2653  REAL(KIND=8) :: x, alphah
2654 
2655  TYPE(mesh_type), TARGET :: mesh
2656  INTEGER, DIMENSION(:,:), POINTER :: jjs
2657  INTEGER, POINTER :: mes
2658 
2659  CALL gauss(mesh)
2660  jjs => mesh%jjs
2661  mes => mesh%mes
2662 
2663  bloc_size = mesh%np
2664 
2665  DO ms = 1, mes
2666  alphah = stab_b*sum(rjs(:,ms))**(exp_b)
2667  m =mesh%neighs(ms)
2668 
2669  DO ls = 1, l_gs
2670 
2671  DO ni = 1, n_ws; i_b = jjs(ni, ms)
2672  DO k = 1, 2
2673  i = i_b + (k-1)*bloc_size
2674  k1 = modulo(k,2) + 1
2675 
2676  DO nj = 1, n_ws; j_b = jjs(nj,ms)
2677  DO h = 1, 2
2678  j = j_b + (h-1)*bloc_size
2679  h1 = modulo(h,2) + 1
2680 
2681  x = alphah*(-1)**(k1+h1)*wws(ni,ls)*rnorms(k1,ls,ms)*wws(nj,ls)*rnorms(h1,ls,ms)
2682  x = x * rjs(ls,ms)
2683 
2684  DO p = ia(i), ia(i+1) - 1
2685  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
2686  ENDIF
2687  ENDDO
2688 
2689  ENDDO
2690  ENDDO
2691 
2692  ENDDO
2693  ENDDO
2694 
2695  DO ni = 1, n_ws; i_b = jjs(ni, ms)
2696  DO k = 1, 2
2697  i = i_b + (k-1)*bloc_size
2698  k1 = modulo(k,2) + 1
2699 
2700  DO nj = 1, n_w; j_b = mesh%jj(nj,m)
2701  DO h = 1, 2
2702  j = j_b + (h-1)*bloc_size
2703  h1 = modulo(h,2) + 1
2704 
2705  x =consist*(-1)**(k1+h) *wws(ni,ls)*rnorms(k1,ls,ms)*dw_s(h1,nj,ls,ms)
2706  x = x * rjs(ls,ms)
2707 
2708  DO p = ia(i), ia(i+1) - 1
2709  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
2710  ENDIF
2711  ENDDO
2712 
2713  ENDDO
2714  ENDDO
2715 
2716  ENDDO
2717  ENDDO
2718 
2719  DO ni = 1, n_w; i_b = mesh%jj(ni,m)
2720  DO k = 1, 2
2721  i = i_b + (k-1)*bloc_size
2722  k1 = modulo(k,2) + 1
2723 
2724  DO nj = 1, n_ws; j_b = jjs(nj,ms)
2725  DO h = 1, 2
2726  j = j_b + (h-1)*bloc_size
2727  h1 = modulo(h,2) + 1
2728 
2729  x = adj*(-1)**(k+h1) *wws(nj,ls)*rnorms(h1,ls,ms)*dw_s(k1,ni,ls,ms)
2730  x = x * rjs(ls,ms)
2731 
2732  DO p = ia(i), ia(i+1) - 1
2733  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
2734  ENDIF
2735  ENDDO
2736 
2737  ENDDO
2738  ENDDO
2739 
2740  ENDDO
2741  ENDDO
2742 
2743  ENDDO
2744  END DO
2745 
2746  END SUBROUTINE curl_surf_2d_m
2747 
2748  SUBROUTINE qs_1x1x_m (mesh, alpha, ia, ja, a0)
2749  !=================================================
2750 
2751  ! alpha << (Dw), (D_) >> ===> A0 incremental accumulation
2752 
2753  USE gauss_points
2754 
2755  IMPLICIT NONE
2756  REAL(KIND=8), INTENT(IN) :: alpha
2757  INTEGER, DIMENSION(:), INTENT(IN) :: ia
2758  INTEGER, DIMENSION(:), INTENT(IN) :: ja
2759  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: a0
2760 
2761  INTEGER :: m, l, ni, nj, i, j, p
2762  REAL(KIND=8) :: al, x
2763 
2764  TYPE(mesh_type), TARGET :: mesh
2765  INTEGER, DIMENSION(:,:), POINTER :: jj
2766  INTEGER, POINTER :: me
2767 
2768  CALL gauss(mesh)
2769  jj => mesh%jj
2770  me => mesh%me
2771 
2772  DO m = 1, me
2773  DO l = 1, l_g
2774 
2775  al = alpha * rj(l,m)
2776 
2777  DO nj = 1, n_w; j = jj(nj, m)
2778  DO ni = 1, n_w; i = jj(ni, m)
2779 
2780  x = al * dw(1,nj,l,m) * dw(1,ni,l,m)
2781  DO p = ia(i), ia(i+1) - 1
2782  IF (ja(p) == j) then; a0(p) = a0(p) + x; exit;
2783  ENDIF
2784  ENDDO
2785 
2786  ENDDO
2787  ENDDO
2788 
2789  ENDDO
2790  ENDDO
2791 
2792  END SUBROUTINE qs_1x1x_m
2793 
2794  SUBROUTINE edge_stab_m(mesh, coeff_visc, ia, ja, aa)
2795  USE def_type_mesh
2796  !USE solve_sp
2797  IMPLICIT NONE
2798  TYPE(mesh_type) :: mesh
2799  !TYPE(csr_matrix) :: mat
2800  INTEGER, DIMENSION(:), INTENT(IN) :: ia
2801  INTEGER, DIMENSION(:), INTENT(IN) :: ja
2802  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: aa
2803  REAL(KIND=8) :: coeff_visc
2804  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%l_Gs, 2) :: dwni_loc
2805  INTEGER, DIMENSION(mesh%gauss%n_w,2) :: jji_loc
2806  INTEGER :: ls, ms, cotei, cotej, p, ni, nj, i, j, type_fe
2807  REAL(KIND=8) :: x, h2, coeff != 0.02! 0.003! .003
2808  IF (mesh%gauss%n_w==3) THEN
2809  type_fe = 1
2810  ELSE
2811  type_fe = 2
2812  END IF
2813  coeff = coeff_visc/type_fe**2
2814  DO ms = 1, mesh%mi
2815  dwni_loc = mesh%gauss%dwni(:,:,:,ms)
2816  jji_loc = mesh%jji(:,:,ms)
2817  h2 = coeff*sum(mesh%gauss%rji(:,ms))**2
2818  DO cotei = 1, 2
2819  DO ni = 1, mesh%gauss%n_w
2820  i = jji_loc(ni,cotei)
2821  DO cotej = 1, 2
2822  DO nj = 1, mesh%gauss%n_w
2823  j = jji_loc(nj,cotej)
2824  x = 0.d0
2825  DO ls = 1, mesh%gauss%l_Gs
2826  x = x + dwni_loc(ni,ls,cotei)*dwni_loc(nj,ls,cotej)*mesh%gauss%rji(ls,ms)
2827  END DO
2828  DO p = ia(i), ia(i+1) - 1
2829  IF (ja(p) == j) then;
2830  aa(p) = aa(p) + x*h2;
2831  EXIT
2832  ENDIF
2833  ENDDO
2834  END DO
2835  END DO
2836  END DO
2837  END DO
2838  END DO
2839  END SUBROUTINE edge_stab_m
2840 
2841 END MODULE fem_s_m
2842 
2843 
2844 
2845 
2846 
2847 
2848 
2849 
fem_s_m::qs_00_s_m
subroutine qs_00_s_m(mesh, fs, ia, ja, a0)
Definition: fem_sparsekit.f90:1653
fem_s_m::qs_varden_adv_m
subroutine qs_varden_adv_m(mesh, mass, ff, gg, ia, ja, a0, a_skew)
Definition: fem_sparsekit.f90:2140
fem_s_m::qs_mass_adv_m
subroutine qs_mass_adv_m(mesh, alpha, gg, ia, ja, a0, a_skew)
Definition: fem_sparsekit.f90:2228
fem_s_m::qs_h_v_grad_v_grad_w
subroutine qs_h_v_grad_v_grad_w(mesh, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:492
fem_s_m::edge_stab_m
subroutine edge_stab_m(mesh, coeff_visc, ia, ja, aa)
Definition: fem_sparsekit.f90:2794
fem_s_m::qs_adv_m
subroutine qs_adv_m(mesh, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:976
fem_s_m::qs_diff_mass_adv_m
subroutine qs_diff_mass_adv_m(mesh, visco, alpha, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:660
fem_s_m::cv_11_m
subroutine cv_11_m(mesh, alpha, ia, ja, a0)
Definition: fem_sparsekit.f90:1703
fem_s_m::curl_surf_2d_m
subroutine curl_surf_2d_m(mesh, stab_b, exp_b, consist, adj, ia, ja, a0)
Definition: fem_sparsekit.f90:2638
gauss_points
Definition: associate_gauss.f90:4
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:73
fem_s_m::qs_001_m
subroutine qs_001_m(mesh, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:1059
fem_s_m::qs_dif_mass_adv_m
subroutine qs_dif_mass_adv_m(mesh, visco, alpha, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:569
fem_s_m::qs_ls_mass_adv
subroutine qs_ls_mass_adv(mesh, alpha, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:341
fem_s_m::qs_adv_m_bis
subroutine qs_adv_m_bis(mesh, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:901
fem_s_m::qs_v_grad_v_grad_w
subroutine qs_v_grad_v_grad_w(mesh, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:271
fem_s_m::qs_dif_massvar_m
subroutine qs_dif_massvar_m(mesh, visco, ff, ia, ja, a0)
Definition: fem_sparsekit.f90:1981
fem_s_m::qs_massvar_adv_m
subroutine qs_massvar_adv_m(mesh, ff, gg, ia, ja, a0, a_skew)
Definition: fem_sparsekit.f90:2050
fem_s_m
Definition: fem_sparsekit.f90:1
fem_s_m::qs_100_m
subroutine qs_100_m(mesh, vv, ia, ja, a0)
Definition: fem_sparsekit.f90:55
f
section doc_intro_frame_work_num_app Numerical approximation subsection doc_intro_fram_work_num_app_Fourier_FEM Fourier Finite element representation The SFEMaNS code uses a hybrid Fourier Finite element formulation The Fourier decomposition allows to approximate the problem’s solutions for each Fourier mode modulo nonlinear terms that are made explicit The variables are then approximated on a meridian section of the domain with a finite element method The numerical approximation of a function f $f f is written in the following generic f with f(r,\theta, z)\f $the cylindrical coordinates
fem_s_m::cc_101_m
subroutine cc_101_m(mesh, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:1768
def_type_mesh
Definition: def_type_mesh.f90:4
fem_s_m::bs_101_m
subroutine bs_101_m(mesh, ff, ia, ja, a0)
Definition: fem_sparsekit.f90:1601
fem_s_m::qs_11_bb_p1_m
subroutine qs_11_bb_p1_m(mesh, alpha, bcd)
Definition: fem_sparsekit.f90:221
fem_s_m::qs_plap_m
subroutine qs_plap_m(mesh, p, ff, ia, ja, a0)
Definition: fem_sparsekit.f90:1390
fem_s_m::qs_h_100_m
subroutine qs_h_100_m(mesh, alpha, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:1122
fem_s_m::qs_gals_mass_adv_m
subroutine qs_gals_mass_adv_m(mesh, param, alpha, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:414
fem_s_m::qs_11_m
subroutine qs_11_m(mesh, alpha, ia, ja, a0)
Definition: fem_sparsekit.f90:171
gauss_points::gauss
subroutine gauss(mesh)
Definition: associate_gauss.f90:23
fem_s_m::qs_v_plap_m
subroutine qs_v_plap_m(mesh, p, vstar, ff, ia, ja, a0)
Definition: fem_sparsekit.f90:1318
fem_s_m::elast_m
subroutine elast_m(mesh, alpha, ia, ja, a0)
Definition: fem_sparsekit.f90:2395
fem_s_m::qs_1x1x_m
subroutine qs_1x1x_m(mesh, alpha, ia, ja, a0)
Definition: fem_sparsekit.f90:2748
fem_s_m::cv_00_s_m
subroutine cv_00_s_m(mesh, fs, ia, ja, a0)
Definition: fem_sparsekit.f90:1841
fem_s_m::curl_div_2d_m
subroutine curl_div_2d_m(mesh, alpha, stabh, expdiv, ia, ja, a0)
Definition: fem_sparsekit.f90:2462
fem_s_m::qs_00_m
subroutine qs_00_m(mesh, alpha, ia, ja, a0)
Definition: fem_sparsekit.f90:6
fem_s_m::cc_00_s_m
subroutine cc_00_s_m(mesh, alpha, ia, ja, a0)
Definition: fem_sparsekit.f90:1912
fem_s_m::qs_mass_div_m
subroutine qs_mass_div_m(mesh, alpha, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:2315
fem_s_m::qs_000_m
subroutine qs_000_m(mesh, ff, ia, ja, a0)
Definition: fem_sparsekit.f90:115
fem_s_m::curl_grad_2d_m
subroutine curl_grad_2d_m(mesh, alpha, stab, exp_sta, ia, ja, a0)
Definition: fem_sparsekit.f90:2529
fem_s_m::qs_dif_mass_adv_skew_m
subroutine qs_dif_mass_adv_skew_m(mesh, visco, alpha, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:752
fem_s_m::qs_dif_mass_m
subroutine qs_dif_mass_m(mesh, visco, alpha, ia, ja, a0)
Definition: fem_sparsekit.f90:836
fem_s_m::qs_101_m
subroutine qs_101_m(mesh, ff, ia, ja, a0)
Definition: fem_sparsekit.f90:1262
fem_s_m::qs_1_sgs_1_m
subroutine qs_1_sgs_1_m(mesh, ff, ia, ja, a0)
Definition: fem_sparsekit.f90:1192
fem_s_m::qs_adv_stab_m
subroutine qs_adv_stab_m(mesh, gg, stab, istab, ia, ja, a0)
Definition: fem_sparsekit.f90:1459