From: Magdalena Schreter <schreter.magdalena@gmail.com> 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 <peter.muench@uni-a.de> --- 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. +<br> +(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<spacedim, Number> unit_to_real(const Point<spacedim, Number> &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<spacedim, Number> &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<spacedim, Number> &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 <deal.II/base/config.h> +#include <deal.II/base/bounding_box.h> #include <deal.II/base/function.h> #include <array> @@ -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 <int dim> + class Rectangle : public Function<dim> + { + 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<dim> &bottom_left, const Point<dim> &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<dim> & p, + const unsigned int component = 0) const override; + + private: + const BoundingBox<dim> 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 <deal.II/base/geometry_info.h> #include <limits> +#include <numeric> DEAL_II_NAMESPACE_OPEN @@ -335,6 +336,61 @@ BoundingBox<spacedim, Number>::unit_to_real( +template <int spacedim, typename Number> +Number +BoundingBox<spacedim, Number>::signed_distance( + const Point<spacedim, Number> &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 <int spacedim, typename Number> +Number +BoundingBox<spacedim, Number>::signed_distance( + const Point<spacedim, Number> &point) const +{ + // calculate vector of orthogonal signed distances + std::array<Number, spacedim> 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 <int dim, typename Number> BoundingBox<dim, Number> 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 <int dim> + Rectangle<dim>::Rectangle(const Point<dim> &bottom_left, + const Point<dim> &top_right) + : bounding_box({bottom_left, top_right}) + {} + + + + template <int dim> + double + Rectangle<dim>::value(const Point<dim> & 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<deal_II_dimension>; template class Functions::SignedDistance::Plane<deal_II_dimension>; template class Functions::SignedDistance::Ellipsoid<deal_II_dimension>; + template class Functions::SignedDistance::Rectangle<deal_II_dimension>; } 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 <deal.II/base/function_signed_distance.h> + +#include "../tests.h" + +namespace +{ + template <int dim> + void + print_value_at_point(const Function<dim> &function, const Point<dim> &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<dim> bl(0); + const Point<dim> tr(1); + const Point<dim> center(0.5); + + const Functions::SignedDistance::Rectangle<dim> rectangle(bl, tr); + + deallog << "center" << std::endl; + print_value_at_point(rectangle, center); + deallog << "inside" << std::endl; + print_value_at_point(rectangle, Point<dim>(0.25)); + print_value_at_point(rectangle, Point<dim>(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<dim>(-0.5)); + print_value_at_point(rectangle, Point<dim>(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<dim> bl(0, 0); + const Point<dim> tr(1, 1); + const Point<dim> center(0.5, 0.5); + + const Functions::SignedDistance::Rectangle<dim> rectangle(bl, tr); + + deallog << "center" << std::endl; + print_value_at_point(rectangle, center); + deallog << "inside" << std::endl; + print_value_at_point(rectangle, Point<dim>(0.25, 0.25)); + print_value_at_point(rectangle, Point<dim>(0.75, 0.25)); + print_value_at_point(rectangle, Point<dim>(0.75, 0.75)); + print_value_at_point(rectangle, Point<dim>(0.25, 0.75)); + deallog << "on rectangle" << std::endl; + print_value_at_point(rectangle, Point<dim>(0., 0.5)); + print_value_at_point(rectangle, Point<dim>(1.0, 0.5)); + print_value_at_point(rectangle, Point<dim>(1.0, 0.5123)); + print_value_at_point(rectangle, Point<dim>(0.0, 0.5123)); + deallog << "outside" << std::endl; + + for (double step = 0; step < 2.1; step += 0.5) + { + print_value_at_point(rectangle, Point<dim>(step, -1)); + print_value_at_point(rectangle, Point<dim>(step, 2)); + } + for (double step = 0; step < 2.1; step += 0.5) + { + print_value_at_point(rectangle, Point<dim>(-1, step)); + print_value_at_point(rectangle, Point<dim>(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<dim> bl(0, 0, 0); + const Point<dim> tr(1, 1, 1); + const Point<dim> center(0.5, 0.5, 0.5); + + const Functions::SignedDistance::Rectangle<dim> rectangle(bl, tr); + + deallog << "center" << std::endl; + print_value_at_point(rectangle, center); + deallog << "inside" << std::endl; + print_value_at_point(rectangle, Point<dim>(0.25, 0.25, 0.25)); + print_value_at_point(rectangle, Point<dim>(0.75, 0.25, 0.25)); + print_value_at_point(rectangle, Point<dim>(0.75, 0.75, 0.25)); + print_value_at_point(rectangle, Point<dim>(0.25, 0.75, 0.25)); + deallog << "on rectangle" << std::endl; + print_value_at_point(rectangle, Point<dim>(0., 0.5, 0.5)); + print_value_at_point(rectangle, Point<dim>(0.5, 0., 0.5)); + print_value_at_point(rectangle, Point<dim>(0.5, 0.5, 0.)); + print_value_at_point(rectangle, Point<dim>(1., 0.5, 0.5)); + print_value_at_point(rectangle, Point<dim>(0.5, 1., 0.5)); + print_value_at_point(rectangle, Point<dim>(0.5, 0.5, 1.)); + deallog << "outside" << std::endl; + print_value_at_point(rectangle, Point<dim>(-1, 0, 0)); + print_value_at_point(rectangle, Point<dim>(2, 0, 0)); + print_value_at_point(rectangle, Point<dim>(0, -1, 0)); + print_value_at_point(rectangle, Point<dim>(0, 2, 0)); + print_value_at_point(rectangle, Point<dim>(0, 0, -1)); + print_value_at_point(rectangle, Point<dim>(0, 0, 2)); + print_value_at_point(rectangle, Point<dim>(0, -1, 2)); + print_value_at_point(rectangle, Point<dim>(-1, 2, 0)); + print_value_at_point(rectangle, Point<dim>(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::