From c9a92e6312d45e3a642bd614f3d58e9b728ca77e Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Mon, 4 Dec 2017 17:52:50 -0700 Subject: [PATCH] Avoid more memory allocation --- include/deal.II/base/array_view.h | 48 ++++++++++++++++ include/deal.II/grid/manifold_lib.h | 19 +++++++ source/grid/manifold_lib.cc | 88 ++++++++++++++++++----------- 3 files changed, 121 insertions(+), 34 deletions(-) diff --git a/include/deal.II/base/array_view.h b/include/deal.II/base/array_view.h index b3f9f6c02f..013e86a034 100644 --- a/include/deal.II/base/array_view.h +++ b/include/deal.II/base/array_view.h @@ -614,6 +614,54 @@ make_array_view (Table<2,ElementType> &table, +/** + * Create a view to an entire Table<2> object. This is equivalent to + * initializing an ArrayView object with a pointer to the first element of the + * given table, and the number of table entries as the length of the view. + * + * This function is used for non-@p const references to objects of Table type. + * Such objects contain elements that can be written to. Consequently, the + * return type of this function is a view to a set of writable objects. + * + * @param[in] table The Table for which we want to have an array view object. + * The array view corresponds to the entire table. + * + * @relates ArrayView + */ +template +inline +ArrayView +make_array_view (Table<2,ElementType> &table) +{ + return ArrayView (&table[0][0], table.n_elements()); +} + + + +/** + * Create a view to an entire Table<2> object. This is equivalent to + * initializing an ArrayView object with a pointer to the first element of the + * given table, and the number of table entries as the length of the view. + * + * This function is used for @p const references to objects of Table type + * because they contain immutable elements. Consequently, the return type of + * this function is a view to a set of @p const objects. + * + * @param[in] table The Table for which we want to have an array view object. + * The array view corresponds to the entire table. + * + * @relates ArrayView + */ +template +inline +ArrayView +make_array_view (const Table<2,ElementType> &table) +{ + return ArrayView (&table[0][0], table.n_elements()); +} + + + /** * Create a view to an entire row of a Table<2> object. This is equivalent to * initializing an ArrayView object with a pointer to the first element of the diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h index 2d25d7abfb..8958348370 100644 --- a/include/deal.II/grid/manifold_lib.h +++ b/include/deal.II/grid/manifold_lib.h @@ -301,6 +301,25 @@ private: const ArrayView &weights, const Point &candidate_point) const; + /** + * Compute a new set of points that interpolate between the given points @p + * surrounding_points. @p weights is an array view with as many entries as @p + * surrounding_points.size() times @p new_points.size(). + * + * This function is optimized to perform on a collection + * of new points, by collecting operations that are not dependent on the + * weights outside of the loop over all new points. + * + * The implementation does not allow for @p surrounding_points and + * @p new_points to point to the same array, so make sure to pass different + * objects into the function. + */ + virtual + void + get_new_points (const ArrayView> &surrounding_points, + const ArrayView &weights, + ArrayView> new_points) const; + /** * A manifold description to be used for get_new_point in 2D. **/ diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index c061714c06..3bd2c3c25a 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -376,15 +376,53 @@ void SphericalManifold:: get_new_points (const ArrayView> &surrounding_points, const Table<2,double> &weights, - ArrayView> new_points) const + ArrayView> new_points) const { AssertDimension(new_points.size(), weights.size(0)); AssertDimension(surrounding_points.size(), weights.size(1)); + get_new_points(surrounding_points, + make_array_view(weights), + new_points); + + return; +} + + + +template +Point +SphericalManifold:: +get_new_point (const ArrayView> &vertices, + const ArrayView &weights) const +{ + // To avoid duplicating all of the logic in get_new_points, simply call it + // for one position. + Point new_point; + get_new_points(vertices, + weights, + make_array_view(&new_point,&new_point+1)); + + return new_point; +} + + + +template +void +SphericalManifold:: +get_new_points (const ArrayView> &surrounding_points, + const ArrayView &weights, + ArrayView> new_points) const +{ + AssertDimension(weights.size(), new_points.size() * surrounding_points.size()); + const unsigned int weight_rows = new_points.size(); + const unsigned int weight_columns = surrounding_points.size(); + if (surrounding_points.size() == 2) { - for (unsigned int row=0; row> &surrounding_points, boost::container::small_vector accurate_point_was_found(new_points.size(),false); const ArrayView> array_directions = make_array_view(directions.begin(),directions.end()); const ArrayView array_distances = make_array_view(distances.begin(),distances.end()); - for (unsigned int row=0; row(&weights[row*weight_columns],weight_columns)); // If the candidate is the center, mark it as found to avoid entering // the Newton iteration in step 2, which would crash. @@ -437,7 +475,7 @@ get_new_points (const ArrayView> &surrounding_points, // after we verified that the candidate is not the center. if (spacedim<3) { - new_points[row] = polar_manifold.get_new_point(surrounding_points, make_array_view(weights,row)); + new_points[row] = polar_manifold.get_new_point(surrounding_points, ArrayView(&weights[row*weight_columns],weight_columns)); accurate_point_was_found[row] = true; continue; } @@ -450,7 +488,7 @@ get_new_points (const ArrayView> &surrounding_points, // doesn't already succeed. if (max_distance < 2e-2 || spacedim<3) { - for (unsigned int row=0; row> &surrounding_points, // Search for duplicate directions and merge them to minimize the cost of // the get_new_point function call below. - boost::container::small_vector merged_weights(weights.size(0)*weights.size(1)); + boost::container::small_vector merged_weights(weights.size()); boost::container::small_vector, 100> merged_directions(surrounding_points.size(),Point()); boost::container::small_vector merged_distances(surrounding_points.size(),0.0); @@ -478,8 +516,8 @@ get_new_points (const ArrayView> &surrounding_points, if (!found_duplicate && squared_distance < 1e-28) { found_duplicate = true; - for (unsigned int row = 0; row < weights.size(0); ++row) - merged_weights[row*weights.size(1) + j] += weights[row][i]; + for (unsigned int row = 0; row < weight_rows; ++row) + merged_weights[row*weight_columns + j] += weights[row*weight_columns + i]; } } @@ -487,8 +525,8 @@ get_new_points (const ArrayView> &surrounding_points, { merged_directions[n_unique_directions] = directions[i]; merged_distances[n_unique_directions] = distances[i]; - for (unsigned int row = 0; row < weights.size(0); ++row) - merged_weights[row*weights.size(1) + n_unique_directions] = weights[row][i]; + for (unsigned int row = 0; row < weight_rows; ++row) + merged_weights[row*weight_columns + n_unique_directions] = weights[row*weight_columns + i]; ++n_unique_directions; } @@ -498,7 +536,7 @@ get_new_points (const ArrayView> &surrounding_points, // the get_new_point function call below. boost::container::small_vector merged_weights_index(new_points.size(),numbers::invalid_unsigned_int); unsigned int unique_weight_rows = 0; - for (unsigned int row = 0; row < new_points.size(); ++row) + for (unsigned int row = 0; row < weight_rows; ++row) { bool found_identical_row = false; for (unsigned int existing_row = 0; existing_row < row; ++existing_row) @@ -506,7 +544,7 @@ get_new_points (const ArrayView> &surrounding_points, bool identical_weights = true; for (unsigned int weight_index = 0; weight_index < n_unique_directions; ++weight_index) - if (merged_weights[row*weights.size(1) + weight_index] != merged_weights[existing_row*weights.size(1) + weight_index]) + if (std::abs(merged_weights[row*weight_columns + weight_index] - merged_weights[existing_row*weight_columns + weight_index]) > tolerance) { identical_weights = false; break; @@ -527,12 +565,12 @@ get_new_points (const ArrayView> &surrounding_points, const ArrayView array_merged_distances = make_array_view(merged_distances.begin(), merged_distances.begin()+n_unique_directions); - for (unsigned int row=0; row array_merged_weights (&merged_weights[row*weights.size(1)],n_unique_directions); + const ArrayView array_merged_weights (&merged_weights[row*weight_columns],n_unique_directions); new_candidates[row].second = get_new_point(array_merged_directions, array_merged_distances, array_merged_weights, @@ -547,24 +585,6 @@ get_new_points (const ArrayView> &surrounding_points, -template -Point -SphericalManifold:: -get_new_point (const ArrayView> &vertices, - const ArrayView &weights) const -{ - // To avoid duplicating all of the logic in get_new_points, simply call it - // for one position. - const Table<2,double> table_weights(1,weights.size(),weights.begin()); - Point new_point; - ArrayView> new_points(&new_point,1); - - get_new_points(vertices,table_weights,new_points); - return new_points[0]; -} - - - template std::pair > SphericalManifold:: -- 2.39.5