From cb932f3b6f46a3f6630872b53df2bc760587f026 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 17 May 2021 19:41:57 +0200 Subject: [PATCH] Optimize FEPointEvaluation::integrate() with tensor product path --- include/deal.II/fe/fe_point_evaluation.h | 88 +++++++++++++------ include/deal.II/fe/mapping_q_generic.h | 69 ++++++++++++--- .../matrix_free/tensor_product_kernels.h | 71 +++++++++++---- 3 files changed, 172 insertions(+), 56 deletions(-) diff --git a/include/deal.II/fe/fe_point_evaluation.h b/include/deal.II/fe/fe_point_evaluation.h index 552e3bf7e1..75213d9640 100644 --- a/include/deal.II/fe/fe_point_evaluation.h +++ b/include/deal.II/fe/fe_point_evaluation.h @@ -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 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>` format. + */ + AlignedVector>::value_type> + solution_renumbered_vectorized; + /** * Temporary array to store the values at the points. */ @@ -606,7 +631,7 @@ FEPointEvaluation::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::evaluate( mapping_q_generic->transform_variable( cell, mapping_covariant, + /* apply_from_left */ true, unit_points, ArrayView(gradients.data(), gradients.size()), ArrayView(gradients.data(), gradients.size())); @@ -789,17 +815,13 @@ FEPointEvaluation::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::integrate( mapping_q_generic->transform_variable( cell, mapping_covariant, + /* apply_from_left */ false, unit_points, ArrayView(gradients.data(), gradients.size()), ArrayView(gradients.data(), gradients.size())); @@ -817,6 +840,15 @@ FEPointEvaluation::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>::value_type()); + const std::size_t n_points = unit_points.size(); const std::size_t n_lanes = VectorizedArray::size(); for (unsigned int i = 0; i < n_points; i += n_lanes) @@ -829,7 +861,7 @@ FEPointEvaluation::integrate( typename internal::ProductTypeNoPoint>::type - value; + value = {}; Tensor<1, dim, typename internal::ProductTypeNoPoint< @@ -837,7 +869,6 @@ FEPointEvaluation::integrate( VectorizedArray>::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::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::write_value( - solution_values[renumber[comp * dofs_per_component + i]], - comp, - solution_renumbered[i]); + { + VectorizedArray result; + internal::FEPointEvaluation:: + EvaluatorTypeTraits>:: + 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)) diff --git a/include/deal.II/fe/mapping_q_generic.h b/include/deal.II/fe/mapping_q_generic.h index 3a24ed6d7b..c941236ea9 100644 --- a/include/deal.II/fe/mapping_q_generic.h +++ b/include/deal.II/fe/mapping_q_generic.h @@ -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 void transform_variable( const typename Triangulation::cell_iterator &cell, const MappingKind kind, + const bool apply_from_left, const ArrayView> & unit_points, const ArrayView & input, const ArrayView & output) const; @@ -937,6 +960,7 @@ inline void MappingQGeneric::transform_variable( const typename Triangulation::cell_iterator &cell, const MappingKind kind, + const bool apply_from_left, const ArrayView> & unit_points, const ArrayView & input, const ArrayView & output) const @@ -976,15 +1000,31 @@ MappingQGeneric::transform_variable( renumber_lexicographic_to_hierarchic) .second; - const DerivativeForm<1, spacedim, dim, VectorizedArray> 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> + 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> + 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::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]); } } diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 153b9199cb..ae7eae100c 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -2524,25 +2524,60 @@ namespace internal */ template inline void - integrate_tensor_product_value_and_gradient( - const std::vector> & poly, - const std::vector & values, - const typename ProductTypeNoPoint::type &value, - const Tensor<1, dim, typename ProductTypeNoPoint::type> - & gradient, - const Point & p, - const bool d_linear = false, - const std::vector &renumber = {}) + integrate_add_tensor_product_value_and_gradient( + const std::vector> &poly, + const Number2 & value, + const Tensor<1, dim, Number2> & gradient, + const Point & p, + AlignedVector & values, + const std::vector & 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 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; + } + } } -- 2.39.5