]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Optimize FEPointEvaluation::integrate() with tensor product path
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 17 May 2021 17:41:57 +0000 (19:41 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 18 May 2021 07:30:06 +0000 (09:30 +0200)
include/deal.II/fe/fe_point_evaluation.h
include/deal.II/fe/mapping_q_generic.h
include/deal.II/matrix_free/tensor_product_kernels.h

index 552e3bf7e1ccd1698857f14abcdd5f851d49bdfb..75213d9640c545242125d00aa4890a39a8352731 100644 (file)
@@ -408,8 +408,9 @@ public:
                     const unsigned int        first_selected_component = 0);
 
   /**
-   * This is one of the main functions in this class and the one that does the
-   * heavy lifting.
+   * This function interpolates the finite element solution, represented by
+   * `solution_values` on the given cell, to the `unit_points` passed to the
+   * function.
    *
    * @param[in] cell An iterator to the current cell in question
    *
@@ -430,16 +431,28 @@ public:
            const EvaluationFlags::EvaluationFlags &evaluation_flags);
 
   /**
-   * This is one of the main functions in this class and the one that does the
-   * heavy lifting.
+   * This function multiplies the quantities passed in by previous
+   * submit_value() or submit_gradient() calls by the value or gradient of the
+   * test functions, and performs summation over all given points. This is
+   * similar to the integration of a bilinear form in terms of the test
+   * function, with the difference that this formula does not include a `JxW`
+   * factor. This allows the class to naturally embed point information
+   * (e.g. particles) into a finite element formulation. Of course, by
+   * multiplication of a `JxW` information of the data given to
+   * submit_value(), the integration can also be represented by this class.
    *
    * @param[in] cell An iterator to the current cell in question
    *
    * @param[in] unit_points List of points in the reference locations of the
    * current cell where the FiniteElement object should be evaluated
    *
-   * @param[out] solution_values This array is filled and can be used to
-   * during `cell->set_dof_values(global_vector, solution_values)`.
+   * @param[out] solution_values This array will contain the result of the
+   * integral, which can be used to during
+   * `cell->set_dof_values(solution_values, global_vector)` or
+   * `cell->distribute_local_to_global(solution_values, global_vector)`. Note
+   * that for multi-component systems where only some of the components are
+   * selected by the present class, the entries not touched by this class will
+   * be zeroed out.
    *
    * @param[in] integration_flags Flags specifying which quantities should be
    * integrated at the points.
@@ -537,6 +550,18 @@ private:
    */
   std::vector<value_type> solution_renumbered;
 
+  /**
+   * Temporary array to store a vectorized version of the `solution_values`
+   * computed during `integrate()` in a format compatible with the tensor
+   * product evaluators. For vector-valued setups, this array uses a
+   * `Tensor<1, n_components, VectorizedArray<Number>>` format.
+   */
+  AlignedVector<typename internal::FEPointEvaluation::EvaluatorTypeTraits<
+    dim,
+    n_components,
+    VectorizedArray<Number>>::value_type>
+    solution_renumbered_vectorized;
+
   /**
    * Temporary array to store the values at the points.
    */
@@ -606,7 +631,7 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::FEPointEvaluation(
          poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
          poly[1].value(1.) == 1.);
     }
