From 5db8fdad335a2fc35fac338ae5e422cae3943302 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sun, 9 Apr 2023 20:44:53 +0200 Subject: [PATCH] Optimize FlatManifold::get_new_points --- source/grid/manifold.cc | 88 +++++++++++++++++++++++++++++------------ 1 file changed, 62 insertions(+), 26 deletions(-) diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc index b2f85d0e18..88b2bd776f 100644 --- a/source/grid/manifold.cc +++ b/source/grid/manifold.cc @@ -630,6 +630,35 @@ FlatManifold::get_new_point( +namespace +{ + template + std::array, n_eval_points> + interpolate_points(const unsigned int n_points, + const Point * surrounding_points, + const double * weights, + const Tensor<1, spacedim> &periodicity) + { + std::array, n_eval_points> new_point; + for (unsigned int p = 0; p < n_points; ++p) + { + const Point surrounding_point = surrounding_points[p]; + for (unsigned int i = 0; i < n_eval_points; ++i) + new_point[i] += surrounding_point * weights[p + i * n_points]; + } + + for (unsigned int d = 0; d < spacedim; ++d) + if (periodicity[d] > 0) + for (unsigned int i = 0; i < n_eval_points; ++i) + if (new_point[i][d] < 0) + new_point[i][d] += periodicity[d]; + + return new_point; + } +} // namespace + + + template void FlatManifold::get_new_points( @@ -671,10 +700,9 @@ FlatManifold::get_new_points( } if (adjust_periodicity == true) { - modified_points.resize(surrounding_points.size()); - std::copy(surrounding_points.begin(), - surrounding_points.end(), - modified_points.begin()); + modified_points.insert(modified_points.end(), + surrounding_points.begin(), + surrounding_points.end()); for (unsigned int d = 0; d < spacedim; ++d) if (periodicity[d] > 0) for (unsigned int i = 0; i < n_points; ++i) @@ -683,30 +711,38 @@ FlatManifold::get_new_points( surrounding_points_start = modified_points.data(); } + for (unsigned int r = 0; r < weights.size(0); ++r) + Assert(std::abs( + std::accumulate(&weights(r, 0), &weights(r, 0) + n_points, 0.0) - + 1.0) < 1e-10, + ExcMessage("The weights each point should sum to 1!")); + // Now perform the interpolation - for (unsigned int row = 0; row < weights.size(0); ++row) + for (unsigned int row = 0; row < weights.size(0);) { - Assert( - std::abs( - std::accumulate(&weights(row, 0), &weights(row, 0) + n_points, 0.0) - - 1.0) < 1e-10, - ExcMessage("The weights for the individual points should sum to 1!")); - Point new_point; - for (unsigned int p = 0; p < n_points; ++p) - new_point += surrounding_points_start[p] * weights(row, p); - - // if necessary, also adjust the weighted point by the periodicity - for (unsigned int d = 0; d < spacedim; ++d) - if (periodicity[d] > 0) - if (new_point[d] < 0) - new_point[d] += periodicity[d]; - - // TODO should this use surrounding_points_start or surrounding_points? - // The older version used surrounding_points - new_points[row] = - project_to_manifold(make_array_view(surrounding_points.begin(), - surrounding_points.end()), - new_point); + // For performance reasons, try to compute three points at once + if (row + 3 <= weights.size(0)) + { + const std::array, 3> points = interpolate_points<3>( + n_points, surrounding_points_start, &weights(row, 0), periodicity); + for (unsigned int i = 0; i < 3; ++i) + new_points[row + i] = + project_to_manifold(make_array_view(surrounding_points.begin(), + surrounding_points.end()), + points[i]); + row += 3; + } + else + { + new_points[row] = + project_to_manifold(make_array_view(surrounding_points.begin(), + surrounding_points.end()), + interpolate_points<1>(n_points, + surrounding_points_start, + &weights(row, 0), + periodicity)[0]); + ++row; + } } } -- 2.39.5