FEPointEvaluation: Work towards dim < spacedim

diff --git a/include/deal.II/base/derivative_form.h b/include/deal.II/base/derivative_form.h
index 624f1781b19ab1dfd8227b88bb56175c970fa3ce..455bc111888d36e880a1c46451f4c611d96e7993 100644 (file)
--- a/include/deal.II/base/derivative_form.h
+++ b/include/deal.II/base/derivative_form.h
@@ -609,9 +609,7 @@ template <int dim, int spacedim, typename Number>
 inline DerivativeForm<1, spacedim, dim, Number>
 transpose(const DerivativeForm<1, dim, spacedim, Number> &DF)
 {
-  DerivativeForm<1, spacedim, dim, Number> tt;
-  tt = DF.transpose();
-  return tt;
+  return DF.transpose();
 }
 
 
@@ -691,6 +689,34 @@ apply_diagonal_transformation(
 }
 
 
+
+/**
+ * Similar to the previous apply_transformation().
+ * Each row of the result corresponds to one of the rows of @p D_X transformed
+ * by @p grad_F, equivalent to $\mathrm{D\_X} \, \mathrm{grad\_F}^T$ in matrix
+ * notation.
+ *
+ * @relatesalso DerivativeForm
+ */
+// rank=2
+template <int spacedim, int dim, typename Number1, typename Number2>
+inline DerivativeForm<1,
+                      spacedim,
+                      dim,
+                      typename ProductType<Number1, Number2>::type>
+apply_diagonal_transformation(
+  const DerivativeForm<1, dim, spacedim, Number1> &grad_F,
+  const Tensor<2, dim, Number2>                   &D_X)
+{
+  DerivativeForm<1, spacedim, dim, typename ProductType<Number1, Number2>::type>
+    dest;
+  for (unsigned int i = 0; i < dim; ++i)
+    dest[i] = apply_diagonal_transformation(grad_F, D_X[i]);
+
+  return dest;
+}
+
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif

diff --git a/include/deal.II/base/floating_point_comparator.h b/include/deal.II/base/floating_point_comparator.h
index 9a1d8774570c93a503d70bd933b3b67e60af8a88..c3a468784ac26dcd746dd372e980b4be090f08d0 100644 (file)
--- a/include/deal.II/base/floating_point_comparator.h
+++ b/include/deal.II/base/floating_point_comparator.h
@@ -232,7 +232,7 @@ FloatingPointComparator<Number>::compare(
   const DerivativeForm<rank, dim, spacedim, T> &t1,
   const DerivativeForm<rank, dim, spacedim, T> &t2) const
 {
-  for (unsigned int i = 0; i < dim; ++i)
+  for (unsigned int i = 0; i < spacedim; ++i)
     {
       const ComparisonResult result = compare(t1[i], t2[i]);
       if (result != ComparisonResult::equal)

diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h
index ca0c20726a9dfc14121021b250aee2e7837e06ca..357d1922b5b69eded50c86ab36c5337c0218f556 100644 (file)
--- a/include/deal.II/matrix_free/fe_point_evaluation.h
+++ b/include/deal.II/matrix_free/fe_point_evaluation.h
@@ -56,7 +56,11 @@ namespace internal
      * Struct to distinguish between the value and gradient types of different
      * numbers of components used by the FlexibleEvaluator class.
      */
-    template <int dim, int n_components, typename Number>
+    template <int dim,
+              int spacedim,
+              int n_components,
+              typename Number,
+              typename Enable = void>
     struct EvaluatorTypeTraits
     {
       using ScalarNumber =
@@ -69,6 +73,10 @@ namespace internal
       using vectorized_value_type =
         Tensor<1, n_components, VectorizedArrayType>;
       using gradient_type = Tensor<1, n_components, Tensor<1, dim, Number>>;
+      using real_gradient_type = std::conditional_t<
+        n_components == spacedim,
+        Tensor<2, spacedim, Number>,
+        Tensor<1, n_components, Tensor<1, spacedim, Number>>>;
       using scalar_gradient_type =
         Tensor<1, n_components, Tensor<1, dim, ScalarNumber>>;
       using vectorized_gradient_type =
@@ -134,10 +142,22 @@ namespace internal
       }
 
       static void
-      set_zero_gradient(gradient_type &value, const unsigned int vector_lane)
+      get_gradient(interface_vectorized_gradient_type &value,
+                   const unsigned int                  vector_lane,
+                   const DerivativeForm<1, dim, n_components, Number> &result)
       {
         for (unsigned int i = 0; i < n_components; ++i)
           for (unsigned int d = 0; d < dim; ++d)
+            internal::VectorizedArrayTrait<Number>::get_from_vectorized(
+              value[d][i], vector_lane) = result[i][d];
+      }
+
+      static void
+      set_zero_gradient(real_gradient_type &value,
+                        const unsigned int  vector_lane)
+      {
+        for (unsigned int i = 0; i < n_components; ++i)
+          for (unsigned int d = 0; d < spacedim; ++d)
             internal::VectorizedArrayTrait<Number>::get(value[i][d],
                                                         vector_lane) = 0.;
       }
@@ -204,24 +224,24 @@ namespace internal
       }
 
       static void
-      access(gradient_type                      &value,
-             const unsigned int                  vector_lane,
-             const unsigned int                  component,
-             const Tensor<1, dim, ScalarNumber> &shape_gradient)
+      access(real_gradient_type                      &value,
+             const unsigned int                       vector_lane,
+             const unsigned int                       component,
+             const Tensor<1, spacedim, ScalarNumber> &shape_gradient)
       {
-        for (unsigned int d = 0; d < dim; ++d)
+        for (unsigned int d = 0; d < spacedim; ++d)
           internal::VectorizedArrayTrait<Number>::get(value[component][d],
                                                       vector_lane) +=
             shape_gradient[d];
       }
 
-      static Tensor<1, dim, ScalarNumber>
-      access(const gradient_type &value,
-             const unsigned int   vector_lane,
-             const unsigned int   component)
+      static Tensor<1, spacedim, ScalarNumber>
+      access(const real_gradient_type &value,
+             const unsigned int        vector_lane,
+             const unsigned int        component)
       {
-        Tensor<1, dim, ScalarNumber> result;
-        for (unsigned int d = 0; d < dim; ++d)
+        Tensor<1, spacedim, ScalarNumber> result;
+        for (unsigned int d = 0; d < spacedim; ++d)
           result[d] =
             internal::VectorizedArrayTrait<Number>::get(value[component][d],
                                                         vector_lane);
@@ -229,8 +249,8 @@ namespace internal
       }
     };
 
