add signed distance for rectangle

Co-authored-by: Peter Munch <peter.muench@uni-a.de>

diff --git a/doc/doxygen/images/signed_distance_hyper_rectangle.png b/doc/doxygen/images/signed_distance_hyper_rectangle.png
new file mode 100644 (file)
index 0000000..c983c22
Binary files /dev/null and b/doc/doxygen/images/signed_distance_hyper_rectangle.png differ

diff --git a/doc/news/changes/minor/20230331Schreter b/doc/news/changes/minor/20230331Schreter
new file mode 100644 (file)
index 0000000..433e089
--- /dev/null
+++ b/doc/news/changes/minor/20230331Schreter
@@ -0,0 +1,4 @@
+New: Added Functions::SignedDistance::Rectangle() and BoundingBox::signed_distance()
+to compute the signed distance function of a rectangle.
+<br>
+(Magdalena Schreter, Peter Munch, 2023/03/31)

diff --git a/include/deal.II/base/bounding_box.h b/include/deal.II/base/bounding_box.h
index 24d5a7befa1896d2ad3830d6d8feb7b8f082eba0..610776612290a3b388df7e0fa037f5fa6c6268b2 100644 (file)
--- a/include/deal.II/base/bounding_box.h
+++ b/include/deal.II/base/bounding_box.h
@@ -337,6 +337,24 @@ public:
    */
   Point<spacedim, Number>
   unit_to_real(const Point<spacedim, Number> &point) const;
+  /**
+   * Returns the signed distance from a @p point orthogonal to the bounds of the
+   * box in @p direction. The signed distance is negative for points inside the
+   * interval described by the bounds of the rectangle in the respective
+   * direction, zero for points on the interval boundary and positive for points
+   * outside.
+   */
+  Number
+  signed_distance(const Point<spacedim, Number> &point,
+                  const unsigned int             direction) const;
+
+  /**
+   * Returns the signed distance from a @p point to the bounds of the box. The
+   * signed distance is negative for points inside the rectangle, zero for
+   * points on the rectangle and positive for points outside the rectangle.
+   */
+  Number
+  signed_distance(const Point<spacedim, Number> &point) const;
 
   /**
    * Write or read the data of this object to or from a stream for the

diff --git a/include/deal.II/base/function_signed_distance.h b/include/deal.II/base/function_signed_distance.h
index 06274083e9fbac4bccba834d1735595688f7b879..62cff816562225373564e763671c7aed84ba7073 100644 (file)
--- a/include/deal.II/base/function_signed_distance.h
+++ b/include/deal.II/base/function_signed_distance.h
@@ -18,6 +18,7 @@
 
 #include <deal.II/base/config.h>
 
+#include <deal.II/base/bounding_box.h>
 #include <deal.II/base/function.h>
 
 #include <array>
@@ -212,6 +213,47 @@ namespace Functions
       const double                  tolerance;
       const unsigned int            max_iter;
     };
+
+
+    /**
+     * Signed-distance level set function of a rectangle.
+     *
+     * This function is zero on the rectangle, negative "inside" and positive
+     * in the rest of $\mathbb{R}^{dim}$.
+     *
+     * Contour surfaces of the signed distance function of a 3D rectangle are
+     * illustrated below:
+     *
+     * @image html signed_distance_hyper_rectangle.png
+     *
+     * @ingroup functions
+     */
+    template <int dim>
+    class Rectangle : public Function<dim>
+    {
+    public:
+      /**
+       * Constructor, takes the bottom left point and the top right
+       * point of the rectangle.
+       *
+       * @param bottom_left Bottom left point of the rectangle.
+       * @param top_right Top right point of the rectangle.
+       */
+      Rectangle(const Point<dim> &bottom_left, const Point<dim> &top_right);
+
+      /**
+       * Calculates the signed distance from a given point @p p to the rectangle.
+       * The signed distance is negative for points inside the rectangle, zero
+       * for points on the rectangle and positive for points outside the
+       * rectangle.
+       */
+      double
+      value(const Point<dim> & p,
+            const unsigned int component = 0) const override;
+
+    private:
+      const BoundingBox<dim> bounding_box;
+    };
   } // namespace SignedDistance
 } // namespace Functions
 

