From: Martin Kronbichler <martin.kronbichler@it.uu.se>
Date: Tue, 30 Nov 2021 09:43:11 +0000 (+0100)
Subject: Move implementation for ShapeInfo to templates.h file
X-Git-Tag: v9.4.0-rc1~781^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=634d503791869fb2cd8ae89ce6d37a4aee7b59eb;p=dealii.git

Move implementation for ShapeInfo to templates.h file
---

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index 75e7f55842..4fc5c6cc2a 100644
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -232,7 +232,6 @@ public:
   std::array<unsigned int, VectorizedArrayType::size()>
   get_cell_or_face_ids() const;
 
-
   /**
    * Return the numbering of local degrees of freedom within the evaluation
    * routines of FEEvaluation in terms of the standard numbering on finite
diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h
index 59fe6a9bc1..241cbe517f 100644
--- a/include/deal.II/matrix_free/shape_info.h
+++ b/include/deal.II/matrix_free/shape_info.h
@@ -22,14 +22,19 @@
 
 #include <deal.II/base/aligned_vector.h>
 #include <deal.II/base/exceptions.h>
-#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/table.h>
 #include <deal.II/base/vectorization.h>
 
-#include <deal.II/fe/fe.h>
-
 
 DEAL_II_NAMESPACE_OPEN
 
+
+// forward declaration
+template <int dim, int spacedim>
+class FiniteElement;
+
+
 namespace internal
 {
   namespace MatrixFreeFunctions
@@ -132,7 +137,7 @@ namespace internal
       ElementType element_type;
 
       /**
-       * Stores the shape values of the 1D finite element evaluated on all 1D
+       * Stores the shape values of the 1D finite element evaluated at all 1D
        * quadrature points. The length of
        * this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature
        * points are the index running fastest.
@@ -140,7 +145,7 @@ namespace internal
       AlignedVector<Number> shape_values;
 
       /**
-       * Stores the shape gradients of the 1D finite element evaluated on all
+       * Stores the shape gradients of the 1D finite element evaluated at all
        * 1D quadrature points. The length of
        * this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature
        * points are the index running fastest.
@@ -148,7 +153,7 @@ namespace internal
       AlignedVector<Number> shape_gradients;
 
       /**
-       * Stores the shape Hessians of the 1D finite element evaluated on all
+       * Stores the shape Hessians of the 1D finite element evaluated at all
        * 1D quadrature points. The length of
        * this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature
        * points are the index running fastest.
@@ -300,16 +305,16 @@ namespace internal
       bool nodal_at_cell_boundaries;
 
       /**
-       * Stores the shape values of the finite element evaluated on all
+       * Stores the shape values of the finite element evaluated at all
        * quadrature points for all faces and orientations (no tensor-product
        * structure exploited).
        */
       Table<3, Number> shape_values_face;
 
       /**
-       * Stores the shape gradients of the finite element evaluated on all
+       * Stores the shape gradients of the finite element evaluated at all
        * quadrature points for all faces, orientations, and directions
-       * (no tensor-product structure  exploited).
+       * (no tensor-product structure exploited).
        */
       Table<4, Number> shape_gradients_face;
     };
@@ -344,10 +349,10 @@ namespace internal
       /**
        * Constructor that initializes the data fields using the reinit method.
        */
-      template <int dim, int dim_q>
-      ShapeInfo(const Quadrature<dim_q> & quad,
-                const FiniteElement<dim> &fe,
-                const unsigned int        base_element = 0);
+      template <int dim, int spacedim, int dim_q>
+      ShapeInfo(const Quadrature<dim_q> &           quad,
+                const FiniteElement<dim, spacedim> &fe,
+                const unsigned int                  base_element = 0);
 
       /**
        * Initializes the data fields. Takes a one-dimensional quadrature
@@ -357,11 +362,11 @@ namespace internal
        * dimensional element by a tensor product and that the zeroth shape
        * function in zero evaluates to one.
        */
-      template <int dim, int dim_q>
+      template <int dim, int spacedim, int dim_q>
       void