-    template <int dim, typename Number>
-    struct EvaluatorTypeTraits<dim, 1, Number>
+    template <int dim, int spacedim, typename Number>
+    struct EvaluatorTypeTraits<dim, spacedim, 1, Number>
     {
       using ScalarNumber =
         typename internal::VectorizedArrayTrait<Number>::value_type;
@@ -241,6 +261,7 @@ namespace internal
       using scalar_value_type        = ScalarNumber;
       using vectorized_value_type    = VectorizedArrayType;
       using gradient_type            = Tensor<1, dim, Number>;
+      using real_gradient_type       = Tensor<1, spacedim, Number>;
       using scalar_gradient_type     = Tensor<1, dim, ScalarNumber>;
       using vectorized_gradient_type = Tensor<1, dim, VectorizedArrayType>;
       using interface_vectorized_gradient_type = vectorized_gradient_type;
@@ -306,9 +327,10 @@ namespace internal
       }
 
       static void
-      set_zero_gradient(gradient_type &value, const unsigned int vector_lane)
+      set_zero_gradient(real_gradient_type &value,
+                        const unsigned int  vector_lane)
       {
-        for (unsigned int d = 0; d < dim; ++d)
+        for (unsigned int d = 0; d < spacedim; ++d)
           internal::VectorizedArrayTrait<Number>::get(value[d], vector_lane) =
             0.;
       }
@@ -370,23 +392,23 @@ namespace internal
       }
 
       static void
-      access(gradient_type     &value,
-             const unsigned int vector_lane,
+      access(real_gradient_type &value,
+             const unsigned int  vector_lane,
              const unsigned int,
-             const scalar_gradient_type &shape_gradient)
+             const Tensor<1, spacedim, ScalarNumber> &shape_gradient)
       {
-        for (unsigned int d = 0; d < dim; ++d)
+        for (unsigned int d = 0; d < spacedim; ++d)
           internal::VectorizedArrayTrait<Number>::get(value[d], vector_lane) +=
             shape_gradient[d];
       }
 
-      static scalar_gradient_type
-      access(const gradient_type &value,
-             const unsigned int   vector_lane,
+      static Tensor<1, spacedim, ScalarNumber>
+      access(const real_gradient_type &value,
+             const unsigned int        vector_lane,
              const unsigned int)
       {
-        scalar_gradient_type result;
-        for (unsigned int d = 0; d < dim; ++d)
+        Tensor<1, spacedim, ScalarNumber> result;
+        for (unsigned int d = 0; d < spacedim; ++d)
           result[d] =
             internal::VectorizedArrayTrait<Number>::get(value[d], vector_lane);
         return result;
@@ -394,7 +416,11 @@ namespace internal
     };
 
     template <int dim, typename Number>
-    struct EvaluatorTypeTraits<dim, dim, Number>
+    struct EvaluatorTypeTraits<dim,
+                               dim,
+                               dim,
+                               Number,
+                               std::enable_if_t<dim != 1>>
     {
       using ScalarNumber =
         typename internal::VectorizedArrayTrait<Number>::value_type;
@@ -405,6 +431,7 @@ namespace internal
       using scalar_value_type        = Tensor<1, dim, ScalarNumber>;
       using vectorized_value_type    = Tensor<1, dim, VectorizedArrayType>;
       using gradient_type            = Tensor<2, dim, Number>;
+      using real_gradient_type       = gradient_type;
       using scalar_gradient_type     = Tensor<2, dim, ScalarNumber>;
       using vectorized_gradient_type = Tensor<2, dim, VectorizedArrayType>;
       using interface_vectorized_gradient_type =
@@ -538,7 +565,7 @@ namespace internal
       }
 
       static void
-      access(gradient_type                      &value,
+      access(real_gradient_type                 &value,
              const unsigned int                  vector_lane,
              const unsigned int                  component,
              const Tensor<1, dim, ScalarNumber> &shape_gradient)
@@ -550,9 +577,9 @@ namespace internal
       }
 
       static Tensor<1, dim, ScalarNumber>
-      access(const gradient_type &value,
-             const unsigned int   vector_lane,
-             const unsigned int   component)
+      access(const real_gradient_type &value,
+             const unsigned int        vector_lane,
+             const unsigned int        component)
       {
         Tensor<1, dim, ScalarNumber> result;
         for (unsigned int d = 0; d < dim; ++d)
@@ -563,164 +590,6 @@ namespace internal
       }
     };
 