diff --git a/source/base/bounding_box.cc b/source/base/bounding_box.cc
index 79eda17a7af91ca8fadc660e1491696e42958f48..7cd3f31114bcd557fcce45d83c422300806dc34c 100644 (file)
--- a/source/base/bounding_box.cc
+++ b/source/base/bounding_box.cc
@@ -17,6 +17,7 @@
 #include <deal.II/base/geometry_info.h>
 
 #include <limits>
+#include <numeric>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -335,6 +336,61 @@ BoundingBox<spacedim, Number>::unit_to_real(
 
 
 
+template <int spacedim, typename Number>
+Number
+BoundingBox<spacedim, Number>::signed_distance(
+  const Point<spacedim, Number> &point,
+  const unsigned int             direction) const
+{
+  const Number p1 = lower_bound(direction);
+  const Number p2 = upper_bound(direction);
+
+  if (point[direction] > p2)
+    return point[direction] - p2;
+  else if (point[direction] < p1)
+    return p1 - point[direction];
+  else
+    return -std::min(point[direction] - p1, p2 - point[direction]);
+}
+
+
+
+template <int spacedim, typename Number>
+Number
+BoundingBox<spacedim, Number>::signed_distance(
+  const Point<spacedim, Number> &point) const
+{
+  // calculate vector of orthogonal signed distances
+  std::array<Number, spacedim> distances;
+  for (unsigned int d = 0; d < spacedim; ++d)
+    distances[d] = signed_distance(point, d);
+
+  // determine the number of positive signed distances
+  const unsigned int n_positive_signed_distances =
+    std::count_if(distances.begin(), distances.end(), [](const auto &a) {
+      return a > 0.0;
+    });
+
+  if (n_positive_signed_distances <= 1)
+    // point is inside of bounding box (0: all signed distances are
+    // negative; find the index with the smallest absolute value)
+    // or next to a face (1: all signed distances are negative
+    // but one; find this index)
+    return *std::max_element(distances.begin(), distances.end());
+  else
+    // point is next to a corner (2D/3D: all signed distances are
+    // positive) or a line (3D: all signed distances are positive
+    // but one) -> take the l2-norm of all positive signed distances
+    return std::sqrt(std::accumulate(distances.begin(),
+                                     distances.end(),
+                                     0.0,
+                                     [](const auto &a, const auto &b) {
+                                       return a + (b > 0 ? b * b : 0.0);
+                                     }));
+}
+
+
+
 template <int dim, typename Number>
 BoundingBox<dim, Number>
 create_unit_bounding_box()

diff --git a/source/base/function_signed_distance.cc b/source/base/function_signed_distance.cc
index 8d17b5b407156f76852830b02179fa2fbb495a8b..48d094dc27b4c5e4ef53d13fb402e5fa3db621c0 100644 (file)
--- a/source/base/function_signed_distance.cc
+++ b/source/base/function_signed_distance.cc
@@ -332,6 +332,27 @@ namespace Functions
 
       return evaluate_ellipsoid(point) < 0.0 ? -distance : distance;
     }
+
+
+
+    template <int dim>
+    Rectangle<dim>::Rectangle(const Point<dim> &bottom_left,
+                              const Point<dim> &top_right)
+      : bounding_box({bottom_left, top_right})
+    {}
+
+
+
+    template <int dim>
+    double
+    Rectangle<dim>::value(const Point<dim> & p,
+                          const unsigned int component) const
+    {
+      AssertDimension(component, 0);
+      (void)component;
+
+      return bounding_box.signed_distance(p);
+    }
   } // namespace SignedDistance
 } // namespace Functions
 

