From: Magdalena Schreter Date: Fri, 31 Mar 2023 11:00:14 +0000 (+0200) Subject: add signed distance for rectangle X-Git-Tag: v9.5.0-rc1~373^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F15008%2Fhead;p=dealii.git add signed distance for rectangle Co-authored-by: Peter Munch --- diff --git a/doc/doxygen/images/signed_distance_hyper_rectangle.png b/doc/doxygen/images/signed_distance_hyper_rectangle.png new file mode 100644 index 0000000000..c983c222e5 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 index 0000000000..433e0893dd --- /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. +
+(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 24d5a7befa..6107766122 100644 --- a/include/deal.II/base/bounding_box.h +++ b/include/deal.II/base/bounding_box.h @@ -337,6 +337,24 @@ public: */ Point unit_to_real(const Point &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 &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 &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 06274083e9..62cff81656 100644 --- a/include/deal.II/base/function_signed_distance.h +++ b/include/deal.II/base/function_signed_distance.h @@ -18,6 +18,7 @@ #include +#include #include #include @@ -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 + class Rectangle : public Function + { + 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 &bottom_left, const Point &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 & p, + const unsigned int component = 0) const override; + + private: + const BoundingBox bounding_box; + }; } // namespace SignedDistance } // namespace Functions diff --git a/source/base/bounding_box.cc b/source/base/bounding_box.cc index 79eda17a7a..7cd3f31114 100644 --- a/source/base/bounding_box.cc +++ b/source/base/bounding_box.cc @@ -17,6 +17,7 @@ #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -335,6 +336,61 @@ BoundingBox::unit_to_real( +template +Number +BoundingBox::signed_distance( + const Point &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 +Number +BoundingBox::signed_distance( + const Point &point) const +{ + // calculate vector of orthogonal signed distances + std::array 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 BoundingBox create_unit_bounding_box() diff --git a/source/base/function_signed_distance.cc b/source/base/function_signed_distance.cc index 8d17b5b407..48d094dc27 100644 --- 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 + Rectangle::Rectangle(const Point &bottom_left, + const Point &top_right) + : bounding_box({bottom_left, top_right}) + {} + + + + template + double + Rectangle::value(const Point & 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 8c71662932..869e599494 100644 --- 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; template class Functions::SignedDistance::Plane; template class Functions::SignedDistance::Ellipsoid; + template class Functions::SignedDistance::Rectangle; } diff --git a/tests/base/function_signed_distance_03.cc b/tests/base/function_signed_distance_03.cc new file mode 100644 index 0000000000..e134577737 --- /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 + +#include "../tests.h" + +namespace +{ + template + void + print_value_at_point(const Function &function, const Point &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 bl(0); + const Point tr(1); + const Point center(0.5); + + const Functions::SignedDistance::Rectangle rectangle(bl, tr); + + deallog << "center" << std::endl; + print_value_at_point(rectangle, center); + deallog << "inside" << std::endl; + print_value_at_point(rectangle, Point(0.25)); + print_value_at_point(rectangle, Point(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(-0.5)); + print_value_at_point(rectangle, Point(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 bl(0, 0); + const Point tr(1, 1); + const Point center(0.5, 0.5); + + const Functions::SignedDistance::Rectangle rectangle(bl, tr); + + deallog << "center" << std::endl; + print_value_at_point(rectangle, center); + deallog << "inside" << std::endl; + print_value_at_point(rectangle, Point(0.25, 0.25)); + print_value_at_point(rectangle, Point(0.75, 0.25)); + print_value_at_point(rectangle, Point(0.75, 0.75)); + print_value_at_point(rectangle, Point(0.25, 0.75)); + deallog << "on rectangle" << std::endl; + print_value_at_point(rectangle, Point(0., 0.5)); + print_value_at_point(rectangle, Point(1.0, 0.5)); + print_value_at_point(rectangle, Point(1.0, 0.5123)); + print_value_at_point(rectangle, Point(0.0, 0.5123)); + deallog << "outside" << std::endl; + + for (double step = 0; step < 2.1; step += 0.5) + { + print_value_at_point(rectangle, Point(step, -1)); + print_value_at_point(rectangle, Point(step, 2)); + } + for (double step = 0; step < 2.1; step += 0.5) + { + print_value_at_point(rectangle, Point(-1, step)); + print_value_at_point(rectangle, Point(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 bl(0, 0, 0); + const Point tr(1, 1, 1); + const Point center(0.5, 0.5, 0.5); + + const Functions::SignedDistance::Rectangle rectangle(bl, tr); + + deallog << "center" << std::endl; + print_value_at_point(rectangle, center); + deallog << "inside" << std::endl; + print_value_at_point(rectangle, Point(0.25, 0.25, 0.25)); + print_value_at_point(rectangle, Point(0.75, 0.25, 0.25)); + print_value_at_point(rectangle, Point(0.75, 0.75, 0.25)); + print_value_at_point(rectangle, Point(0.25, 0.75, 0.25)); + deallog << "on rectangle" << std::endl; + print_value_at_point(rectangle, Point(0., 0.5, 0.5)); + print_value_at_point(rectangle, Point(0.5, 0., 0.5)); + print_value_at_point(rectangle, Point(0.5, 0.5, 0.)); + print_value_at_point(rectangle, Point(1., 0.5, 0.5)); + print_value_at_point(rectangle, Point(0.5, 1., 0.5)); + print_value_at_point(rectangle, Point(0.5, 0.5, 1.)); + deallog << "outside" << std::endl; + print_value_at_point(rectangle, Point(-1, 0, 0)); + print_value_at_point(rectangle, Point(2, 0, 0)); + print_value_at_point(rectangle, Point(0, -1, 0)); + print_value_at_point(rectangle, Point(0, 2, 0)); + print_value_at_point(rectangle, Point(0, 0, -1)); + print_value_at_point(rectangle, Point(0, 0, 2)); + print_value_at_point(rectangle, Point(0, -1, 2)); + print_value_at_point(rectangle, Point(-1, 2, 0)); + print_value_at_point(rectangle, Point(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 index 0000000000..068ba26338 --- /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::