From 0e8b9a03786983cb304f1e5f2eabceec1c829baa Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Wed, 3 Apr 2024 22:46:41 +0200 Subject: [PATCH] MappingQ: Use get_vertices() instead of compute_mapping_support_points when possible --- include/deal.II/fe/mapping_q_internal.h | 4 +- source/fe/mapping_q.cc | 182 +++++++++++++++++++----- 2 files changed, 147 insertions(+), 39 deletions(-) diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h index d2e8cd9053..bc4e2e4fe9 100644 --- a/include/deal.II/fe/mapping_q_internal.h +++ b/include/deal.II/fe/mapping_q_internal.h @@ -833,8 +833,8 @@ namespace internal * points in real space by a polynomial map. */ InverseQuadraticApproximation( - const std::vector> &real_support_points, - const std::vector> &unit_support_points) + const ArrayView> &real_support_points, + const std::vector> &unit_support_points) : normalization_shift(real_support_points[0]) , normalization_length( 1. / real_support_points[0].distance(real_support_points[1])) diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 38f5f2d7ae..5b77e81603 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -290,12 +290,19 @@ MappingQ::transform_unit_to_real_cell( const typename Triangulation::cell_iterator &cell, const Point &p) const { - return Point(internal::evaluate_tensor_product_value( - polynomials_1d, - make_const_array_view(this->compute_mapping_support_points(cell)), - p, - polynomials_1d.size() == 2, - renumber_lexicographic_to_hierarchic)); + if (polynomial_degree == 1) + { + const auto vertices = this->get_vertices(cell); + return Point( + internal::evaluate_tensor_product_value_linear(vertices.data(), p)); + } + else + return Point(internal::evaluate_tensor_product_value( + polynomials_1d, + make_const_array_view(this->compute_mapping_support_points(cell)), + p, + polynomials_1d.size() == 2, + renumber_lexicographic_to_hierarchic)); } @@ -341,13 +348,25 @@ MappingQ<1, 1>::transform_real_to_unit_cell_internal( { // dispatch to the various specializations for spacedim=dim, // spacedim=dim+1, etc - return internal::MappingQImplementation:: - do_transform_real_to_unit_cell_internal<1>( - p, - initial_p_unit, - make_const_array_view(this->compute_mapping_support_points(cell)), - polynomials_1d, - renumber_lexicographic_to_hierarchic); + if (polynomial_degree == 1) + { + const auto vertices = this->get_vertices(cell); + return internal::MappingQImplementation:: + do_transform_real_to_unit_cell_internal<1>( + p, + initial_p_unit, + ArrayView>(vertices.data(), vertices.size()), + polynomials_1d, + renumber_lexicographic_to_hierarchic); + } + else + return internal::MappingQImplementation:: + do_transform_real_to_unit_cell_internal<1>( + p, + initial_p_unit, + make_const_array_view(this->compute_mapping_support_points(cell)), + polynomials_1d, + renumber_lexicographic_to_hierarchic); } @@ -359,13 +378,25 @@ MappingQ<2, 2>::transform_real_to_unit_cell_internal( const Point<2> &p, const Point<2> &initial_p_unit) const { - return internal::MappingQImplementation:: - do_transform_real_to_unit_cell_internal<2>( - p, - initial_p_unit, - make_const_array_view(this->compute_mapping_support_points(cell)), - polynomials_1d, - renumber_lexicographic_to_hierarchic); + if (polynomial_degree == 1) + { + const auto vertices = this->get_vertices(cell); + return internal::MappingQImplementation:: + do_transform_real_to_unit_cell_internal<2>( + p, + initial_p_unit, + ArrayView>(vertices.data(), vertices.size()), + polynomials_1d, + renumber_lexicographic_to_hierarchic); + } + else + return internal::MappingQImplementation:: + do_transform_real_to_unit_cell_internal<2>( + p, + initial_p_unit, + make_const_array_view(this->compute_mapping_support_points(cell)), + polynomials_1d, + renumber_lexicographic_to_hierarchic); } @@ -377,13 +408,25 @@ MappingQ<3, 3>::transform_real_to_unit_cell_internal( const Point<3> &p, const Point<3> &initial_p_unit) const { - return internal::MappingQImplementation:: - do_transform_real_to_unit_cell_internal<3>( - p, - initial_p_unit, - make_const_array_view(this->compute_mapping_support_points(cell)), - polynomials_1d, - renumber_lexicographic_to_hierarchic); + if (polynomial_degree == 1) + { + const auto vertices = this->get_vertices(cell); + return internal::MappingQImplementation:: + do_transform_real_to_unit_cell_internal<3>( + p, + initial_p_unit, + ArrayView>(vertices.data(), vertices.size()), + polynomials_1d, + renumber_lexicographic_to_hierarchic); + } + else + return internal::MappingQImplementation:: + do_transform_real_to_unit_cell_internal<3>( + p, + initial_p_unit, + make_const_array_view(this->compute_mapping_support_points(cell)), + polynomials_1d, + renumber_lexicographic_to_hierarchic); } @@ -475,7 +518,7 @@ MappingQ::transform_real_to_unit_cell( { // Use an exact formula if one is available. this is only the case // for Q1 mappings in 1d, and in 2d if dim==spacedim - if (this->preserves_vertex_locations() && (polynomial_degree == 1) && + if ((polynomial_degree == 1) && ((dim == 1) || ((dim == 2) && (dim == spacedim)))) { // The dimension-dependent algorithms are much faster (about 25-45x in @@ -605,8 +648,18 @@ MappingQ::transform_points_real_to_unit_cell( } AssertDimension(real_points.size(), unit_points.size()); - const std::vector> support_points = - this->compute_mapping_support_points(cell); + std::vector> support_points_higher_order; + boost::container::small_vector, + GeometryInfo::vertices_per_cell> + vertices; + if (polynomial_degree == 1) + vertices = this->get_vertices(cell); + else + support_points_higher_order = this->compute_mapping_support_points(cell); + const ArrayView> support_points( + polynomial_degree == 1 ? vertices.data() : + support_points_higher_order.data(), + Utilities::pow(polynomial_degree + 1, dim)); // From the given (high-order) support points, now only pick the first // 2^dim points and construct an affine approximation from those. @@ -808,7 +861,16 @@ MappingQ::fill_fe_values( // object attached to the cell and all of its bounding faces/edges, // etc. to reliably test that the "cell" we are on is, therefore, // not easily done - data.mapping_support_points = this->compute_mapping_support_points(cell); + if (polynomial_degree == 1) + { + data.mapping_support_points.resize(GeometryInfo::vertices_per_cell); + const auto vertices = this->get_vertices(cell); + for (unsigned int i = 0; i < GeometryInfo::vertices_per_cell; ++i) + data.mapping_support_points[i] = vertices[i]; + } + else + data.mapping_support_points = this->compute_mapping_support_points(cell); + data.cell_of_current_support_points = cell; // if the order of the mapping is greater than 1, then do not reuse any cell @@ -1026,7 +1088,18 @@ MappingQ::fill_fe_face_values( &data.cell_of_current_support_points->get_triangulation()) || (cell != data.cell_of_current_support_points)) { - data.mapping_support_points = this->compute_mapping_support_points(cell); + if (polynomial_degree == 1) + { + data.mapping_support_points.resize( + GeometryInfo::vertices_per_cell); + const auto vertices = this->get_vertices(cell); + for (unsigned int i = 0; i < GeometryInfo::vertices_per_cell; + ++i) + data.mapping_support_points[i] = vertices[i]; + } + else + data.mapping_support_points = + this->compute_mapping_support_points(cell); data.cell_of_current_support_points = cell; } @@ -1077,7 +1150,18 @@ MappingQ::fill_fe_subface_values( &data.cell_of_current_support_points->get_triangulation()) || (cell != data.cell_of_current_support_points)) { - data.mapping_support_points = this->compute_mapping_support_points(cell); + if (polynomial_degree == 1) + { + data.mapping_support_points.resize( + GeometryInfo::vertices_per_cell); + const auto vertices = this->get_vertices(cell); + for (unsigned int i = 0; i < GeometryInfo::vertices_per_cell; + ++i) + data.mapping_support_points[i] = vertices[i]; + } + else + data.mapping_support_points = + this->compute_mapping_support_points(cell); data.cell_of_current_support_points = cell; } @@ -1121,7 +1205,15 @@ MappingQ::fill_fe_immersed_surface_values( const unsigned int n_q_points = quadrature.size(); - data.mapping_support_points = this->compute_mapping_support_points(cell); + if (polynomial_degree == 1) + { + data.mapping_support_points.resize(GeometryInfo::vertices_per_cell); + const auto vertices = this->get_vertices(cell); + for (unsigned int i = 0; i < GeometryInfo::vertices_per_cell; ++i) + data.mapping_support_points[i] = vertices[i]; + } + else + data.mapping_support_points = this->compute_mapping_support_points(cell); data.cell_of_current_support_points = cell; internal::MappingQImplementation::maybe_update_q_points_Jacobians_generic( @@ -1258,7 +1350,15 @@ MappingQ::fill_mapping_data_for_generic_points( unit_points.end()))); const InternalData &data = static_cast(*internal_data); data.output_data = &output_data; - data.mapping_support_points = this->compute_mapping_support_points(cell); + if (polynomial_degree == 1) + { + data.mapping_support_points.resize(GeometryInfo::vertices_per_cell); + const auto vertices = this->get_vertices(cell); + for (unsigned int i = 0; i < GeometryInfo::vertices_per_cell; ++i) + data.mapping_support_points[i] = vertices[i]; + } + else + data.mapping_support_points = this->compute_mapping_support_points(cell); internal::MappingQImplementation::maybe_update_q_points_Jacobians_generic( CellSimilarity::none, @@ -1291,8 +1391,16 @@ MappingQ::fill_mapping_data_for_face_quadrature( ExcInternalError()); const InternalData &data = static_cast(internal_data); - data.mapping_support_points = this->compute_mapping_support_points(cell); - data.output_data = &output_data; + if (polynomial_degree == 1) + { + data.mapping_support_points.resize(GeometryInfo::vertices_per_cell); + const auto vertices = this->get_vertices(cell); + for (unsigned int i = 0; i < GeometryInfo::vertices_per_cell; ++i) + data.mapping_support_points[i] = vertices[i]; + } + else + data.mapping_support_points = this->compute_mapping_support_points(cell); + data.output_data = &output_data; internal::MappingQImplementation::do_fill_fe_face_values( *this, -- 2.39.5