Add gradient to SignedDistance::Ellipsoid

diff --git a/include/deal.II/base/function_signed_distance.h b/include/deal.II/base/function_signed_distance.h
index f1727d6c05ab41f67a057c3c4e9c4b37cde1fef5..4851396ab314a94cf383715bcb33bba80f0f0e3b 100644 (file)
--- 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<dim> & 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<dim> &,
+               const unsigned int component = 0) const override;
+
     private:
       /**
        * Evaluates the ellipsoid function:
@@ -172,6 +180,27 @@ namespace Functions
       double
       evaluate_ellipsoid(const Point<dim> &point) const;
 
+      /**
+       * Computes the closest point on the given ellipse for an arbitrary point.
+       */
+      Point<dim>
+      compute_closest_point_ellipse(const Point<dim> &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<dim> &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 e1f28544e88fff595f1eba727b091cc1aa63c2ac..d941e3d99dabafca82f08751164f218596a2a925 100644 (file)
--- a/source/base/function_signed_distance.cc
+++ b/source/base/function_signed_distance.cc
@@ -166,6 +166,35 @@ namespace Functions
 
 
 
+    template <int dim>
+    Tensor<1, dim>
+    Ellipsoid<dim>::gradient(const Point<dim> & 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<dim> point_in_centered_coordinate_system =
+            Point<dim>(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 <int dim>
     double
     Ellipsoid<dim>::evaluate_ellipsoid(const Point<dim> &point) const
@@ -179,16 +208,11 @@ namespace Functions
 
 
     template <int dim>
-    double
-    Ellipsoid<dim>::compute_signed_distance_ellipse(
-      const Point<dim> &point) const
+    Point<dim>
+    Ellipsoid<dim>::compute_closest_point_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
@@ -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<dim>(sign_px * x, sign_py * y);
+    }
+
+
+
+    template <int dim>
+    Tensor<1, dim, double>
+    Ellipsoid<dim>::compute_analyical_normal_vector_on_ellipse(
+      const Point<dim> &) 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 <int dim>
+    double
+    Ellipsoid<dim>::compute_signed_distance_ellipse(const Point<dim> &) 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 433a48e75e70fde242180b54e3fe214fda114074..563d0c66a3f089ec6d0661edca5b98bd6433ba11 100644 (file)
--- a/tests/base/function_signed_distance_02.cc
+++ b/tests/base/function_signed_distance_02.cc
@@ -27,14 +27,15 @@ namespace
 {
   template <int dim>
   void
-  print_value_at_point(const Function<dim> &function, const Point<dim> &point)
+  print_value_and_gradient_at_point(const Function<dim> &function,
+                                    const Point<dim> &   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<dim> 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<dim>(0.));
-    print_value_at_point(ellipsoid, Point<dim>(1.5));
+    print_value_and_gradient_at_point(ellipsoid, Point<dim>(0.));
+    print_value_and_gradient_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.));
+    print_value_and_gradient_at_point(ellipsoid, Point<dim>(-1.));
+    print_value_and_gradient_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.));
+    print_value_and_gradient_at_point(ellipsoid, Point<dim>(-2.));
+    print_value_and_gradient_at_point(ellipsoid, Point<dim>(6.));
 
     deallog << std::endl;
   }
@@ -81,20 +82,20 @@ namespace
     const Functions::SignedDistance::Ellipsoid<dim> 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 6b56eda45b933bca10facd0d35256d9cdbeb5bbc..ac98384e97f904b2af986ff78a76c32c5ed97d01 100644 (file)
--- 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::