-    template <typename Number>
-    struct EvaluatorTypeTraits<1, 1, Number>
-    {
-      using ScalarNumber =
-        typename internal::VectorizedArrayTrait<Number>::value_type;
-      using VectorizedArrayType =
-        typename dealii::internal::VectorizedArrayTrait<
-          Number>::vectorized_value_type;
-      using value_type               = Number;
-      using scalar_value_type        = ScalarNumber;
-      using vectorized_value_type    = VectorizedArrayType;
-      using gradient_type            = Tensor<1, 1, Number>;
-      using scalar_gradient_type     = Tensor<1, 1, ScalarNumber>;
-      using vectorized_gradient_type = Tensor<1, 1, VectorizedArrayType>;
-      using interface_vectorized_gradient_type = vectorized_gradient_type;
-
-      static void
-      read_value(const ScalarNumber vector_entry,
-                 const unsigned int,
-                 scalar_value_type &result)
-      {
-        result = vector_entry;
-      }
-
-      static scalar_value_type
-      sum_value(const scalar_value_type &result)
-      {
-        return result;
-      }
-
-      static scalar_value_type
-      sum_value(const vectorized_value_type &result)
-      {
-        return result.sum();
-      }
-
-      static ScalarNumber
-      sum_value(const unsigned int, const vectorized_value_type &result)
-      {
-        return result.sum();
-      }
-
-      static void
-      set_gradient(const vectorized_gradient_type &value,
-                   const unsigned int              vector_lane,
-                   scalar_gradient_type           &result)
-      {
-        result[0] = value[0][vector_lane];
-      }
-
-      static void
-      set_gradient(const vectorized_gradient_type &value,
-                   const unsigned int,
-                   vectorized_gradient_type &result)
-      {
-        result = value;
-      }
-
-      static void
-      get_gradient(vectorized_gradient_type   &value,
-                   const unsigned int          vector_lane,
-                   const scalar_gradient_type &result)
-      {
-        value[0][vector_lane] = result[0];
-      }
-
-      static void
-      get_gradient(vectorized_gradient_type &value,
-                   const unsigned int,
-                   const vectorized_gradient_type &result)
-      {
-        value = result;
-      }
-
-      static void
-      set_zero_gradient(gradient_type &value, const unsigned int vector_lane)
-      {
-        internal::VectorizedArrayTrait<Number>::get(value[0], vector_lane) = 0.;
-      }
-
-      static void
-      set_value(const vectorized_value_type &value,
-                const unsigned int           vector_lane,
-                scalar_value_type           &result)
-      {
-        result = value[vector_lane];
-      }
-
-      static void
-      set_value(const vectorized_value_type &value,
-                const unsigned int,
-                vectorized_value_type &result)
-      {
-        result = value;
-      }
-
-      static void
-      get_value(vectorized_value_type   &value,
-                const unsigned int       vector_lane,
-                const scalar_value_type &result)
-      {
-        value[vector_lane] = result;
-      }
-
-      static void
-      get_value(vectorized_value_type &value,
-                const unsigned int,
-                const vectorized_value_type &result)
-      {
-        value = result;
-      }
-
-      static void
-      set_zero_value(value_type &value, const unsigned int vector_lane)
-      {
-        internal::VectorizedArrayTrait<Number>::get(value, vector_lane) = 0.;
-      }
-
-      static void
-      access(value_type        &value,
-             const unsigned int vector_lane,
-             const unsigned int,
-             const ScalarNumber &shape_value)
-      {
-        internal::VectorizedArrayTrait<Number>::get(value, vector_lane) +=
-          shape_value;
-      }
-
-      static ScalarNumber
-      access(const value_type  &value,
-             const unsigned int vector_lane,
-             const unsigned int)
-      {
-        return internal::VectorizedArrayTrait<Number>::get(value, vector_lane);
-      }
-
-      static void
-      access(gradient_type     &value,
-             const unsigned int vector_lane,
-             const unsigned int,
-             const scalar_gradient_type &shape_gradient)
-      {
-        internal::VectorizedArrayTrait<Number>::get(value[0], vector_lane) +=
-          shape_gradient[0];
-      }
-
-      static scalar_gradient_type
-      access(const gradient_type &value,
-             const unsigned int   vector_lane,
-             const unsigned int)
-      {
-        scalar_gradient_type result;
-        result[0] =
-          internal::VectorizedArrayTrait<Number>::get(value[0], vector_lane);
-        return result;
-      }
-    };
-
     template <int dim, int spacedim>
     bool
     is_fast_path_supported(const FiniteElement<dim, spacedim> &fe,
@@ -761,11 +630,11 @@ public:
   using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait<
     Number>::vectorized_value_type;
   using ETT = typename internal::FEPointEvaluation::
-    EvaluatorTypeTraits<dim, n_components, Number>;
+    EvaluatorTypeTraits<dim, spacedim, n_components, Number>;
   using value_type            = typename ETT::value_type;
   using scalar_value_type     = typename ETT::scalar_value_type;
   using vectorized_value_type = typename ETT::vectorized_value_type;
-  using gradient_type         = typename ETT::gradient_type;
+  using gradient_type         = typename ETT::real_gradient_type;
   using interface_vectorized_gradient_type =
     typename ETT::interface_vectorized_gradient_type;
 
@@ -788,10 +657,10 @@ protected:
    * objects, this parameter allows to select a range of `n_components`
    * components starting from this parameter.
    */
-  FEPointEvaluationBase(const Mapping<dim>       &mapping,
-                        const FiniteElement<dim> &fe,
-                        const UpdateFlags         update_flags,
-                        const unsigned int        first_selected_component = 0);
+  FEPointEvaluationBase(const Mapping<dim, spacedim>       &mapping,
+                        const FiniteElement<dim, spacedim> &fe,
+                        const UpdateFlags                   update_flags,
+                        const unsigned int first_selected_component = 0);
 
   /**
    * Constructor to make the present class able to re-use the geometry
@@ -814,7 +683,7 @@ protected:
    */
   FEPointEvaluationBase(
     NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
-    const FiniteElement<dim>                        &fe,
+    const FiniteElement<dim, spacedim>              &fe,
     const unsigned int first_selected_component = 0,
     const bool         is_interior              = true);
 
@@ -995,7 +864,7 @@ protected:
   /**
    * Pointer to the FiniteElement object passed to the constructor.
    */
-  SmartPointer<const FiniteElement<dim>> fe;
+  SmartPointer<const FiniteElement<dim, spacedim>> fe;
 
   /**
    * Description of the 1d polynomial basis for tensor product elements used
@@ -1237,11 +1106,12 @@ public:
   using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait<
     Number>::vectorized_value_type;
   using ETT = typename internal::FEPointEvaluation::
-    EvaluatorTypeTraits<dim, n_components, Number>;
+    EvaluatorTypeTraits<dim, spacedim, n_components, Number>;
   using value_type            = typename ETT::value_type;
   using scalar_value_type     = typename ETT::scalar_value_type;
   using vectorized_value_type = typename ETT::vectorized_value_type;
-  using gradient_type         = typename ETT::gradient_type;
+  using unit_gradient_type    = typename ETT::gradient_type;
+  using gradient_type         = typename ETT::real_gradient_type;
   using interface_vectorized_gradient_type =
     typename ETT::interface_vectorized_gradient_type;
 
@@ -1263,10 +1133,10 @@ public:
    * objects, this parameter allows to select a range of `n_components`
    * components starting from this parameter.
    */
-  FEPointEvaluation(const Mapping<dim>       &mapping,
-                    const FiniteElement<dim> &fe,
-                    const UpdateFlags         update_flags,
-                    const unsigned int        first_selected_component = 0);
+  FEPointEvaluation(const Mapping<dim, spacedim>       &mapping,
+                    const FiniteElement<dim, spacedim> &fe,
+                    const UpdateFlags                   update_flags,
+                    const unsigned int first_selected_component = 0);
 
   /**
    * Constructor to make the present class able to re-use the geometry
@@ -1286,7 +1156,7 @@ public:
    */
   FEPointEvaluation(
     NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
-    const FiniteElement<dim>                        &fe,
+    const FiniteElement<dim, spacedim>              &fe,
     const unsigned int first_selected_component = 0);
 
   /**
@@ -1622,11 +1492,12 @@ public:
   using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait<
     Number>::vectorized_value_type;
   using ETT = typename internal::FEPointEvaluation::
-    EvaluatorTypeTraits<dim, n_components, Number>;
+    EvaluatorTypeTraits<dim, spacedim, n_components, Number>;
   using value_type            = typename ETT::value_type;
   using scalar_value_type     = typename ETT::scalar_value_type;
   using vectorized_value_type = typename ETT::vectorized_value_type;
-  using gradient_type         = typename ETT::gradient_type;
+  using unit_gradient_type    = typename ETT::gradient_type;
+  using gradient_type         = typename ETT::real_gradient_type;
   using interface_vectorized_gradient_type =
     typename ETT::interface_vectorized_gradient_type;
 
@@ -1635,7 +1506,7 @@ public:
    */
   FEFacePointEvaluation(
     NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
-    const FiniteElement<dim>                        &fe,
+    const FiniteElement<dim, spacedim>              &fe,
     const bool                                       is_interior = true,
     const unsigned int first_selected_component                  = 0);
 
@@ -1869,15 +1740,17 @@ private:
     const bool                              sum_into_values);
 };
 
