--- /dev/null
+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)
#include <deal.II/base/function.h>
+#include <array>
+
DEAL_II_NAMESPACE_OPEN
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
#include <deal.II/base/exceptions.h>
#include <deal.II/base/function_signed_distance.h>
+#include <algorithm>
+
DEAL_II_NAMESPACE_OPEN
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
{
template class Functions::SignedDistance::Sphere<deal_II_dimension>;
template class Functions::SignedDistance::Plane<deal_II_dimension>;
+ template class Functions::SignedDistance::Ellipsoid<deal_II_dimension>;
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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();
+}
--- /dev/null
+
+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::
projects / dealii.git / commitdiff