From d8eaa2e77ca8e0e1959a0d591dfe94fd39671ff5 Mon Sep 17 00:00:00 2001 From: Magdalena Schreter Date: Thu, 6 Apr 2023 15:01:13 +0200 Subject: [PATCH] add signed distance function for Zalesak's disk --- doc/doxygen/references.bib | 11 ++ doc/news/changes/minor/20230406Schreter | 4 + .../deal.II/base/function_signed_distance.h | 63 ++++++++++ source/base/function_signed_distance.cc | 61 ++++++++++ source/base/function_signed_distance.inst.in | 1 + tests/base/function_signed_distance_04.cc | 108 ++++++++++++++++++ tests/base/function_signed_distance_04.output | 85 ++++++++++++++ 7 files changed, 333 insertions(+) create mode 100644 doc/news/changes/minor/20230406Schreter create mode 100644 tests/base/function_signed_distance_04.cc create mode 100644 tests/base/function_signed_distance_04.output diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index f4f11a0151..1433d28ad1 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -2033,3 +2033,14 @@ month = {apr} 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} +} diff --git a/doc/news/changes/minor/20230406Schreter b/doc/news/changes/minor/20230406Schreter new file mode 100644 index 0000000000..bd224af935 --- /dev/null +++ b/doc/news/changes/minor/20230406Schreter @@ -0,0 +1,4 @@ +New: Added Functions::SignedDistance::ZalesakDisk() +to compute the signed distance function of Zalesak's disk. +
+(Magdalena Schreter, Simon Sticko, Peter Munch, 2023/04/06) diff --git a/include/deal.II/base/function_signed_distance.h b/include/deal.II/base/function_signed_distance.h index 62cff81656..5d5d64dc91 100644 --- a/include/deal.II/base/function_signed_distance.h +++ b/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 &bottom_left, const Point &top_right); + /** + * Constructor, takes a bounding box. + * + * @param bounding_box Bounding box. + */ + Rectangle(const BoundingBox bounding_box); /** * Calculates the signed distance from a given point @p p to the rectangle. @@ -254,6 +260,63 @@ namespace Functions private: const BoundingBox bounding_box; }; + + + /** + * Signed-distance level set function of 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 https://www.dealii.org/images/grids/signed_distance_zalesak_disk.png + * + * @ingroup functions + */ + template + class ZalesakDisk : public Function + { + 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 ¢er, + 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 & p, + const unsigned int component = 0) const override; + + private: + /** + * Disk described by a sphere. + */ + const Functions::SignedDistance::Sphere sphere; + + /** + * Notch described by a rectangle. + */ + const Functions::SignedDistance::Rectangle notch; + }; } // namespace SignedDistance } // namespace Functions diff --git a/source/base/function_signed_distance.cc b/source/base/function_signed_distance.cc index 48d094dc27..3583b126bd 100644 --- a/source/base/function_signed_distance.cc +++ b/source/base/function_signed_distance.cc @@ -343,6 +343,13 @@ namespace Functions + template + Rectangle::Rectangle(const BoundingBox bounding_box) + : bounding_box(bounding_box) + {} + + + template double Rectangle::value(const Point & p, @@ -353,6 +360,60 @@ namespace Functions return bounding_box.signed_distance(p); } + + + + template + ZalesakDisk::ZalesakDisk(const Point ¢er, + const double radius, + const double notch_width, + const double notch_height) + : sphere(center, radius) + , notch(// bottom left + (dim == 1) ? + Point(center[0] - 0.5 * notch_width) : + (dim == 2) ? + Point(center[0] - 0.5 * notch_width, + std::numeric_limits::lowest()) : /* notch is open in negative y-direction*/ + Point(center[0] - 0.5 * notch_width, + std::numeric_limits< + double>::lowest() /* notch is open in negative y-direction*/, + std::numeric_limits::lowest()), /* notch is open in negative z-direction*/ + // top right + (dim == 1) ? + Point(center[0] + 0.5 * notch_width) : + (dim == 2) ? + Point(center[0] + 0.5 * notch_width, + center[1] + notch_height - radius) : + Point(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 + double + ZalesakDisk::value(const Point & 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 diff --git a/source/base/function_signed_distance.inst.in b/source/base/function_signed_distance.inst.in index 869e599494..cb5305d718 100644 --- a/source/base/function_signed_distance.inst.in +++ b/source/base/function_signed_distance.inst.in @@ -19,4 +19,5 @@ for (deal_II_dimension : SPACE_DIMENSIONS) template class Functions::SignedDistance::Plane; template class Functions::SignedDistance::Ellipsoid; template class Functions::SignedDistance::Rectangle; + template class Functions::SignedDistance::ZalesakDisk; } diff --git a/tests/base/function_signed_distance_04.cc b/tests/base/function_signed_distance_04.cc new file mode 100644 index 0000000000..293e3c0042 --- /dev/null +++ b/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 + +#include "../tests.h" + + +namespace +{ + template + void + print_values_at_point(const Function &function, const Point &point) + { + deallog << "point = " << point << std::endl; + deallog << "value = " << function.value(point) << std::endl; + } + + + + template + void + test_disk_level_set() + { + deallog << "test_zalesak_disk_level_set" << std::endl; + + const Point center; + const double radius = 1; + + const Functions::SignedDistance::ZalesakDisk level_set(center, + radius, + radius * 0.5, + 0.75 * radius); + + { + Point 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 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 + void + run_test() + { + deallog << "dim = " << dim << std::endl; + deallog << std::endl; + test_disk_level_set(); + deallog << std::endl; + } +} // namespace + + + +int +main() +{ + initlog(); + run_test<1>(); + run_test<2>(); + run_test<3>(); +} diff --git a/tests/base/function_signed_distance_04.output b/tests/base/function_signed_distance_04.output new file mode 100644 index 0000000000..608fcfca70 --- /dev/null +++ b/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:: -- 2.39.5