From e21855dd0b71b9473eab2ca9b3fe39ff91829abd Mon Sep 17 00:00:00 2001 From: Maximilian Bergbauer Date: Mon, 4 Apr 2022 18:02:33 +0200 Subject: [PATCH] Add gradient to SignedDistance::Ellipsoid --- .../deal.II/base/function_signed_distance.h | 29 +++++ source/base/function_signed_distance.cc | 100 ++++++++++++++++-- tests/base/function_signed_distance_02.cc | 41 +++---- tests/base/function_signed_distance_02.output | 82 ++++++++------ 4 files changed, 190 insertions(+), 62 deletions(-) diff --git a/include/deal.II/base/function_signed_distance.h b/include/deal.II/base/function_signed_distance.h index f1727d6c05..4851396ab3 100644 --- a/include/deal.II/base/function_signed_distance.h +++ b/include/deal.II/base/function_signed_distance.h @@ -161,6 +161,14 @@ namespace Functions value(const Point & point, const unsigned int component = 0) const override; + /** + * Calculates the gradient of the signed distance function via the normal + * of the closest point on the ellipsoid. + */ + Tensor<1, dim> + gradient(const Point &, + const unsigned int component = 0) const override; + private: /** * Evaluates the ellipsoid function: @@ -172,6 +180,27 @@ namespace Functions double evaluate_ellipsoid(const Point &point) const; + /** + * Computes the closest point on the given ellipse for an arbitrary point. + */ + Point + compute_closest_point_ellipse(const Point &point) const; + + /** + * Computes the analytical normal vector of a origin centred ellipse at a + * given @p point according to: + * + * / \ / b x a y \ + * | n_x , n_y | = | ----- , ----- | + * \ / \ a b / + * + * + * @param point [x, y] + * @return [n_x , n_y] + */ + Tensor<1, dim, double> + compute_analyical_normal_vector_on_ellipse(const Point &point) const; + /** * Compute the signed distance to a 2D ellipsoid i.e. ellipse. */ diff --git a/source/base/function_signed_distance.cc b/source/base/function_signed_distance.cc index e1f28544e8..d941e3d99d 100644 --- a/source/base/function_signed_distance.cc +++ b/source/base/function_signed_distance.cc @@ -166,6 +166,35 @@ namespace Functions + template + Tensor<1, dim> + Ellipsoid::gradient(const Point & point, + const unsigned int component) const + { + AssertIndexRange(component, this->n_components); + (void)component; + + Tensor<1, dim> grad; + if (dim == 1) + grad = point - center; + else if (dim == 2) + { + const Point point_in_centered_coordinate_system = + Point(compute_closest_point_ellipse(point) - center); + grad = compute_analyical_normal_vector_on_ellipse( + point_in_centered_coordinate_system); + } + else + AssertThrow(false, ExcNotImplemented()); + + if (grad.norm() > 1e-12) + return grad / grad.norm(); + else + return grad; + } + + + template double Ellipsoid::evaluate_ellipsoid(const Point &point) const @@ -179,16 +208,11 @@ namespace Functions template - double - Ellipsoid::compute_signed_distance_ellipse( - const Point &point) const + Point + Ellipsoid::compute_closest_point_ellipse(const Point &point) const { AssertDimension(dim, 2); - // point corresponds to center - if (point.distance(center) < tolerance) - return *std::min_element(radii.begin(), radii.end()) * -1.; - /* * Function to compute the closest point on an ellipse (adopted from * https://wet-robots.ghost.io/simple-method-for-distance-to-ellipse/ and @@ -204,8 +228,10 @@ namespace Functions * 4. Repeat 2.-4. until convergence. */ // get equivalent point in first quadrant of centered ellipse - const double px = std::abs(point[0] - center[0]); - const double py = std::abs(point[1] - center[1]); + const double px = std::abs(point[0] - center[0]); + const double py = std::abs(point[1] - center[1]); + const double sign_px = std::copysign(1.0, point[0] - center[0]); + const double sign_py = std::copysign(1.0, point[1] - center[1]); // get semi axes radii const double &a = radii[0]; const double &b = radii[1]; @@ -245,7 +271,61 @@ namespace Functions while (std::abs(delta_t) > tolerance && iter < max_iter); AssertIndexRange(iter, max_iter); - const double distance = std::hypot(x - px, y - py); + AssertIsFinite(x); + AssertIsFinite(y); + + return center + Point(sign_px * x, sign_py * y); + } + + + + template + Tensor<1, dim, double> + Ellipsoid::compute_analyical_normal_vector_on_ellipse( + const Point &) const + { + AssertThrow(false, ExcNotImplemented()); + return Tensor<1, dim, double>(); + } + + + + template <> + Tensor<1, 2, double> + Ellipsoid<2>::compute_analyical_normal_vector_on_ellipse( + const Point<2> &point) const + { + const auto &a = radii[0]; + const auto &b = radii[1]; + const auto &x = point[0]; + const auto &y = point[1]; + return Tensor<1, 2, double>({b * x / a, a * y / b}); + } + + + + template + double + Ellipsoid::compute_signed_distance_ellipse(const Point &) const + { + AssertThrow(false, ExcNotImplemented()); + return 0; + } + + + + template <> + double + Ellipsoid<2>::compute_signed_distance_ellipse(const Point<2> &point) const + { + // point corresponds to center + if (point.distance(center) < tolerance) + return *std::min_element(radii.begin(), radii.end()) * -1.; + + const Point<2> &closest_point = compute_closest_point_ellipse(point); + + const double distance = + std::hypot(closest_point[0] - point[0], closest_point[1] - point[1]); return evaluate_ellipsoid(point) < 0.0 ? -distance : distance; } diff --git a/tests/base/function_signed_distance_02.cc b/tests/base/function_signed_distance_02.cc index 433a48e75e..563d0c66a3 100644 --- a/tests/base/function_signed_distance_02.cc +++ b/tests/base/function_signed_distance_02.cc @@ -27,14 +27,15 @@ namespace { template void - print_value_at_point(const Function &function, const Point &point) + print_value_and_gradient_at_point(const Function &function, + const Point & point) { deallog << "point = " << point << std::endl; deallog << "value = " << function.value(point) << std::endl; + deallog << "gradient = " << function.gradient(point) << std::endl; } - void test_ellipsoid_signed_distance_1d() { @@ -50,16 +51,16 @@ namespace const Functions::SignedDistance::Ellipsoid ellipsoid(center, radii); deallog << "center" << std::endl; - print_value_at_point(ellipsoid, center); + print_value_and_gradient_at_point(ellipsoid, center); deallog << "inside" << std::endl; - print_value_at_point(ellipsoid, Point(0.)); - print_value_at_point(ellipsoid, Point(1.5)); + print_value_and_gradient_at_point(ellipsoid, Point(0.)); + print_value_and_gradient_at_point(ellipsoid, Point(1.5)); deallog << "on ellipse" << std::endl; - print_value_at_point(ellipsoid, Point(-1.)); - print_value_at_point(ellipsoid, Point(3.)); + print_value_and_gradient_at_point(ellipsoid, Point(-1.)); + print_value_and_gradient_at_point(ellipsoid, Point(3.)); deallog << "outside" << std::endl; - print_value_at_point(ellipsoid, Point(-2.)); - print_value_at_point(ellipsoid, Point(6.)); + print_value_and_gradient_at_point(ellipsoid, Point(-2.)); + print_value_and_gradient_at_point(ellipsoid, Point(6.)); deallog << std::endl; } @@ -81,20 +82,20 @@ namespace const Functions::SignedDistance::Ellipsoid ellipsoid(center, radii); deallog << "center" << std::endl; - print_value_at_point(ellipsoid, center); + print_value_and_gradient_at_point(ellipsoid, center); deallog << "inside" << std::endl; - print_value_at_point(ellipsoid, {3., 1.}); - print_value_at_point(ellipsoid, {3.5, 2.}); - print_value_at_point(ellipsoid, {3.1, 1.9}); + print_value_and_gradient_at_point(ellipsoid, {3., 1.}); + print_value_and_gradient_at_point(ellipsoid, {3.5, 2.}); + print_value_and_gradient_at_point(ellipsoid, {3.1, 1.9}); deallog << "on ellipse" << std::endl; - print_value_at_point(ellipsoid, {2., 2.}); - print_value_at_point(ellipsoid, {3., 4.}); - print_value_at_point(ellipsoid, {2.5, 2. + std::sqrt(3.)}); + print_value_and_gradient_at_point(ellipsoid, {2., 2.}); + print_value_and_gradient_at_point(ellipsoid, {3., 4.}); + print_value_and_gradient_at_point(ellipsoid, {2.5, 2. + std::sqrt(3.)}); deallog << "outside" << std::endl; - print_value_at_point(ellipsoid, {3., -1.}); - print_value_at_point(ellipsoid, {-1., 2.}); - print_value_at_point(ellipsoid, {0., 0.}); - print_value_at_point(ellipsoid, {4., 5.}); + print_value_and_gradient_at_point(ellipsoid, {3., -1.}); + print_value_and_gradient_at_point(ellipsoid, {-1., 2.}); + print_value_and_gradient_at_point(ellipsoid, {0., 0.}); + print_value_and_gradient_at_point(ellipsoid, {4., 5.}); deallog << std::endl; } diff --git a/tests/base/function_signed_distance_02.output b/tests/base/function_signed_distance_02.output index 6b56eda45b..ac98384e97 100644 --- a/tests/base/function_signed_distance_02.output +++ b/tests/base/function_signed_distance_02.output @@ -3,51 +3,69 @@ DEAL::test_ellipsoid_signed_distance DEAL::dim = 1 DEAL:: DEAL::center -DEAL::point = 1.0 -DEAL::value = -2.0 +DEAL::point = 1.00000 +DEAL::value = -2.00000 +DEAL::gradient = 0.00000 DEAL::inside -DEAL::point = 0.0 -DEAL::value = -1.0 -DEAL::point = 1.5 -DEAL::value = -1.5 +DEAL::point = 0.00000 +DEAL::value = -1.00000 +DEAL::gradient = -1.00000 +DEAL::point = 1.50000 +DEAL::value = -1.50000 +DEAL::gradient = 1.00000 DEAL::on ellipse -DEAL::point = -1.0 -DEAL::value = 0.0 -DEAL::point = 3.0 -DEAL::value = 0.0 +DEAL::point = -1.00000 +DEAL::value = 0.00000 +DEAL::gradient = -1.00000 +DEAL::point = 3.00000 +DEAL::value = 0.00000 +DEAL::gradient = 1.00000 DEAL::outside -DEAL::point = -2.0 -DEAL::value = 1.0 -DEAL::point = 6.0 -DEAL::value = 3.0 +DEAL::point = -2.00000 +DEAL::value = 1.00000 +DEAL::gradient = -1.00000 +DEAL::point = 6.00000 +DEAL::value = 3.00000 +DEAL::gradient = 1.00000 DEAL:: DEAL::test_ellipsoid_signed_distance DEAL::dim = 2 DEAL:: DEAL::center -DEAL::point = 3.0 2.0 -DEAL::value = -1.0 +DEAL::point = 3.00000 2.00000 +DEAL::value = -1.00000 +DEAL::gradient = 1.00000 0.00000 DEAL::inside -DEAL::point = 3.0 1.0 +DEAL::point = 3.00000 1.00000 DEAL::value = -0.816497 -DEAL::point = 3.5 2.0 -DEAL::value = -0.50 -DEAL::point = 3.1 1.9 +DEAL::gradient = 0.912871 -0.408248 +DEAL::point = 3.50000 2.00000 +DEAL::value = -0.500000 +DEAL::gradient = 1.00000 0.00000 +DEAL::point = 3.10000 1.90000 DEAL::value = -0.898386 +DEAL::gradient = 0.999478 -0.0323064 DEAL::on ellipse -DEAL::point = 2.0 2.0 -DEAL::value = 0.0 -DEAL::point = 3.0 4.0 -DEAL::value = 0.0 -DEAL::point = 2.5 3.73205 -DEAL::value = 0.0 +DEAL::point = 2.00000 2.00000 +DEAL::value = 0.00000 +DEAL::gradient = -1.00000 0.00000 +DEAL::point = 3.00000 4.00000 +DEAL::value = 0.00000 +DEAL::gradient = 0.00000 1.00000 +DEAL::point = 2.50000 3.73205 +DEAL::value = 0.00000 +DEAL::gradient = -0.755929 0.654654 DEAL::outside -DEAL::point = 3.0 -1.0 -DEAL::value = 1.0 -DEAL::point = -1.0 2.0 -DEAL::value = 3.0 -DEAL::point = 0.0 0.0 +DEAL::point = 3.00000 -1.00000 +DEAL::value = 1.00000 +DEAL::gradient = 0.00000 -1.00000 +DEAL::point = -1.00000 2.00000 +DEAL::value = 3.00000 +DEAL::gradient = -1.00000 0.00000 +DEAL::point = 0.00000 0.00000 DEAL::value = 2.34088 -DEAL::point = 4.0 5.0 +DEAL::gradient = -0.938013 -0.346599 +DEAL::point = 4.00000 5.00000 DEAL::value = 1.29718 +DEAL::gradient = 0.537059 0.843545 DEAL:: -- 2.39.5