Add level set functions for a plane and a sphere

Add two level set functions in a new namespace
dealii::Functions::LevelSet, one being the signed distance function
from an arbitrary point to a sphere, and the other being the
(possibly scaled) signed distance from a point to a plane.

diff --git a/include/deal.II/base/function_level_set.h b/include/deal.II/base/function_level_set.h
new file mode 100644 (file)
index 0000000..fb1219f
--- /dev/null
+++ b/include/deal.II/base/function_level_set.h
@@ -0,0 +1,132 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 - 2020 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_function_level_set_h
+#define dealii_function_level_set_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/function.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Functions
+{
+  namespace LevelSet
+  {
+    /**
+     * Signed-distance level set function of a sphere:
+     * $\psi(x) = \| x - x^c \| - R$.
+     * Here, $x^c$ is the center of the sphere and $R$ is its radius. This
+     * function is thus zero on the sphere, negative "inside" the ball having
+     * the sphere as its boundary, and positive in the rest of
+     * $\mathbb{R}^{dim}$.
+     *
+     * This function has gradient and Hessian equal to
+     * $\partial_i \psi(x) = (x - x^c)/\| x - x^c \|$,
+     * $\partial_i \partial_j \psi =
+     * \delta_{ij}/\| x - x^c \| - (x_i - x_i^c)(x_j - x_j^c)/\| x - x^c \|^3$,
+     * where $\delta_{ij}$ is the Kronecker delta function.
+     *
+     * @ingroup functions
+     */
+    template <int dim>
+    class Sphere : public Function<dim>
+    {
+    public:
+      /**
+       * Constructor, takes the center and radius of the sphere.
+       */
+      Sphere(const Point<dim> &center = Point<dim>(), const double radius = 1);
+
+      double
+      value(const Point<dim> & point,
+            const unsigned int component = 0) const override;
+
+      /**
+       * @copydoc Function::gradient()
+       *
+       * @note The gradient is singular at the center of the sphere. If the
+       * incoming @p point is too close to the center, a floating-point
+       * exception may be thrown or entries in the gradient may be +inf/-inf
+       * or +nan/-nan, depending on how the point is located relative to the
+       * singularity.
+       */
+      Tensor<1, dim>
+      gradient(const Point<dim> & point,
+               const unsigned int component = 0) const override;
+
+      /**
+       * @copydoc Function::hessian()
+       *
+       * @note The Hessian is singular at the center of the sphere. If the
+       * incoming @p point is too close to the center, a floating-point
+       * exception may be thrown or entries in the Hessian may be +inf/-inf
+       * or +nan/-nan, depending on how the point is located relative to the
+       * singularity.
+       */
+      SymmetricTensor<2, dim>
+      hessian(const Point<dim> & point,
+              const unsigned int component = 0) const override;
+
+    private:
+      const Point<dim> center;
+      const double     radius;
+    };
+
+
+    /**
+     * Signed level set function of a plane in $\mathbb{R}^{dim}$:
+     * $\psi(x) = n \cdot (x - x_p)$.
+     * Here, $n$ is the plane normal and $x_p$ is a point in the plane.
+     * Thus, with respect to the direction of the normal, this function is
+     * positive above the plane, zero in the plane, and negative below the
+     * plane. If the normal is normalized, $\psi$ will be the signed distance to
+     * the closest point in the plane.
+     *
+     * @ingroup functions
+     */
+    template <int dim>
+    class Plane : public Function<dim>
+    {
+    public:
+      /**
+       * Constructor, takes a point in the plane and the plane normal.
+       */
+      Plane(const Point<dim> &point, const Tensor<1, dim> &normal);
+
+      double
+      value(const Point<dim> & point,
+            const unsigned int component = 0) const override;
+
+      Tensor<1, dim>
+      gradient(const Point<dim> &,
+               const unsigned int component = 0) const override;
+
+      SymmetricTensor<2, dim>
+      hessian(const Point<dim> &,
+              const unsigned int component = 0) const override;
+
+    private:
+      const Point<dim>     point_in_plane;
+      const Tensor<1, dim> normal;
+    };
+
+  } // namespace LevelSet
+} // namespace Functions
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif /* CE8BBB3E_B726_40A7_B963_561AE7B84973 */

