Introduce QGaussWedge

diff --git a/include/deal.II/simplex/quadrature_lib.h b/include/deal.II/simplex/quadrature_lib.h
index 09e2e97abda841651c06d5b78f7574f427261e12..09037ff897b829431f3f00e15d167e75ab026651 100644 (file)
--- 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 bb60716e6c6ed7e67790fbead015319ba4af3a2f..277633a373b2006899bccae142fd6e1bd492d6b2 100644 (file)
--- 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