From 8da044893752e30bb0d47dd5e60c3ac7d48c95f2 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Wed, 3 Apr 2024 13:39:35 +0200 Subject: [PATCH] Use ArrayView for interface to tensor-product point evaluator --- include/deal.II/fe/mapping_q_internal.h | 34 +++++++++---------- .../tensor_product_point_kernels.h | 12 +++---- source/fe/mapping_q.cc | 12 +++---- source/non_matching/quadrature_generator.cc | 15 ++++---- .../mapping_q_real_to_unit_internal.cc | 2 +- .../matrix_free/tensor_product_evaluate_01.cc | 12 +++++-- .../matrix_free/tensor_product_evaluate_02.cc | 12 +++++-- .../matrix_free/tensor_product_evaluate_03.cc | 12 +++++-- .../matrix_free/tensor_product_evaluate_04.cc | 4 +-- .../matrix_free/tensor_product_evaluate_05.cc | 6 ++-- .../matrix_free/tensor_product_evaluate_06.cc | 5 ++- .../matrix_free/tensor_product_evaluate_07.cc | 2 +- .../matrix_free/tensor_product_evaluate_08.cc | 2 +- 13 files changed, 77 insertions(+), 53 deletions(-) diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h index 90ba787dbc..d2e8cd9053 100644 --- a/include/deal.II/fe/mapping_q_internal.h +++ b/include/deal.II/fe/mapping_q_internal.h @@ -477,7 +477,7 @@ namespace internal do_transform_real_to_unit_cell_internal( const Point &p, const Point &initial_p_unit, - const std::vector> &points, + const ArrayView> &points, const std::vector> &polynomials_1d, const std::vector &renumber, const bool print_iterations_to_deallog = false) @@ -734,7 +734,7 @@ namespace internal do_transform_real_to_unit_cell_internal_codim1( const Point &p, const Point &initial_p_unit, - const std::vector> &points, + const ArrayView> &points, const std::vector> &polynomials_1d, const std::vector &renumber) { @@ -1245,9 +1245,9 @@ namespace internal std::vector> &jacobians, std::vector> &inverse_jacobians) { - const UpdateFlags update_flags = data.update_each; - const std::vector> &support_points = - data.mapping_support_points; + const UpdateFlags update_flags = data.update_each; + const ArrayView> support_points( + data.mapping_support_points); const unsigned int n_points = unit_points.size(); const unsigned int n_lanes = VectorizedArray::size(); @@ -1406,8 +1406,8 @@ namespace internal { if (data.update_each & update_jacobian_grads) { - const std::vector> &support_points = - data.mapping_support_points; + const ArrayView> support_points( + data.mapping_support_points); const unsigned int n_q_points = jacobian_grads.size(); if (cell_similarity != CellSimilarity::translation) @@ -1448,8 +1448,8 @@ namespace internal { if (data.update_each & update_jacobian_pushed_forward_grads) { - const std::vector> &support_points = - data.mapping_support_points; + const ArrayView> support_points( + data.mapping_support_points); const unsigned int n_q_points = jacobian_pushed_forward_grads.size(); if (cell_similarity != CellSimilarity::translation) @@ -1563,8 +1563,8 @@ namespace internal { if (data.update_each & update_jacobian_2nd_derivatives) { - const std::vector> &support_points = - data.mapping_support_points; + const ArrayView> support_points( + data.mapping_support_points); const unsigned int n_q_points = jacobian_2nd_derivatives.size(); if (cell_similarity != CellSimilarity::translation) @@ -1603,8 +1603,8 @@ namespace internal { if (data.update_each & update_jacobian_pushed_forward_2nd_derivatives) { - const std::vector> &support_points = - data.mapping_support_points; + const ArrayView> support_points( + data.mapping_support_points); const unsigned int n_q_points = jacobian_pushed_forward_2nd_derivatives.size(); @@ -1750,8 +1750,8 @@ namespace internal { if (data.update_each & update_jacobian_3rd_derivatives) { - const std::vector> &support_points = - data.mapping_support_points; + const ArrayView> support_points( + data.mapping_support_points); const unsigned int n_q_points = jacobian_3rd_derivatives.size(); if (cell_similarity != CellSimilarity::translation) @@ -1790,8 +1790,8 @@ namespace internal { if (data.update_each & update_jacobian_pushed_forward_3rd_derivatives) { - const std::vector> &support_points = - data.mapping_support_points; + const ArrayView> support_points( + data.mapping_support_points); const unsigned int n_q_points = jacobian_pushed_forward_3rd_derivatives.size(); diff --git a/include/deal.II/matrix_free/tensor_product_point_kernels.h b/include/deal.II/matrix_free/tensor_product_point_kernels.h index d6dbfccaca..0fed110392 100644 --- a/include/deal.II/matrix_free/tensor_product_point_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_point_kernels.h @@ -461,7 +461,7 @@ namespace internal Tensor<1, dim, typename ProductTypeNoPoint::type>> evaluate_tensor_product_value_and_gradient( const std::vector> &poly, - const std::vector &values, + const ArrayView &values, const Point &p, const bool d_linear = false, const std::vector &renumber = {}) @@ -682,7 +682,7 @@ namespace internal inline typename ProductTypeNoPoint::type evaluate_tensor_product_value( const std::vector> &poly, - const std::vector &values, + const ArrayView &values, const Point &p, const bool d_linear = false, const std::vector &renumber = {}) @@ -718,7 +718,7 @@ namespace internal inline Tensor<1, 1, typename ProductTypeNoPoint::type> evaluate_tensor_product_higher_derivatives( const std::vector> &poly, - const std::vector &values, + const ArrayView &values, const Point<1, Number2> &p, const std::vector &renumber = {}) { @@ -758,7 +758,7 @@ namespace internal typename ProductTypeNoPoint::type> evaluate_tensor_product_higher_derivatives( const std::vector> &poly, - const std::vector &values, + const ArrayView &values, const Point<2, Number2> &p, const std::vector &renumber = {}) { @@ -810,7 +810,7 @@ namespace internal typename ProductTypeNoPoint::type> evaluate_tensor_product_higher_derivatives( const std::vector> &poly, - const std::vector &values, + const ArrayView &values, const Point<3, Number2> &p, const std::vector &renumber = {}) { @@ -874,7 +874,7 @@ namespace internal SymmetricTensor<2, dim, typename ProductTypeNoPoint::type> evaluate_tensor_product_hessian( const std::vector> &poly, - const std::vector &values, + const ArrayView &values, const Point &p, const std::vector &renumber = {}) { diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index f1631738ad..38f5f2d7ae 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -292,7 +292,7 @@ MappingQ::transform_unit_to_real_cell( { return Point(internal::evaluate_tensor_product_value( polynomials_1d, - this->compute_mapping_support_points(cell), + make_const_array_view(this->compute_mapping_support_points(cell)), p, polynomials_1d.size() == 2, renumber_lexicographic_to_hierarchic)); @@ -345,7 +345,7 @@ MappingQ<1, 1>::transform_real_to_unit_cell_internal( do_transform_real_to_unit_cell_internal<1>( p, initial_p_unit, - this->compute_mapping_support_points(cell), + make_const_array_view(this->compute_mapping_support_points(cell)), polynomials_1d, renumber_lexicographic_to_hierarchic); } @@ -363,7 +363,7 @@ MappingQ<2, 2>::transform_real_to_unit_cell_internal( do_transform_real_to_unit_cell_internal<2>( p, initial_p_unit, - this->compute_mapping_support_points(cell), + make_const_array_view(this->compute_mapping_support_points(cell)), polynomials_1d, renumber_lexicographic_to_hierarchic); } @@ -381,7 +381,7 @@ MappingQ<3, 3>::transform_real_to_unit_cell_internal( do_transform_real_to_unit_cell_internal<3>( p, initial_p_unit, - this->compute_mapping_support_points(cell), + make_const_array_view(this->compute_mapping_support_points(cell)), polynomials_1d, renumber_lexicographic_to_hierarchic); } @@ -414,7 +414,7 @@ MappingQ<1, 2>::transform_real_to_unit_cell_internal( do_transform_real_to_unit_cell_internal_codim1<1>( p, initial_p_unit, - mdata->mapping_support_points, + make_const_array_view(mdata->mapping_support_points), polynomials_1d, renumber_lexicographic_to_hierarchic); } @@ -447,7 +447,7 @@ MappingQ<2, 3>::transform_real_to_unit_cell_internal( do_transform_real_to_unit_cell_internal_codim1<2>( p, initial_p_unit, - mdata->mapping_support_points, + make_const_array_view(mdata->mapping_support_points), polynomials_1d, renumber_lexicographic_to_hierarchic); } diff --git a/source/non_matching/quadrature_generator.cc b/source/non_matching/quadrature_generator.cc index 3dc9d32c95..0558999688 100644 --- a/source/non_matching/quadrature_generator.cc +++ b/source/non_matching/quadrature_generator.cc @@ -1631,12 +1631,12 @@ namespace NonMatching return this->is_fe_q_iso_q1() ? dealii::internal::evaluate_tensor_product_value( poly, - local_dof_values_subcell, + make_array_view(local_dof_values_subcell), subcell_box.real_to_unit(point), polynomials_are_hat_functions) : dealii::internal::evaluate_tensor_product_value( poly, - local_dof_values, + make_array_view(local_dof_values), point, polynomials_are_hat_functions, renumber); @@ -1670,13 +1670,13 @@ namespace NonMatching dealii::internal:: evaluate_tensor_product_value_and_gradient( poly, - local_dof_values_subcell, + make_array_view(local_dof_values_subcell), subcell_box.real_to_unit(point), polynomials_are_hat_functions) : dealii::internal:: evaluate_tensor_product_value_and_gradient( poly, - local_dof_values, + make_array_view(local_dof_values), point, polynomials_are_hat_functions, renumber)) @@ -1710,10 +1710,13 @@ namespace NonMatching return this->is_fe_q_iso_q1() ? dealii::internal::evaluate_tensor_product_hessian( poly, - local_dof_values_subcell, + make_array_view(local_dof_values_subcell), subcell_box.real_to_unit(point)) : dealii::internal::evaluate_tensor_product_hessian( - poly, local_dof_values, point, renumber); + poly, + make_array_view(local_dof_values), + point, + renumber); } else { diff --git a/tests/mappings/mapping_q_real_to_unit_internal.cc b/tests/mappings/mapping_q_real_to_unit_internal.cc index d1e5688662..5e02266e61 100644 --- a/tests/mappings/mapping_q_real_to_unit_internal.cc +++ b/tests/mappings/mapping_q_real_to_unit_internal.cc @@ -74,7 +74,7 @@ print_result(const unsigned int mapping_degree, do_transform_real_to_unit_cell_internal( p, cell->real_to_unit_cell_affine_approximation(p), - fe_values.get_quadrature_points(), + make_array_view(fe_values.get_quadrature_points()), polynomials, renumber, /* print_iterations = */ true); diff --git a/tests/matrix_free/tensor_product_evaluate_01.cc b/tests/matrix_free/tensor_product_evaluate_01.cc index 952f1e5fec..fcc477d74a 100644 --- a/tests/matrix_free/tensor_product_evaluate_01.cc +++ b/tests/matrix_free/tensor_product_evaluate_01.cc @@ -71,7 +71,11 @@ test(const unsigned int degree) for (const auto &p : evaluation_points) { const auto val = internal::evaluate_tensor_product_value_and_gradient( - polynomials, coefficients, p, false, renumbering); + polynomials, + make_const_array_view(coefficients), + p, + false, + renumbering); deallog << "Value " << val.first << " vs " << transform * (matrix * p + offset) << " ; gradient " << val.second << " vs " << transform * matrix << std::endl; @@ -83,7 +87,11 @@ test(const unsigned int degree) for (const auto &p : evaluation_points) { const auto val = internal::evaluate_tensor_product_value_and_gradient( - polynomials, coefficients, p, true, renumbering); + polynomials, + make_const_array_view(coefficients), + p, + true, + renumbering); deallog << "Value " << val.first << " vs " << transform * (matrix * p + offset) << " ; gradient " << val.second << " vs " << transform * matrix << std::endl; diff --git a/tests/matrix_free/tensor_product_evaluate_02.cc b/tests/matrix_free/tensor_product_evaluate_02.cc index eeb7f37f2c..d41525d234 100644 --- a/tests/matrix_free/tensor_product_evaluate_02.cc +++ b/tests/matrix_free/tensor_product_evaluate_02.cc @@ -66,7 +66,11 @@ test(const unsigned int degree) for (const auto &p : evaluation_points) { const auto val = internal::evaluate_tensor_product_value_and_gradient( - polynomials, coefficients, p, false, renumbering); + polynomials, + make_const_array_view(coefficients), + p, + false, + renumbering); deallog << "Value " << val.first << " vs " << matrix * p + offset << std::endl; deallog << "Gradient " << val.second << " vs " << transpose(matrix) @@ -79,7 +83,11 @@ test(const unsigned int degree) for (const auto &p : evaluation_points) { const auto val = internal::evaluate_tensor_product_value_and_gradient( - polynomials, coefficients, p, true, renumbering); + polynomials, + make_const_array_view(coefficients), + p, + true, + renumbering); deallog << "Value " << val.first << " vs " << matrix * p + offset << std::endl; deallog << "Gradient " << val.second << " vs " << transpose(matrix) diff --git a/tests/matrix_free/tensor_product_evaluate_03.cc b/tests/matrix_free/tensor_product_evaluate_03.cc index 2e97ca301d..7507516f06 100644 --- a/tests/matrix_free/tensor_product_evaluate_03.cc +++ b/tests/matrix_free/tensor_product_evaluate_03.cc @@ -72,7 +72,11 @@ test(const unsigned int degree) p_vec[d][v] = p[d] + 0.01 * v; const auto val = internal::evaluate_tensor_product_value_and_gradient( - polynomials, coefficients, p_vec, false, renumbering); + polynomials, + make_const_array_view(coefficients), + p_vec, + false, + renumbering); const auto error_vec = val.first - matrix * p_vec; double error = 0; @@ -100,7 +104,11 @@ test(const unsigned int degree) p_vec[d][v] = p[d] + 0.01 * v; const auto val = internal::evaluate_tensor_product_value_and_gradient( - polynomials, coefficients, p_vec, true, renumbering); + polynomials, + make_const_array_view(coefficients), + p_vec, + true, + renumbering); const auto error_vec = val.first - matrix * p_vec; double error = 0; diff --git a/tests/matrix_free/tensor_product_evaluate_04.cc b/tests/matrix_free/tensor_product_evaluate_04.cc index b27fca1b43..e2b4c19598 100644 --- a/tests/matrix_free/tensor_product_evaluate_04.cc +++ b/tests/matrix_free/tensor_product_evaluate_04.cc @@ -69,7 +69,7 @@ test(const unsigned int degree) p_vec[d][v] = p[d] + 0.01 * v; const auto val = internal::evaluate_tensor_product_value_and_gradient( - polynomials, coefficients, p_vec, false); + polynomials, make_const_array_view(coefficients), p_vec, false); const auto error_vec = val.first - matrix * p_vec; double error = 0; @@ -97,7 +97,7 @@ test(const unsigned int degree) p_vec[d][v] = p[d] + 0.01 * v; const auto val = internal::evaluate_tensor_product_value_and_gradient( - polynomials, coefficients, p_vec, true); + polynomials, make_const_array_view(coefficients), p_vec, true); const auto error_vec = val.first - matrix * p_vec; double error = 0; diff --git a/tests/matrix_free/tensor_product_evaluate_05.cc b/tests/matrix_free/tensor_product_evaluate_05.cc index a3f784dd3a..d1c4c8685a 100644 --- a/tests/matrix_free/tensor_product_evaluate_05.cc +++ b/tests/matrix_free/tensor_product_evaluate_05.cc @@ -69,10 +69,8 @@ test(const unsigned int degree) for (unsigned int d = 0; d < dim; ++d) p_vec[d][v] = p[d] + 0.01 * v; - const auto hess = internal::evaluate_tensor_product_hessian(polynomials, - coefficients, - p_vec, - renumbering); + const auto hess = internal::evaluate_tensor_product_hessian( + polynomials, make_const_array_view(coefficients), p_vec, renumbering); double error = 0; for (unsigned int v = 0; v < VectorizedArray::size(); ++v) diff --git a/tests/matrix_free/tensor_product_evaluate_06.cc b/tests/matrix_free/tensor_product_evaluate_06.cc index ae8a9aa6fc..754a53d486 100644 --- a/tests/matrix_free/tensor_product_evaluate_06.cc +++ b/tests/matrix_free/tensor_product_evaluate_06.cc @@ -65,9 +65,8 @@ test(const unsigned int degree) for (unsigned int d = 0; d < dim; ++d) p_vec[d][v] = p[d] + 0.01 * v; - const auto hess = internal::evaluate_tensor_product_hessian(polynomials, - coefficients, - p_vec); + const auto hess = internal::evaluate_tensor_product_hessian( + polynomials, make_const_array_view(coefficients), p_vec); std::cout << hess << " " << matrix << std::endl; diff --git a/tests/matrix_free/tensor_product_evaluate_07.cc b/tests/matrix_free/tensor_product_evaluate_07.cc index ab3e09b33d..2dc825b1d0 100644 --- a/tests/matrix_free/tensor_product_evaluate_07.cc +++ b/tests/matrix_free/tensor_product_evaluate_07.cc @@ -75,7 +75,7 @@ test(const unsigned int degree) deallog << "]: "; const auto derivative = internal::evaluate_tensor_product_higher_derivatives( - polynomials, function_values, p, renumbering); + polynomials, make_const_array_view(function_values), p, renumbering); for (unsigned int d = 0; d < derivative.dimension; ++d) deallog << (std::abs(derivative[d]) < 1e-11 ? 0. : derivative[d]) diff --git a/tests/matrix_free/tensor_product_evaluate_08.cc b/tests/matrix_free/tensor_product_evaluate_08.cc index 3b9fad9c21..b561a04cbc 100644 --- a/tests/matrix_free/tensor_product_evaluate_08.cc +++ b/tests/matrix_free/tensor_product_evaluate_08.cc @@ -72,7 +72,7 @@ test(const unsigned int degree) deallog << "]: "; const auto derivative = internal::evaluate_tensor_product_higher_derivatives( - polynomials, function_values, p); + polynomials, make_const_array_view(function_values), p); for (unsigned int d = 0; d < derivative.dimension; ++d) deallog << (std::abs(derivative[d]) < 1e-11 ? 0. : derivative[d]) -- 2.39.5