diff --git a/source/base/CMakeLists.txt b/source/base/CMakeLists.txt
index 6cc778da4bb17374fcb610032d2c4d639103aff0..d4e24e4952134840b83be211176e325cd7214f35 100644 (file)
--- a/source/base/CMakeLists.txt
+++ b/source/base/CMakeLists.txt
@@ -32,6 +32,7 @@ SET(_unity_include_src
   function.cc
   function_cspline.cc
   function_derivative.cc
+  function_level_set.cc
   function_lib.cc
   function_lib_cutoff.cc
   function_parser.cc
@@ -128,6 +129,7 @@ SET(_inst
   bounding_box.inst.in
   data_out_base.inst.in
   function.inst.in
+  function_level_set.inst.in
   function_time.inst.in
   incremental_function.inst.in
   geometric_utilities.inst.in

diff --git a/source/base/function_level_set.cc b/source/base/function_level_set.cc
new file mode 100644 (file)
index 0000000..cb6aabe
--- /dev/null
+++ b/source/base/function_level_set.cc
@@ -0,0 +1,134 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 - 2020 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/function_level_set.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Functions
+{
+  namespace LevelSet
+  {
+    template <int dim>
+    Sphere<dim>::Sphere(const Point<dim> &center, const double radius)
+      : center(center)
+      , radius(radius)
+    {
+      Assert(radius > 0, ExcMessage("Radius must be positive."))
+    }
+
+
+
+    template <int dim>
+    double
+    Sphere<dim>::value(const Point<dim> & point,
+                       const unsigned int component) const
+    {
+      AssertIndexRange(component, this->n_components);
+      (void)component;
+
+      return point.distance(center) - radius;
+    }
+
+
+
+    template <int dim>
+    Tensor<1, dim>
+    Sphere<dim>::gradient(const Point<dim> & point,
+                          const unsigned int component) const
+    {
+      AssertIndexRange(component, this->n_components);
+      (void)component;
+
+      const Tensor<1, dim> center_to_point = point - center;
+      const Tensor<1, dim> grad = center_to_point / center_to_point.norm();
+      return grad;
+    }
+
+
+
+    template <int dim>
+    SymmetricTensor<2, dim>
+    Sphere<dim>::hessian(const Point<dim> & point,
+                         const unsigned int component) const
+    {
+      AssertIndexRange(component, this->n_components);
+      (void)component;
+
+      const Tensor<1, dim> center_to_point = point - center;
+      const double         distance        = center_to_point.norm();
+
+      const SymmetricTensor<2, dim> hess =
+        unit_symmetric_tensor<dim>() / distance -
+        symmetrize(outer_product(center_to_point, center_to_point)) /
+          std::pow(distance, 3);
+
+      return hess;
+    }
+
+
+
+    template <int dim>
+    Plane<dim>::Plane(const Point<dim> &point, const Tensor<1, dim> &normal)
+      : point_in_plane(point)
+      , normal(normal)
+    {
+      Assert(normal.norm() > 0, ExcMessage("Plane normal must not be 0."))
+    }
+
+
+
+    template <int dim>
+    double
+    Plane<dim>::value(const Point<dim> & point,
+                      const unsigned int component) const
+    {
+      AssertIndexRange(component, this->n_components);
+      (void)component;
+
+      return normal * (point - point_in_plane);
+    }
+
+
+
+    template <int dim>
+    Tensor<1, dim>
+    Plane<dim>::gradient(const Point<dim> &, const unsigned int component) const
+    {
+      AssertIndexRange(component, this->n_components);
+      (void)component;
+
+      return normal;
+    }
+
+
+
+    template <int dim>
+    SymmetricTensor<2, dim>
+    Plane<dim>::hessian(const Point<dim> &, const unsigned int component) const
+    {
+      AssertIndexRange(component, this->n_components);
+      (void)component;
+
+      return SymmetricTensor<2, dim>();
+    }
+
+  } // namespace LevelSet
+} // namespace Functions
+
+#include "function_level_set.inst"
+
+DEAL_II_NAMESPACE_CLOSE

diff --git a/source/base/function_level_set.inst.in b/source/base/function_level_set.inst.in
new file mode 100644 (file)
index 0000000..418273b
--- /dev/null
+++ b/source/base/function_level_set.inst.in
@@ -0,0 +1,20 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 - 2020 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.
+//
+// ---------------------------------------------------------------------
+
+for (deal_II_dimension : SPACE_DIMENSIONS)
+  {
+    template class Functions::LevelSet::Sphere<deal_II_dimension>;
+    template class Functions::LevelSet::Plane<deal_II_dimension>;
+  }

diff --git a/tests/base/function_level_set.cc b/tests/base/function_level_set.cc
new file mode 100644 (file)
index 0000000..6e81091
--- /dev/null
+++ b/tests/base/function_level_set.cc
@@ -0,0 +1,104 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2007 - 2018 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::LevelSet::Sphere and Functions::LevelSet::Plane,
+ * by creating an object of each class and calling the value, gradient and
+ * hessian 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_derivatives_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;
+    deallog << "Hessian = " << function.hessian(point) << std::endl;
+  }
+
+
+
+  template <int dim>
+  void
+  test_sphere_level_set()
+  {
+    deallog << "test_sphere_level_set" << std::endl;
+
+    const Point<dim> center;
+    const double     radius = 1;
+
+    const Functions::LevelSet::Sphere<dim> level_set(center, radius);
+
+    Point<dim> point;
+    point[0] = 2 * radius;
+
+    print_derivatives_at_point(level_set, point);
+  }
+
+
+
+  template <int dim>
+  void
+  test_plane_level_set()
+  {
+    deallog << "test_plane_level_set" << std::endl;
+
+    const Point<dim>     point_in_plane;
+    const Tensor<1, dim> normal = Point<dim>::unit_vector(0);
+
+    const Functions::LevelSet::Plane<dim> level_set(point_in_plane, normal);
+
+    Point<dim> point;
+    for (unsigned int i = 0; i < dim; i++)
+      point[i] = 1;
+
+    print_derivatives_at_point(level_set, point);
+  }
+
+
+
+  template <int dim>
+  void
+  run_test()
+  {
+    deallog << "dim = " << dim << std::endl;
+    deallog << std::endl;
+    test_sphere_level_set<dim>();
+    deallog << std::endl;
+    test_plane_level_set<dim>();
+    deallog << std::endl;
+  }
+} // namespace
+
+
+
+int
+main()
+{
+  initlog();
+  run_test<1>();
+  run_test<2>();
+  run_test<3>();
+}

