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 <int dim>
+ 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<dim> & level_set,
+ const BoundingBox<dim> &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<dim - 1> &
+ 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<dim - 1> &
+ 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<dim - 1, dim> &
+ 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<dim - 1> quadrature_generator;
+
+ /**
+ * The same surface quadrature as created by the quadrature_generator,
+ * but having dim-dimensional normals.
+ */
+ ImmersedSurfaceQuadrature<dim - 1, dim> surface_quadrature;
+ };
+
+
namespace internal
{
namespace QuadratureGeneratorImplementation
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/function_level_set.h>
+#include <deal.II/base/geometry_info.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/non_matching/quadrature_generator.h>
+
+#include <vector>
+
+#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 <int dim>
+void
+create_and_print_quadratures(const Function<dim> &levelset,
+ const unsigned int face_index)
+{
+ const hp::QCollection<1> quadratures1D(QGauss<1>(2));
+ const BoundingBox<dim> box = create_unit_bounding_box<dim>();
+
+ NonMatching::FaceQuadratureGenerator<dim> 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 <int dim>
+void
+test_plane_cuts_through_center()
+{
+ deallog << "test_vertical_cuts_through_center" << std::endl;
+ deallog << "dim=" << dim << std::endl;
+
+ Point<dim> 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<dim>::unit_vector(plane_direction);
+ const Functions::LevelSet::Plane<dim> levelset(center, normal);
+
+ // Test all faces that are intersected by the plane.
+ for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; f++)
+ {
+ const int normal_direction =
+ GeometryInfo<dim>::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>();
+}
--- /dev/null
+
+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::