add signed distance function for Zalesak's disk

doc/doxygen/references.bib

@@ -2033,3 +2033,14 @@ @article{zampini2016pcbddc
   year={2016},
   publisher={SIAM}
 }
+
+@article{zalesak1979fully,
+  title={Fully multidimensional flux-corrected transport algorithms for fluids},
+  author={Zalesak, Steven T},
+  journal={Journal of Computational Physics},
+  volume={31},
+  number={3},
+  pages={335--362},
+  year={1979},
+  publisher={Elsevier}
+}

doc/news/changes/minor/20230406Schreter

@@ -0,0 +1,4 @@
+New: Added Functions::SignedDistance::ZalesakDisk()
+to compute the signed distance function of Zalesak's disk.
+<br>
+(Magdalena Schreter, Simon Sticko, Peter Munch, 2023/04/06)

include/deal.II/base/function_signed_distance.h

@@ -240,6 +240,12 @@ namespace Functions
        * @param top_right Top right point of the rectangle.
        */
       Rectangle(const Point<dim> &bottom_left, const Point<dim> &top_right);
+      /**
+       * Constructor, takes a bounding box.
+       *
+       * @param bounding_box Bounding box.
+       */
+      Rectangle(const BoundingBox<dim> bounding_box);
 
       /**
        * Calculates the signed distance from a given point @p p to the rectangle.
@@ -254,6 +260,63 @@ namespace Functions
     private:
       const BoundingBox<dim> bounding_box;
     };
+
+
+    /**
+     * Signed-distance level set function of a Zalesak's disk proposed in
+     * @cite zalesak1979fully.
+     *
+     * It is calculated by the set difference $\psi(x) = \max(\psi_{S}(x),
+     * -\psi_{N}(x))$ of the level set functions of a sphere $\psi_{S}$ and a
+     * rectangle $\psi_{N}$. This function is zero on the surface of the disk,
+     * negative "inside" and positive in the rest of $\mathbb{R}^{dim}$.
+     *
+     * Contour surfaces of the signed distance function of a 3D Zalesak's disk
+     * are illustrated below:
+     *
+     * @image html signed_distance_zalesak_disk.png
+     *
+     * @ingroup functions
+     */
+    template <int dim>
+    class ZalesakDisk : public Function<dim>
+    {
+    public:
+      /**
+       * Constructor, takes the radius and center of the disk together
+       * with the width and the height of the notch.
+       *
+       * @param radius Radius of the disk.
+       * @param center Center of the disk.
+       * @param notch_width Width of the notch of the disk.
+       * @param notch_height Height of the notch of the disk.
+       */
+      ZalesakDisk(const Point<dim> &center,
+                  const double      radius,
+                  const double      notch_width,
+                  const double      notch_height);
+
+      /**
+       * Calculates the signed distance from a given point @p p to the disk.
+       * The signed distance is negative for points inside the disk, zero
+       * for points on the disk and positive for points outside the
+       * disk.
+       */
+      double
+      value(const Point<dim> & p,
+            const unsigned int component = 0) const override;
+
+    private:
+      /**
+       * Disk described by a sphere.
+       */
+      const Functions::SignedDistance::Sphere<dim> sphere;
+
+      /**
+       * Notch described by a rectangle.
+       */
+      const Functions::SignedDistance::Rectangle<dim> notch;
+    };
   } // namespace SignedDistance
 } // namespace Functions
 

source/base/function_signed_distance.cc

@@ -343,6 +343,13 @@ namespace Functions
 
 
 
+    template <int dim>
+    Rectangle<dim>::Rectangle(const BoundingBox<dim> bounding_box)
+      : bounding_box(bounding_box)
+    {}
+
+
+
     template <int dim>
     double
     Rectangle<dim>::value(const Point<dim> & p,
@@ -353,6 +360,60 @@ namespace Functions
 
       return bounding_box.signed_distance(p);
     }