+
+
 // ----------------------- template and inline function ----------------------
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
 FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
-  FEPointEvaluationBase(const Mapping<dim>       &mapping,
-                        const FiniteElement<dim> &fe,
-                        const UpdateFlags         update_flags,
-                        const unsigned int        first_selected_component)
+  FEPointEvaluationBase(const Mapping<dim, spacedim>       &mapping,
+                        const FiniteElement<dim, spacedim> &fe,
+                        const UpdateFlags                   update_flags,
+                        const unsigned int first_selected_component)
   : n_q_batches(numbers::invalid_unsigned_int)
   , n_q_points(numbers::invalid_unsigned_int)
   , n_q_points_scalar(numbers::invalid_unsigned_int)
@@ -1904,7 +1777,7 @@ template <int n_components_, int dim, int spacedim, typename Number>
 FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
   FEPointEvaluationBase(
     NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
-    const FiniteElement<dim>                        &fe,
+    const FiniteElement<dim, spacedim>              &fe,
     const unsigned int                               first_selected_component,
     const bool                                       is_interior)
   : n_q_batches(numbers::invalid_unsigned_int)
@@ -1972,9 +1845,7 @@ FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
 
 template <int n_components_, int dim, int spacedim, typename Number>
 FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
-  FEPointEvaluationBase(
-    FEPointEvaluationBase<n_components_, dim, spacedim, Number>
-      &&other) noexcept
+  FEPointEvaluationBase(FEPointEvaluationBase &&other) noexcept
   : n_q_batches(other.n_q_batches)
   , n_q_points(other.n_q_points)
   , n_q_points_scalar(other.n_q_points_scalar)