-      reinit(const Quadrature<dim_q> & quad,
-             const FiniteElement<dim> &fe_dim,
-             const unsigned int        base_element = 0);
+      reinit(const Quadrature<dim_q> &           quad,
+             const FiniteElement<dim, spacedim> &fe_dim,
+             const unsigned int                  base_element = 0);
 
       /**
        * Return which kinds of elements are supported by MatrixFree.
@@ -556,22 +561,6 @@ namespace internal
 
     // ------------------------------------------ inline functions
 
-    template <typename Number>
-    template <int dim, int dim_q>
-    inline ShapeInfo<Number>::ShapeInfo(const Quadrature<dim_q> & quad,
-                                        const FiniteElement<dim> &fe_in,
-                                        const unsigned int base_element_number)
-      : element_type(tensor_general)
-      , n_dimensions(0)
-      , n_components(0)
-      , n_q_points(0)
-      , dofs_per_component_on_cell(0)
-      , n_q_points_face(0)
-      , dofs_per_component_on_face(0)
-    {
-      reinit(quad, fe_in, base_element_number);
-    }
-
     template <typename Number>
     inline const UnivariateShapeData<Number> &
     ShapeInfo<Number>::get_shape_data(const unsigned int dimension,
diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index c8a482d752..b982760224 100644
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -26,6 +26,7 @@
 #include <deal.II/base/tensor_product_polynomials.h>
 #include <deal.II/base/utilities.h>
 
+#include <deal.II/fe/fe.h>
 #include <deal.II/fe/fe_dgp.h>
 #include <deal.II/fe/fe_dgq.h>
 #include <deal.II/fe/fe_poly.h>
@@ -62,8 +63,6 @@ namespace internal
 
 
 
-    // ----------------- actual ShapeInfo functions --------------------
-
     template <typename Number>
     Number
     get_first_array_element(const Number a)
@@ -196,6 +195,8 @@ namespace internal
 
 
 
+    // ----------------- actual ShapeInfo implementation --------------------
+
     template <typename Number>
     ShapeInfo<Number>::ShapeInfo()
       : element_type(tensor_general)
@@ -210,62 +211,33 @@ namespace internal
 
 
     template <typename Number>
-    template <int dim, int spacedim>
-    bool
-    ShapeInfo<Number>::is_supported(const FiniteElement<dim, spacedim> &fe)
+    template <int dim, int spacedim, int dim_q>
+    inline ShapeInfo<Number>::ShapeInfo(
+      const Quadrature<dim_q> &           quad,
+      const FiniteElement<dim, spacedim> &fe_in,
+      const unsigned int                  base_element_number)
+      : element_type(tensor_general)
+      , n_dimensions(0)
+      , n_components(0)
+      , n_q_points(0)
+      , dofs_per_component_on_cell(0)
+      , n_q_points_face(0)
+      , dofs_per_component_on_face(0)
     {
-      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);
-              // Simplices are a special case since the polynomial family is not
-              // indicative of their support
-              if (dynamic_cast<const FE_SimplexP<dim> *>(fe_poly_ptr) ||
-                  dynamic_cast<const FE_SimplexDGP<dim> *>(fe_poly_ptr) ||
-                  dynamic_cast<const FE_WedgeP<dim> *>(fe_poly_ptr) ||
-                  dynamic_cast<const FE_PyramidP<dim> *>(fe_poly_ptr))
-                return true;
-
-              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;
+      reinit(quad, fe_in, base_element_number);
     }
 
 
 
     template <typename Number>
-    template <int dim, int dim_q>
+    template <int dim, int spacedim, int dim_q>
     void
-    ShapeInfo<Number>::reinit(const Quadrature<dim_q> & quad_in,
-                              const FiniteElement<dim> &fe_in,
-                              const unsigned int        base_element_number)
+    ShapeInfo<Number>::reinit(const Quadrature<dim_q> &           quad_in,
+                              const FiniteElement<dim, spacedim> &fe_in,
+                              const unsigned int base_element_number)
     {
+      static_assert(dim == spacedim,
+                    "Currently, only the case dim=spacedim is implemented");
       if (quad_in.is_tensor_product() == false ||
           dynamic_cast<const FE_SimplexP<dim> *>(
             &fe_in.base_element(base_element_number)) ||
@@ -1135,6 +1107,56 @@ 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);
+              // Simplices are a special case since the polynomial family is not
+              // indicative of their support
+              if (dynamic_cast<const FE_SimplexP<dim> *>(fe_poly_ptr) ||
+                  dynamic_cast<const FE_SimplexDGP<dim> *>(fe_poly_ptr) ||
+                  dynamic_cast<const FE_WedgeP<dim> *>(fe_poly_ptr) ||
+                  dynamic_cast<const FE_PyramidP<dim> *>(fe_poly_ptr))
+                return true;
+
+              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>
     std::size_t
     ShapeInfo<Number>::memory_consumption() const
diff --git a/source/matrix_free/shape_info.inst.in b/source/matrix_free/shape_info.inst.in
index 30be5ee7aa..fe5de6c930 100644
--- a/source/matrix_free/shape_info.inst.in
+++ b/source/matrix_free/shape_info.inst.in
@@ -16,18 +16,28 @@
 
 for (deal_II_dimension : DIMENSIONS; deal_II_scalar : REAL_SCALARS)
   {
-    template void internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar>::
-      reinit<deal_II_dimension>(
-        const Quadrature<deal_II_dimension> &,
-        const FiniteElement<deal_II_dimension, deal_II_dimension> &,
-        const unsigned int);
+    template internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar>::
+      ShapeInfo(const Quadrature<deal_II_dimension> &,
+                const FiniteElement<deal_II_dimension, deal_II_dimension> &,
+                const unsigned int);
+
+    template void
+    internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar>::reinit(
+      const Quadrature<deal_II_dimension> &,
+      const FiniteElement<deal_II_dimension, deal_II_dimension> &,
+      const unsigned int);
 
 #if deal_II_dimension > 1
-    template void internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar>::
-      reinit<deal_II_dimension>(
-        const Quadrature<1> &,
-        const FiniteElement<deal_II_dimension, deal_II_dimension> &,
-        const unsigned int);
+    template internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar>::
+      ShapeInfo(const Quadrature<1> &,
+                const FiniteElement<deal_II_dimension, deal_II_dimension> &,
+                const unsigned int);
+
+    template void
+    internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar>::reinit(
+      const Quadrature<1> &,
+      const FiniteElement<deal_II_dimension, deal_II_dimension> &,
+      const unsigned int);
 #endif
   }
 
@@ -45,18 +55,30 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
 for (deal_II_dimension : DIMENSIONS;
      deal_II_scalar_vectorized : REAL_SCALARS_VECTORIZED)
   {
-    template void internal::MatrixFreeFunctions::
-      ShapeInfo<deal_II_scalar_vectorized>::reinit<deal_II_dimension>(
+    template internal::MatrixFreeFunctions::
+      ShapeInfo<deal_II_scalar_vectorized>::ShapeInfo(
         const Quadrature<deal_II_dimension> &,
         const FiniteElement<deal_II_dimension, deal_II_dimension> &,
         const unsigned int);
 
+    template void
+    internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar_vectorized>::reinit(
+      const Quadrature<deal_II_dimension> &,
+      const FiniteElement<deal_II_dimension, deal_II_dimension> &,
+      const unsigned int);
+
 #if deal_II_dimension > 1
-    template void internal::MatrixFreeFunctions::
-      ShapeInfo<deal_II_scalar_vectorized>::reinit<deal_II_dimension>(
+    template internal::MatrixFreeFunctions::
+      ShapeInfo<deal_II_scalar_vectorized>::ShapeInfo(
         const Quadrature<1> &,
         const FiniteElement<deal_II_dimension, deal_II_dimension> &,
         const unsigned int);
+
+    template void
+    internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar_vectorized>::reinit(
+      const Quadrature<1> &,
+      const FiniteElement<deal_II_dimension, deal_II_dimension> &,
+      const unsigned int);
 #endif
   }