From af7a95cb2e0038adf76ad067b3d0e3b9190b0d73 Mon Sep 17 00:00:00 2001 From: Simon Sticko Date: Sun, 5 Apr 2020 21:14:54 +0200 Subject: [PATCH] Implement functionality from Hypercube on BoundingBox. The two classes are very similar. Implement the functionality from Hypercube on Boundingbox so that the Hypercube class can be removed. --- include/deal.II/base/bounding_box.h | 153 ++++++++++- source/base/bounding_box.cc | 156 ++++++++++++ source/base/bounding_box.inst.in | 2 + tests/base/bounding_box_6.cc | 241 ++++++++++++++++++ tests/base/bounding_box_6.output | 124 +++++++++ .../coordinate_to_one_dim_higher.cc | 6 +- .../coordinate_to_one_dim_higher.output | 0 7 files changed, 666 insertions(+), 16 deletions(-) create mode 100644 tests/base/bounding_box_6.cc create mode 100644 tests/base/bounding_box_6.output rename tests/{non_matching => base}/coordinate_to_one_dim_higher.cc (84%) rename tests/{non_matching => base}/coordinate_to_one_dim_higher.output (100%) diff --git a/include/deal.II/base/bounding_box.h b/include/deal.II/base/bounding_box.h index 60ee0bb4a4..62470ee7b7 100644 --- a/include/deal.II/base/bounding_box.h +++ b/include/deal.II/base/bounding_box.h @@ -69,24 +69,27 @@ enum class NeighborType }; /** - * A class that represents a bounding box in a space with arbitrary dimension - * spacedim. + * A class that represents a box of arbitrary dimension spacedim and + * with sides parallel to the coordinate axes. That is, a region * - * Objects of this class are used to represent bounding boxes. They are, - * among other uses, useful in parallel distributed meshes to give a general - * description of the owners of each portion of the mesh. + * @f[ + * [x_0^L, x_0^U] \times ... \times [x_{spacedim-1}^L, x_{spacedim-1}^U], + * @f] * - * Bounding boxes are represented by two vertices (bottom left and top right). - * Geometrically, a bounding box is: - * - 1 d: a segment (represented by its vertices in the proper order) - * - 2 d: a rectangle (represented by the vertices V at bottom left, top right) + * where $(x_0^L , ..., x_{spacedim-1}^L) and $(x_0^U , ..., x_{spacedim-1}^U) + * denote the two vertices (bottom left and top right) which are used to + * represent the box. + * + * Geometrically, a bounding box is thus: + * - 1D: a segment (represented by its vertices in the proper order) + * - 2D: a rectangle (represented by the vertices V at bottom left, top right) * @code * .--------V * | | * V--------. * @endcode * - * - 3 d: a cuboid (in which case the two vertices V follow the convention and + * - 3D: a cuboid (in which case the two vertices V follow the convention and * are not owned by the same face) * @code * .------V @@ -96,9 +99,30 @@ enum class NeighborType * | |/ * V------. * @endcode - * Notice that the sides are always parallel to the respective axis. * - * @author Giovanni Alzetta, 2017. + * Bounding boxes are, for example, useful in parallel distributed meshes to + * give a general description of the owners of each portion of the mesh. + * + * Taking the cross section of a BoundingBox orthogonal to a given + * direction gives a box in one dimension lower: BoundingBox. + * In 3D, the 2 coordinates of the cross section of BoundingBox<3> can be + * ordered in 2 different ways. That is, if we take the cross section orthogonal + * to the y direction we could either order a 3D-coordinate into a + * 2D-coordinate as $(x,z)$ or as $(z,x)$. This class uses the second + * convention, corresponding to the coordinates being ordered cyclicly + * $x \rightarrow y \rightarrow z \rightarrow x \rightarrow ... $ + * To be precise, if we take a cross section: + * + * | Orthogonal to | Cross section coordinates ordered as | + * |:-------------:|:------------------------------------:| + * | x | (y, z) | + * | y | (z, x) | + * | z | (x, y) | + * + * This is according to the convention set by the function + * coordinate_to_one_dim_higher. + * + * @author Giovanni Alzetta, 2017, Simon Sticko 2020. */ template class BoundingBox @@ -182,6 +206,56 @@ public: double volume() const; + /** + * Returns the point in the center of the box. + */ + Point + center() const; + + /** + * Returns the side length of the box in @p direction. + */ + Number + side_length(const unsigned int direction) const; + + /** + * Return the lower bound of the box in @p direction. + */ + Number + lower_bound(const unsigned int direction) const; + + /** + * Return the upper bound of the box in @p direction. + */ + Number + upper_bound(const unsigned int direction) const; + + /** + * Return the bounds of the box in @p direction, as a one-dimensional box. + */ + BoundingBox<1, Number> + bounds(const unsigned int direction) const; + + /** + * Returns the indexth vertex of the box. + */ + Point + vertex(const unsigned int index) const; + + /** + * Returns the indexth child of the box. Child is meant in the same way as for + * a cell. + */ + BoundingBox + child(const unsigned int index) const; + + /** + * Returns the cross section of the box orthogonal to @p direction. + * This is a box in one dimension lower. + */ + BoundingBox + cross_section(const unsigned int direction) const; + /** * Boost serialization function */ @@ -193,6 +267,61 @@ private: std::pair, Point> boundary_points; }; + +/** + * Returns the unit box [0,1]^dim. + */ +template +BoundingBox +create_unit_box(); + + +/** + * This function defines a convention for how coordinates in dim + * dimensions should translate to the coordinates in dim + 1 dimensions, + * when one of the coordinates in dim + 1 dimensions is locked to a given + * value. + * + * The convention is the following: Starting from the locked coordinate we + * store the lower dimensional coordinates consecutively and wrapping + * around when going over the dimension. This relationship can be + * described by the following tables: + * + * 2D + * ------------------------------------------------ + * | locked in 2D | 1D coordinate | 2D coordinate | + * |--------------------------------------------- | + * | x0 | (a) | (x0, a) | + * | x1 | (a) | (a , x1) | + * ------------------------------------------------ + * + * 3D + * -------------------------------------------------- + * | locked in 3D | 2D coordinates | 3D coordinates | + * |----------------------------------------------- | + * | x0 | (a, b) | (x0, a, b) | + * | x1 | (a, b) | ( b, x1, a) | + * | x2 | (a, b) | ( a, b, x2) | + * -------------------------------------------------- + * + * Given a locked coordinate, this function maps a coordinate index in dim + * dimension to a coordinate index in dim + 1 dimensions. + * + * @param locked_coordinate should be in the range [0, dim+1). + * @param coordiante_in_dim should be in the range [0, dim). + * @return A coordinate index in the range [0, dim+1) + */ +template +inline unsigned int +coordinate_to_one_dim_higher(const unsigned int locked_coordinate, + const unsigned int coordiante_in_dim) +{ + AssertIndexRange(locked_coordinate, dim + 1); + AssertIndexRange(coordiante_in_dim, dim); + return (locked_coordinate + coordiante_in_dim + 1) % (dim + 1); +} + + /*------------------------ Inline functions: BoundingBox --------------------*/ #ifndef DOXYGEN diff --git a/source/base/bounding_box.cc b/source/base/bounding_box.cc index 6720c9e878..978bc5c23c 100644 --- a/source/base/bounding_box.cc +++ b/source/base/bounding_box.cc @@ -13,6 +13,7 @@ // // --------------------------------------------------------------------- #include +#include DEAL_II_NAMESPACE_OPEN @@ -172,5 +173,160 @@ BoundingBox::volume() const return vol; } + + +template +Number +BoundingBox::lower_bound(const unsigned int direction) const +{ + AssertIndexRange(direction, spacedim); + + return boundary_points.first[direction]; +} + + + +template +Number +BoundingBox::upper_bound(const unsigned int direction) const +{ + AssertIndexRange(direction, spacedim); + + return boundary_points.second[direction]; +} + + + +template +Point +BoundingBox::center() const +{ + Point point; + for (unsigned int i = 0; i < spacedim; ++i) + point[i] = .5 * (boundary_points.first[i] + boundary_points.second[i]); + + return point; +} + + + +template +BoundingBox<1, Number> +BoundingBox::bounds(const unsigned int direction) const +{ + AssertIndexRange(direction, spacedim); + + std::pair, Point<1, Number>> lower_upper_bounds; + lower_upper_bounds.first[0] = lower_bound(direction); + lower_upper_bounds.second[0] = upper_bound(direction); + + return BoundingBox<1, Number>(lower_upper_bounds); +} + + + +template +Number +BoundingBox::side_length(const unsigned int direction) const +{ + AssertIndexRange(direction, spacedim); + + return boundary_points.second[direction] - boundary_points.first[direction]; +} + + + +template +Point +BoundingBox::vertex(const unsigned int index) const +{ + AssertIndexRange(index, GeometryInfo::vertices_per_cell); + + const Point unit_cell_vertex = + GeometryInfo::unit_cell_vertex(index); + + Point point; + for (unsigned int i = 0; i < spacedim; ++i) + point[i] = boundary_points.first[i] + side_length(i) * unit_cell_vertex[i]; + + return point; +} + + + +template +BoundingBox +BoundingBox::child(const unsigned int index) const +{ + AssertIndexRange(index, GeometryInfo::max_children_per_cell); + + // Vertex closest to child. + const Point parent_vertex = vertex(index); + const Point parent_center = center(); + + const Point upper_corner_unit_cell = + GeometryInfo::unit_cell_vertex( + GeometryInfo::vertices_per_cell - 1); + + const Point lower_corner_unit_cell = + GeometryInfo::unit_cell_vertex(0); + + std::pair, Point> + child_lower_upper_corner; + for (unsigned int i = 0; i < spacedim; ++i) + { + const double child_side_length = side_length(i) / 2; + + const double child_center = (parent_center[i] + parent_vertex[i]) / 2; + + child_lower_upper_corner.first[i] = + child_center + child_side_length * (lower_corner_unit_cell[i] - .5); + child_lower_upper_corner.second[i] = + child_center + child_side_length * (upper_corner_unit_cell[i] - .5); + } + + return BoundingBox(child_lower_upper_corner); +} + + + +template +BoundingBox +BoundingBox::cross_section(const unsigned int direction) const +{ + AssertIndexRange(direction, spacedim); + + std::pair, Point> + cross_section_lower_upper_corner; + for (unsigned int d = 0; d < spacedim - 1; ++d) + { + const int index_to_write_from = + coordinate_to_one_dim_higher(direction, d); + + cross_section_lower_upper_corner.first[d] = + boundary_points.first[index_to_write_from]; + + cross_section_lower_upper_corner.second[d] = + boundary_points.second[index_to_write_from]; + } + + return BoundingBox(cross_section_lower_upper_corner); +} + + + +template +BoundingBox +create_unit_box() +{ + std::pair, Point> lower_upper_corner; + for (unsigned int i = 0; i < dim; ++i) + { + lower_upper_corner.second[i] = 1; + } + return BoundingBox(lower_upper_corner); +} + + #include "bounding_box.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/source/base/bounding_box.inst.in b/source/base/bounding_box.inst.in index 162673be27..f44e86ffa7 100644 --- a/source/base/bounding_box.inst.in +++ b/source/base/bounding_box.inst.in @@ -17,4 +17,6 @@ for (deal_II_dimension : SPACE_DIMENSIONS; number : REAL_SCALARS) { template class BoundingBox; + + template BoundingBox create_unit_box(); } diff --git a/tests/base/bounding_box_6.cc b/tests/base/bounding_box_6.cc new file mode 100644 index 0000000000..ebeeac7004 --- /dev/null +++ b/tests/base/bounding_box_6.cc @@ -0,0 +1,241 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 - 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 following functions of the BoundingBox class +// cross_section +// center +// side_length +// child +// vertex +// bounds +// and the non-member but related function create_unit_box +// Each function is tested in the function named test_{member_function_name}. + +#include +#include + +#include "../tests.h" + +using namespace dealii; + +// Print the bounds of the incoming box to deallog. +template +void +print_box(const BoundingBox &box) +{ + for (unsigned int d = 0; d < dim; ++d) + { + if (d > 0) + deallog << "x"; + deallog << "[" << box.lower_bound(d) << ", " << box.upper_bound(d) << "]"; + } + deallog << std::endl; +} + + + +// Construct the cross section orthogonal to all axes and print them. +template +void +test_cross_section() +{ + deallog << "test_cross_section" << std::endl; + + std::pair, Point> lower_upper_corners; + for (int d = 0; d < dim; ++d) + { + lower_upper_corners.first[d] = -(d + 1); + lower_upper_corners.second[d] = d + 1; + } + + const BoundingBox box(lower_upper_corners); + + for (unsigned int d = 0; d < dim; ++d) + { + deallog << "orthogonal to " << d << std::endl; + const BoundingBox cross_section = box.cross_section(d); + print_box(cross_section); + } +} + + + +// Compute and print the center of the box. +template +void +test_center() +{ + deallog << "test_center" << std::endl; + + std::pair, Point> lower_upper_corners; + for (int d = 0; d < dim; ++d) + lower_upper_corners.second[d] = 1; + + const BoundingBox box(lower_upper_corners); + + deallog << "center " << box.center() << std::endl; +} + + + +// Print all side lengths of a box. +template +void +test_side_length() +{ + deallog << "test_side_length" << std::endl; + + std::pair, Point> lower_upper_corners; + for (int d = 0; d < dim; ++d) + lower_upper_corners.second[d] = d + 1; + + const BoundingBox box(lower_upper_corners); + + for (unsigned int d = 0; d < dim; ++d) + deallog << box.side_length(d) << " "; + deallog << std::endl; +} + + + +// Construct all possible children of a box and print them. +template +void +test_child() +{ + deallog << "test_child" << std::endl; + + std::pair, Point> lower_upper_corners; + for (int d = 0; d < dim; ++d) + { + lower_upper_corners.first[d] = -(d + 1); + lower_upper_corners.second[d] = (d + 1); + } + + const BoundingBox box(lower_upper_corners); + + for (unsigned int i = 0; i < GeometryInfo::max_children_per_cell; ++i) + { + deallog << "child " << i << std::endl; + const BoundingBox child = box.child(i); + print_box(child); + } +} + + + +// Print all the vertices of a box. +template +void +test_vertex() +{ + deallog << "test_vertex" << std::endl; + + std::pair, Point> lower_upper_corners; + for (int d = 0; d < dim; ++d) + { + lower_upper_corners.first[d] = -(d + 1); + lower_upper_corners.second[d] = d + 1; + } + + const BoundingBox box(lower_upper_corners); + + for (unsigned int i = 0; i < GeometryInfo::vertices_per_cell; ++i) + { + deallog << "vertex " << i << " = " << box.vertex(i) << std::endl; + } +} + + + +// Print the results of the bounds function for all coordinate directions. +template +void +test_bounds() +{ + deallog << "test_bounds" << std::endl; + + std::pair, Point> lower_upper_corners; + for (int d = 0; d < dim; ++d) + { + lower_upper_corners.first[d] = -(d + 1); + lower_upper_corners.second[d] = d + 1; + } + + const BoundingBox box(lower_upper_corners); + + for (unsigned int i = 0; i < dim; ++i) + { + const BoundingBox<1> bounds = box.bounds(i); + deallog << "Bounds in direction " << i << std::endl; + print_box(bounds); + } +} + + + +// Test that we can call the create_unit_box function. +template +void +test_create_unit_box() +{ + deallog << "test_create_unit_box" << std::endl; + BoundingBox box = create_unit_box(); + print_box(box); +} + + + +template +void +run_test() +{ + deallog << "dim = " << dim << std::endl; + + test_cross_section(); + deallog << std::endl; + + test_center(); + deallog << std::endl; + + test_side_length(); + deallog << std::endl; + + test_child(); + deallog << std::endl; + + test_vertex(); + deallog << std::endl; + + test_bounds(); + deallog << std::endl; + + test_create_unit_box(); + deallog << std::endl; + + deallog << std::endl; +} + + + +int +main() +{ + initlog(); + + run_test<1>(); + run_test<2>(); + run_test<3>(); +} diff --git a/tests/base/bounding_box_6.output b/tests/base/bounding_box_6.output new file mode 100644 index 0000000000..c16b85b8da --- /dev/null +++ b/tests/base/bounding_box_6.output @@ -0,0 +1,124 @@ + +DEAL::dim = 1 +DEAL::test_cross_section +DEAL::orthogonal to 0 +DEAL:: +DEAL:: +DEAL::test_center +DEAL::center 0.500000 +DEAL:: +DEAL::test_side_length +DEAL::1.00000 +DEAL:: +DEAL::test_child +DEAL::child 0 +DEAL::[-1.00000, 0.00000] +DEAL::child 1 +DEAL::[0.00000, 1.00000] +DEAL:: +DEAL::test_vertex +DEAL::vertex 0 = -1.00000 +DEAL::vertex 1 = 1.00000 +DEAL:: +DEAL::test_bounds +DEAL::Bounds in direction 0 +DEAL::[-1.00000, 1.00000] +DEAL:: +DEAL::test_create_unit_box +DEAL::[0.00000, 1.00000] +DEAL:: +DEAL:: +DEAL::dim = 2 +DEAL::test_cross_section +DEAL::orthogonal to 0 +DEAL::[-2.00000, 2.00000] +DEAL::orthogonal to 1 +DEAL::[-1.00000, 1.00000] +DEAL:: +DEAL::test_center +DEAL::center 0.500000 0.500000 +DEAL:: +DEAL::test_side_length +DEAL::1.00000 2.00000 +DEAL:: +DEAL::test_child +DEAL::child 0 +DEAL::[-1.00000, 0.00000]x[-2.00000, 0.00000] +DEAL::child 1 +DEAL::[0.00000, 1.00000]x[-2.00000, 0.00000] +DEAL::child 2 +DEAL::[-1.00000, 0.00000]x[0.00000, 2.00000] +DEAL::child 3 +DEAL::[0.00000, 1.00000]x[0.00000, 2.00000] +DEAL:: +DEAL::test_vertex +DEAL::vertex 0 = -1.00000 -2.00000 +DEAL::vertex 1 = 1.00000 -2.00000 +DEAL::vertex 2 = -1.00000 2.00000 +DEAL::vertex 3 = 1.00000 2.00000 +DEAL:: +DEAL::test_bounds +DEAL::Bounds in direction 0 +DEAL::[-1.00000, 1.00000] +DEAL::Bounds in direction 1 +DEAL::[-2.00000, 2.00000] +DEAL:: +DEAL::test_create_unit_box +DEAL::[0.00000, 1.00000]x[0.00000, 1.00000] +DEAL:: +DEAL:: +DEAL::dim = 3 +DEAL::test_cross_section +DEAL::orthogonal to 0 +DEAL::[-2.00000, 2.00000]x[-3.00000, 3.00000] +DEAL::orthogonal to 1 +DEAL::[-3.00000, 3.00000]x[-1.00000, 1.00000] +DEAL::orthogonal to 2 +DEAL::[-1.00000, 1.00000]x[-2.00000, 2.00000] +DEAL:: +DEAL::test_center +DEAL::center 0.500000 0.500000 0.500000 +DEAL:: +DEAL::test_side_length +DEAL::1.00000 2.00000 3.00000 +DEAL:: +DEAL::test_child +DEAL::child 0 +DEAL::[-1.00000, 0.00000]x[-2.00000, 0.00000]x[-3.00000, 0.00000] +DEAL::child 1 +DEAL::[0.00000, 1.00000]x[-2.00000, 0.00000]x[-3.00000, 0.00000] +DEAL::child 2 +DEAL::[-1.00000, 0.00000]x[0.00000, 2.00000]x[-3.00000, 0.00000] +DEAL::child 3 +DEAL::[0.00000, 1.00000]x[0.00000, 2.00000]x[-3.00000, 0.00000] +DEAL::child 4 +DEAL::[-1.00000, 0.00000]x[-2.00000, 0.00000]x[0.00000, 3.00000] +DEAL::child 5 +DEAL::[0.00000, 1.00000]x[-2.00000, 0.00000]x[0.00000, 3.00000] +DEAL::child 6 +DEAL::[-1.00000, 0.00000]x[0.00000, 2.00000]x[0.00000, 3.00000] +DEAL::child 7 +DEAL::[0.00000, 1.00000]x[0.00000, 2.00000]x[0.00000, 3.00000] +DEAL:: +DEAL::test_vertex +DEAL::vertex 0 = -1.00000 -2.00000 -3.00000 +DEAL::vertex 1 = 1.00000 -2.00000 -3.00000 +DEAL::vertex 2 = -1.00000 2.00000 -3.00000 +DEAL::vertex 3 = 1.00000 2.00000 -3.00000 +DEAL::vertex 4 = -1.00000 -2.00000 3.00000 +DEAL::vertex 5 = 1.00000 -2.00000 3.00000 +DEAL::vertex 6 = -1.00000 2.00000 3.00000 +DEAL::vertex 7 = 1.00000 2.00000 3.00000 +DEAL:: +DEAL::test_bounds +DEAL::Bounds in direction 0 +DEAL::[-1.00000, 1.00000] +DEAL::Bounds in direction 1 +DEAL::[-2.00000, 2.00000] +DEAL::Bounds in direction 2 +DEAL::[-3.00000, 3.00000] +DEAL:: +DEAL::test_create_unit_box +DEAL::[0.00000, 1.00000]x[0.00000, 1.00000]x[0.00000, 1.00000] +DEAL:: +DEAL:: diff --git a/tests/non_matching/coordinate_to_one_dim_higher.cc b/tests/base/coordinate_to_one_dim_higher.cc similarity index 84% rename from tests/non_matching/coordinate_to_one_dim_higher.cc rename to tests/base/coordinate_to_one_dim_higher.cc index a5b9bc4ceb..26973e9608 100644 --- a/tests/non_matching/coordinate_to_one_dim_higher.cc +++ b/tests/base/coordinate_to_one_dim_higher.cc @@ -13,7 +13,7 @@ // // --------------------------------------------------------------------- -#include +#include #include "../tests.h" @@ -33,9 +33,7 @@ print_what_each_coordinate_maps_to() deallog << "locked coordinate = " << locked << std::endl; for (unsigned int i = 0; i < dim; ++i) { - deallog << i << " -> " - << NonMatching::internal::QuadratureGeneratorImplementation:: - coordinate_to_one_dim_higher(locked, i) + deallog << i << " -> " << coordinate_to_one_dim_higher(locked, i) << std::endl; } } diff --git a/tests/non_matching/coordinate_to_one_dim_higher.output b/tests/base/coordinate_to_one_dim_higher.output similarity index 100% rename from tests/non_matching/coordinate_to_one_dim_higher.output rename to tests/base/coordinate_to_one_dim_higher.output -- 2.39.5