-  if (true /*TODO: as long as the fast path of integrate() is not working*/)
+  else
     {
       nonzero_shape_function_component.resize(fe.n_dofs_per_cell());
       for (unsigned int d = 0; d < n_components; ++d)
@@ -702,6 +727,7 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate(
           mapping_q_generic->transform_variable(
             cell,
             mapping_covariant,
+            /* apply_from_left */ true,
             unit_points,
             ArrayView<const gradient_type>(gradients.data(), gradients.size()),
             ArrayView<gradient_type>(gradients.data(), gradients.size()));
@@ -789,17 +815,13 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
     }
 
   AssertDimension(solution_values.size(), fe->dofs_per_cell);
-  if (false /*TODO*/ && (((integration_flags & EvaluationFlags::values) ||
-                          (integration_flags & EvaluationFlags::gradients)) &&
-                         !poly.empty()))
+  if (((integration_flags & EvaluationFlags::values) ||
+       (integration_flags & EvaluationFlags::gradients)) &&
+      !poly.empty())
     {
-      Assert(false, ExcNotImplemented());
       // fast path with tensor product integration
 
-      if (solution_renumbered.size() != dofs_per_component)
-        solution_renumbered.resize(dofs_per_component);
-
-      // let mapping compute the transformation
+      // let mapping apply the transformation
       if (integration_flags & EvaluationFlags::gradients)
         {
           Assert(mapping_q_generic != nullptr, ExcInternalError());
@@ -807,6 +829,7 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
           mapping_q_generic->transform_variable(
             cell,
             mapping_covariant,
+            /* apply_from_left */ false,
             unit_points,
             ArrayView<const gradient_type>(gradients.data(), gradients.size()),
             ArrayView<gradient_type>(gradients.data(), gradients.size()));
@@ -817,6 +840,15 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
       if (integration_flags & EvaluationFlags::gradients)
         AssertIndexRange(unit_points.size(), gradients.size() + 1);
 
+      if (solution_renumbered_vectorized.size() != dofs_per_component)
+        solution_renumbered_vectorized.resize(dofs_per_component);
+      // zero content
+      solution_renumbered_vectorized.fill(
+        typename internal::FEPointEvaluation::EvaluatorTypeTraits<
+          dim,
+          n_components,
+          VectorizedArray<Number>>::value_type());
+
       const std::size_t n_points = unit_points.size();
       const std::size_t n_lanes  = VectorizedArray<Number>::size();
       for (unsigned int i = 0; i < n_points; i += n_lanes)
@@ -829,7 +861,7 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
 
           typename internal::ProductTypeNoPoint<value_type,
                                                 VectorizedArray<Number>>::type
-            value;
+            value = {};
           Tensor<1,
                  dim,
                  typename internal::ProductTypeNoPoint<
@@ -837,7 +869,6 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
                    VectorizedArray<Number>>::type>
             gradient;
 
-          // convert back to standard format
           if (integration_flags & EvaluationFlags::values)
             for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
               internal::FEPointEvaluation::
@@ -850,22 +881,29 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
                   gradient, j, gradients[i + j]);
 
           // compute
-          internal::integrate_tensor_product_value_and_gradient(
+          internal::integrate_add_tensor_product_value_and_gradient(
             poly,
-            solution_renumbered,
             value,
             gradient,
             vectorized_points,
-            poly.size() == 2);
+            solution_renumbered_vectorized);
         }
 
+      // add between the lanes and write into the result
+      std::fill(solution_values.begin(), solution_values.end(), Number());
       for (unsigned int comp = 0; comp < n_components; ++comp)
         for (unsigned int i = 0; i < dofs_per_component; ++i)
-          internal::FEPointEvaluation::
-            EvaluatorTypeTraits<dim, n_components, Number>::write_value(
-              solution_values[renumber[comp * dofs_per_component + i]],
-              comp,
-              solution_renumbered[i]);
+          {
+            VectorizedArray<Number> result;
+            internal::FEPointEvaluation::
+              EvaluatorTypeTraits<dim, n_components, VectorizedArray<Number>>::
+                write_value(result, comp, solution_renumbered_vectorized[i]);
+            for (unsigned int lane = n_lanes / 2; lane > 0; lane /= 2)
+              for (unsigned int j = 0; j < lane; ++j)
+                result[j] += result[lane + j];
+            solution_values[renumber[comp * dofs_per_component + i]] =
+              result[0];
+          }
     }
   else if ((integration_flags & EvaluationFlags::values) ||
            (integration_flags & EvaluationFlags::gradients))
