From b1b5ac915c1246d5575f0ec6ff0468edfba2e2a0 Mon Sep 17 00:00:00 2001 From: Simon Sticko Date: Thu, 22 Apr 2021 15:28:56 +0200 Subject: [PATCH] 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. --- include/deal.II/base/function_level_set.h | 132 +++++++++++++++++++++ source/base/CMakeLists.txt | 2 + source/base/function_level_set.cc | 134 ++++++++++++++++++++++ source/base/function_level_set.inst.in | 20 ++++ tests/base/function_level_set.cc | 104 +++++++++++++++++ tests/base/function_level_set.output | 43 +++++++ 6 files changed, 435 insertions(+) create mode 100644 include/deal.II/base/function_level_set.h create mode 100644 source/base/function_level_set.cc create mode 100644 source/base/function_level_set.inst.in create mode 100644 tests/base/function_level_set.cc create mode 100644 tests/base/function_level_set.output diff --git a/include/deal.II/base/function_level_set.h b/include/deal.II/base/function_level_set.h new file mode 100644 index 0000000000..fb1219f68f --- /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 + +#include + +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 + class Sphere : public Function + { + public: + /** + * Constructor, takes the center and radius of the sphere. + */ + Sphere(const Point ¢er = Point(), const double radius = 1); + + double + value(const Point & 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 & 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 & point, + const unsigned int component = 0) const override; + + private: + const Point 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 + class Plane : public Function + { + public: + /** + * Constructor, takes a point in the plane and the plane normal. + */ + Plane(const Point &point, const Tensor<1, dim> &normal); + + double + value(const Point & point, + const unsigned int component = 0) const override; + + Tensor<1, dim> + gradient(const Point &, + const unsigned int component = 0) const override; + + SymmetricTensor<2, dim> + hessian(const Point &, + const unsigned int component = 0) const override; + + private: + const Point 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 6cc778da4b..d4e24e4952 100644 --- 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 index 0000000000..cb6aabeef4 --- /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 +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Functions +{ + namespace LevelSet + { + template + Sphere::Sphere(const Point ¢er, const double radius) + : center(center) + , radius(radius) + { + Assert(radius > 0, ExcMessage("Radius must be positive.")) + } + + + + template + double + Sphere::value(const Point & point, + const unsigned int component) const + { + AssertIndexRange(component, this->n_components); + (void)component; + + return point.distance(center) - radius; + } + + + + template + Tensor<1, dim> + Sphere::gradient(const Point & 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 + SymmetricTensor<2, dim> + Sphere::hessian(const Point & 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() / distance - + symmetrize(outer_product(center_to_point, center_to_point)) / + std::pow(distance, 3); + + return hess; + } + + + + template + Plane::Plane(const Point &point, const Tensor<1, dim> &normal) + : point_in_plane(point) + , normal(normal) + { + Assert(normal.norm() > 0, ExcMessage("Plane normal must not be 0.")) + } + + + + template + double + Plane::value(const Point & point, + const unsigned int component) const + { + AssertIndexRange(component, this->n_components); + (void)component; + + return normal * (point - point_in_plane); + } + + + + template + Tensor<1, dim> + Plane::gradient(const Point &, const unsigned int component) const + { + AssertIndexRange(component, this->n_components); + (void)component; + + return normal; + } + + + + template + SymmetricTensor<2, dim> + Plane::hessian(const Point &, 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 index 0000000000..418273bab5 --- /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; + template class Functions::LevelSet::Plane; + } diff --git a/tests/base/function_level_set.cc b/tests/base/function_level_set.cc new file mode 100644 index 0000000000..6e81091263 --- /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 + +#include "../tests.h" + + +namespace +{ + template + void + print_derivatives_at_point(const Function &function, + const Point & 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 + void + test_sphere_level_set() + { + deallog << "test_sphere_level_set" << std::endl; + + const Point center; + const double radius = 1; + + const Functions::LevelSet::Sphere level_set(center, radius); + + Point point; + point[0] = 2 * radius; + + print_derivatives_at_point(level_set, point); + } + + + + template + void + test_plane_level_set() + { + deallog << "test_plane_level_set" << std::endl; + + const Point point_in_plane; + const Tensor<1, dim> normal = Point::unit_vector(0); + + const Functions::LevelSet::Plane level_set(point_in_plane, normal); + + Point point; + for (unsigned int i = 0; i < dim; i++) + point[i] = 1; + + print_derivatives_at_point(level_set, point); + } + + + + template + void + run_test() + { + deallog << "dim = " << dim << std::endl; + deallog << std::endl; + test_sphere_level_set(); + deallog << std::endl; + test_plane_level_set(); + 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 index 0000000000..e8b578ef4b --- /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:: -- 2.39.5