From 598ee3146cce3aa27ff2df95e75ec40c4156396e Mon Sep 17 00:00:00 2001 From: Simon Sticko Date: Wed, 6 Apr 2022 17:52:58 +0200 Subject: [PATCH] Specialize NonMatching::FaceQuadratureGenerator in 1D This class must be specialized in 1D, because the general FaceQuadratureGenerator class uses the QuadratureGenerator class internally, which does not make sense when dim-1 = 0. --- .../non_matching/quadrature_generator.h | 97 +++++++++++++++++++ source/non_matching/quadrature_generator.cc | 77 +++++++++++++++ .../non_matching/face_quadrature_generator.cc | 29 ++++++ .../face_quadrature_generator.output | 13 +++ 4 files changed, 216 insertions(+) diff --git a/include/deal.II/non_matching/quadrature_generator.h b/include/deal.II/non_matching/quadrature_generator.h index 2fa8214601..97653b097a 100644 --- a/include/deal.II/non_matching/quadrature_generator.h +++ b/include/deal.II/non_matching/quadrature_generator.h @@ -362,6 +362,103 @@ namespace NonMatching }; + /** + * Specialization of the FaceQuadratureGenerator class for the 1-dimensional + * case. + * + * In 1D, a face is only a point. Thus to generate the immersed + * quadrature rules we add a single 0-dimensional quadrature point to the + * inside or outside quadrature rule depending on if the level set function is + * positive or negative at the face. The added quadrature point will have + * weight equal to 1. The immersed surface quadrature over a face corresponds + * to integrating over a dim-1 dimensional curve. Thus, surface quadrature + * generated by this specialized class is always empty. + * + * This class must be specialized in 1D, because the general + * FaceQuadratureGenerator class uses the QuadratureGenerator + * class internally, which does not make sense when dim-1 = 0. + */ + template <> + class FaceQuadratureGenerator<1> + { + public: + using AdditionalData = AdditionalQGeneratorData; + + /** + * Constructor. The incoming hp::QCollection is not used. But this class + * must have the same signature as the non-specialized class. + */ + 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(). + */ + void + generate(const Function<1> & level_set, + const BoundingBox<1> &box, + const unsigned int face_index); + + /** + * @copydoc FaceQuadratureGenerator::get_inside_quadrature() + */ + const Quadrature<0> & + get_inside_quadrature() const; + + /** + * @copydoc FaceQuadratureGenerator::get_outside_quadrature() + */ + const Quadrature<0> & + get_outside_quadrature() const; + + /** + * Return the quadrature rule for the region + * $\{x \in F : \psi(x) = 0 \}$ + * where, $F$ is the face of the BoundingBox passed to generate(). + * + * @note In 1D, this quadrature always contains 0 points. + */ + const ImmersedSurfaceQuadrature<0, 1> & + get_surface_quadrature() const; + + /** + * This function does nothing. It only exist to be compatible with + * FaceQuadratureGenerator. + */ + void + set_1D_quadrature(const unsigned int q_index); + + private: + /** + * Quadrature for the region + * $\{x \in F : \psi(x) < 0 \}$. + * Created in the last call to generate(). + */ + Quadrature<0> inside_quadrature; + + /** + * Quadrature for the region + * $\{x \in F : \psi(x) > 0 \}$. + * Created in the last call to generate(). + */ + Quadrature<0> outside_quadrature; + + /** + * Quadrature for the region + * $\{x \in F : \psi(x) = 0 \}$. + * This quadrature always contains zero points in 1D. + */ + const ImmersedSurfaceQuadrature<0, 1> surface_quadrature; + }; + + namespace internal { namespace QuadratureGeneratorImplementation diff --git a/source/non_matching/quadrature_generator.cc b/source/non_matching/quadrature_generator.cc index 9f446a375d..da551e41a1 100644 --- a/source/non_matching/quadrature_generator.cc +++ b/source/non_matching/quadrature_generator.cc @@ -1448,6 +1448,83 @@ namespace NonMatching { return surface_quadrature; } + + + + FaceQuadratureGenerator<1>::FaceQuadratureGenerator( + const hp::QCollection<1> &quadratures1D, + const AdditionalData & additional_data) + { + (void)quadratures1D; + (void)additional_data; + } + + + + void + FaceQuadratureGenerator<1>::generate(const Function<1> & level_set, + const BoundingBox<1> &box, + const unsigned int face_index) + { + AssertIndexRange(face_index, GeometryInfo<1>::faces_per_cell); + + // The only vertex the 1D-face has. + const Point<1> vertex = + box.vertex(GeometryInfo<1>::face_to_cell_vertices(face_index, 0)); + + const unsigned int n_points = 1; + const double weight = 1; + const std::vector> points(n_points); + const std::vector weights(n_points, weight); + + const double level_set_value = level_set.value(vertex); + if (level_set_value < 0) + { + inside_quadrature = Quadrature<0>(points, weights); + outside_quadrature = Quadrature<0>(); + } + else if (level_set_value > 0) + { + inside_quadrature = Quadrature<0>(); + outside_quadrature = Quadrature<0>(points, weights); + } + else + { + inside_quadrature = Quadrature<0>(); + outside_quadrature = Quadrature<0>(); + } + } + + + + void + FaceQuadratureGenerator<1>::set_1D_quadrature(const unsigned int q_index) + { + (void)q_index; + } + + + + const Quadrature<0> & + FaceQuadratureGenerator<1>::get_inside_quadrature() const + { + return inside_quadrature; + } + + + const Quadrature<0> & + FaceQuadratureGenerator<1>::get_outside_quadrature() const + { + return outside_quadrature; + } + + + + const ImmersedSurfaceQuadrature<0, 1> & + FaceQuadratureGenerator<1>::get_surface_quadrature() const + { + return surface_quadrature; + } } // namespace NonMatching #include "quadrature_generator.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/tests/non_matching/face_quadrature_generator.cc b/tests/non_matching/face_quadrature_generator.cc index 4c0964ee7a..afb24dd8ce 100644 --- a/tests/non_matching/face_quadrature_generator.cc +++ b/tests/non_matching/face_quadrature_generator.cc @@ -105,11 +105,40 @@ test_plane_cuts_through_center() +/* + * Test the 1D-specialization of the FaceQuadratureGenerator class. + * Set up a 1D-box [0,1] and a level set function: psi(x) = x - 0.5, + * so that the level set function is negative at the left face and positive at + * the right. Generate quadrature rules at both faces and check that a single + * quadrature point at the inside/outside quadrature at the left/right face is + * generated. + */ +void +test_1D() +{ + deallog << "test_1D" << std::endl; + + const int dim = 1; + Point center(.5); + const Tensor<1, dim> normal = Point::unit_vector(0); + const Functions::LevelSet::Plane levelset(center, normal); + + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + { + deallog << "face = " << f << std::endl; + create_and_print_quadratures(levelset, f); + deallog << std::endl; + } +} + + + int main() { initlog(); + test_1D(); 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 index c3a17d9259..397a723acf 100644 --- a/tests/non_matching/face_quadrature_generator.output +++ b/tests/non_matching/face_quadrature_generator.output @@ -1,4 +1,17 @@ +DEAL::test_1D +DEAL::face = 0 +DEAL::Inside quadrature +DEAL::1.00000 +DEAL::Outside quadrature +DEAL::Surface quadrature +DEAL:: +DEAL::face = 1 +DEAL::Inside quadrature +DEAL::Outside quadrature +DEAL::1.00000 +DEAL::Surface quadrature +DEAL:: DEAL::test_vertical_cuts_through_center DEAL::dim=2 DEAL::plane direction = 0, face = 2 -- 2.39.5