From 90d10ab79a26dd4a51e5cbbf9b35bbaec731bce1 Mon Sep 17 00:00:00 2001 From: Simon Sticko Date: Mon, 14 Jun 2021 15:26:42 +0200 Subject: [PATCH] Add a spacedim template parameter on ImmersedSurfaceQuadrature To handle the case when we have a line integral immersed in a face, the normals on ImmersedSurfaceQuadrature need to be vectors in spacedim. The reason is that the conormal can not be computed without knowing the normal to the surface. Add a spacedim template parameter and change the normals to Tensor<1, spacedim>. Rewrite the documentation and modify the test to a account for this. --- .../immersed_surface_quadrature.h | 95 ++++++++++++------- .../immersed_surface_quadrature.cc | 40 ++++---- .../immersed_surface_quadrature.cc | 52 +++++----- .../immersed_surface_quadrature.output | 8 ++ 4 files changed, 120 insertions(+), 75 deletions(-) diff --git a/include/deal.II/non_matching/immersed_surface_quadrature.h b/include/deal.II/non_matching/immersed_surface_quadrature.h index b917cd82f1..a594aab63c 100644 --- a/include/deal.II/non_matching/immersed_surface_quadrature.h +++ b/include/deal.II/non_matching/immersed_surface_quadrature.h @@ -29,48 +29,75 @@ DEAL_II_NAMESPACE_OPEN namespace NonMatching { /** - * This class defines a quadrature formula for integration over an - * oriented surface, $\hat{S}$, immersed in the unit cell. By - * immersed it is meant that the surface may intersect the unit cell - * in an arbitrary way. The quadrature formula is described by a set - * of quadrature points, $\hat{x}_q$, weights, $w_q$, and normalized - * surface normals, $\hat{n}_q$. + * This class defines a quadrature formula to integrate over the intersection + * between an oriented surface, $\hat{S}$, and a cell or face. The word + * "immersed" in the class name reflects that the surface may intersect the + * cell/face in an arbitrary way. * - * We typically want to compute surface integrals in real space. - * A surface $S$ intersecting a cell $K$ in real space, can be mapped onto a - * surface $\hat{S}$ intersecting the unit cell $\hat{K}$. - * Thus a surface integral over $S\cap K$ in real space can be transformed to - * a surface integral over $\hat{S} \cap \hat{K}$ according to + * The spacedim template parameter of this class is the dimension that the + * (spacedim-1)-dimensional surface is embedded in: + * $\hat{S} \subset \mathbb{R}^{\text{spacedim}}$. The dim parameter + * describes the dimension of the "object" that the surface intersects. That + * is, dim = spacedim corresponds to the surface intersecting a cell and + * dim = spacedim - 1 corresponds to the surface intersecting a face. The + * quadrature formula is described by a set of quadrature points, + * $\hat{x}_q \in \mathbb{R}^{\text{dim}}$, weights, $w_q$, and normalized + * surface normals, $\hat{n}_q \in \mathbb{R}^{\text{spacedim}}$. + * + * Consider first the case dim = spacedim. We typically want to compute + * integrals in real space. A surface, $S$, intersecting a cell, $K$, in + * real space can be mapped onto a surface, $\hat{S}$, intersecting the unit + * cell, $\hat{K}$. Thus an integral over $S\cap K$ in real space can be + * transformed to an integral over $\hat{S} \cap \hat{K}$ according to * @f[ - * \int_{S\cap K} f(x) dS = - * \int_{S\cap K} f(x) |d\bar{S}| = - * \int_{\hat{S}\cap\hat{K}} f(F_{K}(\hat{x})) \det(J) |\left( J^{-1} \right + * \int_{S\cap K} f dS = + * \int_{S\cap K} f |d\bar{S}| = + * \int_{\hat{S}\cap\hat{K}} f \circ F_{K} \det(J) |\left( J^{-1} \right * )^T d\hat{S}|, * @f] * where $F_K$ is the mapping from reference to real space and $J$ is its - * Jacobian. This transformation is possible since the continuous surface - * elements are vectors: $d\bar{S}, d\hat{S} \in \mathbb{R}^{dim}$ which are - * parallel to the normals of $S$ and $\hat{S}$. So in order to compute the - * integral in real space one needs information about the normal to do the - * transformation. - * - * Thus, in addition to storing points and weights, this quadrature stores - * also the normalized normal for each quadrature point. This can be viewed - * as storing a discrete surface element, + * Jacobian matrix. This transformation is possible since the continuous + * surface elements are vectors: + * $d\bar{S}, d\hat{S} \in \mathbb{R}^{spacedim}$, which are parallel to the + * normals of $S$ and $\hat{S}$. That is, the normal is needed to do the + * transformation. Thus, in addition to storing points and weights, this + * quadrature stores also the normalized normal for each quadrature point. + * This can be viewed as storing a discrete surface element, * @f[ * \Delta \hat{S}_q \dealcoloneq w_q \hat{n}_q \approx d\hat{S}(\hat{x}_q), * @f] * for each quadrature point. The surface integral in real space would then be * approximated as * @f[ - * \int_{S\cap K} f(x) dS \approx + * \int_{S\cap K} f dS \approx * \sum_{q} f \left(F_{K}(\hat{x}_{q}) \right) \det(J_q) * |\left( J_q^{-1} \right)^T \hat{n}_q| w_q. * @f] * * @image html immersed_surface_quadrature.svg + * + * When dim = spacedim - 1, this class represents a (spacedim-2)-dimensional + * integral. That is, if spacedim = 3 we have a line integral immersed in a + * face. Here, the transformation between the face, $F$, and reference face, + * $\hat{F}$, reads + * @f[ + * \int_{S\cap F} f dr = + * \int_{S\cap F} f |d\bar{r}| = + * \int_{\hat{S}\cap\hat{F}} f \circ F_{K} | J d\hat{r}| + * \approx \sum_{q} f \left(F_{K}(\hat{x}_{q}) \right) |J_q \hat{t}_q| w_q, + * @f] + * where $\hat{t}_q = \hat{n}_q \times \hat{n}_F$ is the tangent to the curve + * at $\hat{x}_q$ and $\hat{n}_F$ is the face normal. It would be possible to + * compute the tangent by only knowing the normal to the curve in the face + * plane (i.e. the dim-dimensional normal). However, when these quadratures + * are used, the weak form typically involves the so-called conormal, which + * can not be computed without knowing the surface normal in + * $\mathbb{R}^{\text{spacedim}}$. The conormal is the unit vector parallel to + * the projection of the face normal into the surface plane. This is + * essentially the same thing as the normalized @ref GlossBoundaryForm + * "boundary form". */ - template + template class ImmersedSurfaceQuadrature : public Quadrature { public: @@ -85,9 +112,9 @@ namespace NonMatching * surface normals. The points, weights and normals should be with respect * to reference space, and the normals should be normalized. */ - ImmersedSurfaceQuadrature(const std::vector> & points, - const std::vector & weights, - const std::vector> &normals); + ImmersedSurfaceQuadrature(const std::vector> & points, + const std::vector & weights, + const std::vector> &normals); /** * Extend the given formula by an additional quadrature point. @@ -105,27 +132,27 @@ namespace NonMatching * quadrature formula. */ void - push_back(const Point & point, - const double weight, - const Tensor<1, dim> &normal); + push_back(const Point & point, + const double weight, + const Tensor<1, spacedim> &normal); /** * Return a reference to the ith surface normal. */ - const Tensor<1, dim> & + const Tensor<1, spacedim> & normal_vector(const unsigned int i) const; /** * Return a reference to the whole %vector of normals. */ - const std::vector> & + const std::vector> & get_normal_vectors() const; protected: /** * %Vector of surface normals at each quadrature point. */ - std::vector> normals; + std::vector> normals; }; } // namespace NonMatching diff --git a/source/non_matching/immersed_surface_quadrature.cc b/source/non_matching/immersed_surface_quadrature.cc index 117b71c06d..b5529864af 100644 --- a/source/non_matching/immersed_surface_quadrature.cc +++ b/source/non_matching/immersed_surface_quadrature.cc @@ -18,11 +18,11 @@ DEAL_II_NAMESPACE_OPEN namespace NonMatching { - template - ImmersedSurfaceQuadrature::ImmersedSurfaceQuadrature( - const std::vector> & points, - const std::vector & weights, - const std::vector> &normals) + template + ImmersedSurfaceQuadrature::ImmersedSurfaceQuadrature( + const std::vector> & points, + const std::vector & weights, + const std::vector> &normals) : Quadrature(points, weights) , normals(normals) { @@ -38,11 +38,12 @@ namespace NonMatching - template + template void - ImmersedSurfaceQuadrature::push_back(const Point & point, - const double weight, - const Tensor<1, dim> &normal) + ImmersedSurfaceQuadrature::push_back( + const Point & point, + const double weight, + const Tensor<1, spacedim> &normal) { this->quadrature_points.push_back(point); this->weights.push_back(weight); @@ -53,9 +54,10 @@ namespace NonMatching - template - const Tensor<1, dim> & - ImmersedSurfaceQuadrature::normal_vector(const unsigned int i) const + template + const Tensor<1, spacedim> & + ImmersedSurfaceQuadrature::normal_vector( + const unsigned int i) const { AssertIndexRange(i, this->size()); return normals[i]; @@ -63,18 +65,20 @@ namespace NonMatching - template - const std::vector> & - ImmersedSurfaceQuadrature::get_normal_vectors() const + template + const std::vector> & + ImmersedSurfaceQuadrature::get_normal_vectors() const { return normals; } - template class ImmersedSurfaceQuadrature<1>; - template class ImmersedSurfaceQuadrature<2>; - template class ImmersedSurfaceQuadrature<3>; + template class ImmersedSurfaceQuadrature<1, 1>; + template class ImmersedSurfaceQuadrature<2, 2>; + template class ImmersedSurfaceQuadrature<3, 3>; + template class ImmersedSurfaceQuadrature<1, 2>; + template class ImmersedSurfaceQuadrature<2, 3>; } // namespace NonMatching DEAL_II_NAMESPACE_CLOSE diff --git a/tests/non_matching/immersed_surface_quadrature.cc b/tests/non_matching/immersed_surface_quadrature.cc index ba769d12c0..7c6e8a9311 100644 --- a/tests/non_matching/immersed_surface_quadrature.cc +++ b/tests/non_matching/immersed_surface_quadrature.cc @@ -26,9 +26,10 @@ -template +template void -print_quadrature(const NonMatching::ImmersedSurfaceQuadrature &quadrature) +print_quadrature( + const NonMatching::ImmersedSurfaceQuadrature &quadrature) { for (unsigned int i = 0; i < quadrature.size(); ++i) { @@ -41,9 +42,10 @@ print_quadrature(const NonMatching::ImmersedSurfaceQuadrature &quadrature) // Check that get_normals() are callable and are of the same size as // points and weights. -template +template void -check_get_normals(const NonMatching::ImmersedSurfaceQuadrature &quadrature) +check_get_normals( + const NonMatching::ImmersedSurfaceQuadrature &quadrature) { const std::vector> & points = quadrature.get_points(); const std::vector> &normals = quadrature.get_normal_vectors(); @@ -52,34 +54,34 @@ check_get_normals(const NonMatching::ImmersedSurfaceQuadrature &quadrature) -template +template void test_non_default_constructor() { deallog << "Using constructor" << std::endl; - std::vector> points(1); - std::vector weights(1, 1); - std::vector> normals; - normals.push_back(Point::unit_vector(dim - 1)); - NonMatching::ImmersedSurfaceQuadrature quadrature(points, - weights, - normals); + std::vector> points(1); + std::vector weights(1, 1); + std::vector> normals; + normals.push_back(Point::unit_vector(spacedim - 1)); + NonMatching::ImmersedSurfaceQuadrature quadrature(points, + weights, + normals); print_quadrature(quadrature); } -template +template void test_push_back() { deallog << "Using push_back" << std::endl; - const Point point; - const double weight = 1; - const Tensor<1, dim> normal = Point::unit_vector(dim - 1); + const Point point; + const double weight = 1; + const Tensor<1, spacedim> normal = Point::unit_vector(spacedim - 1); - NonMatching::ImmersedSurfaceQuadrature quadrature; + NonMatching::ImmersedSurfaceQuadrature quadrature; quadrature.push_back(point, weight, normal); print_quadrature(quadrature); @@ -87,12 +89,12 @@ test_push_back() -template +template void construct_quadrature_and_print_points() { - test_push_back(); - test_non_default_constructor(); + test_push_back(); + test_non_default_constructor(); } @@ -101,7 +103,11 @@ int main() { initlog(); - construct_quadrature_and_print_points<1>(); - construct_quadrature_and_print_points<2>(); - construct_quadrature_and_print_points<3>(); + construct_quadrature_and_print_points<1, 1>(); + construct_quadrature_and_print_points<2, 2>(); + construct_quadrature_and_print_points<3, 3>(); + + // Face quadrature + construct_quadrature_and_print_points<1, 2>(); + construct_quadrature_and_print_points<2, 3>(); } diff --git a/tests/non_matching/immersed_surface_quadrature.output b/tests/non_matching/immersed_surface_quadrature.output index 7143d449ee..efb667c09a 100644 --- a/tests/non_matching/immersed_surface_quadrature.output +++ b/tests/non_matching/immersed_surface_quadrature.output @@ -11,3 +11,11 @@ DEAL::Using push_back DEAL::0.00000 0.00000 0.00000, 1.00000, 0.00000 0.00000 1.00000 DEAL::Using constructor DEAL::0.00000 0.00000 0.00000, 1.00000, 0.00000 0.00000 1.00000 +DEAL::Using push_back +DEAL::0.00000, 1.00000, 0.00000 1.00000 +DEAL::Using constructor +DEAL::0.00000, 1.00000, 0.00000 1.00000 +DEAL::Using push_back +DEAL::0.00000 0.00000, 1.00000, 0.00000 0.00000 1.00000 +DEAL::Using constructor +DEAL::0.00000 0.00000, 1.00000, 0.00000 0.00000 1.00000 -- 2.39.5