compute_affine_transformation in codim1 and codim2 case

diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h
index 87b2718a06c41fab5bd06a27a468c2933b7cc041..6f1ad61c531337b548b1b8386aa5ab64a456ae08 100644 (file)
--- a/include/deal.II/base/quadrature_lib.h
+++ b/include/deal.II/base/quadrature_lib.h
@@ -609,6 +609,7 @@ public:
   QSimplex(const Quadrature<dim> &quad);
 
   /**
+   *
    * Return an affine transformation of this quadrature, that can be used to
    * integrate on the simplex identified by `vertices`.
    *
@@ -626,14 +627,19 @@ public:
    * that is $J \dealcoloneq |\text{det}(B)|$. If $J$ is zero, an empty
    * quadrature is returned. This may happen, in two dimensions, if the three
    * vertices are aligned, or in three dimensions if the four vertices are on
-   * the same plane.
+   * the same plane. The present function works also in the codimension one and
+   * codimension two case. For instance, when `dim=2` and `spacedim=3`, we can
+   * map the quadrature points so that they live on the physical triangle
+   * embedded in the three dimensional space. In such a case, the matrix $B$ is
+   * not square anymore.
    *
    * @param[in] vertices The vertices of the simplex you wish to integrate on
    * @return A quadrature object that can be used to integrate on the simplex
    */
