From: Simon Sticko Date: Thu, 6 Feb 2020 16:26:04 +0000 (+0100) Subject: Add header quadrature_generator.h X-Git-Tag: v9.2.0-rc1~321^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fddd4db6f0ca434c16ef01cac8823ab176b0d708;p=dealii.git Add header quadrature_generator.h Add a header to be used for implementing a QuadratureGenerator class, which should generate immersed quadrature rules for the different regions of a cell intersected by the zero contour of a level set function. Add a hypercube class in a internal namespace QuadratureGeneratorImplementation, which is needed to implement the underlying algorithm. --- diff --git a/include/deal.II/non_matching/quadrature_generator.h b/include/deal.II/non_matching/quadrature_generator.h new file mode 100644 index 0000000000..ab6693f9d5 --- /dev/null +++ b/include/deal.II/non_matching/quadrature_generator.h @@ -0,0 +1,231 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii_non_matching_quadrature_generator +#define dealii_non_matching_quadrature_generator + +#include + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace NonMatching +{ + namespace internal + { + namespace QuadratureGeneratorImplementation + { + /** + * Returns the center $(0.5, ..., 0.5)$ of the unit hypercube: + * $[0,1]^{dim}$. + */ + template + Point + unit_hypercube_center(); + + + /** + * Describes a dim-dimensional hypercube aligned with the coordinate + * axes. That is, a region + * + * @f[ + * [x_0-L/2, x_0+L/2] \times ... \times [x_{dim-1}-L/2, x_{dim-1}+L/2] + * @f] + * + * where L is the side length and $(x_0, ..., x_{dim-1})$ are the + * coordinates of the center of the hypercube. + * + * The purpose of this class is to help us keep track of where we are on + * the reference cell when generating immersed quadrature rules. There + * are two main things we want to be able to do. + * + * 1. Split a hypercube into its $2^{dim}$ children. + * 2. Take a cross section orthogonal to a given direction and get this + * as a Hypercube. + * + * These are needed because the algorithm is allowed to recurse both + * over children and over dimensions. + * + * In 3D, the 2 coordinates of the cross section of Hypercube<3> can be + * ordered in 2 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. + */ + template + class Hypercube + { + public: + /** + * Constructor, takes the center of the hypercube and its side length. + * Passing no inarguments creates a hypercube between 0 and 1 in all + * directions: $[0, 1]^{dim}$. + */ + Hypercube( + const Point ¢er_coordinate = unit_hypercube_center(), + const double length_of_side = 1); + + /** + * Return the lower bound of the hypercube in @p direction. + */ + double + lower_bound(const unsigned int direction) const; + + /** + * Return the upper bound of the hypercube in @p direction. + */ + double + upper_bound(const unsigned int direction) const; + + /** + * Return the bounds of the hypercube in @p direction, as a one-dimensional + * hypercube. + */ + Hypercube<1> + bounds(const unsigned int direction) const; + + /** + * Returns the cross section of the hypercube orthogonal to @p direction. + * This is a hypercube in one dimension lower. + */ + Hypercube + cross_section(const unsigned int direction) const; + + /** + * Returns the point in the center of the hypercube. + */ + const Point & + center() const; + + /** + * Returns the length of the side of the hypercube. + */ + double + side_length() const; + + /** + * Returns the volume of the hypercube. + */ + double + volume() const; + + /** + * Returns the indexth vertex of the hypercube. + */ + Point + vertex(const unsigned int index) const; + + /** + * Returns the indexth child of the hypercube. Child is meant in the + * same way as for a cell. + */ + Hypercube + child(const unsigned int index) const; + + private: + /** + * Coordinate in the center of the hypercube. + */ + const Point center_coordinate; + + /** + * Side length of the hypercube. + */ + const double length_of_side; + }; + + + /** + * 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); + } + + + + /* -------------- declaration of explicit specializations ------------- */ + + template <> + Hypercube<-1> + Hypercube<0>::cross_section(const unsigned int) const; + + template <> + Point<0> + Hypercube<0>::vertex(const unsigned int) const; + + template <> + Hypercube<0> + Hypercube<0>::child(const unsigned int) const; + + } // namespace QuadratureGeneratorImplementation + } // namespace internal + +} // namespace NonMatching +DEAL_II_NAMESPACE_CLOSE + +#endif /* dealii_non_matching_quadrature_generator */ diff --git a/source/non_matching/CMakeLists.txt b/source/non_matching/CMakeLists.txt index 6d9ccbb32a..eecc07cec5 100644 --- a/source/non_matching/CMakeLists.txt +++ b/source/non_matching/CMakeLists.txt @@ -18,6 +18,7 @@ INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) SET(_src coupling.cc immersed_surface_quadrature.cc + quadrature_generator.cc ) SET(_inst diff --git a/source/non_matching/quadrature_generator.cc b/source/non_matching/quadrature_generator.cc new file mode 100644 index 0000000000..905754aec7 --- /dev/null +++ b/source/non_matching/quadrature_generator.cc @@ -0,0 +1,194 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 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. +// +// --------------------------------------------------------------------- + +#include +#include + +#include + + +DEAL_II_NAMESPACE_OPEN + +namespace NonMatching +{ + namespace internal + { + namespace QuadratureGeneratorImplementation + { + template + Point + unit_hypercube_center() + { + Point center; + for (unsigned int i = 0; i < dim; ++i) + center[i] = .5; + + return center; + } + + + + template + Hypercube::Hypercube(const Point ¢er_coordinate, + const double length_of_side) + : center_coordinate(center_coordinate) + , length_of_side(length_of_side) + { + Assert(0 < dim, ExcImpossibleInDim(dim)); + } + + + + template + double + Hypercube::lower_bound(const unsigned int direction) const + { + AssertIndexRange(direction, dim); + + return center_coordinate[direction] - length_of_side / 2; + } + + + + template + double + Hypercube::upper_bound(const unsigned int direction) const + { + AssertIndexRange(direction, dim); + + return center_coordinate[direction] + length_of_side / 2; + } + + + + template + Hypercube<1> + Hypercube::bounds(const unsigned int direction) const + { + AssertIndexRange(direction, dim); + + const Point<1> center(center_coordinate[direction]); + return Hypercube<1>(center, side_length()); + } + + + + template + const Point & + Hypercube::center() const + { + return center_coordinate; + } + + + + template + double + Hypercube::side_length() const + { + return length_of_side; + } + + + + template + double + Hypercube::volume() const + { + return std::pow(length_of_side, dim); + } + + + + template + Hypercube + Hypercube::cross_section(const unsigned int direction) const + { + AssertIndexRange(direction, dim); + + Point center_of_cross_section; + for (unsigned int d = 0; d < dim - 1; ++d) + { + const int index_to_write_from = + coordinate_to_one_dim_higher(direction, d); + center_of_cross_section[d] = center_coordinate[index_to_write_from]; + } + + return Hypercube(center_of_cross_section, length_of_side); + } + + + + template + Point + Hypercube::vertex(const unsigned int index) const + { + AssertIndexRange(index, GeometryInfo::vertices_per_cell); + + Point bottom_corner; + for (unsigned int i = 0; i < dim; ++i) + bottom_corner[i] = center_coordinate[i] - length_of_side / 2; + + const Point point = + bottom_corner + + length_of_side * GeometryInfo::unit_cell_vertex(index); + + return point; + } + + + + template + Hypercube + Hypercube::child(const unsigned int index) const + { + AssertIndexRange(index, GeometryInfo::max_children_per_cell); + + const double child_side_length = length_of_side / 2; + + const Point unit_cell_vertex = + GeometryInfo::unit_cell_vertex(index); + + Point child_center; + for (unsigned int i = 0; i < dim; ++i) + child_center[i] = center_coordinate[i] + + child_side_length * (unit_cell_vertex[i] - .5); + + return Hypercube(child_center, child_side_length); + } + + + + template Point<0> + unit_hypercube_center<0>(); + template Point<1> + unit_hypercube_center<1>(); + template Point<2> + unit_hypercube_center<2>(); + template Point<3> + unit_hypercube_center<3>(); + + template class Hypercube<0>; + template class Hypercube<1>; + template class Hypercube<2>; + template class Hypercube<3>; + + } // namespace QuadratureGeneratorImplementation + } // namespace internal + + + +} // namespace NonMatching +DEAL_II_NAMESPACE_CLOSE diff --git a/tests/non_matching/coordinate_to_one_dim_higher.cc b/tests/non_matching/coordinate_to_one_dim_higher.cc new file mode 100644 index 0000000000..a5b9bc4ceb --- /dev/null +++ b/tests/non_matching/coordinate_to_one_dim_higher.cc @@ -0,0 +1,56 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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. +// +// --------------------------------------------------------------------- + +#include + +#include "../tests.h" + +/* + * Lock all possible coordinates in dim + 1. For each locked coordinate, print + * how the coordinates in dim dimensions map to the coordinates in dim + 1 + * dimensions. + */ +template +void +print_what_each_coordinate_maps_to() +{ + deallog << "dim = " << dim << std::endl; + + for (unsigned int locked = 0; locked < dim + 1; ++locked) + { + 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) + << std::endl; + } + } +} + + + +int +main() +{ + initlog(); + + print_what_each_coordinate_maps_to<1>(); + deallog << std::endl; + print_what_each_coordinate_maps_to<2>(); + + return 0; +} diff --git a/tests/non_matching/coordinate_to_one_dim_higher.output b/tests/non_matching/coordinate_to_one_dim_higher.output new file mode 100644 index 0000000000..cbb82ceaa5 --- /dev/null +++ b/tests/non_matching/coordinate_to_one_dim_higher.output @@ -0,0 +1,17 @@ + +DEAL::dim = 1 +DEAL::locked coordinate = 0 +DEAL::0 -> 1 +DEAL::locked coordinate = 1 +DEAL::0 -> 0 +DEAL:: +DEAL::dim = 2 +DEAL::locked coordinate = 0 +DEAL::0 -> 1 +DEAL::1 -> 2 +DEAL::locked coordinate = 1 +DEAL::0 -> 2 +DEAL::1 -> 0 +DEAL::locked coordinate = 2 +DEAL::0 -> 0 +DEAL::1 -> 1 diff --git a/tests/non_matching/hypercube.cc b/tests/non_matching/hypercube.cc new file mode 100644 index 0000000000..16585651f1 --- /dev/null +++ b/tests/non_matching/hypercube.cc @@ -0,0 +1,217 @@ +// --------------------------------------------------------------------- +// +// 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 all member functions of the HyperCube class. Each function is +// tested in the function named test_{member_function_name}. + +#include + +#include + +#include "../tests.h" + + +using namespace dealii; +using NonMatching::internal::QuadratureGeneratorImplementation::Hypercube; + +// Print the center and the side length of the incoming hypercube, which are +// its degrees of freedom. +template +void +print_cube(const Hypercube &hyper_cube) +{ + deallog << "center = " << hyper_cube.center() << ", "; + deallog << "side = " << hyper_cube.side_length() << std::endl; +} + + + +// Construct the cross section orthogonal to all axes and print them. +template +void +test_cross_section() +{ + deallog << "test_cross_section" << std::endl; + + const double width = 1; + Point center; + for (unsigned int d = 0; d < dim; ++d) + center[d] = d; + + const Hypercube hypercube(center, width); + + for (unsigned int d = 0; d < dim; ++d) + { + deallog << "orthogonal to " << d << std::endl; + const Hypercube down = hypercube.cross_section(d); + print_cube(down); + } +} + + + +// 1D-specialization of the above test. We can't take a cross section in 1D, so +// we do nothing. +template <> +void +test_cross_section<1>() +{} + + + +// Construct all possible children of a hypercube and print them. +template +void +test_child() +{ + deallog << "test_child" << std::endl; + + const double width = 4; + Point center; + + const Hypercube hypercube(center, width); + + // Take the cross section in all directions. + for (unsigned int i = 0; i < GeometryInfo::max_children_per_cell; ++i) + { + deallog << "child " << i << std::endl; + const Hypercube child = hypercube.child(i); + print_cube(child); + } +} + + + +// Comptue and print the volume of a hypercube. +template +void +test_volume() +{ + deallog << "test_volume" << std::endl; + + const double width = 2; + Point center; + + const Hypercube hypercube(center, width); + deallog << "volume = " << hypercube.volume() << std::endl; +} + + + +// Print the result of the lower_bound and upper_bound functions in all +// coordinate directions. +template +void +test_lower_upper_bound() +{ + deallog << "test_lower_upper_bound" << std::endl; + + const double width = 2; + Point center; + for (unsigned int d = 0; d < dim; ++d) + center[d] = d + 1; + + const Hypercube hypercube(center, width); + for (int d = 0; d < dim; ++d) + { + deallog << "[" << hypercube.lower_bound(d) << ", " + << hypercube.upper_bound(d) << "]"; + if (dim - 1 != d) + { + deallog << "x"; + } + } +} + + + +// Extract and print all the vertices of a hypercube. +template +void +test_vertex() +{ + deallog << "test_vertex" << std::endl; + + const Hypercube hypercube; + for (unsigned int i = 0; i < GeometryInfo::vertices_per_cell; ++i) + { + deallog << "vertex " << i << " = " << hypercube.vertex(i) << std::endl; + } +} + + + +// Print the results of the bounds function for all coordinate directions. +template +void +test_bounds() +{ + deallog << "test_bounds" << std::endl; + + const double width = 2; + Point center; + for (unsigned int d = 0; d < dim; ++d) + center[d] = d + 1; + + const Hypercube hypercube(center, width); + + for (unsigned int i = 0; i < dim; ++i) + { + const Hypercube<1> bounds = hypercube.bounds(i); + deallog << "Bounds in direction " << i << std::endl; + print_cube(bounds); + } +} + + + +template +void +run_test() +{ + deallog << "dim = " << dim << std::endl; + + test_cross_section(); + deallog << std::endl; + + test_child(); + deallog << std::endl; + + test_volume(); + deallog << std::endl; + + test_lower_upper_bound(); + deallog << std::endl; + + test_vertex(); + deallog << std::endl; + + test_bounds(); + deallog << std::endl; + + deallog << std::endl; +} + + + +int +main() +{ + initlog(); + + run_test<1>(); + run_test<2>(); + run_test<3>(); +} diff --git a/tests/non_matching/hypercube.output b/tests/non_matching/hypercube.output new file mode 100644 index 0000000000..abcd206164 --- /dev/null +++ b/tests/non_matching/hypercube.output @@ -0,0 +1,109 @@ + +DEAL::dim = 1 +DEAL:: +DEAL::test_child +DEAL::child 0 +DEAL::center = -1.00000, side = 2.00000 +DEAL::child 1 +DEAL::center = 1.00000, side = 2.00000 +DEAL:: +DEAL::test_volume +DEAL::volume = 2.00000 +DEAL:: +DEAL::test_lower_upper_bound +DEAL::[0.00000, 2.00000] +DEAL::test_vertex +DEAL::vertex 0 = 0.00000 +DEAL::vertex 1 = 1.00000 +DEAL:: +DEAL::test_bounds +DEAL::Bounds in direction 0 +DEAL::center = 1.00000, side = 2.00000 +DEAL:: +DEAL:: +DEAL::dim = 2 +DEAL::test_cross_section +DEAL::orthogonal to 0 +DEAL::center = 1.00000, side = 1.00000 +DEAL::orthogonal to 1 +DEAL::center = 0.00000, side = 1.00000 +DEAL:: +DEAL::test_child +DEAL::child 0 +DEAL::center = -1.00000 -1.00000, side = 2.00000 +DEAL::child 1 +DEAL::center = 1.00000 -1.00000, side = 2.00000 +DEAL::child 2 +DEAL::center = -1.00000 1.00000, side = 2.00000 +DEAL::child 3 +DEAL::center = 1.00000 1.00000, side = 2.00000 +DEAL:: +DEAL::test_volume +DEAL::volume = 4.00000 +DEAL:: +DEAL::test_lower_upper_bound +DEAL::[0.00000, 2.00000]x[1.00000, 3.00000] +DEAL::test_vertex +DEAL::vertex 0 = 0.00000 0.00000 +DEAL::vertex 1 = 1.00000 0.00000 +DEAL::vertex 2 = 0.00000 1.00000 +DEAL::vertex 3 = 1.00000 1.00000 +DEAL:: +DEAL::test_bounds +DEAL::Bounds in direction 0 +DEAL::center = 1.00000, side = 2.00000 +DEAL::Bounds in direction 1 +DEAL::center = 2.00000, side = 2.00000 +DEAL:: +DEAL:: +DEAL::dim = 3 +DEAL::test_cross_section +DEAL::orthogonal to 0 +DEAL::center = 1.00000 2.00000, side = 1.00000 +DEAL::orthogonal to 1 +DEAL::center = 2.00000 0.00000, side = 1.00000 +DEAL::orthogonal to 2 +DEAL::center = 0.00000 1.00000, side = 1.00000 +DEAL:: +DEAL::test_child +DEAL::child 0 +DEAL::center = -1.00000 -1.00000 -1.00000, side = 2.00000 +DEAL::child 1 +DEAL::center = 1.00000 -1.00000 -1.00000, side = 2.00000 +DEAL::child 2 +DEAL::center = -1.00000 1.00000 -1.00000, side = 2.00000 +DEAL::child 3 +DEAL::center = 1.00000 1.00000 -1.00000, side = 2.00000 +DEAL::child 4 +DEAL::center = -1.00000 -1.00000 1.00000, side = 2.00000 +DEAL::child 5 +DEAL::center = 1.00000 -1.00000 1.00000, side = 2.00000 +DEAL::child 6 +DEAL::center = -1.00000 1.00000 1.00000, side = 2.00000 +DEAL::child 7 +DEAL::center = 1.00000 1.00000 1.00000, side = 2.00000 +DEAL:: +DEAL::test_volume +DEAL::volume = 8.00000 +DEAL:: +DEAL::test_lower_upper_bound +DEAL::[0.00000, 2.00000]x[1.00000, 3.00000]x[2.00000, 4.00000] +DEAL::test_vertex +DEAL::vertex 0 = 0.00000 0.00000 0.00000 +DEAL::vertex 1 = 1.00000 0.00000 0.00000 +DEAL::vertex 2 = 0.00000 1.00000 0.00000 +DEAL::vertex 3 = 1.00000 1.00000 0.00000 +DEAL::vertex 4 = 0.00000 0.00000 1.00000 +DEAL::vertex 5 = 1.00000 0.00000 1.00000 +DEAL::vertex 6 = 0.00000 1.00000 1.00000 +DEAL::vertex 7 = 1.00000 1.00000 1.00000 +DEAL:: +DEAL::test_bounds +DEAL::Bounds in direction 0 +DEAL::center = 1.00000, side = 2.00000 +DEAL::Bounds in direction 1 +DEAL::center = 2.00000, side = 2.00000 +DEAL::Bounds in direction 2 +DEAL::center = 3.00000, side = 2.00000 +DEAL:: +DEAL::