]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add a spacedim template parameter on ImmersedSurfaceQuadrature 12459/head
authorSimon Sticko <simon@sticko.se>
Mon, 14 Jun 2021 13:26:42 +0000 (15:26 +0200)
committerSimon Sticko <simon@sticko.se>
Mon, 14 Jun 2021 13:26:42 +0000 (15:26 +0200)
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.

include/deal.II/non_matching/immersed_surface_quadrature.h
source/non_matching/immersed_surface_quadrature.cc
tests/non_matching/immersed_surface_quadrature.cc
tests/non_matching/immersed_surface_quadrature.output

index b917cd82f1eca5eaeb759e20d7c804b87b7e8c46..a594aab63c60cbd37bf754e676776795107e7e51 100644 (file)
@@ -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 <int dim>
+  template <int dim, int spacedim = dim>
   class ImmersedSurfaceQuadrature : public Quadrature<dim>
   {
   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<Point<dim>> &    points,
-                              const std::vector<double> &        weights,
-                              const std::vector<Tensor<1, dim>> &normals);
+    ImmersedSurfaceQuadrature(const std::vector<Point<dim>> &         points,
+                              const std::vector<double> &             weights,
+                              const std::vector<Tensor<1, spacedim>> &normals);
 
     /**
      * Extend the given formula by an additional quadrature point.
@@ -105,27 +132,27 @@ namespace NonMatching
      * quadrature formula.
      */
     void
-    push_back(const Point<dim> &    point,
-              const double          weight,
-              const Tensor<1, dim> &normal);
+    push_back(const Point<dim> &         point,
+              const double               weight,
+              const Tensor<1, spacedim> &normal);
 
     /**
      * Return a reference to the <tt>i</tt>th 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<Tensor<1, dim>> &
+    const std::vector<Tensor<1, spacedim>> &
     get_normal_vectors() const;
 
   protected:
     /**
      * %Vector of surface normals at each quadrature point.
      */
-    std::vector<Tensor<1, dim>> normals;
+    std::vector<Tensor<1, spacedim>> normals;
   };
 
 } // namespace NonMatching