index 3a24ed6d7b37e7d07d9ae2fc194736711f47400d..c941236ea91451a8fcbd9cb858fbbea83def4058 100644 (file)
@@ -245,15 +245,38 @@ public:
    * As compared to the other transform functions that rely on pre-computed
    * information of InternalDataBase, this function chooses the flexible
    * evaluation path on the cell and points passed in to the current
-   * function. The types `Number` and `Number2` of the input and output arrays
-   * must be such that `Number2 = apply_transformation(DerivativeForm<1,
-   * spacedim, dim>, Number)`.
+   * function.
+   *
+   * @param cell The cell where to evaluate the mapping
+   *
+   * @param kind Select the kind of the mapping to be applied; currently, this
+   * class only implements `mapping_covariant`.
+   *
+   * @param apply_from_left If true, the mapping is applied to a input
+   * vector from the left, the usual choice for the transformation of field
+   * variables. If false, the mapping is applied to the vector from the right,
+   * representing the transpose of the initial operation; this is the choice
+   * to be used for implementing the action by a test function in an
+   * integration step.
+   *
+   * @param unit_points The points in reference coordinates where the
+   * transformation should be computed and applied to the vector.
+   *
+   * @param input The array of vectors (e.g., gradients for
+   * `mapping_covariant`) to be transformed.
+   *
+   * @param output The array where the result will be stored. The types
+   * `Number` and `Number2` of the input and output arrays must be such that
+   * `Number2 = apply_transformation(DerivativeForm<1, spacedim, dim>,
+   * Number)`. In case the two number types match, this array can be the same
+   * as input array.
    */
   template <typename Number, typename Number2>
   void
   transform_variable(
     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
     const MappingKind                                           kind,
+    const bool                                                  apply_from_left,
     const ArrayView<const Point<dim>> &                         unit_points,
     const ArrayView<const Number> &                             input,
     const ArrayView<Number2> &                                  output) const;
@@ -937,6 +960,7 @@ inline void
 MappingQGeneric<dim, spacedim>::transform_variable(
   const typename Triangulation<dim, spacedim>::cell_iterator &cell,
   const MappingKind                                           kind,
+  const bool                                                  apply_from_left,
   const ArrayView<const Point<dim>> &                         unit_points,
   const ArrayView<const Number> &                             input,
   const ArrayView<Number2> &                                  output) const
@@ -976,15 +1000,31 @@ MappingQGeneric<dim, spacedim>::transform_variable(
             renumber_lexicographic_to_hierarchic)
             .second;
 
-        const DerivativeForm<1, spacedim, dim, VectorizedArray<double>> jac =
-          grad.transpose().covariant_form();
-        for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+        if (apply_from_left)
           {
-            DerivativeForm<1, spacedim, dim> jac_j;
-            for (unsigned int d = 0; d < spacedim; ++d)
-              for (unsigned int e = 0; e < dim; ++e)
-                jac_j[d][e] = jac[d][e][j];
-            output[i + j] = apply_transformation(jac_j, input[i + j]);
+            const DerivativeForm<1, spacedim, dim, VectorizedArray<double>>
+              jac = grad.transpose().covariant_form();
+            for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+              {
+                DerivativeForm<1, spacedim, dim> jac_j;
+                for (unsigned int d = 0; d < spacedim; ++d)
+                  for (unsigned int e = 0; e < dim; ++e)
+                    jac_j[d][e] = jac[d][e][j];
+                output[i + j] = apply_transformation(jac_j, input[i + j]);
+              }
+          }
+        else
+          {
+            const DerivativeForm<1, dim, spacedim, VectorizedArray<double>>
+              jac = grad.covariant_form();
+            for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+              {
+                DerivativeForm<1, dim, spacedim> jac_j;
+                for (unsigned int d = 0; d < dim; ++d)
+                  for (unsigned int e = 0; e < spacedim; ++e)
+                    jac_j[d][e] = jac[d][e][j];
+                output[i + j] = apply_transformation(jac_j, input[i + j]);
+              }
           }
       }
     else
