From 7efece5350563b44fb9a1034a1e6e90ab5466ae2 Mon Sep 17 00:00:00 2001 From: Simon Sticko Date: Mon, 14 Jun 2021 20:26:11 +0200 Subject: [PATCH] Add class generating immersed quadrature rules over a box face Add a class NonMatching::FaceQuadratureGenerator, that in the same way as QuadratureGenerator generates immersed quadrature rules, but over a given face of a BoundingBox. --- .../non_matching/quadrature_generator.h | 112 +++++++++ source/non_matching/quadrature_generator.cc | 96 ++++++++ .../non_matching/quadrature_generator.inst.in | 4 + .../non_matching/face_quadrature_generator.cc | 115 +++++++++ .../face_quadrature_generator.output | 225 ++++++++++++++++++ tests/non_matching/quadrature_printing.h | 11 +- 6 files changed, 558 insertions(+), 5 deletions(-) create mode 100644 tests/non_matching/face_quadrature_generator.cc create mode 100644 tests/non_matching/face_quadrature_generator.output diff --git a/include/deal.II/non_matching/quadrature_generator.h b/include/deal.II/non_matching/quadrature_generator.h index 6d105f2ec7..2fa8214601 100644 --- a/include/deal.II/non_matching/quadrature_generator.h +++ b/include/deal.II/non_matching/quadrature_generator.h @@ -250,6 +250,118 @@ namespace NonMatching q_generator; }; + + /** + * This class creates immersed quadrature rules over a face, $F$, of a + * BoundingBox, when the domain is described by a level set function, $\psi$. + * + * In the same way as in the QuadratureGenerator class, this class generates + * quadrature rules to integrate over 3 different regions of the face, + * $F \subset \mathbb{R}^{dim}$: + * @f[ + * N = \{x \in F : \psi(x) < 0 \}, \\ + * P = \{x \in F : \psi(x) > 0 \}, \\ + * S = \{x \in F : \psi(x) = 0 \}, + * @f] + * which are again referred to as the "inside", $N$, "outside", $P$, + * and "surface" region, $S$. These type of quadrature rules are in general + * needed in immersed discontinuous Galerkin methods. + * + * Under the hood, this class uses the QuadratureGenerator class to build + * these rules. This is done by restricting the dim-dimensional level set + * function to the face, thus creating a (dim-1)-dimensional level set + * function, $\phi$. It then creates the (dim-1)-dimensional quadratures by + * calling QuadratureGenerator with $\phi$. This means that what holds for the + * QuadratureGenerator class in general also holds for this class. In + * particular, if the 1D-quadrature that is used as base contains $n$ points, + * the number of points will be proportional to $n^{dim-1}$ in the in the + * inside/outside quadratures and to $n^{dim-2}$ in the surface quadrature. + */ + template + class FaceQuadratureGenerator + { + public: + using AdditionalData = AdditionalQGeneratorData; + + /** + * Constructor. Each Quadrature<1> in @p quadratures1D can be chosen as base + * for generating the immersed quadrature rules. + * + * @note It is important that each 1D-quadrature rule in the + * hp::QCollection does not contain the points 0 and 1. + */ + FaceQuadratureGenerator( + const hp::QCollection<1> &quadratures1D, + const AdditionalData & additional_data = AdditionalData()); + + /** + * Construct immersed quadratures rules for the incoming level set + * function on a given face of the BoundingBox. + * + * To get the constructed quadratures, use the functions + * get_inside_quadrature(), + * get_outside_quadrature(), + * get_surface_quadrature(). + * + * @note Both value, gradient and hessian need to be implemented on the + * incoming function. + */ + void + generate(const Function & level_set, + const BoundingBox &box, + const unsigned int face_index); + + /** + * Return the quadrature rule for the region + * $\{x \in F : \psi(x) < 0 \}$ + * created in the previous call to generate(). + * Here, $F$ is the face of the BoundingBox passed to generate(). + */ + const Quadrature & + get_inside_quadrature() const; + + /** + * Return the quadrature rule for the region + * $\{x \in F : \psi(x) > 0 \}$ + * created in the previous call to generate(). + * Here, $F$ is the face of the BoundingBox passed to generate(). + */ + const Quadrature & + get_outside_quadrature() const; + + /** + * Return the quadrature rule for the region + * $\{x \in F : \psi(x) = 0 \}$ + * created in the previous call to generate(). + * Here, $F$ is the face of the BoundingBox passed to generate(). + * + * @note The normal at the quadrature points will be parallel to $\nabla \psi$. + */ + const ImmersedSurfaceQuadrature & + get_surface_quadrature() const; + + /** + * Set which 1D-quadrature in the collection passed to the constructor + * should be used to create the immersed quadratures. + */ + void + set_1D_quadrature(const unsigned int q_index); + + private: + /** + * Lower-dimensional quadrature generator used to build the quadratures over + * the face. + */ + QuadratureGenerator quadrature_generator; + + /** + * The same surface quadrature as created by the quadrature_generator, + * but having dim-dimensional normals. + */ + ImmersedSurfaceQuadrature surface_quadrature; + }; + + namespace internal { namespace QuadratureGeneratorImplementation diff --git a/source/non_matching/quadrature_generator.cc b/source/non_matching/quadrature_generator.cc index ea0245e069..4680dd8536 100644 --- a/source/non_matching/quadrature_generator.cc +++ b/source/non_matching/quadrature_generator.cc @@ -1353,6 +1353,102 @@ namespace NonMatching q_generator.set_1D_quadrature(q_index); } + + + template + FaceQuadratureGenerator::FaceQuadratureGenerator( + const hp::QCollection<1> &quadratures1D, + const AdditionalData & additional_data) + : quadrature_generator(quadratures1D, additional_data) + {} + + + + template + void + FaceQuadratureGenerator::generate(const Function & level_set, + const BoundingBox &box, + const unsigned int face_index) + { + AssertIndexRange(face_index, GeometryInfo::faces_per_cell); + + // We restrict the level set function to the face, by locking the coordinate + // that is constant over the face. This will be the same as the direction of + // the face normal. + const unsigned int face_normal_direction = + GeometryInfo::unit_normal_direction[face_index]; + + const Point vertex0 = + box.vertex(GeometryInfo::face_to_cell_vertices(face_index, 0)); + const double coordinate_value = vertex0(face_normal_direction); + + const Functions::CoordinateRestriction face_restriction( + level_set, face_normal_direction, coordinate_value); + + // Reuse the lower dimensional QuadratureGenerator on the face. + const BoundingBox cross_section = + box.cross_section(face_normal_direction); + quadrature_generator.generate(face_restriction, cross_section); + + // We need the dim-dimensional normals of the zero-contour. + // Recompute these. + const ImmersedSurfaceQuadrature + &surface_quadrature_wrong_normal = + quadrature_generator.get_surface_quadrature(); + + std::vector> normals; + normals.reserve(surface_quadrature_wrong_normal.size()); + for (unsigned int i = 0; i < surface_quadrature_wrong_normal.size(); i++) + { + const Point point = dealii::internal::create_higher_dim_point( + surface_quadrature_wrong_normal.point(i), + face_normal_direction, + coordinate_value); + + Tensor<1, dim> normal = level_set.gradient(point); + normal /= normal.norm(); + normals.push_back(normal); + } + surface_quadrature = ImmersedSurfaceQuadrature( + surface_quadrature_wrong_normal.get_points(), + surface_quadrature_wrong_normal.get_weights(), + normals); + } + + + + template + void + FaceQuadratureGenerator::set_1D_quadrature(const unsigned int q_index) + { + quadrature_generator.set_1D_quadrature(q_index); + } + + + + template + const Quadrature & + FaceQuadratureGenerator::get_inside_quadrature() const + { + return quadrature_generator.get_inside_quadrature(); + } + + + template + const Quadrature & + FaceQuadratureGenerator::get_outside_quadrature() const + { + return quadrature_generator.get_outside_quadrature(); + } + + + + template + const ImmersedSurfaceQuadrature & + FaceQuadratureGenerator::get_surface_quadrature() const + { + return surface_quadrature; + } } // namespace NonMatching #include "quadrature_generator.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/source/non_matching/quadrature_generator.inst.in b/source/non_matching/quadrature_generator.inst.in index 736ff97ae9..3b6d03c08a 100644 --- a/source/non_matching/quadrature_generator.inst.in +++ b/source/non_matching/quadrature_generator.inst.in @@ -20,6 +20,10 @@ for (deal_II_dimension : DIMENSIONS) \{ template class QuadratureGenerator; +#if 1 < deal_II_dimension + template class FaceQuadratureGenerator; +#endif + namespace internal \{ namespace QuadratureGeneratorImplementation diff --git a/tests/non_matching/face_quadrature_generator.cc b/tests/non_matching/face_quadrature_generator.cc new file mode 100644 index 0000000000..c9e8824441 --- /dev/null +++ b/tests/non_matching/face_quadrature_generator.cc @@ -0,0 +1,115 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 - 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 FaceQuadratureGenerator class, by setting up simple cuts over a face + * of the unit box and writing the generated quadrature rules to the output + * file. + * + * Each function beginning with "test_" sets up a level set function and then + * calls the function create_and_print_quadratures() to generate the + * quadratures. + */ + +#include +#include +#include + +#include + +#include + +#include "../tests.h" + +#include "quadrature_printing.h" + + +/* + * Create immersed quadrature rules over the face_index:th unit-box-face for the + * incoming level set function. Print the constructed quadrature rules to + * deallog. + */ +template +void +create_and_print_quadratures(const Function &levelset, + const unsigned int face_index) +{ + const hp::QCollection<1> quadratures1D(QGauss<1>(2)); + const BoundingBox box = create_unit_bounding_box(); + + NonMatching::FaceQuadratureGenerator quadrature_generator(quadratures1D); + quadrature_generator.generate(levelset, box, face_index); + + deallog << "Inside quadrature" << std::endl; + print_quadrature(quadrature_generator.get_inside_quadrature()); + deallog << "Outside quadrature" << std::endl; + print_quadrature(quadrature_generator.get_outside_quadrature()); + deallog << "Surface quadrature" << std::endl; + print_surface_quadrature(quadrature_generator.get_surface_quadrature()); +} + + + +/* + * For each coordinate direction, set up a level set function corresponding to a + * plane cutting through the unit cell, having the normal aligned with the + * coordinate direction. Create and print the constructed quadratures for the + * two faces that are intersected. We expect that the constructed quadrature has + * the same number of points in the inside/outside region and that they are + * tensor products. + */ +template +void +test_plane_cuts_through_center() +{ + deallog << "test_vertical_cuts_through_center" << std::endl; + deallog << "dim=" << dim << std::endl; + + Point center; + for (int d = 0; d < dim; ++d) + center(d) = .5; + + // For each coordinate direction set up a plane throught the center. + for (int plane_direction = 0; plane_direction < dim; ++plane_direction) + { + const Tensor<1, dim> normal = Point::unit_vector(plane_direction); + const Functions::LevelSet::Plane levelset(center, normal); + + // Test all faces that are intersected by the plane. + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; f++) + { + const int normal_direction = + GeometryInfo::unit_normal_direction[f]; + if (plane_direction != normal_direction) + { + deallog << "plane direction = " << plane_direction << ", "; + deallog << "face = " << f << std::endl; + create_and_print_quadratures(levelset, f); + deallog << std::endl; + } + } + } +} + + + +int +main() +{ + initlog(); + + test_plane_cuts_through_center<2>(); + test_plane_cuts_through_center<3>(); +} diff --git a/tests/non_matching/face_quadrature_generator.output b/tests/non_matching/face_quadrature_generator.output new file mode 100644 index 0000000000..c3a17d9259 --- /dev/null +++ b/tests/non_matching/face_quadrature_generator.output @@ -0,0 +1,225 @@ + +DEAL::test_vertical_cuts_through_center +DEAL::dim=2 +DEAL::plane direction = 0, face = 2 +DEAL::Inside quadrature +DEAL::0.105662, 0.250000 +DEAL::0.394338, 0.250000 +DEAL::Outside quadrature +DEAL::0.605662, 0.250000 +DEAL::0.894338, 0.250000 +DEAL::Surface quadrature +DEAL::0.500000, 1.00000, 1.00000, 0.00000 +DEAL:: +DEAL::plane direction = 0, face = 3 +DEAL::Inside quadrature +DEAL::0.105662, 0.250000 +DEAL::0.394338, 0.250000 +DEAL::Outside quadrature +DEAL::0.605662, 0.250000 +DEAL::0.894338, 0.250000 +DEAL::Surface quadrature +DEAL::0.500000, 1.00000, 1.00000, 0.00000 +DEAL:: +DEAL::plane direction = 1, face = 0 +DEAL::Inside quadrature +DEAL::0.105662, 0.250000 +DEAL::0.394338, 0.250000 +DEAL::Outside quadrature +DEAL::0.605662, 0.250000 +DEAL::0.894338, 0.250000 +DEAL::Surface quadrature +DEAL::0.500000, 1.00000, 0.00000, 1.00000 +DEAL:: +DEAL::plane direction = 1, face = 1 +DEAL::Inside quadrature +DEAL::0.105662, 0.250000 +DEAL::0.394338, 0.250000 +DEAL::Outside quadrature +DEAL::0.605662, 0.250000 +DEAL::0.894338, 0.250000 +DEAL::Surface quadrature +DEAL::0.500000, 1.00000, 0.00000, 1.00000 +DEAL:: +DEAL::test_vertical_cuts_through_center +DEAL::dim=3 +DEAL::plane direction = 0, face = 2 +DEAL::Inside quadrature +DEAL::0.211325, 0.105662, 0.125000 +DEAL::0.211325, 0.394338, 0.125000 +DEAL::0.788675, 0.105662, 0.125000 +DEAL::0.788675, 0.394338, 0.125000 +DEAL::Outside quadrature +DEAL::0.211325, 0.605662, 0.125000 +DEAL::0.211325, 0.894338, 0.125000 +DEAL::0.788675, 0.605662, 0.125000 +DEAL::0.788675, 0.894338, 0.125000 +DEAL::Surface quadrature +DEAL::0.211325, 0.500000, 0.500000, 1.00000, 0.00000, 0.00000 +DEAL::0.788675, 0.500000, 0.500000, 1.00000, 0.00000, 0.00000 +DEAL:: +DEAL::plane direction = 0, face = 3 +DEAL::Inside quadrature +DEAL::0.211325, 0.105662, 0.125000 +DEAL::0.211325, 0.394338, 0.125000 +DEAL::0.788675, 0.105662, 0.125000 +DEAL::0.788675, 0.394338, 0.125000 +DEAL::Outside quadrature +DEAL::0.211325, 0.605662, 0.125000 +DEAL::0.211325, 0.894338, 0.125000 +DEAL::0.788675, 0.605662, 0.125000 +DEAL::0.788675, 0.894338, 0.125000 +DEAL::Surface quadrature +DEAL::0.211325, 0.500000, 0.500000, 1.00000, 0.00000, 0.00000 +DEAL::0.788675, 0.500000, 0.500000, 1.00000, 0.00000, 0.00000 +DEAL:: +DEAL::plane direction = 0, face = 4 +DEAL::Inside quadrature +DEAL::0.105662, 0.211325, 0.125000 +DEAL::0.394338, 0.211325, 0.125000 +DEAL::0.105662, 0.788675, 0.125000 +DEAL::0.394338, 0.788675, 0.125000 +DEAL::Outside quadrature +DEAL::0.605662, 0.211325, 0.125000 +DEAL::0.894338, 0.211325, 0.125000 +DEAL::0.605662, 0.788675, 0.125000 +DEAL::0.894338, 0.788675, 0.125000 +DEAL::Surface quadrature +DEAL::0.500000, 0.211325, 0.500000, 1.00000, 0.00000, 0.00000 +DEAL::0.500000, 0.788675, 0.500000, 1.00000, 0.00000, 0.00000 +DEAL:: +DEAL::plane direction = 0, face = 5 +DEAL::Inside quadrature +DEAL::0.105662, 0.211325, 0.125000 +DEAL::0.394338, 0.211325, 0.125000 +DEAL::0.105662, 0.788675, 0.125000 +DEAL::0.394338, 0.788675, 0.125000 +DEAL::Outside quadrature +DEAL::0.605662, 0.211325, 0.125000 +DEAL::0.894338, 0.211325, 0.125000 +DEAL::0.605662, 0.788675, 0.125000 +DEAL::0.894338, 0.788675, 0.125000 +DEAL::Surface quadrature +DEAL::0.500000, 0.211325, 0.500000, 1.00000, 0.00000, 0.00000 +DEAL::0.500000, 0.788675, 0.500000, 1.00000, 0.00000, 0.00000 +DEAL:: +DEAL::plane direction = 1, face = 0 +DEAL::Inside quadrature +DEAL::0.105662, 0.211325, 0.125000 +DEAL::0.394338, 0.211325, 0.125000 +DEAL::0.105662, 0.788675, 0.125000 +DEAL::0.394338, 0.788675, 0.125000 +DEAL::Outside quadrature +DEAL::0.605662, 0.211325, 0.125000 +DEAL::0.894338, 0.211325, 0.125000 +DEAL::0.605662, 0.788675, 0.125000 +DEAL::0.894338, 0.788675, 0.125000 +DEAL::Surface quadrature +DEAL::0.500000, 0.211325, 0.500000, 0.00000, 1.00000, 0.00000 +DEAL::0.500000, 0.788675, 0.500000, 0.00000, 1.00000, 0.00000 +DEAL:: +DEAL::plane direction = 1, face = 1 +DEAL::Inside quadrature +DEAL::0.105662, 0.211325, 0.125000 +DEAL::0.394338, 0.211325, 0.125000 +DEAL::0.105662, 0.788675, 0.125000 +DEAL::0.394338, 0.788675, 0.125000 +DEAL::Outside quadrature +DEAL::0.605662, 0.211325, 0.125000 +DEAL::0.894338, 0.211325, 0.125000 +DEAL::0.605662, 0.788675, 0.125000 +DEAL::0.894338, 0.788675, 0.125000 +DEAL::Surface quadrature +DEAL::0.500000, 0.211325, 0.500000, 0.00000, 1.00000, 0.00000 +DEAL::0.500000, 0.788675, 0.500000, 0.00000, 1.00000, 0.00000 +DEAL:: +DEAL::plane direction = 1, face = 4 +DEAL::Inside quadrature +DEAL::0.211325, 0.105662, 0.125000 +DEAL::0.211325, 0.394338, 0.125000 +DEAL::0.788675, 0.105662, 0.125000 +DEAL::0.788675, 0.394338, 0.125000 +DEAL::Outside quadrature +DEAL::0.211325, 0.605662, 0.125000 +DEAL::0.211325, 0.894338, 0.125000 +DEAL::0.788675, 0.605662, 0.125000 +DEAL::0.788675, 0.894338, 0.125000 +DEAL::Surface quadrature +DEAL::0.211325, 0.500000, 0.500000, 0.00000, 1.00000, 0.00000 +DEAL::0.788675, 0.500000, 0.500000, 0.00000, 1.00000, 0.00000 +DEAL:: +DEAL::plane direction = 1, face = 5 +DEAL::Inside quadrature +DEAL::0.211325, 0.105662, 0.125000 +DEAL::0.211325, 0.394338, 0.125000 +DEAL::0.788675, 0.105662, 0.125000 +DEAL::0.788675, 0.394338, 0.125000 +DEAL::Outside quadrature +DEAL::0.211325, 0.605662, 0.125000 +DEAL::0.211325, 0.894338, 0.125000 +DEAL::0.788675, 0.605662, 0.125000 +DEAL::0.788675, 0.894338, 0.125000 +DEAL::Surface quadrature +DEAL::0.211325, 0.500000, 0.500000, 0.00000, 1.00000, 0.00000 +DEAL::0.788675, 0.500000, 0.500000, 0.00000, 1.00000, 0.00000 +DEAL:: +DEAL::plane direction = 2, face = 0 +DEAL::Inside quadrature +DEAL::0.211325, 0.105662, 0.125000 +DEAL::0.211325, 0.394338, 0.125000 +DEAL::0.788675, 0.105662, 0.125000 +DEAL::0.788675, 0.394338, 0.125000 +DEAL::Outside quadrature +DEAL::0.211325, 0.605662, 0.125000 +DEAL::0.211325, 0.894338, 0.125000 +DEAL::0.788675, 0.605662, 0.125000 +DEAL::0.788675, 0.894338, 0.125000 +DEAL::Surface quadrature +DEAL::0.211325, 0.500000, 0.500000, 0.00000, 0.00000, 1.00000 +DEAL::0.788675, 0.500000, 0.500000, 0.00000, 0.00000, 1.00000 +DEAL:: +DEAL::plane direction = 2, face = 1 +DEAL::Inside quadrature +DEAL::0.211325, 0.105662, 0.125000 +DEAL::0.211325, 0.394338, 0.125000 +DEAL::0.788675, 0.105662, 0.125000 +DEAL::0.788675, 0.394338, 0.125000 +DEAL::Outside quadrature +DEAL::0.211325, 0.605662, 0.125000 +DEAL::0.211325, 0.894338, 0.125000 +DEAL::0.788675, 0.605662, 0.125000 +DEAL::0.788675, 0.894338, 0.125000 +DEAL::Surface quadrature +DEAL::0.211325, 0.500000, 0.500000, 0.00000, 0.00000, 1.00000 +DEAL::0.788675, 0.500000, 0.500000, 0.00000, 0.00000, 1.00000 +DEAL:: +DEAL::plane direction = 2, face = 2 +DEAL::Inside quadrature +DEAL::0.105662, 0.211325, 0.125000 +DEAL::0.394338, 0.211325, 0.125000 +DEAL::0.105662, 0.788675, 0.125000 +DEAL::0.394338, 0.788675, 0.125000 +DEAL::Outside quadrature +DEAL::0.605662, 0.211325, 0.125000 +DEAL::0.894338, 0.211325, 0.125000 +DEAL::0.605662, 0.788675, 0.125000 +DEAL::0.894338, 0.788675, 0.125000 +DEAL::Surface quadrature +DEAL::0.500000, 0.211325, 0.500000, 0.00000, 0.00000, 1.00000 +DEAL::0.500000, 0.788675, 0.500000, 0.00000, 0.00000, 1.00000 +DEAL:: +DEAL::plane direction = 2, face = 3 +DEAL::Inside quadrature +DEAL::0.105662, 0.211325, 0.125000 +DEAL::0.394338, 0.211325, 0.125000 +DEAL::0.105662, 0.788675, 0.125000 +DEAL::0.394338, 0.788675, 0.125000 +DEAL::Outside quadrature +DEAL::0.605662, 0.211325, 0.125000 +DEAL::0.894338, 0.211325, 0.125000 +DEAL::0.605662, 0.788675, 0.125000 +DEAL::0.894338, 0.788675, 0.125000 +DEAL::Surface quadrature +DEAL::0.500000, 0.211325, 0.500000, 0.00000, 0.00000, 1.00000 +DEAL::0.500000, 0.788675, 0.500000, 0.00000, 0.00000, 1.00000 +DEAL:: diff --git a/tests/non_matching/quadrature_printing.h b/tests/non_matching/quadrature_printing.h index 004c8fccbf..827d745776 100644 --- a/tests/non_matching/quadrature_printing.h +++ b/tests/non_matching/quadrature_printing.h @@ -45,12 +45,12 @@ print_quadrature(const Quadrature &quadrature) /* * Print the incoming surface quadrature to deallog as comma separated values: - * p[0], ..., p[dim-1], weight, normal[0], ..., normal[dim-1] + * p[0], ..., p[dim-1], weight, normal[0], ..., normal[spacedim-1] */ -template +template void print_surface_quadrature( - const NonMatching::ImmersedSurfaceQuadrature &quadrature) + const NonMatching::ImmersedSurfaceQuadrature &quadrature) { for (unsigned int i = 0; i < quadrature.size(); ++i) { @@ -60,11 +60,12 @@ print_surface_quadrature( deallog << quadrature.weight(i); - const Tensor<1, dim> &normal = quadrature.normal_vector(i); - for (int d = 0; d < dim; d++) + const Tensor<1, spacedim> &normal = quadrature.normal_vector(i); + for (int d = 0; d < spacedim; d++) deallog << ", " << normal[d]; deallog << std::endl; } } + #endif -- 2.39.5