-  Quadrature<dim>
+  template <int spacedim = dim>
+  Quadrature<spacedim>
   compute_affine_transformation(
-    const std::array<Point<dim>, dim + 1> &vertices) const;
+    const std::array<Point<spacedim>, dim + 1> &vertices) const;
 };
 
 /**

diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc
index 368cd980feee46a4e7afe66c741f415bea794854..bfe1d52b236df51479b971888d7c80e76175652d 100644 (file)
--- a/source/base/quadrature_lib.cc
+++ b/source/base/quadrature_lib.cc
@@ -1227,31 +1227,35 @@ QSimplex<dim>::QSimplex(const Quadrature<dim> &quad)
 
 
 template <int dim>
-Quadrature<dim>
+template <int spacedim>
+Quadrature<spacedim>
 QSimplex<dim>::compute_affine_transformation(
-  const std::array<Point<dim>, dim + 1> &vertices) const
+  const std::array<Point<spacedim>, dim + 1> &vertices) const
 {
-  Tensor<2, dim> B;
+  Assert(dim <= spacedim,
+         ExcMessage("Invalid combination of dim and spacedim ."));
+  DerivativeForm<1, spacedim, dim> Bt;
   for (unsigned int d = 0; d < dim; ++d)
-    B[d] = vertices[d + 1] - vertices[0];
+    Bt[d] = vertices[d + 1] - vertices[0];
 
-  B              = transpose(B);
-  const double J = std::abs(determinant(B));
+  const auto   B = Bt.transpose();
+  const double J = std::abs(B.determinant());
 
   // if the determinant is zero, we return an empty quadrature
   if (J < 1e-12)
-    return Quadrature<dim>();
+    return Quadrature<spacedim>();
 
-  std::vector<Point<dim>> qp(this->size());
-  std::vector<double>     w(this->size());
+  std::vector<Point<spacedim>> qp(this->size());
+  std::vector<double>          w(this->size());
 
   for (unsigned int i = 0; i < this->size(); ++i)
     {
-      qp[i] = Point<dim>(vertices[0] + B * this->point(i));
-      w[i]  = J * this->weight(i);
+      qp[i] =
+        Point<spacedim>(vertices[0] + apply_transformation(B, this->point(i)));
+      w[i] = J * this->weight(i);
     }
 
-  return Quadrature<dim>(qp, w);
+  return Quadrature<spacedim>(qp, w);
 }
 
 
@@ -2208,3 +2212,18 @@ template class QWitherdenVincentSimplex<2>;
 template class QWitherdenVincentSimplex<3>;
 
 DEAL_II_NAMESPACE_CLOSE
+
+namespace dealii
+{
+  template Quadrature<2>
+  QSimplex<1>::compute_affine_transformation(
+    const std::array<Point<2>, 1 + 1> &vertices) const;
+
+  template Quadrature<3>
+  QSimplex<1>::compute_affine_transformation(
+    const std::array<Point<3>, 1 + 1> &vertices) const;
+
+  template Quadrature<3>
+  QSimplex<2>::compute_affine_transformation(
+    const std::array<Point<3>, 2 + 1> &vertices) const;
+} // namespace dealii

diff --git a/tests/base/quadrature_simplex_01_codim1.cc b/tests/base/quadrature_simplex_01_codim1.cc
new file mode 100644 (file)
index 0000000..e50ff08
--- /dev/null
+++ b/tests/base/quadrature_simplex_01_codim1.cc
@@ -0,0 +1,59 @@
+// ---------------------------------------------------------------------
+
+// Copyright (C) 2022 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.
+
+// ---------------------------------------------------------------------
+
+// construct a simplex quadrature, and check that the affine transformation
+// gives the correct results in the codimension one case
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <numeric>
+
+#include "../tests.h"
+
+// #include "simplex.h"
+
+template <int dim, int spacedim>
+void
+test(const int degree, const std::array<Point<spacedim>, spacedim> &vertices)
+{
+  QGaussSimplex<dim> quad(degree);
+
+  deallog << "# dim = " << dim << std::endl;
+  deallog << "# spacedim = " << spacedim << std::endl;
+
+  auto quad2 = quad.compute_affine_transformation(vertices);
+
+  for (auto p : quad2.get_points())
+    deallog << p << std::endl;
+
+  deallog << std::endl
+          << "# Area: " << std::setprecision(15)
+          << std::accumulate(quad2.get_weights().begin(),
+                             quad2.get_weights().end(),
+                             0.0)
+          << std::endl
+          << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+
+  test<1, 2>(1, {{Point<2>(0., 2.), Point<2>(1., 0.)}});
+  test<2, 3>(
+    1, {{Point<3>(0., 2., 0.), Point<3>(1., 2., 4.), Point<3>(1., 1., 1.)}});
+}

diff --git a/tests/base/quadrature_simplex_01_codim1.output b/tests/base/quadrature_simplex_01_codim1.output
new file mode 100644 (file)
index 0000000..268045a
--- /dev/null
+++ b/tests/base/quadrature_simplex_01_codim1.output
@@ -0,0 +1,13 @@
+
+DEAL::# dim = 1
+DEAL::# spacedim = 2
+DEAL::0.500000 1.00000
+DEAL::
+DEAL::# Area: 2.23606797749979
+DEAL::
+DEAL::# dim = 2
+DEAL::# spacedim = 3
+DEAL::0.666666666666667 1.66666666666667 1.66666666666667
+DEAL::
+DEAL::# Area: 2.54950975679639
+DEAL::

diff --git a/tests/base/quadrature_simplex_02_codim2.cc b/tests/base/quadrature_simplex_02_codim2.cc
new file mode 100644 (file)
index 0000000..be07845
--- /dev/null
+++ b/tests/base/quadrature_simplex_02_codim2.cc
@@ -0,0 +1,56 @@
+// ---------------------------------------------------------------------
+
+// Copyright (C) 2022 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.
+
+// ---------------------------------------------------------------------
+
+// construct a simplex quadrature, and check that the affine transformation
+// gives the correct results in the codimension one case
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <numeric>
+
+#include "../tests.h"
+
+// #include "simplex.h"
+
+template <int dim, int spacedim>
+void
+test(const int degree, const std::array<Point<spacedim>, dim + 1> &vertices)
+{
+  QGaussSimplex<dim> quad(degree);
+
+  deallog << "# dim = " << dim << std::endl;
+  deallog << "# spacedim = " << spacedim << std::endl;
+
+  auto quad2 = quad.compute_affine_transformation(vertices);
+
+  for (auto p : quad2.get_points())
+    deallog << p << std::endl;
+
+  deallog << std::endl
+          << "# Area: " << std::setprecision(15)
+          << std::accumulate(quad2.get_weights().begin(),
+                             quad2.get_weights().end(),
+                             0.0)
+          << std::endl
+          << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+  test<1, 3>(1, {{Point<3>(0., 1., 0.), Point<3>(0., 0., 1.)}});
+}

diff --git a/tests/base/quadrature_simplex_02_codim2.output b/tests/base/quadrature_simplex_02_codim2.output
new file mode 100644 (file)
index 0000000..afcbea4
--- /dev/null
+++ b/tests/base/quadrature_simplex_02_codim2.output
@@ -0,0 +1,7 @@
+
+DEAL::# dim = 1
+DEAL::# spacedim = 3
+DEAL::0.00000 0.500000 0.500000
+DEAL::
+DEAL::# Area: 1.41421356237310
+DEAL::