index 117b71c06da383000bf1ba4a4b2b866f008b796d..b5529864afe00153e650f567941065a4c7afa74f 100644 (file)
 DEAL_II_NAMESPACE_OPEN
 namespace NonMatching
 {
-  template <int dim>
-  ImmersedSurfaceQuadrature<dim>::ImmersedSurfaceQuadrature(
-    const std::vector<Point<dim>> &    points,
-    const std::vector<double> &        weights,
-    const std::vector<Tensor<1, dim>> &normals)
+  template <int dim, int spacedim>
+  ImmersedSurfaceQuadrature<dim, spacedim>::ImmersedSurfaceQuadrature(
+    const std::vector<Point<dim>> &         points,
+    const std::vector<double> &             weights,
+    const std::vector<Tensor<1, spacedim>> &normals)
     : Quadrature<dim>(points, weights)
     , normals(normals)
   {
@@ -38,11 +38,12 @@ namespace NonMatching
 
 
 
-  template <int dim>
+  template <int dim, int spacedim>
   void
-  ImmersedSurfaceQuadrature<dim>::push_back(const Point<dim> &    point,
-                                            const double          weight,
-                                            const Tensor<1, dim> &normal)
+  ImmersedSurfaceQuadrature<dim, spacedim>::push_back(
+    const Point<dim> &         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 <int dim>
-  const Tensor<1, dim> &
-  ImmersedSurfaceQuadrature<dim>::normal_vector(const unsigned int i) const
+  template <int dim, int spacedim>
+  const Tensor<1, spacedim> &
+  ImmersedSurfaceQuadrature<dim, spacedim>::normal_vector(
+    const unsigned int i) const
   {
     AssertIndexRange(i, this->size());
     return normals[i];
@@ -63,18 +65,20 @@ namespace NonMatching
 
 
 
-  template <int dim>
-  const std::vector<Tensor<1, dim>> &
-  ImmersedSurfaceQuadrature<dim>::get_normal_vectors() const
+  template <int dim, int spacedim>
+  const std::vector<Tensor<1, spacedim>> &
+  ImmersedSurfaceQuadrature<dim, spacedim>::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
index ba769d12c0a8dc236c65af94f24c2cbef4cd56ea..7c6e8a9311a1550d22d9d41c9ec5982287b5dc9c 100644 (file)
 
 
 
-template <int dim>
+template <int dim, int spacedim>
 void
-print_quadrature(const NonMatching::ImmersedSurfaceQuadrature<dim> &quadrature)
+print_quadrature(
+  const NonMatching::ImmersedSurfaceQuadrature<dim, spacedim> &quadrature)
 {
   for (unsigned int i = 0; i < quadrature.size(); ++i)
     {
@@ -41,9 +42,10 @@ print_quadrature(const NonMatching::ImmersedSurfaceQuadrature<dim> &quadrature)
 
 // Check that get_normals() are callable and are of the same size as
 // points and weights.
-template <int dim>
+template <int dim, int spacedim>
 void
-check_get_normals(const NonMatching::ImmersedSurfaceQuadrature<dim> &quadrature)
+check_get_normals(
+  const NonMatching::ImmersedSurfaceQuadrature<dim, spacedim> &quadrature)
 {
   const std::vector<Point<dim>> &    points  = quadrature.get_points();
   const std::vector<Tensor<1, dim>> &normals = quadrature.get_normal_vectors();
@@ -52,34 +54,34 @@ check_get_normals(const NonMatching::ImmersedSurfaceQuadrature<dim> &quadrature)
 
 
 
-template <int dim>
+template <int dim, int spacedim>
 void
 test_non_default_constructor()
 {
   deallog << "Using constructor" << std::endl;
-  std::vector<Point<dim>>     points(1);
-  std::vector<double>         weights(1, 1);
-  std::vector<Tensor<1, dim>> normals;
-  normals.push_back(Point<dim>::unit_vector(dim - 1));
-  NonMatching::ImmersedSurfaceQuadrature<dim> quadrature(points,
-                                                         weights,
-                                                         normals);
+  std::vector<Point<dim>>          points(1);
+  std::vector<double>              weights(1, 1);
+  std::vector<Tensor<1, spacedim>> normals;
+  normals.push_back(Point<spacedim>::unit_vector(spacedim - 1));
+  NonMatching::ImmersedSurfaceQuadrature<dim, spacedim> quadrature(points,
+                                                                   weights,
+                                                                   normals);
 
   print_quadrature(quadrature);
 }
 
 
 
-template <int dim>
+template <int dim, int spacedim>
 void
 test_push_back()
 {
   deallog << "Using push_back" << std::endl;
-  const Point<dim>     point;
-  const double         weight = 1;
-  const Tensor<1, dim> normal = Point<dim>::unit_vector(dim - 1);
+  const Point<dim>          point;
+  const double              weight = 1;
+  const Tensor<1, spacedim> normal = Point<spacedim>::unit_vector(spacedim - 1);
 
-  NonMatching::ImmersedSurfaceQuadrature<dim> quadrature;
+  NonMatching::ImmersedSurfaceQuadrature<dim, spacedim> quadrature;
   quadrature.push_back(point, weight, normal);
 
   print_quadrature(quadrature);
@@ -87,12 +89,12 @@ test_push_back()
 
 
 
-template <int dim>
+template <int dim, int spacedim>
 void
 construct_quadrature_and_print_points()
 {
-  test_push_back<dim>();
-  test_non_default_constructor<dim>();
+  test_push_back<dim, spacedim>();
+  test_non_default_constructor<dim, spacedim>();
 }
 
 
@@ -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>();
 }
index 7143d449ee29c61f6fcd2ceb2bc0f2fef3e3e04f..efb667c09afb9fec94987b5ce09b8a17742d6657 100644 (file)
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.