SFEMaNS  version 5.3
Reference documentation for SFEMaNS
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 
gauss_points::k_d
integer, public k_d
Definition: associate_gauss.f90:9
fem_s_m::cv_11_m
subroutine cv_11_m(mesh, alpha, ia, ja, a0)
Definition: fem_sparsekit.f90:1704
gauss_points
Definition: associate_gauss.f90:4
fem_s_m::cc_00_s_m
subroutine cc_00_s_m(mesh, alpha, ia, ja, a0)
Definition: fem_sparsekit.f90:1913
gauss_points::l_g
integer, public l_g
Definition: associate_gauss.f90:9
fem_s_m::qs_100_m
subroutine qs_100_m(mesh, vv, ia, ja, a0)
Definition: fem_sparsekit.f90:56
gauss_points::l_gs
integer, public l_gs
Definition: associate_gauss.f90:9
gauss_points::n_ws
integer, public n_ws
Definition: associate_gauss.f90:9
fem_s_m::elast_m
subroutine elast_m(mesh, alpha, ia, ja, a0)
Definition: fem_sparsekit.f90:2396
fem_s_m::qs_mass_div_m
subroutine qs_mass_div_m(mesh, alpha, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:2316
fem_s_m::qs_dif_mass_m
subroutine qs_dif_mass_m(mesh, visco, alpha, ia, ja, a0)
Definition: fem_sparsekit.f90:837
fem_s_m::qs_plap_m
subroutine qs_plap_m(mesh, p, ff, ia, ja, a0)
Definition: fem_sparsekit.f90:1391
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:74
fem_s_m::edge_stab_m
subroutine edge_stab_m(mesh, coeff_visc, ia, ja, aa)
Definition: fem_sparsekit.f90:2795
gauss_points::rj
real(kind=8), dimension(:,:), pointer rj
Definition: associate_gauss.f90:14
fem_s_m::qs_1x1x_m
subroutine qs_1x1x_m(mesh, alpha, ia, ja, a0)
Definition: fem_sparsekit.f90:2749
fem_s_m::qs_ls_mass_adv
subroutine qs_ls_mass_adv(mesh, alpha, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:342
gauss_points::rnorms
real(kind=8), dimension(:,:,:), pointer rnorms
Definition: associate_gauss.f90:13
fem_s_m::qs_101_m
subroutine qs_101_m(mesh, ff, ia, ja, a0)
Definition: fem_sparsekit.f90:1263
fem_s_m::qs_001_m
subroutine qs_001_m(mesh, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:1060
fem_s_m::curl_div_2d_m
subroutine curl_div_2d_m(mesh, alpha, stabh, expdiv, ia, ja, a0)
Definition: fem_sparsekit.f90:2463
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:661
fem_s_m::qs_massvar_adv_m
subroutine qs_massvar_adv_m(mesh, ff, gg, ia, ja, a0, a_skew)
Definition: fem_sparsekit.f90:2051
fem_s_m::qs_00_s_m
subroutine qs_00_s_m(mesh, fs, ia, ja, a0)
Definition: fem_sparsekit.f90:1654
fem_s_m::qs_mass_adv_m
subroutine qs_mass_adv_m(mesh, alpha, gg, ia, ja, a0, a_skew)
Definition: fem_sparsekit.f90:2229
fem_s_m::bs_101_m
subroutine bs_101_m(mesh, ff, ia, ja, a0)
Definition: fem_sparsekit.f90:1602
gauss_points::n_w
integer, public n_w
Definition: associate_gauss.f90:9
fem_s_m::qs_adv_stab_m
subroutine qs_adv_stab_m(mesh, gg, stab, istab, ia, ja, a0)
Definition: fem_sparsekit.f90:1460
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:493
fem_s_m::qs_11_bb_p1_m
subroutine qs_11_bb_p1_m(mesh, alpha, bcd)
Definition: fem_sparsekit.f90:222
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:415
fem_s_m
Definition: fem_sparsekit.f90:1
fem_s_m::qs_v_plap_m
subroutine qs_v_plap_m(mesh, p, vstar, ff, ia, ja, a0)
Definition: fem_sparsekit.f90:1319
fem_s_m::qs_00_m
subroutine qs_00_m(mesh, alpha, ia, ja, a0)
Definition: fem_sparsekit.f90:7
fem_s_m::qs_h_100_m
subroutine qs_h_100_m(mesh, alpha, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:1123
fem_s_m::qs_000_m
subroutine qs_000_m(mesh, ff, ia, ja, a0)
Definition: fem_sparsekit.f90:116
fem_s_m::qs_11_m
subroutine qs_11_m(mesh, alpha, ia, ja, a0)
Definition: fem_sparsekit.f90:172
fem_s_m::qs_adv_m_bis
subroutine qs_adv_m_bis(mesh, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:902
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:753
fem_s_m::qs_adv_m
subroutine qs_adv_m(mesh, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:977
fem_s_m::qs_dif_massvar_m
subroutine qs_dif_massvar_m(mesh, visco, ff, ia, ja, a0)
Definition: fem_sparsekit.f90:1982
gauss_points::rjs
real(kind=8), dimension(:,:), pointer rjs
Definition: associate_gauss.f90:15
gauss_points::wws
real(kind=8), dimension(:,:), pointer wws
Definition: associate_gauss.f90:11
gauss_points::gauss
subroutine gauss(mesh)
Definition: associate_gauss.f90:24
fem_s_m::cv_00_s_m
subroutine cv_00_s_m(mesh, fs, ia, ja, a0)
Definition: fem_sparsekit.f90:1842
fem_s_m::curl_grad_2d_m
subroutine curl_grad_2d_m(mesh, alpha, stab, exp_sta, ia, ja, a0)
Definition: fem_sparsekit.f90:2530
fem_s_m::qs_1_sgs_1_m
subroutine qs_1_sgs_1_m(mesh, ff, ia, ja, a0)
Definition: fem_sparsekit.f90:1193
def_type_mesh
Definition: def_type_mesh.f90:4
fem_s_m::cc_101_m
subroutine cc_101_m(mesh, gg, ia, ja, a0)
Definition: fem_sparsekit.f90:1769
gauss_points::dw
real(kind=8), dimension(:,:,:,:), pointer dw
Definition: associate_gauss.f90:12
gauss_points::ww
real(kind=8), dimension(:,:), pointer ww
Definition: associate_gauss.f90:10
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:2639
gauss_points::dw_s
real(kind=8), dimension(:,:,:,:), pointer dw_s
Definition: associate_gauss.f90:16
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:570
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:2141
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:272