+
+
+
+    template <int dim>
+    ZalesakDisk<dim>::ZalesakDisk(const Point<dim> &center,
+                                  const double      radius,
+                                  const double      notch_width,
+                                  const double      notch_height)
+      : sphere(center, radius)
+      , notch(// bottom left
+           (dim == 1) ?
+             Point<dim>(center[0] - 0.5 * notch_width) :
+           (dim == 2) ?
+             Point<dim>(center[0] - 0.5 * notch_width,
+                        std::numeric_limits<double>::lowest()) : /* notch is open in negative y-direction*/
+             Point<dim>(center[0] - 0.5 * notch_width,
+                        std::numeric_limits<
+                          double>::lowest() /* notch is open in negative y-direction*/,
+                        std::numeric_limits<double>::lowest()), /* notch is open in negative z-direction*/
+           // top right
+           (dim == 1) ?
+             Point<dim>(center[0] + 0.5 * notch_width) :
+           (dim == 2) ?
+             Point<dim>(center[0] + 0.5 * notch_width,
+                        center[1] + notch_height - radius) :
+             Point<dim>(center[0] + 0.5 * notch_width,
+                        std::numeric_limits<
+                          double>::max() /* notch is open in y-direction*/,
+                        center[2] + notch_height - radius))
+    {
+      Assert(
+        notch_width <= 2 * radius,
+        ExcMessage(
+          "The width of the notch must be less than the circle diameter."));
+      Assert(
+        notch_height <= 2 * radius,
+        ExcMessage(
+          "The height of the notch must be less than the circle diameter."));
+    }
+
+
+
+    template <int dim>
+    double
+    ZalesakDisk<dim>::value(const Point<dim> & p,
+                            const unsigned int component) const
+    {
+      (void)component;
+      AssertDimension(component, 0);
+
+      // calculate the set difference between the level set functions of the
+      // sphere and the notch
+      return std::max(sphere.value(p), -notch.value(p));
+    }
   } // namespace SignedDistance
 } // namespace Functions
 

source/base/function_signed_distance.inst.in

@@ -19,4 +19,5 @@ for (deal_II_dimension : SPACE_DIMENSIONS)
     template class Functions::SignedDistance::Plane<deal_II_dimension>;
     template class Functions::SignedDistance::Ellipsoid<deal_II_dimension>;
     template class Functions::SignedDistance::Rectangle<deal_II_dimension>;
+    template class Functions::SignedDistance::ZalesakDisk<deal_II_dimension>;
   }

tests/base/function_signed_distance_04.cc

@@ -0,0 +1,108 @@
+// ---------------------------------------------------------------------
+//
+// 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 classes Functions::SignedDistance::ZalesakDisk by creating an
+ * object of each class and calling the value functions at a point where we know
+ * what the function values should be.
+ */
+
+#include <deal.II/base/function_level_set.h>
+
+#include "../tests.h"
+
+
+namespace
+{
+  template <int dim>
+  void
+  print_values_at_point(const Function<dim> &function, const Point<dim> &point)
+  {
+    deallog << "point = " << point << std::endl;
+    deallog << "value = " << function.value(point) << std::endl;
+  }
+
+
+
+  template <int dim>
+  void
+  test_disk_level_set()
+  {
+    deallog << "test_zalesak_disk_level_set" << std::endl;
+
+    const Point<dim> center;
+    const double     radius = 1;
+
+    const Functions::SignedDistance::ZalesakDisk<dim> level_set(center,
+                                                                radius,
+                                                                radius * 0.5,
+                                                                0.75 * radius);
+
+    {
+      Point<dim> p;
+
+      deallog << "center" << std::endl;
+      print_values_at_point(level_set, p);
+      deallog << "inside" << std::endl;
+      p[0] += 0.1;
+      print_values_at_point(level_set, p);
+      p[0] = 0.3;
+      print_values_at_point(level_set, p);
+    }
+    deallog << "notch" << std::endl;
+    Point<dim> p;
+    p[dim - 1]              = -0.45;
+    p[std::max(0, dim - 2)] = -0.1;
+    print_values_at_point(level_set, p);
+    deallog << "outside" << std::endl;
+    p[dim - 1] = -2;
+    print_values_at_point(level_set, p);
+    p[dim - 1]              = -1.2;
+    p[0]                    = -0.3;
+    p[std::max(0, dim - 2)] = -0.3;
+    print_values_at_point(level_set, p);
+    p[dim - 1]              = -2;
+    p[0]                    = 0;
+    p[std::max(0, dim - 2)] = 0;
+    print_values_at_point(level_set, p);
+    p[dim - 1] = 2;
+    print_values_at_point(level_set, p);
+    p[0] = -2;
+    print_values_at_point(level_set, p);
+    p[0] = 2;
+    print_values_at_point(level_set, p);
+  }
+
+  template <int dim>
+  void
+  run_test()
+  {
+    deallog << "dim = " << dim << std::endl;
+    deallog << std::endl;
+    test_disk_level_set<dim>();
+    deallog << std::endl;
+  }
+} // namespace
+
+
+
+int
+main()
+{
+  initlog();
+  run_test<1>();
+  run_test<2>();
+  run_test<3>();
+}

