From d2b14a3162ce6e54d6eacc6ce25ccd6e54564342 Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Sat, 28 Nov 2020 09:31:24 +0100
Subject: [PATCH] Introduce QGaussWedge

---
 include/deal.II/simplex/quadrature_lib.h | 16 ++++++++++++++
 source/simplex/quadrature_lib.cc         | 28 ++++++++++++++++++++++++
 2 files changed, 44 insertions(+)

diff --git a/include/deal.II/simplex/quadrature_lib.h b/include/deal.II/simplex/quadrature_lib.h
index 09e2e97abd..09037ff897 100644
--- a/include/deal.II/simplex/quadrature_lib.h
+++ b/include/deal.II/simplex/quadrature_lib.h
@@ -50,6 +50,22 @@ namespace Simplex
      */
     explicit QGauss(const unsigned int n_points_1D);
   };
+
+  /**
+   * Integration rule for wedge entities.
+   */
+  template <int dim>
+  class QGaussWedge : public Quadrature<dim>
+  {
+  public:
+    /**
+     * Users specify a number `n_points_1D` as an indication of what polynomial
+     * degree to be integrated exactly. For details, see the comments of
+     * Simplex::QGauss.
+     */
+    explicit QGaussWedge(const unsigned int n_points_1D);
+  };
+
 } // namespace Simplex
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/source/simplex/quadrature_lib.cc b/source/simplex/quadrature_lib.cc
index bb60716e6c..277633a373 100644
--- a/source/simplex/quadrature_lib.cc
+++ b/source/simplex/quadrature_lib.cc
@@ -143,11 +143,39 @@ namespace Simplex
            ExcMessage("No valid quadrature points!"));
   }
 
+
+
+  template <int dim>
+  QGaussWedge<dim>::QGaussWedge(const unsigned int n_points)
+    : Quadrature<dim>()
+  {
+    AssertDimension(dim, 3);
+
+    Simplex::QGauss<2> quad_tri(n_points);
+    QGauss<1>          quad_line(n_points);
+
+    for (unsigned int i = 0; i < quad_line.size(); ++i)
+      for (unsigned int j = 0; j < quad_tri.size(); ++j)
+        {
+          this->quadrature_points.emplace_back(quad_tri.point(j)[0],
+                                               quad_tri.point(j)[1],
+                                               quad_line.point(i)[0]);
+          this->weights.emplace_back(quad_tri.weight(j) * quad_line.weight(i));
+        }
+
+    AssertDimension(this->quadrature_points.size(), this->weights.size());
+    Assert(this->quadrature_points.size() > 0,
+           ExcMessage("No valid quadrature points!"));
+  }
+
 } // namespace Simplex
 
 
 template class Simplex::QGauss<1>;
 template class Simplex::QGauss<2>;
 template class Simplex::QGauss<3>;
+template class Simplex::QGaussWedge<1>;
+template class Simplex::QGaussWedge<2>;
+template class Simplex::QGaussWedge<3>;
 
 DEAL_II_NAMESPACE_CLOSE
-- 
2.39.5