From cc4513cd2d2cdfca87082e97f9b51eebbace8022 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sun, 19 Mar 2017 10:16:53 +0100 Subject: [PATCH] Performance improvements. Augment documentation. --- include/deal.II/grid/manifold.h | 40 ++++++++++++++++------------ source/grid/manifold.cc | 46 +++++++++++++++++++-------------- 2 files changed, 50 insertions(+), 36 deletions(-) diff --git a/include/deal.II/grid/manifold.h b/include/deal.II/grid/manifold.h index 424b331bf4..fd355a490f 100644 --- a/include/deal.II/grid/manifold.h +++ b/include/deal.II/grid/manifold.h @@ -400,7 +400,7 @@ public: const std::vector &weights) const; /** - * Compute a new set of points around the given points + * Compute a new set of points that interpolate between the given points * @p surrounding_points. @p weights is a table with as many columns as * @p surrounding_points.size(). The number of rows in @p weights determines * how many new points will be computed and appended to the last input @@ -408,16 +408,18 @@ public: * @p new_points equals the size at entry plus the number of rows in * @p weights. * - * In its default implementation, this function simply calls - * get_new_point() on each row of @weights and appends those points to the - * output vector @p new_points. However, this function is more efficient if - * multiple new points need to be generated like in MappingQGeneric and the - * manifold does expensive transformations between a chart space and the - * physical space, such as ChartManifold. If efficiency is not important, - * you may get away by implementing only the get_new_point() function. + * In its default implementation, this function simply calls get_new_point() + * on each row of @weights and appends those points to the output vector + * @p new_points. However, this function is more efficient if multiple new + * points need to be generated like in MappingQGeneric and the manifold does + * expensive transformations between a chart space and the physical space, + * such as ChartManifold. For this function, the surrounding points need to + * be transformed back to the chart sparse only once, rather than for every + * call to get_new_point(). If efficiency is not important, you may get away + * by implementing only the get_new_point() function. * * The implementation does not allow for @p surrounding_points and - * @p new_points to point to the same vector, so make sure pass different + * @p new_points to point to the same vector, so make sure to pass different * objects into the function. */ virtual @@ -754,7 +756,7 @@ public: const std::vector &weights) const; /** - * Compute a new set of points around the given points + * Compute a new set of points that interpolate between the given points * @p surrounding_points. @p weights is a table with as many columns as * @p surrounding_points.size(). The number of rows in @p weights determines * how many new points will be computed and appended to the last input @@ -762,8 +764,9 @@ public: * @p new_points equals the size at entry plus the number of rows in * @p weights. * - * For this particular implementation, an interpolation of the - * @p surrounding_points according to the @p weights is performed. + * For this particular implementation, the interpolation of the + * @p surrounding_points according to the @p weights is simply performed in + * Cartesian space. */ virtual void @@ -982,7 +985,7 @@ public: const std::vector &weights) const; /** - * Compute a new set of points around the given points + * Compute a new set of points that interpolate between the given points * @p surrounding_points. @p weights is a table with as many columns as * @p surrounding_points.size(). The number of rows in @p weights determines * how many new points will be computed and appended to the last input @@ -991,16 +994,19 @@ public: * @p weights. * * The implementation of this function first transforms the - * @p surrounding_points to the chart by calling pull_back(). Then, new + * @p surrounding_points to the chart space by calling pull_back(). Then, new * points are computed on the chart by usual interpolation according to the * given @p weights, which are finally transformed to the image space by * push_forward(). * * This implementation can be much more efficient for computing multiple new * points from the same surrounding points than separate calls to - * get_new_point() in case the pull_back() operation is expensive. Often, - * this is indeed the case because pull_back() might involve Newton - * iterations or something similar in non-trivial manifolds. + * get_new_point() in case the pull_back() operation is expensive. This is + * because pull_back() is only called once for the surrounding points and + * the interpolation is done for all given weights using this set of + * points. Often, pull_back() is also more expensive than push_forward() + * because the former might involve some kind of Newton iteration in + * non-trivial manifolds. */ virtual void diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc index 93a949cd00..eb88fe9e35 100644 --- a/source/grid/manifold.cc +++ b/source/grid/manifold.cc @@ -33,16 +33,16 @@ struct CompareWeights public: CompareWeights(const std::vector &weights) : - compare_weights(&weights) + compare_weights(weights) {} bool operator() (unsigned int a, unsigned int b) const { - return (*compare_weights)[a] < (*compare_weights)[b]; + return compare_weights[a] < compare_weights[b]; } private: - const std::vector *compare_weights; + const std::vector &compare_weights; }; /* -------------------------- Manifold --------------------- */ @@ -111,12 +111,22 @@ get_new_point (const std::vector > &surrounding_points, // First sort points in the order of their weights. This is done to // produce unique points even if get_intermediate_points is not // associative (as for the SphericalManifold). - std::vector permutation(n_points); + unsigned int permutation_short[30]; + std::vector permutation_long; + unsigned int *permutation; + if (n_points > 30) + { + permutation_long.resize(n_points); + permutation = &permutation_long[0]; + } + else + permutation = &permutation_short[0]; + for (unsigned int i=0; i > &surrounding_points, else weight = w/(weights[permutation[i]] + w); - p = get_intermediate_point(p, surrounding_points[permutation[i]],1.0 - weight ); + if (std::abs(weight) > 1e-14) + p = get_intermediate_point(p, surrounding_points[permutation[i]],1.0 - weight ); w += weights[permutation[i]]; } @@ -152,11 +163,12 @@ add_new_points (const std::vector > &surrounding_points, ExcMessage("surrounding_points and new_points cannot be the same " "array")); - std::vector local_weights(surrounding_points.size()); + const unsigned int n_points = surrounding_points.size(); + std::vector local_weights(n_points); for (unsigned int row=0; row > &surrounding_points, const unsigned int n_points = surrounding_points.size(); std::vector > chart_points(n_points); - std::vector local_weights(n_points); - for (unsigned int i=0; i > new_points_on_chart; + new_points_on_chart.reserve(weights.size(0)); + sub_manifold.add_new_points(chart_points, weights, new_points_on_chart); + for (unsigned int row=0; row p_chart = sub_manifold.get_new_point(chart_points, - local_weights); - new_points.push_back(push_forward(p_chart)); - } + new_points.push_back(push_forward(new_points_on_chart[row])); } -- 2.39.5