From 1c730e86d1f542eb8b889372c1cd5c55d6be85af Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 28 Jan 2022 14:52:08 +0100 Subject: [PATCH] Speed up evaluation of RefSpaceFEFieldFunction --- .../matrix_free/tensor_product_kernels.h | 99 ++++++++++++ source/non_matching/fe_values.cc | 148 ++++++++++++++---- source/non_matching/quadrature_generator.cc | 2 +- 3 files changed, 216 insertions(+), 33 deletions(-) diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 5948e7c355..c09b652bb3 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -2523,6 +2523,105 @@ namespace internal + template + SymmetricTensor<2, dim, typename ProductTypeNoPoint::type> + evaluate_tensor_product_hessian( + const std::vector> &poly, + const std::vector & values, + const Point & p, + const std::vector & renumber = {}) + { + static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented"); + + using Number3 = typename ProductTypeNoPoint::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 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. */ diff --git a/source/non_matching/fe_values.cc b/source/non_matching/fe_values.cc index 2b26e6aa85..d4e0bcd40a 100644 --- a/source/non_matching/fe_values.cc +++ b/source/non_matching/fe_values.cc @@ -28,6 +28,8 @@ #include #include +#include + #include DEAL_II_NAMESPACE_OPEN @@ -163,6 +165,23 @@ namespace NonMatching * set_active_cell() */ std::vector 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> poly; + + /** + * Renumbering for the tensor-product evaluator in the fast path. + */ + std::vector 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 + 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 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 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], diff --git a/source/non_matching/quadrature_generator.cc b/source/non_matching/quadrature_generator.cc index 0775c86910..9f446a375d 100644 --- a/source/non_matching/quadrature_generator.cc +++ b/source/non_matching/quadrature_generator.cc @@ -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; -- 2.39.5