]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Speed up evaluation of RefSpaceFEFieldFunction
authorMartin Kronbichler <martin.kronbichler@it.uu.se>
Fri, 28 Jan 2022 13:52:08 +0000 (14:52 +0100)
committerMartin Kronbichler <martin.kronbichler@it.uu.se>
Fri, 28 Jan 2022 21:25:55 +0000 (22:25 +0100)
include/deal.II/matrix_free/tensor_product_kernels.h
source/non_matching/fe_values.cc
source/non_matching/quadrature_generator.cc

index 5948e7c355704903c27b88ed3d039d1587552d1d..c09b652bb3f1a7b4ec8f1fb8abe155db9d01a5fc 100644 (file)
@@ -2523,6 +2523,105 @@ namespace internal
 
 
 
+  template <int dim, typename Number, typename Number2>
+  SymmetricTensor<2, dim, typename ProductTypeNoPoint<Number, Number2>::type>
+  evaluate_tensor_product_hessian(
+    const std::vector<Polynomials::Polynomial<double>> &poly,
+    const std::vector<Number> &                         values,
+    const Point<dim, Number2> &                         p,
+    const std::vector<unsigned int> &                   renumber = {})
+  {
+    static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
+
+    using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
+
+    // use `int` type for this variable and the loops below to inform the
+    // compiler that the loops below will never overflow, which allows it to
+    // generate more optimized code for the variable loop bounds in the
+    // present context
+    const int n_shapes = poly.size();
+    AssertDimension(Utilities::pow(n_shapes, dim), values.size());
+    Assert(renumber.empty() || renumber.size() == values.size(),
+           ExcDimensionMismatch(renumber.size(), values.size()));
+
+    AssertIndexRange(n_shapes, 200);
+    std::array<Number2, 3 * dim * 200> shapes;
+
+    // Evaluate 1D polynomials and their derivatives
+    for (unsigned int d = 0; d < dim; ++d)
+      for (int i = 0; i < n_shapes; ++i)
+        poly[i].value(p[d], 2, shapes.data() + 3 * (d * n_shapes + i));
+
+    // Go through the tensor product of shape functions and interpolate
+    // with optimal algorithm
+    SymmetricTensor<2, dim, Number3> result;
+    for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
+      {
+        Number3 value_y = {}, deriv_x = {}, deriv_y = {}, deriv_xx = {},
+                deriv_xy = {}, deriv_yy = {};
+        for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
+          {
+            // Interpolation + derivative x direction
+            Number3 value = {}, deriv_1 = {}, deriv_2 = {};
+
+            // Distinguish the inner loop based on whether we have a
+            // renumbering or not
+            if (renumber.empty())
+              for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
+                {
+                  value += shapes[3 * i0] * values[i];
+                  deriv_1 += shapes[3 * i0 + 1] * values[i];
+                  deriv_2 += shapes[3 * i0 + 2] * values[i];
+                }
+            else
+              for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
+                {
+                  value += shapes[3 * i0] * values[renumber[i]];
+                  deriv_1 += shapes[3 * i0 + 1] * values[renumber[i]];
+                  deriv_2 += shapes[3 * i0 + 2] * values[renumber[i]];
+                }
+
+            // Interpolation + derivative in y direction
+            if (dim > 1)
+              {
+                if (dim > 2)
+                  {
+                    value_y += value * shapes[3 * n_shapes + 3 * i1];
+                    deriv_x += deriv_1 * shapes[3 * n_shapes + 3 * i1];
+                    deriv_y += value * shapes[3 * n_shapes + 3 * i1 + 1];
+                  }
+                deriv_xx += deriv_2 * shapes[3 * n_shapes + 3 * i1];
+                deriv_xy += deriv_1 * shapes[3 * n_shapes + 3 * i1 + 1];
+                deriv_yy += value * shapes[3 * n_shapes + 3 * i1 + 2];
+              }
+            else
+              {
+                result[0][0] = deriv_2;
+              }
+          }
+        if (dim == 3)
+          {
+            // Interpolation + derivative in z direction
+            result[0][0] += deriv_xx * shapes[6 * n_shapes + 3 * i2];
+            result[0][1] += deriv_xy * shapes[6 * n_shapes + 3 * i2];
+            result[0][2] += deriv_x * shapes[6 * n_shapes + 3 * i2 + 1];
+            result[1][1] += deriv_yy * shapes[6 * n_shapes + 3 * i2];
+            result[1][2] += deriv_y * shapes[6 * n_shapes + 3 * i2 + 1];
+            result[2][2] += value_y * shapes[6 * n_shapes + 3 * i2 + 2];
+          }
+        else if (dim == 2)
+          {
+            result[0][0] = deriv_xx;
+            result[1][0] = deriv_xy;
+            result[1][1] = deriv_yy;
+          }
+      }
+
+    return result;
+  }
+
+
+
   /**
    * Same as evaluate_tensor_product_value_and_gradient() but for integration.
    */
