From 07fff1974aec7069fb8d9c1c28d60b92db5c5861 Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Tue, 11 May 2021 07:56:31 +0200
Subject: [PATCH] Add and use
 internal::FEPointEvaluation::is_fast_path_supported

---
 include/deal.II/fe/fe_point_evaluation.h |   6 +-
 source/fe/CMakeLists.txt                 |   2 +
 source/fe/fe_point_evaluation.cc         | 101 +++++++++++++++++++++++
 source/fe/fe_point_evaluation.inst.in    |  24 ++++++
 4 files changed, 132 insertions(+), 1 deletion(-)
 create mode 100644 source/fe/fe_point_evaluation.cc
 create mode 100644 source/fe/fe_point_evaluation.inst.in

diff --git a/include/deal.II/fe/fe_point_evaluation.h b/include/deal.II/fe/fe_point_evaluation.h
index 53d87bad3d..2ee5b1165f 100644
--- a/include/deal.II/fe/fe_point_evaluation.h
+++ b/include/deal.II/fe/fe_point_evaluation.h
@@ -344,6 +344,10 @@ namespace internal
         return value;
       }
     };
+
+    template <int dim, int spacedim>
+    bool
+    is_fast_path_supported(const FiniteElement<dim, spacedim> &fe);
   } // namespace FEPointEvaluation
 } // namespace internal
 
@@ -558,7 +562,7 @@ FEPointEvaluation<n_components, dim, spacedim>::FEPointEvaluation(
   , fe(&fe)
 {
   if (mapping_q_generic != nullptr &&
-      internal::MatrixFreeFunctions::ShapeInfo<double>::is_supported(fe))
+      internal::FEPointEvaluation::is_fast_path_supported(fe))
     {
       internal::MatrixFreeFunctions::ShapeInfo<double> shape_info;
       unsigned int                                     base_element_number = 0;
diff --git a/source/fe/CMakeLists.txt b/source/fe/CMakeLists.txt
index 7cbd87bb87..9a75812e3b 100644
--- a/source/fe/CMakeLists.txt
+++ b/source/fe/CMakeLists.txt
@@ -34,6 +34,7 @@ SET(_unity_include_src
   fe_nedelec.cc
   fe_nedelec_sz.cc
   fe_nothing.cc
+  fe_point_evaluation.cc
   fe_poly.cc
   fe_poly_tensor.cc
   fe_pyramid_p.cc
@@ -107,6 +108,7 @@ SET(_inst
   fe_nedelec.inst.in
   fe_nedelec_sz.inst.in
   fe_nothing.inst.in
+  fe_point_evaluation.inst.in
   fe_poly.inst.in
   fe_poly_tensor.inst.in
   fe_pyramid_p.inst.in
diff --git a/source/fe/fe_point_evaluation.cc b/source/fe/fe_point_evaluation.cc
new file mode 100644
index 0000000000..b566432131
--- /dev/null
+++ b/source/fe/fe_point_evaluation.cc
@@ -0,0 +1,101 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/polynomials_piecewise.h>
+
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_point_evaluation.h>
+#include <deal.II/fe/fe_poly.h>
+#include <deal.II/fe/fe_q_dg0.h>
+
+#include <deal.II/matrix_free/shape_info.h>
+
+#include <memory>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace internal
+{
+  namespace FEPointEvaluation
+  {
+    /**
+     * Similar to internal::MatrixFreeFunctions::ShapeInfo but not supporting
+     * non-tensor product elements and FE_DGQHermite
+     */
+    template <int dim, int spacedim>
+    bool
+    is_fast_path_supported(const FiniteElement<dim, spacedim> &fe)
+    {
+      // check if supported
+      const bool flag = [&]() {
+        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 &&
+                    dynamic_cast<const FE_DGQHermite<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;
+      }();
+
+      // make sure that if supported also ShapeInfo is supporting it
+      if (flag)
+        Assert(internal::MatrixFreeFunctions::ShapeInfo<double>::is_supported(
+                 fe),
+               ExcInternalError());
+
+      return flag;
+    }
+  } // namespace FEPointEvaluation
+} // namespace internal
+
+
+// explicit instantiations
+#include "fe_point_evaluation.inst"
+
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/fe/fe_point_evaluation.inst.in b/source/fe/fe_point_evaluation.inst.in
new file mode 100644
index 0000000000..4f6f7d8095
--- /dev/null
+++ b/source/fe/fe_point_evaluation.inst.in
@@ -0,0 +1,24 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+for (deal_II_dimension, deal_II_space_dimension : DIMENSIONS)
+  {
+#if deal_II_dimension <= deal_II_space_dimension
+    template bool internal::FEPointEvaluation::is_fast_path_supported(
+      const FiniteElement<deal_II_dimension, deal_II_space_dimension> &fe);
+#endif
+  }
-- 
2.39.5