diff --git a/tests/base/function_level_set.output b/tests/base/function_level_set.output
new file mode 100644 (file)
index 0000000..e8b578e
--- /dev/null
+++ b/tests/base/function_level_set.output
@@ -0,0 +1,43 @@
+
+DEAL::dim = 1
+DEAL::
+DEAL::test_sphere_level_set
+DEAL::point = 2.00000
+DEAL::value = 1.00000
+DEAL::gradient = 1.00000
+DEAL::Hessian = 0.00000
+DEAL::
+DEAL::test_plane_level_set
+DEAL::point = 1.00000
+DEAL::value = 1.00000
+DEAL::gradient = 1.00000
+DEAL::Hessian = 0.00000
+DEAL::
+DEAL::dim = 2
+DEAL::
+DEAL::test_sphere_level_set
+DEAL::point = 2.00000 0.00000
+DEAL::value = 1.00000
+DEAL::gradient = 1.00000 0.00000
+DEAL::Hessian = 0.00000 0.00000 0.00000 0.500000
+DEAL::
+DEAL::test_plane_level_set
+DEAL::point = 1.00000 1.00000
+DEAL::value = 1.00000
+DEAL::gradient = 1.00000 0.00000
+DEAL::Hessian = 0.00000 0.00000 0.00000 0.00000
+DEAL::
+DEAL::dim = 3
+DEAL::
+DEAL::test_sphere_level_set
+DEAL::point = 2.00000 0.00000 0.00000
+DEAL::value = 1.00000
+DEAL::gradient = 1.00000 0.00000 0.00000
+DEAL::Hessian = 0.00000 0.00000 0.00000 0.00000 0.500000 0.00000 0.00000 0.00000 0.500000
+DEAL::
+DEAL::test_plane_level_set
+DEAL::point = 1.00000 1.00000 1.00000
+DEAL::value = 1.00000
+DEAL::gradient = 1.00000 0.00000 0.00000
+DEAL::Hessian = 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL::