From eb6a00ded2d9a6589ffe0e141f73072d04b680cc Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Sun, 3 Dec 2017 15:46:30 -0700 Subject: [PATCH] Merge horizontally --- include/deal.II/grid/manifold_lib.h | 2 +- source/grid/manifold_lib.cc | 102 ++++++++++++++++++---------- 2 files changed, 69 insertions(+), 35 deletions(-) diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h index 980323d031..2d25d7abfb 100644 --- a/include/deal.II/grid/manifold_lib.h +++ b/include/deal.II/grid/manifold_lib.h @@ -299,7 +299,7 @@ private: get_new_point (const ArrayView> &directions, const ArrayView &distances, const ArrayView &weights, - const Point candidate_point) const; + const Point &candidate_point) 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 b66bb88043..c061714c06 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -378,6 +378,7 @@ get_new_points (const ArrayView> &surrounding_points, const Table<2,double> &weights, ArrayView> new_points) const { + AssertDimension(new_points.size(), weights.size(0)); AssertDimension(surrounding_points.size(), weights.size(1)); if (surrounding_points.size() == 2) @@ -387,6 +388,7 @@ get_new_points (const ArrayView> &surrounding_points, return; } + boost::container::small_vector>, 100> new_candidates(new_points.size()); boost::container::small_vector, 100> directions(surrounding_points.size(),Point()); boost::container::small_vector distances(surrounding_points.size(),0.0); double max_distance = 0.; @@ -418,13 +420,13 @@ get_new_points (const ArrayView> &surrounding_points, const ArrayView array_distances = make_array_view(distances.begin(),distances.end()); for (unsigned int row=0; row> candidate = guess_new_point(array_directions, - array_distances, - make_array_view(weights,row)); + new_candidates[row] = guess_new_point(array_directions, + array_distances, + make_array_view(weights,row)); // If the candidate is the center, mark it as found to avoid entering // the Newton iteration in step 2, which would crash. - if (candidate.first == 0.0) + if (new_candidates[row].first == 0.0) { new_points[row] = center; accurate_point_was_found[row] = true; @@ -440,15 +442,19 @@ get_new_points (const ArrayView> &surrounding_points, continue; } - new_points[row] = center + candidate.first * candidate.second; } // If all the points are close to each other, we expect the estimate to // be good enough. This tolerance was chosen such that the first iteration // for a at least three time refined HyperShell mesh with radii .5 and 1. // doesn't already succeed. - if (max_distance < 1e-2 || spacedim<3) - return; + if (max_distance < 2e-2 || spacedim<3) + { + for (unsigned int row=0; row> &surrounding_points, } } + // Search for duplicate weight rows and merge them to minimize the cost of + // 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) + { + bool found_identical_row = false; + for (unsigned int existing_row = 0; existing_row < row; ++existing_row) + { + 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]) + { + identical_weights = false; + break; + } + + if (identical_weights) + { + found_identical_row = true; + merged_weights_index[row] = existing_row; + break; + } + } + } + // Note that we only use the n_unique_directions first entries in the ArrayView const ArrayView> array_merged_directions = make_array_view(merged_directions.begin(), merged_directions.begin()+n_unique_directions); @@ -497,11 +530,18 @@ get_new_points (const ArrayView> &surrounding_points, for (unsigned int row=0; row array_merged_weights (&merged_weights[row*weights.size(1)],n_unique_directions); - new_points[row] = get_new_point(array_merged_directions, - array_merged_distances, - array_merged_weights, - new_points[row]); + if (merged_weights_index[row] == numbers::invalid_unsigned_int) + { + const ArrayView array_merged_weights (&merged_weights[row*weights.size(1)],n_unique_directions); + new_candidates[row].second = get_new_point(array_merged_directions, + array_merged_distances, + array_merged_weights, + Point(new_candidates[row].second)); + } + else + new_candidates[row].second = new_candidates[merged_weights_index[row]].second; + + new_points[row] = center + new_candidates[row].first * new_candidates[row].second; } } @@ -567,10 +607,11 @@ namespace do_get_new_point(const ArrayView> &directions, const ArrayView &distances, const ArrayView &weights, - const Point candidate_point, - const Point center) + const Point &candidate_point, + const Point ¢er) { Assert(false,ExcNotImplemented()); + return Point(); } template <> @@ -578,21 +619,17 @@ namespace do_get_new_point(const ArrayView> &directions, const ArrayView &distances, const ArrayView &weights, - const Point<3> candidate_point, - const Point<3> center) + const Point<3> &candidate_point, + const Point<3> ¢er) { AssertDimension(directions.size(), distances.size()); AssertDimension(directions.size(), weights.size()); + Point<3> candidate = candidate_point; const unsigned int n_merged_points = directions.size(); const double tolerance = 1e-10; const int max_iterations = 10; - // Recover radius and normalized direction from candidate point - Tensor <1,3> candidate = candidate_point - center; - const double rho = candidate.norm(); - candidate /= rho; - { // If the candidate happens to coincide with a normalized // direction, we return it. Otherwise, the Hessian would be singular. @@ -600,7 +637,7 @@ namespace { const double squared_distance = (candidate - directions[i]).norm_square(); if (squared_distance < tolerance*tolerance) - return center + rho * candidate; + return candidate; } // check if we only have two points now, in which case we can use the @@ -612,7 +649,7 @@ namespace ExcMessage("Weights do not sum up to 1")); Point<3> intermediate = unit_manifold.get_intermediate_point (Point<3>(directions[0]), Point<3>(directions[1]), weights[1]); - return center + rho * intermediate; + return intermediate; } Tensor<1,3> vPerp; @@ -679,18 +716,15 @@ namespace const Tensor<1,3> xDisp = xDisplocal[0]*Clocalx + xDisplocal[1]*Clocaly; // Step 2b: rotate candidate in direction xDisp for a new candidate. - const Tensor<1,3> candidateOld = candidate; - candidate = internal::apply_exponential_map(candidate, xDisp); + const Point<3> candidateOld = candidate; + candidate = Point<3>(internal::apply_exponential_map(candidate, xDisp)); - // Step 3c: return the new candidate if we didn't move + // Step 2c: return the new candidate if we didn't move if ((candidate-candidateOld).norm_square() < tolerance*tolerance) break; } - - Assert (std::abs(candidate[2]) < tolerance, - ExcInternalError()); } - return center + rho*candidate; + return candidate; } } @@ -702,7 +736,7 @@ SphericalManifold:: get_new_point (const ArrayView> &directions, const ArrayView &distances, const ArrayView &weights, - const Point candidate_point) const + const Point &candidate_point) const { Assert (false, ExcNotImplemented()); return Point(); @@ -716,7 +750,7 @@ SphericalManifold<1,3>:: get_new_point (const ArrayView> &directions, const ArrayView &distances, const ArrayView &weights, - const Point<3> candidate_point) const + const Point<3> &candidate_point) const { return do_get_new_point(directions,distances,weights,candidate_point,center); } @@ -729,7 +763,7 @@ SphericalManifold<2,3>:: get_new_point (const ArrayView> &directions, const ArrayView &distances, const ArrayView &weights, - const Point<3> candidate_point) const + const Point<3> &candidate_point) const { return do_get_new_point(directions,distances,weights,candidate_point,center); } @@ -742,7 +776,7 @@ SphericalManifold<3,3>:: get_new_point (const ArrayView> &directions, const ArrayView &distances, const ArrayView &weights, - const Point<3> candidate_point) const + const Point<3> &candidate_point) const { return do_get_new_point(directions,distances,weights,candidate_point,center); } -- 2.39.5