add n_q_points_1d template argument to MatrixFreeOperators::MassOperator

the rationale is support different usage cases. Say, when projecting
a function one would use fe_degree+2. However when projecting the
MatrixFree quadrature data given by Table<>, most often this data
corresponds to fe_degree+1. So while being given a MatrixFree data
object, we use it even though the mass matrix could be slightly
under-integrated.

diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h
index 3b178053382e6d285aa9bc8c919df516aac39d41..4000fb34675d74a2c25eccb142e9321defa7d089 100644 (file)
--- a/include/deal.II/matrix_free/operators.h
+++ b/include/deal.II/matrix_free/operators.h
@@ -301,7 +301,7 @@ namespace MatrixFreeOperators
    *
    * @author Daniel Arndt, 2016
    */
-  template <int dim, int fe_degree, int n_components = 1, typename Number = double>
+  template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1, int n_components = 1, typename Number = double>
   class MassOperator : public Base<dim, Number>
   {
   public:
@@ -901,8 +901,8 @@ namespace MatrixFreeOperators
 
   //-----------------------------MassOperator----------------------------------
 
-  template <int dim, int fe_degree, int n_components, typename Number>
-  MassOperator<dim, fe_degree, n_components, Number>::
+  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
+  MassOperator<dim, fe_degree, n_q_points_1d, n_components, Number>::
   MassOperator ()
     :
     Base<dim, Number>()
@@ -910,9 +910,9 @@ namespace MatrixFreeOperators
 
 
 
-  template <int dim, int fe_degree, int n_components, typename Number>
+  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
   void
-  MassOperator<dim, fe_degree, n_components, Number>::
+  MassOperator<dim, fe_degree, n_q_points_1d, n_components, Number>::
   compute_diagonal()
   {
     Assert((Base<dim, Number>::data != NULL), ExcNotInitialized());
@@ -938,9 +938,9 @@ namespace MatrixFreeOperators
 
 
 
-  template <int dim, int fe_degree, int n_components, typename Number>
+  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
   void
-  MassOperator<dim, fe_degree, n_components, Number>::
+  MassOperator<dim, fe_degree, n_q_points_1d, n_components, Number>::
   apply_add (LinearAlgebra::distributed::Vector<Number>       &dst,
              const LinearAlgebra::distributed::Vector<Number> &src) const
   {
@@ -950,15 +950,15 @@ namespace MatrixFreeOperators
 
 
 
-  template <int dim, int fe_degree, int n_components, typename Number>
+  template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
   void
-  MassOperator<dim, fe_degree, n_components, Number>::
+  MassOperator<dim, fe_degree, n_q_points_1d, n_components, Number>::
   local_apply_cell (const MatrixFree<dim,Number>                     &data,
                     LinearAlgebra::distributed::Vector<Number>       &dst,
                     const LinearAlgebra::distributed::Vector<Number> &src,
                     const std::pair<unsigned int,unsigned int>  &cell_range) const
   {
-    FEEvaluation<dim, fe_degree, fe_degree+1, n_components, Number> phi(data);
+    FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number> phi(data);
     for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
       {
         phi.reinit (cell);

diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h
index 1decd4bbfd413168972b69ce42352c2d20050ba3..e277a25b89a0734315b449f6069867d2a5c0db42 100644 (file)
--- a/include/deal.II/numerics/vector_tools.templates.h
+++ b/include/deal.II/numerics/vector_tools.templates.h
@@ -917,8 +917,8 @@ namespace VectorTools
       additional_data.mapping_update_flags = (update_values | update_JxW_values);
       MatrixFree<dim, Number> matrix_free;
       matrix_free.reinit (mapping, dof, constraints,
-                          QGauss<1>(fe_degree+1), additional_data);
-      typedef MatrixFreeOperators::MassOperator<dim, fe_degree, 1, Number> MatrixType;
+                          QGauss<1>(fe_degree+2), additional_data);
+      typedef MatrixFreeOperators::MassOperator<dim, fe_degree, fe_degree+2, 1, Number> MatrixType;
       MatrixType mass_matrix;
       mass_matrix.initialize(matrix_free);
       mass_matrix.compute_diagonal();
@@ -1004,7 +1004,7 @@ namespace VectorTools
       Assert (dof.get_fe().degree == fe_degree,
               ExcDimensionMismatch(fe_degree, dof.get_fe().degree));
 
-      typedef MatrixFreeOperators::MassOperator<dim, fe_degree, 1, Number> MatrixType;
+      typedef MatrixFreeOperators::MassOperator<dim, fe_degree, n_q_points_1d, 1, Number> MatrixType;
       MatrixType mass_matrix;
       mass_matrix.initialize(matrix_free);
       mass_matrix.compute_diagonal();

diff --git a/tests/matrix_free/mass_operator_01.cc b/tests/matrix_free/mass_operator_01.cc
index 575b0bc596eb667aedee3eafa2890069e28692a8..2e95bd51a95e85a78119ff1ebb8d33a407150e60 100644 (file)
--- a/tests/matrix_free/mass_operator_01.cc
+++ b/tests/matrix_free/mass_operator_01.cc
@@ -116,7 +116,7 @@ void test ()
     mf_data.reinit (dof, constraints, quad, data);
   }
 
-  MatrixFreeOperators::MassOperator<dim,fe_degree, 1, number> mf;
+  MatrixFreeOperators::MassOperator<dim,fe_degree, fe_degree+1, 1, number> mf;
   mf.initialize(mf_data);
   mf.compute_diagonal();
   LinearAlgebra::distributed::Vector<number> in, out, ref;

diff --git a/tests/matrix_free/mass_operator_02.cc b/tests/matrix_free/mass_operator_02.cc
index 4bcc628f4e67c6d0660de92671fd35f5d2d5931a..b7db2771e1f1f3401e23d0eb6daf2725a2258155 100644 (file)
--- a/tests/matrix_free/mass_operator_02.cc
+++ b/tests/matrix_free/mass_operator_02.cc
@@ -80,7 +80,7 @@ void test ()
     mf_data.reinit (dof, constraints, quad, data);
   }
 
-  MatrixFreeOperators::MassOperator<dim,fe_degree, 1, number> mf;
+  MatrixFreeOperators::MassOperator<dim,fe_degree, fe_degree+1, 1, number> mf;
   mf.initialize(mf_data);
   mf.compute_diagonal();
   const LinearAlgebra::distributed::Vector<double> &diagonal