Skip to content

Commit

Permalink
Merge pull request #16847 from kronbichler/avoid_allocations
Browse files Browse the repository at this point in the history 
MappingQ: Favor get_vertices() over compute_mapping_support_points() if possible
  • Loading branch information
@kronbichler
kronbichler committed Apr 4, 2024
2 parents 809ce96 + 0e8b9a0 commit f844528
Show file tree
Hide file tree
Showing 2 changed files with 147 additions and 39 deletions.
4 changes: 2 additions & 2 deletions include/deal.II/fe/mapping_q_internal.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -833,8 +833,8 @@ namespace internal
* points in real space by a polynomial map.
*/
InverseQuadraticApproximation(
const std::vector<Point<spacedim>> &real_support_points,
const std::vector<Point<dim>> &unit_support_points)
const ArrayView<const Point<spacedim>> &real_support_points,
const std::vector<Point<dim>> &unit_support_points)
: normalization_shift(real_support_points[0])
, normalization_length(
1. / real_support_points[0].distance(real_support_points[1]))
Expand Down
182 changes: 145 additions & 37 deletions source/fe/mapping_q.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -290,12 +290,19 @@ MappingQ<dim, spacedim>::transform_unit_to_real_cell(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const Point<dim> &p) const
{
return Point<spacedim>(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<spacedim>(
internal::evaluate_tensor_product_value_linear(vertices.data(), p));
}
else
return Point<spacedim>(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));
}


Expand Down Expand Up @@ -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<const Point<1>>(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);
}


Expand All @@ -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<const Point<2>>(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);
}


Expand All @@ -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<const Point<3>>(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);
}


Expand Down Expand Up @@ -475,7 +518,7 @@ MappingQ<dim, spacedim>::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
Expand Down Expand Up @@ -605,8 +648,18 @@ MappingQ<dim, spacedim>::transform_points_real_to_unit_cell(
}

AssertDimension(real_points.size(), unit_points.size());
const std::vector<Point<spacedim>> support_points =
this->compute_mapping_support_points(cell);
std::vector<Point<spacedim>> support_points_higher_order;
boost::container::small_vector<Point<spacedim>,
GeometryInfo<dim>::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<const Point<spacedim>> 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.
Expand Down Expand Up @@ -808,7 +861,16 @@ MappingQ<dim, spacedim>::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<dim>::vertices_per_cell);
const auto vertices = this->get_vertices(cell);
for (unsigned int i = 0; i < GeometryInfo<dim>::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
Expand Down Expand Up @@ -1026,7 +1088,18 @@ MappingQ<dim, spacedim>::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<dim>::vertices_per_cell);
const auto vertices = this->get_vertices(cell);
for (unsigned int i = 0; i < GeometryInfo<dim>::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;
}

Expand Down Expand Up @@ -1077,7 +1150,18 @@ MappingQ<dim, spacedim>::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<dim>::vertices_per_cell);
const auto vertices = this->get_vertices(cell);
for (unsigned int i = 0; i < GeometryInfo<dim>::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;
}

Expand Down Expand Up @@ -1121,7 +1205,15 @@ MappingQ<dim, spacedim>::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<dim>::vertices_per_cell);
const auto vertices = this->get_vertices(cell);
for (unsigned int i = 0; i < GeometryInfo<dim>::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(
Expand Down Expand Up @@ -1258,7 +1350,15 @@ MappingQ<dim, spacedim>::fill_mapping_data_for_generic_points(
unit_points.end())));
const InternalData &data = static_cast<const InternalData &>(*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<dim>::vertices_per_cell);
const auto vertices = this->get_vertices(cell);
for (unsigned int i = 0; i < GeometryInfo<dim>::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,
Expand Down Expand Up @@ -1291,8 +1391,16 @@ MappingQ<dim, spacedim>::fill_mapping_data_for_face_quadrature(
ExcInternalError());
const InternalData &data = static_cast<const InternalData &>(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<dim>::vertices_per_cell);
const auto vertices = this->get_vertices(cell);
for (unsigned int i = 0; i < GeometryInfo<dim>::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,
Expand Down

0 comments on commit f844528

Please sign in to comment.