From eab5da59a0bd7bea90b1213d827e0d1ce0c39c9f Mon Sep 17 00:00:00 2001 From: Nils Much Date: Thu, 4 Nov 2021 16:10:09 +0100 Subject: [PATCH] Add signed distance function of 1 and 2D ellipsoid --- doc/news/changes/minor/20211105Much | 5 + .../deal.II/base/function_signed_distance.h | 59 +++++++++ source/base/function_signed_distance.cc | 123 ++++++++++++++++++ source/base/function_signed_distance.inst.in | 1 + tests/base/function_signed_distance.cc | 111 ++++++++++++++++ tests/base/function_signed_distance.output | 53 ++++++++ 6 files changed, 352 insertions(+) create mode 100644 doc/news/changes/minor/20211105Much create mode 100644 tests/base/function_signed_distance.cc create mode 100644 tests/base/function_signed_distance.output diff --git a/doc/news/changes/minor/20211105Much b/doc/news/changes/minor/20211105Much new file mode 100644 index 0000000000..0af78948d1 --- /dev/null +++ b/doc/news/changes/minor/20211105Much @@ -0,0 +1,5 @@ +New: The new class Functions::SignedDistance::Ellipsoid computes +a signed distance function for ellipsoids. It is implemented for 1D and +2D ellipsoids (ellipses). +
+(Nils Much, 2021/11/05) diff --git a/include/deal.II/base/function_signed_distance.h b/include/deal.II/base/function_signed_distance.h index 253aaecb3f..f1727d6c05 100644 --- a/include/deal.II/base/function_signed_distance.h +++ b/include/deal.II/base/function_signed_distance.h @@ -20,6 +20,8 @@ #include +#include + DEAL_II_NAMESPACE_OPEN namespace Functions @@ -124,6 +126,63 @@ namespace Functions const Tensor<1, dim> normal; }; + + /** + * Signed-distance level set function to an ellipsoid defined by: + * + * @f[ + * \sum_{i=1}^{dim} \frac{(x_i - c_i)^2}{R_i^2} = 1 + * @f] + * + * Here, $c_i$ are the coordinates of the center of the ellipsoid and $R_i$ + * are the elliptic radii. This function is zero on the ellipsoid, negative + * inside the ellipsoid and positive outside the ellipsoid. + * + * @ingroup functions + */ + template + class Ellipsoid : public Function + { + public: + /** + * Constructor, takes the center and radii of the ellipsoid. + * + * @param center Center of the ellipsoid. + * @param radii Array of radii of the ellipsoid. + * @param tolerance Tolerance of the distance computation. + * @param max_iter Max. number of iteration of the distance computation algorithm. + */ + Ellipsoid(const Point & center, + const std::array &radii, + const double tolerance = 1e-14, + const unsigned int max_iter = 10); + + double + value(const Point & point, + const unsigned int component = 0) const override; + + private: + /** + * Evaluates the ellipsoid function: + * + * @f[ + * f(\vec{x}) = \sum_{i=1}^{dim} \frac{(x_i - c_i)^2}{R_i^2} - 1 + * @f] + */ + double + evaluate_ellipsoid(const Point &point) const; + + /** + * Compute the signed distance to a 2D ellipsoid i.e. ellipse. + */ + double + compute_signed_distance_ellipse(const Point &point) const; + + const Point center; + const std::array radii; + const double tolerance; + const unsigned int max_iter; + }; } // namespace SignedDistance } // namespace Functions diff --git a/source/base/function_signed_distance.cc b/source/base/function_signed_distance.cc index 114b358e17..af14810e22 100644 --- a/source/base/function_signed_distance.cc +++ b/source/base/function_signed_distance.cc @@ -16,6 +16,8 @@ #include #include +#include + DEAL_II_NAMESPACE_OPEN namespace Functions @@ -126,6 +128,127 @@ namespace Functions return SymmetricTensor<2, dim>(); } + + + template + Ellipsoid::Ellipsoid(const Point & center, + const std::array &radii, + const double tolerance, + const unsigned int max_iter) + : center(center) + , radii(radii) + , tolerance(tolerance) + , max_iter(max_iter) + { + for (int d = 0; d < dim; ++d) + Assert(radii[d] > 0, ExcMessage("All radii must be positive.")) + } + + + + template + double + Ellipsoid::value(const Point & point, + const unsigned int component) const + { + AssertIndexRange(component, this->n_components); + (void)component; + + if (dim == 1) + return point.distance(center) - radii[0]; + else if (dim == 2) + return compute_signed_distance_ellipse(point); + else + Assert(false, ExcNotImplemented()); + + return 0.0; + } + + + + template + double + Ellipsoid::evaluate_ellipsoid(const Point &point) const + { + double val = 0.0; + for (int d = 0; d < dim; ++d) + val += std::pow((point[d] - center[d]) / radii[d], 2); + return val - 1.0; + } + + + + template + double + Ellipsoid::compute_signed_distance_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 + * https://github.com/0xfaded/ellipse_demo): + * + * Since the ellipse is symmetric to the two major axes through its + * center, the point is moved so the center coincides with the origin and + * into the first quadrant. + * 1. Choose a point on the ellipse (x), here x = a*cos(pi/4) and y = + * b*sin(pi/4). + * 2. Find second point on the ellipse, that has the same distance. + * 3. Find midpoint on the ellipse (must be closer). + * 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]); + // get semi axes radii + const double &a = radii[0]; + const double &b = radii[1]; + // initial guess (t = angle from x-axis) + double t = numbers::PI_4; + double x = a * std::cos(t); + double y = b * std::sin(t); + + unsigned int iter = 0; + double delta_t; + do + { + // compute the ellipse evolute (center of curvature) for the current t + const double ex = (a * a - b * b) * std::pow(std::cos(t), 3) / a; + const double ey = (b * b - a * a) * std::pow(std::sin(t), 3) / b; + // compute distances from current point on ellipse to its evolute + const double rx = x - ex; + const double ry = y - ey; + // compute distances from point to the current evolute + const double qx = px - ex; + const double qy = py - ey; + // compute the curvature radius at the current point on the ellipse + const double r = std::hypot(rx, ry); + // compute the distance from evolute to the point + const double q = std::hypot(qx, qy); + // compute step size on ellipse + const double delta_c = r * std::asin((rx * qy - ry * qx) / (r * q)); + // compute approximate angle step + delta_t = delta_c / std::sqrt(a * a + b * b - x * x - y * y); + t += delta_t; + // make sure the angle stays in first quadrant + t = std::min(numbers::PI_2, std::max(0.0, t)); + x = a * std::cos(t); + y = b * std::sin(t); + iter++; + } + while (std::abs(delta_t) > tolerance && iter < max_iter); + AssertIndexRange(iter, max_iter); + + const double distance = std::hypot(x - px, y - py); + + return evaluate_ellipsoid(point) < 0.0 ? -distance : distance; + } } // namespace SignedDistance } // namespace Functions diff --git a/source/base/function_signed_distance.inst.in b/source/base/function_signed_distance.inst.in index 9227b43001..8c71662932 100644 --- a/source/base/function_signed_distance.inst.in +++ b/source/base/function_signed_distance.inst.in @@ -17,4 +17,5 @@ for (deal_II_dimension : SPACE_DIMENSIONS) { template class Functions::SignedDistance::Sphere; template class Functions::SignedDistance::Plane; + template class Functions::SignedDistance::Ellipsoid; } diff --git a/tests/base/function_signed_distance.cc b/tests/base/function_signed_distance.cc new file mode 100644 index 0000000000..433a48e75e --- /dev/null +++ b/tests/base/function_signed_distance.cc @@ -0,0 +1,111 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2007 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +/** + * Test the class Functions::SignedDistance::Ellipsoid by creating an object of + * the class and calling the value function at points where we know what the + * function values should be. + */ + +#include + +#include "../tests.h" + +namespace +{ + template + void + print_value_at_point(const Function &function, const Point &point) + { + deallog << "point = " << point << std::endl; + deallog << "value = " << function.value(point) << std::endl; + } + + + + void + test_ellipsoid_signed_distance_1d() + { + static constexpr int dim = 1; + + deallog << "test_ellipsoid_signed_distance" << std::endl; + deallog << "dim = " << dim << std::endl; + deallog << std::endl; + + const Point center(1.); + const std::array radii{{2.}}; + + const Functions::SignedDistance::Ellipsoid ellipsoid(center, radii); + + deallog << "center" << std::endl; + print_value_at_point(ellipsoid, center); + deallog << "inside" << std::endl; + print_value_at_point(ellipsoid, Point(0.)); + print_value_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.)); + deallog << "outside" << std::endl; + print_value_at_point(ellipsoid, Point(-2.)); + print_value_at_point(ellipsoid, Point(6.)); + + deallog << std::endl; + } + + + + void + test_ellipsoid_signed_distance_2d() + { + static constexpr int dim = 2; + + deallog << "test_ellipsoid_signed_distance" << std::endl; + deallog << "dim = " << dim << std::endl; + deallog << std::endl; + + const Point center(3., 2.); + const std::array radii{{1., 2.}}; + + const Functions::SignedDistance::Ellipsoid ellipsoid(center, radii); + + deallog << "center" << std::endl; + print_value_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}); + 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.)}); + 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.}); + + deallog << std::endl; + } +} // namespace + + + +int +main() +{ + initlog(); + test_ellipsoid_signed_distance_1d(); + test_ellipsoid_signed_distance_2d(); +} diff --git a/tests/base/function_signed_distance.output b/tests/base/function_signed_distance.output new file mode 100644 index 0000000000..6b56eda45b --- /dev/null +++ b/tests/base/function_signed_distance.output @@ -0,0 +1,53 @@ + +DEAL::test_ellipsoid_signed_distance +DEAL::dim = 1 +DEAL:: +DEAL::center +DEAL::point = 1.0 +DEAL::value = -2.0 +DEAL::inside +DEAL::point = 0.0 +DEAL::value = -1.0 +DEAL::point = 1.5 +DEAL::value = -1.5 +DEAL::on ellipse +DEAL::point = -1.0 +DEAL::value = 0.0 +DEAL::point = 3.0 +DEAL::value = 0.0 +DEAL::outside +DEAL::point = -2.0 +DEAL::value = 1.0 +DEAL::point = 6.0 +DEAL::value = 3.0 +DEAL:: +DEAL::test_ellipsoid_signed_distance +DEAL::dim = 2 +DEAL:: +DEAL::center +DEAL::point = 3.0 2.0 +DEAL::value = -1.0 +DEAL::inside +DEAL::point = 3.0 1.0 +DEAL::value = -0.816497 +DEAL::point = 3.5 2.0 +DEAL::value = -0.50 +DEAL::point = 3.1 1.9 +DEAL::value = -0.898386 +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::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::value = 2.34088 +DEAL::point = 4.0 5.0 +DEAL::value = 1.29718 +DEAL:: -- 2.39.5