diff --git a/source/base/function_signed_distance.inst.in b/source/base/function_signed_distance.inst.in
index 8c71662932959627306e6e20629fe340655c21ad..869e599494a399124762e4beb63834f008bf5c9f 100644 (file)
--- a/source/base/function_signed_distance.inst.in
+++ b/source/base/function_signed_distance.inst.in
@@ -18,4 +18,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>;
+    template class Functions::SignedDistance::Rectangle<deal_II_dimension>;
   }

diff --git a/tests/base/function_signed_distance_03.cc b/tests/base/function_signed_distance_03.cc
new file mode 100644 (file)
index 0000000..e134577
--- /dev/null
+++ b/tests/base/function_signed_distance_03.cc
@@ -0,0 +1,164 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2007 - 2022 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::Rectangle 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_rectangle_signed_distance_1d()
+  {
+    static constexpr int dim = 1;
+
+    deallog << "test_rectangle_signed_distance" << std::endl;
+    deallog << "dim = " << dim << std::endl;
+    deallog << std::endl;
+
+    const Point<dim> bl(0);
+    const Point<dim> tr(1);
+    const Point<dim> center(0.5);
+
+    const Functions::SignedDistance::Rectangle<dim> rectangle(bl, tr);
+
+    deallog << "center" << std::endl;
+    print_value_at_point(rectangle, center);
+    deallog << "inside" << std::endl;
+    print_value_at_point(rectangle, Point<dim>(0.25));
+    print_value_at_point(rectangle, Point<dim>(0.75));
+    deallog << "on rectangle" << std::endl;
+    print_value_at_point(rectangle, bl);
+    print_value_at_point(rectangle, tr);
+    deallog << "outside" << std::endl;
+    print_value_at_point(rectangle, Point<dim>(-0.5));
+    print_value_at_point(rectangle, Point<dim>(1.5));
+
+    deallog << std::endl;
+  }
+
+
+
+  void
+  test_rectangle_signed_distance_2d()
+  {
+    static constexpr int dim = 2;
+
+    deallog << "test_rectangle_signed_distance" << std::endl;
+    deallog << "dim = " << dim << std::endl;
+    deallog << std::endl;
+
+    const Point<dim> bl(0, 0);
+    const Point<dim> tr(1, 1);
+    const Point<dim> center(0.5, 0.5);
+
+    const Functions::SignedDistance::Rectangle<dim> rectangle(bl, tr);
+
+    deallog << "center" << std::endl;
+    print_value_at_point(rectangle, center);
+    deallog << "inside" << std::endl;
+    print_value_at_point(rectangle, Point<dim>(0.25, 0.25));
+    print_value_at_point(rectangle, Point<dim>(0.75, 0.25));
+    print_value_at_point(rectangle, Point<dim>(0.75, 0.75));
+    print_value_at_point(rectangle, Point<dim>(0.25, 0.75));
+    deallog << "on rectangle" << std::endl;
+    print_value_at_point(rectangle, Point<dim>(0., 0.5));
+    print_value_at_point(rectangle, Point<dim>(1.0, 0.5));
+    print_value_at_point(rectangle, Point<dim>(1.0, 0.5123));
+    print_value_at_point(rectangle, Point<dim>(0.0, 0.5123));
+    deallog << "outside" << std::endl;
+
+    for (double step = 0; step < 2.1; step += 0.5)
+      {
+        print_value_at_point(rectangle, Point<dim>(step, -1));
+        print_value_at_point(rectangle, Point<dim>(step, 2));
+      }
+    for (double step = 0; step < 2.1; step += 0.5)
+      {
+        print_value_at_point(rectangle, Point<dim>(-1, step));
+        print_value_at_point(rectangle, Point<dim>(2, step));
+      }
+
+    deallog << std::endl;
+  }
+
+  void
+  test_rectangle_signed_distance_3d()
+  {
+    static constexpr int dim = 3;
+
+    deallog << "test_rectangle_signed_distance" << std::endl;
+    deallog << "dim = " << dim << std::endl;
+    deallog << std::endl;
+
+    const Point<dim> bl(0, 0, 0);
+    const Point<dim> tr(1, 1, 1);
+    const Point<dim> center(0.5, 0.5, 0.5);
+
+    const Functions::SignedDistance::Rectangle<dim> rectangle(bl, tr);
+
+    deallog << "center" << std::endl;
+    print_value_at_point(rectangle, center);
+    deallog << "inside" << std::endl;
+    print_value_at_point(rectangle, Point<dim>(0.25, 0.25, 0.25));
+    print_value_at_point(rectangle, Point<dim>(0.75, 0.25, 0.25));
+    print_value_at_point(rectangle, Point<dim>(0.75, 0.75, 0.25));
+    print_value_at_point(rectangle, Point<dim>(0.25, 0.75, 0.25));
+    deallog << "on rectangle" << std::endl;
+    print_value_at_point(rectangle, Point<dim>(0., 0.5, 0.5));
+    print_value_at_point(rectangle, Point<dim>(0.5, 0., 0.5));
+    print_value_at_point(rectangle, Point<dim>(0.5, 0.5, 0.));
+    print_value_at_point(rectangle, Point<dim>(1., 0.5, 0.5));
+    print_value_at_point(rectangle, Point<dim>(0.5, 1., 0.5));
+    print_value_at_point(rectangle, Point<dim>(0.5, 0.5, 1.));
+    deallog << "outside" << std::endl;
+    print_value_at_point(rectangle, Point<dim>(-1, 0, 0));
+    print_value_at_point(rectangle, Point<dim>(2, 0, 0));
+    print_value_at_point(rectangle, Point<dim>(0, -1, 0));
+    print_value_at_point(rectangle, Point<dim>(0, 2, 0));
+    print_value_at_point(rectangle, Point<dim>(0, 0, -1));
+    print_value_at_point(rectangle, Point<dim>(0, 0, 2));
+    print_value_at_point(rectangle, Point<dim>(0, -1, 2));
+    print_value_at_point(rectangle, Point<dim>(-1, 2, 0));
+    print_value_at_point(rectangle, Point<dim>(2, 0, -1));
+    deallog << std::endl;
+  }
+} // namespace
+
+
+
+int
+main()
+{
+  initlog();
+  test_rectangle_signed_distance_1d();
+  test_rectangle_signed_distance_2d();
+  test_rectangle_signed_distance_3d();
+}

