From 1b7744c2a161d8bd5af3c4e5465e39b81ff3063f Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Thu, 22 Oct 2020 18:31:13 +0200
Subject: [PATCH] Improve MatrixFree::is_supported()

---
 .../matrix_free/matrix_free.templates.h       | 33 +------------
 include/deal.II/matrix_free/shape_info.h      |  7 +++
 .../matrix_free/shape_info.templates.h        | 46 +++++++++++++++++++
 .../numerics/vector_tools_project.templates.h |  3 +-
 source/matrix_free/shape_info.inst.in         |  9 ++++
 5 files changed, 65 insertions(+), 33 deletions(-)

diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h
index b6152e0076..64c08905f0 100644
--- a/include/deal.II/matrix_free/matrix_free.templates.h
+++ b/include/deal.II/matrix_free/matrix_free.templates.h
@@ -464,38 +464,7 @@ bool
 MatrixFree<dim, Number, VectorizedArrayType>::is_supported(
   const FiniteElement<dim, spacedim> &fe)
 {
-  if (dim != spacedim)
-    return false;
-
-  // first check for degree, number of base_elemnt and number of its components
-  if (fe.degree == 0 || fe.n_base_elements() != 1)
-    return false;
-
-  const FiniteElement<dim, spacedim> *fe_ptr = &(fe.base_element(0));
-  if (fe_ptr->n_components() != 1)
-    return false;
-
-  // then check of the base element is supported
-  if (dynamic_cast<const FE_Poly<dim, spacedim> *>(fe_ptr) != nullptr)
-    {
-      const FE_Poly<dim, spacedim> *fe_poly_ptr =
-        dynamic_cast<const FE_Poly<dim, spacedim> *>(fe_ptr);
-      if (dynamic_cast<const TensorProductPolynomials<dim> *>(
-            &fe_poly_ptr->get_poly_space()) != nullptr)
-        return true;
-      if (dynamic_cast<const TensorProductPolynomials<
-            dim,
-            Polynomials::PiecewisePolynomial<double>> *>(
-            &fe_poly_ptr->get_poly_space()) != nullptr)
-        return true;
-    }
-  if (dynamic_cast<const FE_DGP<dim, spacedim> *>(fe_ptr) != nullptr)
-    return true;
-  if (dynamic_cast<const FE_Q_DG0<dim, spacedim> *>(fe_ptr) != nullptr)
-    return true;
-
-  // if the base element is not in the above list it is not supported
-  return false;
+  return internal::MatrixFreeFunctions::ShapeInfo<double>::is_supported(fe);
 }
 
 
diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h
index 0637a5e775..64daf31106 100644
--- a/include/deal.II/matrix_free/shape_info.h
+++ b/include/deal.II/matrix_free/shape_info.h
@@ -331,6 +331,13 @@ namespace internal
              const FiniteElement<dim> &fe_dim,
              const unsigned int        base_element = 0);
 
+      /**
+       * Return which kinds of elements are supported by MatrixFree.
+       */
+      template <int dim, int spacedim>
+      static bool
+      is_supported(const FiniteElement<dim, spacedim> &fe);
+
       /**
        * Return data of univariate shape functions which defines the
        * dimension @p dimension of tensor product shape functions
diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index 41a90beb54..f41881d869 100644
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -61,6 +61,8 @@ namespace internal
       return a;
     }
 
+
+
     template <typename Number, std::size_t width>
     Number
     get_first_array_element(const VectorizedArray<Number, width> a)
@@ -68,6 +70,8 @@ namespace internal
       return a[0];
     }
 
+
+
     template <typename Number>
     ShapeInfo<Number>::ShapeInfo()
       : element_type(tensor_general)
@@ -81,6 +85,48 @@ namespace internal
 
 
 
+    template <typename Number>
+    template <int dim, int spacedim>
+    bool
+    ShapeInfo<Number>::is_supported(const FiniteElement<dim, spacedim> &fe)
+    {
+      if (dim != spacedim)
+        return false;
+
+      for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
+        {
+          const FiniteElement<dim, spacedim> *fe_ptr = &(fe.base_element(base));
+          if (fe_ptr->n_components() != 1)
+            return false;
+
+          // then check if the base element is supported or not
+          if (dynamic_cast<const FE_Poly<dim, spacedim> *>(fe_ptr) != nullptr)
+            {
+              const FE_Poly<dim, spacedim> *fe_poly_ptr =
+                dynamic_cast<const FE_Poly<dim, spacedim> *>(fe_ptr);
+              if (dynamic_cast<const TensorProductPolynomials<dim> *>(
+                    &fe_poly_ptr->get_poly_space()) == nullptr &&
+                  dynamic_cast<const TensorProductPolynomials<
+                      dim,
+                      Polynomials::PiecewisePolynomial<double>> *>(
+                    &fe_poly_ptr->get_poly_space()) == nullptr &&
+                  dynamic_cast<const FE_DGP<dim, spacedim> *>(fe_ptr) ==
+                    nullptr &&
+                  dynamic_cast<const FE_Q_DG0<dim, spacedim> *>(fe_ptr) ==
+                    nullptr)
+                return false;
+            }
+          else
+            return false;
+        }
+
+      // if we arrived here, all base elements were supported so we can
+      // support the present element
+      return true;
+    }
+
+
+
     template <typename Number>
     template <int dim>
     void
diff --git a/include/deal.II/numerics/vector_tools_project.templates.h b/include/deal.II/numerics/vector_tools_project.templates.h
index 1988f3decd..69083a0bdc 100644
--- a/include/deal.II/numerics/vector_tools_project.templates.h
+++ b/include/deal.II/numerics/vector_tools_project.templates.h
@@ -706,7 +706,8 @@ namespace VectorTools
       // If we can, use the matrix-free implementation
       bool use_matrix_free =
         MatrixFree<dim, typename VectorType::value_type>::is_supported(
-          dof.get_fe());
+          dof.get_fe()) &&
+        dof.get_fe().n_base_elements() == 1;
 
       // enforce_zero_boundary and project_to_boundary_first
       // are not yet supported.
diff --git a/source/matrix_free/shape_info.inst.in b/source/matrix_free/shape_info.inst.in
index f2c859da1d..293099d241 100644
--- a/source/matrix_free/shape_info.inst.in
+++ b/source/matrix_free/shape_info.inst.in
@@ -23,6 +23,15 @@ for (deal_II_dimension : DIMENSIONS; deal_II_scalar : REAL_SCALARS)
         const unsigned int);
   }
 
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
+  {
+#if deal_II_dimension <= deal_II_space_dimension
+    template bool
+    internal::MatrixFreeFunctions::ShapeInfo<double>::is_supported(
+      const FiniteElement<deal_II_dimension, deal_II_space_dimension> &);
+#endif
+  }
+
 
 
 for (deal_II_dimension : DIMENSIONS;
-- 
2.39.5