index 2b26e6aa85b13a95421eb0abc261d7d5ff496577..d4e0bcd40a4fabd01dc828cd6fc4cc8d4c64c90d 100644 (file)
@@ -28,6 +28,8 @@
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/vector.h>
 
+#include <deal.II/matrix_free/fe_point_evaluation.h>
+
 #include <deal.II/non_matching/fe_values.h>
 
 DEAL_II_NAMESPACE_OPEN
@@ -163,6 +165,23 @@ namespace NonMatching
          * set_active_cell()
          */
         std::vector<typename VectorType::value_type> local_dof_values;
+
+        /**
+         * Description of the 1D polynomial basis for tensor product elements
+         * used for the fast path of this class using tensor product
+         * evaluators.
+         */
+        std::vector<Polynomials::Polynomial<double>> poly;
+
+        /**
+         * Renumbering for the tensor-product evaluator in the fast path.
+         */
+        std::vector<unsigned int> renumber;
+
+        /**
+         * Check whether the shape functions are linear.
+         */
+        bool polynomials_are_hat_functions;
       };
 
 
@@ -199,7 +218,36 @@ namespace NonMatching
 
         // Save the element and the local dof values, since this is what we need
         // to evaluate the function.
-        element = &dof_handler_cell->get_fe();
+
+        // Check if we can use the fast path. In case we have a different
+        // element from the one used before we need to set up the data
+        // structures again.
+        if (element != &dof_handler_cell->get_fe())
+          {
+            poly.clear();
+            element = &dof_handler_cell->get_fe();
+
+            if (element->n_base_elements() == 1 &&
+                dealii::internal::FEPointEvaluation::is_fast_path_supported(
+                  *element, 0))
+              {
+                dealii::internal::MatrixFreeFunctions::ShapeInfo<double>
+                  shape_info;
+
+                shape_info.reinit(QMidpoint<1>(), *element, 0);
+                renumber = shape_info.lexicographic_numbering;
+                poly =
+                  dealii::internal::FEPointEvaluation::get_polynomial_space(
+                    element->base_element(0));
+
+                polynomials_are_hat_functions =
+                  (poly.size() == 2 && poly[0].value(0.) == 1. &&
+                   poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
+                   poly[1].value(1.) == 1.);
+              }
+          }
+        else
+          element = &dof_handler_cell->get_fe();
 
         local_dof_indices.resize(element->dofs_per_cell);
         dof_handler_cell->get_dof_indices(local_dof_indices);
@@ -234,12 +282,26 @@ namespace NonMatching
         AssertIndexRange(component, this->n_components);
         Assert(cell_is_set(), ExcCellNotSet());
 
-        double value = 0;
-        for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
-          value += local_dof_values[i] *
-                   element->shape_value_component(i, point, component);
+        if (!poly.empty() && component == 0)
+          {
+            // TODO: this could be extended to a component that is not zero
+            return dealii::internal::evaluate_tensor_product_value_and_gradient(
+                     poly,
+                     local_dof_values,
+                     point,
+                     polynomials_are_hat_functions,
+                     renumber)
+              .first;
+          }
+        else
+          {
+            double value = 0;
+            for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
+              value += local_dof_values[i] *
+                       element->shape_value_component(i, point, component);
 
-        return value;
+            return value;
+          }
       }
 
 
@@ -253,12 +315,26 @@ namespace NonMatching
         AssertIndexRange(component, this->n_components);
         Assert(cell_is_set(), ExcCellNotSet());
 