@@ -997,8 +1037,11 @@ MappingQGeneric<dim, spacedim>::transform_variable(
             polynomial_degree == 1,
             renumber_lexicographic_to_hierarchic)
             .second;
-        output[i] =
-          apply_transformation(grad.transpose().covariant_form(), input[i]);
+        if (apply_from_left)
+          output[i] =
+            apply_transformation(grad.transpose().covariant_form(), input[i]);
+        else
+          output[i] = apply_transformation(grad.covariant_form(), input[i]);
       }
 }
 
index 153b9199cb8b21d1893f3bdc6b3932a1c4e4698b..ae7eae100cccbf095501b37d534a8bc51e2dcfef 100644 (file)
@@ -2524,25 +2524,60 @@ namespace internal
    */
   template <int dim, typename Number, typename Number2>
   inline void
-  integrate_tensor_product_value_and_gradient(
-    const std::vector<Polynomials::Polynomial<double>> &      poly,
-    const std::vector<Number> &                               values,
-    const typename ProductTypeNoPoint<Number, Number2>::type &value,
-    const Tensor<1, dim, typename ProductTypeNoPoint<Number, Number2>::type>
-      &                              gradient,
-    const Point<dim, Number2> &      p,
-    const bool                       d_linear = false,
-    const std::vector<unsigned int> &renumber = {})
+  integrate_add_tensor_product_value_and_gradient(
+    const std::vector<Polynomials::Polynomial<double>> &poly,
+    const Number2 &                                     value,
+    const Tensor<1, dim, Number2> &                     gradient,
+    const Point<dim, Number> &                          p,
+    AlignedVector<Number2> &                            values,
+    const std::vector<unsigned int> &                   renumber = {})
   {
-    Assert(false, ExcNotImplemented());
-
-    (void)poly;
-    (void)values;
-    (void)value;
-    (void)gradient;
-    (void)p;
-    (void)d_linear;
-    (void)renumber;
+    static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
+
+    const unsigned 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<Number, 2 * dim * 200> shapes;
+
+    // Evaluate 1D polynomials and their derivatives
+    for (unsigned int d = 0; d < dim; ++d)
+      for (unsigned int i = 0; i < n_shapes; ++i)
+        poly[i].value(p[d], 1, shapes.data() + 2 * (d * n_shapes + i));
+
+    // Implement the transpose of the function above
+    for (unsigned int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
+      {
+        const Number2 test_value_z =
+          dim > 2 ? (value * shapes[4 * n_shapes + 2 * i2] +
+                     gradient[2] * shapes[4 * n_shapes + 2 * i2 + 1]) :
+                    value;
+        const Number2 test_grad_x =
+          dim > 2 ? gradient[0] * shapes[4 * n_shapes + 2 * i2] : gradient[0];
+        const Number2 test_grad_y =
+          dim > 2 ? gradient[1] * shapes[4 * n_shapes + 2 * i2] :
+                    (dim > 1 ? gradient[1] : Number2());
+        for (unsigned int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
+          {
+            const Number2 test_value_y =
+              dim > 1 ? (test_value_z * shapes[2 * n_shapes + 2 * i1] +
+                         test_grad_y * shapes[2 * n_shapes + 2 * i1 + 1]) :
+                        test_value_z;
+            const Number2 test_grad_xy =
+              dim > 1 ? test_grad_x * shapes[2 * n_shapes + 2 * i1] :
+                        test_grad_x;
+            if (renumber.empty())
+              for (unsigned int i0 = 0; i0 < n_shapes; ++i0, ++i)
+                values[i] += shapes[2 * i0] * test_value_y +
+                             shapes[2 * i0 + 1] * test_grad_xy;
+            else
+              for (unsigned int i0 = 0; i0 < n_shapes; ++i0, ++i)
+                values[renumber[i]] += shapes[2 * i0] * test_value_y +
+                                       shapes[2 * i0 + 1] * test_grad_xy;
+          }
+      }
   }
 
 

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.