@@ -2373,7 +2244,7 @@ FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
 template <int n_components_, int dim, int spacedim, typename Number>
 FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
   NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
-  const FiniteElement<dim>                        &fe,
+  const FiniteElement<dim, spacedim>              &fe,
   const unsigned int                               first_selected_component)
   : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
       mapping_info,
@@ -2385,10 +2256,10 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
 
 template <int n_components_, int dim, int spacedim, typename Number>
 FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
-  const Mapping<dim>       &mapping,
-  const FiniteElement<dim> &fe,
-  const UpdateFlags         update_flags,
-  const unsigned int        first_selected_component)
+  const Mapping<dim, spacedim>       &mapping,
+  const FiniteElement<dim, spacedim> &fe,
+  const UpdateFlags                   update_flags,
+  const unsigned int                  first_selected_component)
   : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
       mapping,
       fe,
@@ -2743,13 +2614,13 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate_fast(
                v < stride && (stride == 1 || offset < this->n_q_points_scalar);
                ++v, ++offset)
             {
-              gradient_type unit_gradient;
+              unit_gradient_type unit_gradient;
               ETT::set_gradient(gradient, v, unit_gradient);
               this->gradients[offset] =
                 this->cell_type <=
                     internal::MatrixFreeFunctions::GeometryType::cartesian ?
-                  apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
-                                                unit_gradient) :
+                  apply_diagonal_transformation(
+                    this->inverse_jacobian_ptr[0].transpose(), unit_gradient) :
                   apply_transformation(
                     this
                       ->inverse_jacobian_ptr[this->cell_type <=
@@ -2997,7 +2868,7 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate_fast(
              v < stride && (stride == 1 || offset < this->n_q_points_scalar);
              ++v, ++offset)
           {
-            const auto grad_w =
+            const gradient_type grad_w =
               do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] :
                        this->gradients[offset];
             ETT::get_gradient(
@@ -3178,7 +3049,7 @@ template <int n_components_, int dim, int spacedim, typename Number>
 FEFacePointEvaluation<n_components_, dim, spacedim, Number>::
   FEFacePointEvaluation(
     NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
-    const FiniteElement<dim>                        &fe,
+    const FiniteElement<dim, spacedim>              &fe,
     const bool                                       is_interior,
     const unsigned int                               first_selected_component)
   : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
@@ -3665,13 +3536,13 @@ FEFacePointEvaluation<n_components_, dim, spacedim, Number>::
                v < stride && (stride == 1 || offset < this->n_q_points_scalar);
                ++v, ++offset)
             {
-              gradient_type unit_gradient;
+              unit_gradient_type unit_gradient;
               ETT::set_gradient(gradient, v, unit_gradient);
               this->gradients[offset] =
                 this->cell_type <=
                     internal::MatrixFreeFunctions::GeometryType::cartesian ?
-                  apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
-                                                unit_gradient) :
+                  apply_diagonal_transformation(
+                    this->inverse_jacobian_ptr[0].transpose(), unit_gradient) :
                   apply_transformation(
                     this
                       ->inverse_jacobian_ptr[this->cell_type <=

diff --git a/include/deal.II/matrix_free/fe_remote_evaluation.h b/include/deal.II/matrix_free/fe_remote_evaluation.h
index 930dd9ef17718cef20e310bf9f4604eb5c6631fb..c3e512070e2664485e7f8599fdd70fb4f6707430 100644 (file)
--- a/include/deal.II/matrix_free/fe_remote_evaluation.h
+++ b/include/deal.II/matrix_free/fe_remote_evaluation.h
@@ -42,10 +42,10 @@ namespace internal
   struct PrecomputedEvaluationData
   {
     using value_type = typename internal::FEPointEvaluation::
-      EvaluatorTypeTraits<dim, n_components, value_type_>::value_type;
+      EvaluatorTypeTraits<dim, dim, n_components, value_type_>::value_type;
 
     using gradient_type = typename internal::FEPointEvaluation::
-      EvaluatorTypeTraits<dim, n_components, value_type_>::gradient_type;
+      EvaluatorTypeTraits<dim, dim, n_components, value_type_>::gradient_type;
 
     /**
      * Values at quadrature points.

diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index dec73fa9920cb26076865419814bf0d26c051af1..2eafca6dbf1344d8a414bfca3229eaa3dd86ac46 100644 (file)
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -84,23 +84,23 @@ namespace internal
 
 
 
-    template <int dim>
+    template <int dim, int spacedim>
     void
     get_element_type_specific_information(
-      const FiniteElement<dim, dim> &fe_in,
-      const FiniteElement<dim, dim> &fe,
-      const unsigned int             base_element_number,
-      ElementType                   &element_type,
-      std::vector<unsigned int>     &scalar_lexicographic,
-      std::vector<unsigned int>     &lexicographic_numbering)
+      const FiniteElement<dim, spacedim> &fe_in,
+      const FiniteElement<dim, spacedim> &fe,
+      const unsigned int                  base_element_number,
+      ElementType                        &element_type,
+      std::vector<unsigned int>          &scalar_lexicographic,
+      std::vector<unsigned int>          &lexicographic_numbering)
     {
       element_type = tensor_general;
 
-      const auto fe_poly = dynamic_cast<const FE_Poly<dim, dim> *>(&fe);
+      const auto fe_poly = dynamic_cast<const FE_Poly<dim, spacedim> *>(&fe);
 
-      if (dynamic_cast<const FE_SimplexPoly<dim, dim> *>(&fe) != nullptr ||
-          dynamic_cast<const FE_WedgePoly<dim, dim> *>(&fe) != nullptr ||
-          dynamic_cast<const FE_PyramidPoly<dim, dim> *>(&fe) != nullptr)
+      if (dynamic_cast<const FE_SimplexPoly<dim, spacedim> *>(&fe) != nullptr ||
+          dynamic_cast<const FE_WedgePoly<dim, spacedim> *>(&fe) != nullptr ||
+          dynamic_cast<const FE_PyramidPoly<dim, spacedim> *>(&fe) != nullptr)
         {
           scalar_lexicographic.resize(fe.n_dofs_per_cell());
           for (unsigned int i = 0; i < scalar_lexicographic.size(); ++i)
@@ -115,14 +115,16 @@ namespace internal
                     Polynomials::PiecewisePolynomial<double>> *>(
                   &fe_poly->get_poly_space()) != nullptr))
         scalar_lexicographic = fe_poly->get_poly_space_numbering_inverse();
-      else if (const auto fe_dgp = dynamic_cast<const FE_DGP<dim> *>(&fe))
+      else if (const auto fe_dgp =
+                 dynamic_cast<const FE_DGP<dim, spacedim> *>(&fe))
         {
           scalar_lexicographic.resize(fe_dgp->n_dofs_per_cell());
           for (unsigned int i = 0; i < fe_dgp->n_dofs_per_cell(); ++i)
             scalar_lexicographic[i] = i;
           element_type = truncated_tensor;
         }
-      else if (const auto fe_q_dg0 = dynamic_cast<const FE_Q_DG0<dim> *>(&fe))
+      else if (const auto fe_q_dg0 =
+                 dynamic_cast<const FE_Q_DG0<dim, spacedim> *>(&fe))
         {
           scalar_lexicographic = fe_q_dg0->get_poly_space_numbering_inverse();
           element_type         = tensor_symmetric_plus_dg0;
@@ -190,7 +192,7 @@ namespace internal
         Assert(fe_name[template_starts + 1] ==
                  (dim == 1 ? '1' : (dim == 2 ? '2' : '3')),
                ExcInternalError());
-        fe_name[template_starts + 1] = std::to_string(dim_to).c_str()[0];
+        fe_name[template_starts + 1] = std::to_string(dim_to)[0];
       }
       return FETools::get_fe_by_name<dim_to, dim_to>(fe_name);
     }
@@ -237,9 +239,6 @@ namespace internal
                               const FiniteElement<dim, spacedim> &fe_in,
                               const unsigned int base_element_number)
     {
-      static_assert(dim == spacedim,
-                    "Currently, only the case dim=spacedim is implemented");
-
       // ShapeInfo for RT elements. Here, data is of size 2 instead of 1.
       // data[0] is univariate_shape_data in normal direction and
       // data[1] is univariate_shape_data in tangential direction
@@ -251,7 +250,7 @@ namespace internal
 
           const auto quad = quad_in.get_tensor_basis()[0];
 
-          const FiniteElement<dim> &fe =
+          const FiniteElement<dim, spacedim> &fe =
             fe_in.base_element(base_element_number);
           n_dimensions = dim;
           n_components = fe_in.n_components();
@@ -312,11 +311,11 @@ namespace internal
           return;
         }
       else if (quad_in.is_tensor_product() == false ||
-               dynamic_cast<const FE_SimplexPoly<dim, dim> *>(
+               dynamic_cast<const FE_SimplexPoly<dim, spacedim> *>(
                  &fe_in.base_element(base_element_number)) != nullptr ||
-               dynamic_cast<const FE_WedgePoly<dim, dim> *>(
+               dynamic_cast<const FE_WedgePoly<dim, spacedim> *>(
                  &fe_in.base_element(base_element_number)) != nullptr ||
-               dynamic_cast<const FE_PyramidPoly<dim, dim> *>(
+               dynamic_cast<const FE_PyramidPoly<dim, spacedim> *>(
                  &fe_in.base_element(base_element_number)) != nullptr)
         {
           // specialization for arbitrary finite elements and quadrature rules
@@ -504,9 +503,10 @@ namespace internal
 
       const auto quad = quad_in.get_tensor_basis()[0];
 
-      const FiniteElement<dim> &fe = fe_in.base_element(base_element_number);
-      n_dimensions                 = dim;
-      n_components                 = fe_in.n_components();
+      const FiniteElement<dim, spacedim> &fe =
+        fe_in.base_element(base_element_number);
+      n_dimensions = dim;
+      n_components = fe_in.n_components();
 
       Assert(fe.n_components() == 1,
              ExcMessage("FEEvaluation only works for scalar finite elements."));
@@ -553,8 +553,8 @@ namespace internal
                                                      scalar_lexicographic,
                                                      0);
 
-      if (dim > 1 && (dynamic_cast<const FE_Q<dim> *>(&fe) ||
-                      dynamic_cast<const FE_Q_iso_Q1<dim> *>(&fe)))
+      if (dim > 1 && (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe) ||
+                      dynamic_cast<const FE_Q_iso_Q1<dim, spacedim> *>(&fe)))
         {
           auto &subface_interpolation_matrix_0 =
             univariate_shape_data.subface_interpolation_matrices[0];
@@ -571,7 +571,8 @@ namespace internal
           subface_interpolation_matrix_scalar_0.resize(nn * nn);
           subface_interpolation_matrix_scalar_1.resize(nn * nn);
 
-          const bool is_feq = dynamic_cast<const FE_Q<dim> *>(&fe) != nullptr;
+          const bool is_feq =
+            dynamic_cast<const FE_Q<dim, spacedim> *>(&fe) != nullptr;
 
           std::vector<Point<1>> fe_q_points =
             is_feq ? QGaussLobatto<1>(nn).get_points() :
@@ -613,7 +614,7 @@ namespace internal
       if (element_type == tensor_general &&
           univariate_shape_data.check_and_set_shapes_symmetric())
         {
-          if (dynamic_cast<const FE_Q_iso_Q1<dim> *>(&fe) &&
+          if (dynamic_cast<const FE_Q_iso_Q1<dim, spacedim> *>(&fe) &&
               fe.tensor_degree() > 1)
             element_type = tensor_symmetric_no_collocation;
           else if (univariate_shape_data.check_shapes_collocation())
@@ -1151,9 +1152,6 @@ namespace internal
     bool
     ShapeInfo<Number>::is_supported(const FiniteElement<dim, spacedim> &fe)
     {
-      if (dim != spacedim)
-        return false;
-
       if (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe))
         return true;
 
@@ -1170,12 +1168,12 @@ namespace internal
                 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_SimplexPoly<dim, dim> *>(fe_poly_ptr) !=
-                    nullptr ||
-                  dynamic_cast<const FE_WedgePoly<dim, dim> *>(fe_poly_ptr) !=
-                    nullptr ||
-                  dynamic_cast<const FE_PyramidPoly<dim, dim> *>(fe_poly_ptr) !=
-                    nullptr)
+              if (dynamic_cast<const FE_SimplexPoly<dim, spacedim> *>(
+                    fe_poly_ptr) != nullptr ||
+                  dynamic_cast<const FE_WedgePoly<dim, spacedim> *>(
+                    fe_poly_ptr) != nullptr ||
+                  dynamic_cast<const FE_PyramidPoly<dim, spacedim> *>(
+                    fe_poly_ptr) != nullptr)
                 return true;
 
               if (dynamic_cast<const TensorProductPolynomials<

diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
index 69cabe4a3b1c9b6de7c4bc6a803ff23947abc2d9..a7fbda27361d58237c9ffc5ae748e5fbefc09ddb 100644 (file)
--- a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
+++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
@@ -4896,8 +4896,8 @@ MGTwoLevelTransferNonNested<dim, VectorType>::prolongate_and_add_internal_comp(
   VectorType       &dst,
   const VectorType &src) const
 {
-  using Traits =
-    internal::FEPointEvaluation::EvaluatorTypeTraits<dim, n_components, Number>;
+  using Traits = internal::FEPointEvaluation::
+    EvaluatorTypeTraits<dim, dim, n_components, Number>;
   using value_type = typename Traits::value_type;
 
   std::vector<value_type> evaluation_point_results;
@@ -5023,8 +5023,8 @@ MGTwoLevelTransferNonNested<dim, VectorType>::restrict_and_add_internal_comp(
   VectorType       &dst,
   const VectorType &src) const
 {
-  using Traits =
-    internal::FEPointEvaluation::EvaluatorTypeTraits<dim, n_components, Number>;
+  using Traits = internal::FEPointEvaluation::
+    EvaluatorTypeTraits<dim, dim, n_components, Number>;
   using value_type = typename Traits::value_type;
 
   std::vector<value_type> evaluation_point_results;

diff --git a/include/deal.II/non_matching/mapping_info.h b/include/deal.II/non_matching/mapping_info.h
index ea51fa843596ebe3761087703cfe8fca6874b68a..1742fd46d345f2da770ceb2158d502677597346b 100644 (file)
--- a/include/deal.II/non_matching/mapping_info.h
+++ b/include/deal.II/non_matching/mapping_info.h
@@ -186,7 +186,7 @@ namespace NonMatching
     dealii::internal::MatrixFreeFunctions::GeometryType
     compute_geometry_type(
       const double diameter,
-      const std::vector<DerivativeForm<1, dim, spacedim, double>>
+      const std::vector<DerivativeForm<1, spacedim, dim, double>>
         &inverse_jacobians)
     {
       const auto   jac_0 = inverse_jacobians[0];
@@ -195,7 +195,7 @@ namespace NonMatching
       bool jacobian_constant = true;
       for (unsigned int q = 1; q < inverse_jacobians.size(); ++q)
         {
-          const DerivativeForm<1, dim, spacedim> &jac = inverse_jacobians[q];
+          const DerivativeForm<1, spacedim, dim> &jac = inverse_jacobians[q];
           for (unsigned int d = 0; d < dim; ++d)
             for (unsigned int e = 0; e < spacedim; ++e)
               if (std::fabs(jac_0[d][e] - jac[d][e]) > zero_tolerance_double)
@@ -206,7 +206,7 @@ namespace NonMatching
 
       // check whether the Jacobian is diagonal to machine
       // accuracy
-      bool cell_cartesian = jacobian_constant;
+      bool cell_cartesian = jacobian_constant && dim == spacedim;
       for (unsigned int d = 0; d < dim; ++d)
         for (unsigned int e = 0; e < dim; ++e)
           if (d != e)
@@ -296,8 +296,8 @@ namespace NonMatching
      * @param additional_data Additional data for the class to specify the
      * behavior during reinitialization.
      */
-    MappingInfo(const Mapping<dim>  &mapping,
-                const UpdateFlags    update_flags,
+    MappingInfo(const Mapping<dim, spacedim> &mapping,
+                const UpdateFlags             update_flags,
                 const AdditionalData additional_data = AdditionalData());
 
     /**
@@ -448,7 +448,7 @@ namespace NonMatching
      * Getter function for real points. The offset can be obtained with
      * compute_data_index_offset().
      */
-    const Point<dim, Number> *
+    const Point<spacedim, Number> *
     get_real_point(const unsigned int offset) const;
 
     /**
@@ -797,9 +797,9 @@ namespace NonMatching
 
   template <int dim, int spacedim, typename Number>
   MappingInfo<dim, spacedim, Number>::MappingInfo(
-    const Mapping<dim>  &mapping,
-    const UpdateFlags    update_flags,
-    const AdditionalData additional_data)
+    const Mapping<dim, spacedim> &mapping,
+    const UpdateFlags             update_flags,
+    const AdditionalData          additional_data)
     : mapping(&mapping)
     , update_flags(update_flags)
     , update_flags_mapping(update_default)
@@ -1960,8 +1960,8 @@ namespace NonMatching
                         v) = mapping_data.jacobians[q * n_lanes + v][d][s];
                 if (update_flags_mapping &
                     UpdateFlags::update_inverse_jacobians)
-                  for (unsigned int d = 0; d < dim; ++d)
-                    for (unsigned int s = 0; s < spacedim; ++s)
+                  for (unsigned int d = 0; d < spacedim; ++d)
+                    for (unsigned int s = 0; s < dim; ++s)
                       dealii::internal::VectorizedArrayTrait<Number>::get(
                         inverse_jacobians[is_interior ? 0 : 1]
                                          [compressed_offset][s][d],
@@ -2076,7 +2076,7 @@ namespace NonMatching
 
 
   template <int dim, int spacedim, typename Number>
-  inline const Point<dim, Number> *
+  inline const Point<spacedim, Number> *
   MappingInfo<dim, spacedim, Number>::get_real_point(
     const unsigned int offset) const
   {

diff --git a/source/matrix_free/fe_point_evaluation.cc b/source/matrix_free/fe_point_evaluation.cc
index 958e1ef2e15d16f9ef54e1f634f260c5792a0181..149569b8fa1cf605a277873e37cd19f1c9b0b0f5 100644 (file)
--- a/source/matrix_free/fe_point_evaluation.cc
+++ b/source/matrix_free/fe_point_evaluation.cc
@@ -43,9 +43,6 @@ namespace internal
     {
       // check if supported
       const bool flag = [&]() {
-        if (dim != spacedim)
-          return false;
-
         const FiniteElement<dim, spacedim> *fe_ptr =
           &(fe.base_element(base_element_number));
         if (fe_ptr->n_components() != 1)

diff --git a/source/matrix_free/shape_info.inst.in b/source/matrix_free/shape_info.inst.in
index c07dcea0e8d4bd5d3556fd61f95e9c033b3b14c0..d013a949e0abaf90e0445396c8c0f0841ff7ffd1 100644 (file)
--- a/source/matrix_free/shape_info.inst.in
+++ b/source/matrix_free/shape_info.inst.in
@@ -13,30 +13,35 @@
 // ------------------------------------------------------------------------
 
 
-for (deal_II_dimension : DIMENSIONS; deal_II_scalar : REAL_SCALARS)
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS;
+     deal_II_scalar : REAL_SCALARS)
   {
-    template internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar>::
-      ShapeInfo(const Quadrature<deal_II_dimension> &,
-                const FiniteElement<deal_II_dimension, deal_II_dimension> &,
-                const unsigned int);
+#if deal_II_dimension <= deal_II_space_dimension
+    template internal::MatrixFreeFunctions::ShapeInfo<
+      deal_II_scalar>::ShapeInfo(const Quadrature<deal_II_dimension> &,
+                                 const FiniteElement<deal_II_dimension,
+                                                     deal_II_space_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 FiniteElement<deal_II_dimension, deal_II_space_dimension> &,
       const unsigned int);
 
-#if deal_II_dimension > 1
-    template internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar>::
-      ShapeInfo(const Quadrature<1> &,
-                const FiniteElement<deal_II_dimension, deal_II_dimension> &,
-                const unsigned int);
+#  if deal_II_dimension > 1
+    template internal::MatrixFreeFunctions::ShapeInfo<
+      deal_II_scalar>::ShapeInfo(const Quadrature<1> &,
+                                 const FiniteElement<deal_II_dimension,
+                                                     deal_II_space_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 FiniteElement<deal_II_dimension, deal_II_space_dimension> &,
       const unsigned int);
+#  endif
 #endif
   }
 
@@ -51,33 +56,35 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
 
 
 
-for (deal_II_dimension : DIMENSIONS;
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS;
      deal_II_scalar_vectorized : REAL_SCALARS_VECTORIZED)
   {
+#if deal_II_dimension <= deal_II_space_dimension
     template internal::MatrixFreeFunctions::
       ShapeInfo<deal_II_scalar_vectorized>::ShapeInfo(
         const Quadrature<deal_II_dimension> &,
-        const FiniteElement<deal_II_dimension, deal_II_dimension> &,
+        const FiniteElement<deal_II_dimension, deal_II_space_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 FiniteElement<deal_II_dimension, deal_II_space_dimension> &,
       const unsigned int);
 
-#if deal_II_dimension > 1
+#  if deal_II_dimension > 1
     template internal::MatrixFreeFunctions::
       ShapeInfo<deal_II_scalar_vectorized>::ShapeInfo(
         const Quadrature<1> &,
-        const FiniteElement<deal_II_dimension, deal_II_dimension> &,
+        const FiniteElement<deal_II_dimension, deal_II_space_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 FiniteElement<deal_II_dimension, deal_II_space_dimension> &,
       const unsigned int);
+#  endif
 #endif
   }