-        Tensor<1, dim> gradient;
-        for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
-          gradient += local_dof_values[i] *
-                      element->shape_grad_component(i, point, component);
+        if (!poly.empty() && component == 0)
+          {
+            // TODO: this could be extended to a component that is not zero
+            return dealii::internal::evaluate_tensor_product_value_and_gradient(
+                     poly,
+                     local_dof_values,
+                     point,
+                     polynomials_are_hat_functions,
+                     renumber)
+              .second;
+          }
+        else
+          {
+            Tensor<1, dim> gradient;
+            for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
+              gradient += local_dof_values[i] *
+                          element->shape_grad_component(i, point, component);
 
-        return gradient;
+            return gradient;
+          }
       }
 
 
@@ -272,12 +348,22 @@ namespace NonMatching
         AssertIndexRange(component, this->n_components);
         Assert(cell_is_set(), ExcCellNotSet());
 
-        Tensor<2, dim> hessian;
-        for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
-          hessian += local_dof_values[i] *
-                     element->shape_grad_grad_component(i, point, component);
+        if (!poly.empty() && component == 0)
+          {
+            // TODO: this could be extended to a component that is not zero
+            return dealii::internal::evaluate_tensor_product_hessian(
+              poly, local_dof_values, point, renumber);
+          }
+        else
+          {
+            Tensor<2, dim> hessian;
+            for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
+              hessian +=
+                local_dof_values[i] *
+                element->shape_grad_grad_component(i, point, component);
 
-        return symmetrize(hessian);
+            return symmetrize(hessian);
+          }
       }
     } // namespace FEValuesImplementation
   }   // namespace internal
@@ -353,20 +439,18 @@ namespace NonMatching
 
     Assert(fe_collection->size() > 0,
            ExcMessage("Incoming hp::FECollection can not be empty."));
-    Assert(
-      mapping_collection->size() == fe_collection->size() ||
-        mapping_collection->size() == 1,
-      ExcMessage(
-        "Size of hp::MappingCollection must be the same as hp::FECollection or 1."));
-    Assert(
-      q_collection.size() == fe_collection->size() || q_collection.size() == 1,
-      ExcMessage(
-        "Size of hp::QCollection<dim> must be the same as hp::FECollection or 1."));
-    Assert(
-      q_collection_1D.size() == fe_collection->size() ||
-        q_collection_1D.size() == 1,
-      ExcMessage(
-        "Size of hp::QCollection<1> must be the same as hp::FECollection or 1."));
+    Assert(mapping_collection->size() == fe_collection->size() ||
+             mapping_collection->size() == 1,
+           ExcMessage("Size of hp::MappingCollection must be "
+                      "the same as hp::FECollection or 1."));
+    Assert(q_collection.size() == fe_collection->size() ||
+             q_collection.size() == 1,
+           ExcMessage("Size of hp::QCollection<dim> must be the "
+                      "same as hp::FECollection or 1."));
+    Assert(q_collection_1D.size() == fe_collection->size() ||
+             q_collection_1D.size() == 1,
+           ExcMessage("Size of hp::QCollection<1> must be the "
+                      "same as hp::FECollection or 1."));
 
     // For each element in fe_collection, create dealii::FEValues objects to use
     // on the non-intersected cells.
@@ -442,8 +526,8 @@ namespace NonMatching
               quadrature_generator.get_surface_quadrature();
 
             // Even if a cell is formally intersected the number of created
-            // quadrature points can be 0. Avoid creating an FEValues object if
-            // that is the case.
+            // quadrature points can be 0. Avoid creating an FEValues object
+            // if that is the case.
             if (inside_quadrature.size() > 0)
               {
                 fe_values_inside.emplace((*mapping_collection)[mapping_index],
index 0775c86910432f59dc6e668d9edcedc4cfbc3621..9f446a375db6e133655334ca8c67d3bd6b576165 100644 (file)
@@ -597,7 +597,7 @@ namespace NonMatching
             const double function_max =
               std::max(std::max(left_value, right_value), value_bounds.second);
 
-            // If the functions is negative there are no roots.
+            // If the function is negative there are no roots.
             if (function_max < 0)
               return;
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.