tests/base/function_signed_distance_04.output

@@ -0,0 +1,85 @@
+
+DEAL::dim = 1
+DEAL::
+DEAL::test_zalesak_disk_level_set
+DEAL::center
+DEAL::point = 0.00000
+DEAL::value = 0.250000
+DEAL::inside
+DEAL::point = 0.100000
+DEAL::value = 0.150000
+DEAL::point = 0.300000
+DEAL::value = -0.0500000
+DEAL::notch
+DEAL::point = -0.100000
+DEAL::value = 0.150000
+DEAL::outside
+DEAL::point = -2.00000
+DEAL::value = 1.00000
+DEAL::point = -0.300000
+DEAL::value = -0.0500000
+DEAL::point = 0.00000
+DEAL::value = 0.250000
+DEAL::point = 2.00000
+DEAL::value = 1.00000
+DEAL::point = -2.00000
+DEAL::value = 1.00000
+DEAL::point = 2.00000
+DEAL::value = 1.00000
+DEAL::
+DEAL::dim = 2
+DEAL::
+DEAL::test_zalesak_disk_level_set
+DEAL::center
+DEAL::point = 0.00000 0.00000
+DEAL::value = -0.250000
+DEAL::inside
+DEAL::point = 0.100000 0.00000
+DEAL::value = -0.250000
+DEAL::point = 0.300000 0.00000
+DEAL::value = -0.254951
+DEAL::notch
+DEAL::point = -0.100000 -0.450000
+DEAL::value = 0.150000
+DEAL::outside
+DEAL::point = -0.100000 -2.00000
+DEAL::value = 1.00250
+DEAL::point = -0.300000 -1.20000
+DEAL::value = 0.236932
+DEAL::point = 0.00000 -2.00000
+DEAL::value = 1.00000
+DEAL::point = 0.00000 2.00000
+DEAL::value = 1.00000
+DEAL::point = -2.00000 2.00000
+DEAL::value = 1.82843
+DEAL::point = 2.00000 2.00000
+DEAL::value = 1.82843
+DEAL::
+DEAL::dim = 3
+DEAL::
+DEAL::test_zalesak_disk_level_set
+DEAL::center
+DEAL::point = 0.00000 0.00000 0.00000
+DEAL::value = -0.250000
+DEAL::inside
+DEAL::point = 0.100000 0.00000 0.00000
+DEAL::value = -0.250000
+DEAL::point = 0.300000 0.00000 0.00000
+DEAL::value = -0.254951
+DEAL::notch
+DEAL::point = 0.00000 -0.100000 -0.450000
+DEAL::value = 0.200000
+DEAL::outside
+DEAL::point = 0.00000 -0.100000 -2.00000
+DEAL::value = 1.00250
+DEAL::point = -0.300000 -0.300000 -1.20000
+DEAL::value = 0.272792
+DEAL::point = 0.00000 0.00000 -2.00000
+DEAL::value = 1.00000
+DEAL::point = 0.00000 0.00000 2.00000
+DEAL::value = 1.00000
+DEAL::point = -2.00000 0.00000 2.00000
+DEAL::value = 1.82843
+DEAL::point = 2.00000 0.00000 2.00000
+DEAL::value = 1.82843
+DEAL::