From c50dd3feb3b44f94bcd1297db0d3e5cbb7a48590 Mon Sep 17 00:00:00 2001
From: kanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sun, 2 Feb 2014 22:41:39 +0000
Subject: [PATCH] add weighted mass matrix

git-svn-id: https://svn.dealii.org/trunk@32381 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/include/deal.II/integrators/l2.h | 42 ++++++++++++++++++++++++
 1 file changed, 42 insertions(+)

diff --git a/deal.II/include/deal.II/integrators/l2.h b/deal.II/include/deal.II/integrators/l2.h
index 2313997c04..c3004be3e0 100644
--- a/deal.II/include/deal.II/integrators/l2.h
+++ b/deal.II/include/deal.II/integrators/l2.h
@@ -75,6 +75,48 @@ namespace LocalIntegrators
         }
     }
 
+    /**
+     * The weighted mass matrix for scalar or vector values finite elements.
+     * \f[
+     * \int_Z \omega(x) uv\,dx \quad \text{or} \quad \int_Z \omega(x) \mathbf u\cdot \mathbf v\,dx
+     * \f]
+     *
+     * Likewise, this term can be used on faces, where it computes  the integrals
+     * \f[
+     * \int_F \omega(x) uv\,ds \quad \text{or} \quad \int_F \omega(x) \mathbf u\cdot \mathbf v\,ds
+     * \f]
+     *
+     * The size of the vector <tt>weights</tt> must be equal to the
+     * number of quadrature points in the finite element.
+     *
+     * @author Guido Kanschat
+     * @date 2014
+     */
+    template <int dim>
+    void weighted_mass_matrix (
+      FullMatrix<double> &M,
+      const FEValuesBase<dim> &fe,
+      const std::vector<double> weights)
+    {
+      const unsigned int n_dofs = fe.dofs_per_cell;
+      const unsigned int n_components = fe.get_fe().n_components();
+      AssertDimension(M.m(), n_dofs);
+      AssertDimension(M.n(), n_dofs);
+      AssertDimension(weights.size(), fe.n_quadrature_points);
+      
+      for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
+        {
+          const double dx = fe.JxW(k) * weights[k];
+
+          for (unsigned int i=0; i<n_dofs; ++i)
+            for (unsigned int j=0; j<n_dofs; ++j)
+              for (unsigned int d=0; d<n_components; ++d)
+                M(i,j) += dx
+                          * fe.shape_value_component(j,k,d)
+                          * fe.shape_value_component(i,k,d);
+        }
+    }
+
     /**
      * <i>L<sup>2</sup></i>-inner product for scalar functions.
      *
-- 
2.39.5