]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add signed distance function of 1 and 2D ellipsoid 12920/head
authorNils Much <nils.much@tum.de>
Thu, 4 Nov 2021 15:10:09 +0000 (16:10 +0100)
committerNils Much <nils.much@tum.de>
Sat, 6 Nov 2021 12:43:23 +0000 (13:43 +0100)
doc/news/changes/minor/20211105Much [new file with mode: 0644]
include/deal.II/base/function_signed_distance.h
source/base/function_signed_distance.cc
source/base/function_signed_distance.inst.in
tests/base/function_signed_distance.cc [new file with mode: 0644]
tests/base/function_signed_distance.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20211105Much b/doc/news/changes/minor/20211105Much
new file mode 100644 (file)
index 0000000..0af7894
--- /dev/null
@@ -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).
+<br>
+(Nils Much, 2021/11/05)
index 253aaecb3f85eee7ab51519aba9f450c6150b9ba..f1727d6c05ab41f67a057c3c4e9c4b37cde1fef5 100644 (file)
@@ -20,6 +20,8 @@
 
 #include <deal.II/base/function.h>
 
+#include <array>
+
 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 <int dim>
+    class Ellipsoid : public Function<dim>
+    {
+    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<dim> &             center,
+                const std::array<double, dim> &radii,
+                const double                   tolerance = 1e-14,
+                const unsigned int             max_iter  = 10);
+
+      double
+      value(const Point<dim> & 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<dim> &point) const;
+
+      /**
+       * Compute the signed distance to a 2D ellipsoid i.e. ellipse.
+       */
+      double
+      compute_signed_distance_ellipse(const Point<dim> &point) const;
+
+      const Point<dim>              center;
+      const std::array<double, dim> radii;
+      const double                  tolerance;
+      const unsigned int            max_iter;
+    };
   } // namespace SignedDistance
 } // namespace Functions
 
index 114b358e1725d392b41a059a89a9f405839bc902..af14810e22f0cd3e79df1262de7a3498bf67a237 100644 (file)
@@ -16,6 +16,8 @@
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/function_signed_distance.h>
 
+#include <algorithm>
+
 DEAL_II_NAMESPACE_OPEN
 
 namespace Functions
@@ -126,6 +128,127 @@ namespace Functions
       return SymmetricTensor<2, dim>();
     }
 
+
+
+    template <int dim>
+    Ellipsoid<dim>::Ellipsoid(const Point<dim> &             center,
+                              const std::array<double, dim> &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 <int dim>
+    double
+    Ellipsoid<dim>::value(const Point<dim> & 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 <int dim>
+    double
+    Ellipsoid<dim>::evaluate_ellipsoid(const Point<dim> &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 <int dim>
+    double
+    Ellipsoid<dim>::compute_signed_distance_ellipse(
+      const Point<dim> &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
 
index 9227b430013ff9bb34f3ad5f08bd8a68a95a3d3d..8c71662932959627306e6e20629fe340655c21ad 100644 (file)
@@ -17,4 +17,5 @@ for (deal_II_dimension : SPACE_DIMENSIONS)
   {
     template class Functions::SignedDistance::Sphere<deal_II_dimension>;
     template class Functions::SignedDistance::Plane<deal_II_dimension>;
+    template class Functions::SignedDistance::Ellipsoid<deal_II_dimension>;
   }
diff --git a/tests/base/function_signed_distance.cc b/tests/base/function_signed_distance.cc
new file mode 100644 (file)
index 0000000..433a48e
--- /dev/null
@@ -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 <deal.II/base/function_signed_distance.h>
+
+#include "../tests.h"
+
+namespace
+{
+  template <int dim>
+  void
+  print_value_at_point(const Function<dim> &function, const Point<dim> &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<dim>              center(1.);
+    const std::array<double, dim> radii{{2.}};
+
+    const Functions::SignedDistance::Ellipsoid<dim> ellipsoid(center, radii);
+
+    deallog << "center" << std::endl;
+    print_value_at_point(ellipsoid, center);
+    deallog << "inside" << std::endl;
+    print_value_at_point(ellipsoid, Point<dim>(0.));
+    print_value_at_point(ellipsoid, Point<dim>(1.5));
+    deallog << "on ellipse" << std::endl;
+    print_value_at_point(ellipsoid, Point<dim>(-1.));
+    print_value_at_point(ellipsoid, Point<dim>(3.));
+    deallog << "outside" << std::endl;
+    print_value_at_point(ellipsoid, Point<dim>(-2.));
+    print_value_at_point(ellipsoid, Point<dim>(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<dim>              center(3., 2.);
+    const std::array<double, dim> radii{{1., 2.}};
+
+    const Functions::SignedDistance::Ellipsoid<dim> 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 (file)
index 0000000..6b56eda
--- /dev/null
@@ -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::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.