diff --git a/tests/base/function_signed_distance_03.output b/tests/base/function_signed_distance_03.output
new file mode 100644 (file)
index 0000000..068ba26
--- /dev/null
+++ b/tests/base/function_signed_distance_03.output
@@ -0,0 +1,137 @@
+
+DEAL::test_rectangle_signed_distance
+DEAL::dim = 1
+DEAL::
+DEAL::center
+DEAL::point = 0.500000
+DEAL::value = -0.500000
+DEAL::inside
+DEAL::point = 0.250000
+DEAL::value = -0.250000
+DEAL::point = 0.750000
+DEAL::value = -0.250000
+DEAL::on rectangle
+DEAL::point = 0.00000
+DEAL::value = 0.00000
+DEAL::point = 1.00000
+DEAL::value = 0.00000
+DEAL::outside
+DEAL::point = -0.500000
+DEAL::value = 0.500000
+DEAL::point = 1.50000
+DEAL::value = 0.500000
+DEAL::
+DEAL::test_rectangle_signed_distance
+DEAL::dim = 2
+DEAL::
+DEAL::center
+DEAL::point = 0.500000 0.500000
+DEAL::value = -0.500000
+DEAL::inside
+DEAL::point = 0.250000 0.250000
+DEAL::value = -0.250000
+DEAL::point = 0.750000 0.250000
+DEAL::value = -0.250000
+DEAL::point = 0.750000 0.750000
+DEAL::value = -0.250000
+DEAL::point = 0.250000 0.750000
+DEAL::value = -0.250000
+DEAL::on rectangle
+DEAL::point = 0.00000 0.500000
+DEAL::value = 0.00000
+DEAL::point = 1.00000 0.500000
+DEAL::value = 0.00000
+DEAL::point = 1.00000 0.512300
+DEAL::value = 0.00000
+DEAL::point = 0.00000 0.512300
+DEAL::value = 0.00000
+DEAL::outside
+DEAL::point = 0.00000 -1.00000
+DEAL::value = 1.00000
+DEAL::point = 0.00000 2.00000
+DEAL::value = 1.00000
+DEAL::point = 0.500000 -1.00000
+DEAL::value = 1.00000
+DEAL::point = 0.500000 2.00000
+DEAL::value = 1.00000
+DEAL::point = 1.00000 -1.00000
+DEAL::value = 1.00000
+DEAL::point = 1.00000 2.00000
+DEAL::value = 1.00000
+DEAL::point = 1.50000 -1.00000
+DEAL::value = 1.11803
+DEAL::point = 1.50000 2.00000
+DEAL::value = 1.11803
+DEAL::point = 2.00000 -1.00000
+DEAL::value = 1.41421
+DEAL::point = 2.00000 2.00000
+DEAL::value = 1.41421
+DEAL::point = -1.00000 0.00000
+DEAL::value = 1.00000
+DEAL::point = 2.00000 0.00000
+DEAL::value = 1.00000
+DEAL::point = -1.00000 0.500000
+DEAL::value = 1.00000
+DEAL::point = 2.00000 0.500000
+DEAL::value = 1.00000
+DEAL::point = -1.00000 1.00000
+DEAL::value = 1.00000
+DEAL::point = 2.00000 1.00000
+DEAL::value = 1.00000
+DEAL::point = -1.00000 1.50000
+DEAL::value = 1.11803
+DEAL::point = 2.00000 1.50000
+DEAL::value = 1.11803
+DEAL::point = -1.00000 2.00000
+DEAL::value = 1.41421
+DEAL::point = 2.00000 2.00000
+DEAL::value = 1.41421
+DEAL::
+DEAL::test_rectangle_signed_distance
+DEAL::dim = 3
+DEAL::
+DEAL::center
+DEAL::point = 0.500000 0.500000 0.500000
+DEAL::value = -0.500000
+DEAL::inside
+DEAL::point = 0.250000 0.250000 0.250000
+DEAL::value = -0.250000
+DEAL::point = 0.750000 0.250000 0.250000
+DEAL::value = -0.250000
+DEAL::point = 0.750000 0.750000 0.250000
+DEAL::value = -0.250000
+DEAL::point = 0.250000 0.750000 0.250000
+DEAL::value = -0.250000
+DEAL::on rectangle
+DEAL::point = 0.00000 0.500000 0.500000
+DEAL::value = 0.00000
+DEAL::point = 0.500000 0.00000 0.500000
+DEAL::value = 0.00000
+DEAL::point = 0.500000 0.500000 0.00000
+DEAL::value = 0.00000
+DEAL::point = 1.00000 0.500000 0.500000
+DEAL::value = 0.00000
+DEAL::point = 0.500000 1.00000 0.500000
+DEAL::value = 0.00000
+DEAL::point = 0.500000 0.500000 1.00000
+DEAL::value = 0.00000
+DEAL::outside
+DEAL::point = -1.00000 0.00000 0.00000
+DEAL::value = 1.00000
+DEAL::point = 2.00000 0.00000 0.00000
+DEAL::value = 1.00000
+DEAL::point = 0.00000 -1.00000 0.00000
+DEAL::value = 1.00000
+DEAL::point = 0.00000 2.00000 0.00000
+DEAL::value = 1.00000
+DEAL::point = 0.00000 0.00000 -1.00000
+DEAL::value = 1.00000
+DEAL::point = 0.00000 0.00000 2.00000
+DEAL::value = 1.00000
+DEAL::point = 0.00000 -1.00000 2.00000
+DEAL::value = 1.41421
+DEAL::point = -1.00000 2.00000 0.00000
+DEAL::value = 1.41421
+DEAL::point = 2.00000 0.00000 -1.00000
+DEAL::value = 1.41421
+DEAL::