From 93260d538797dc3d4e4bb69f228e9ea74a17b0ff Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 5 Dec 2020 13:10:38 +0100 Subject: [PATCH] Add QGaussPyramid --- include/deal.II/simplex/quadrature_lib.h | 15 ++++++++ source/simplex/quadrature_lib.cc | 47 ++++++++++++++++++++++++ 2 files changed, 62 insertions(+) diff --git a/include/deal.II/simplex/quadrature_lib.h b/include/deal.II/simplex/quadrature_lib.h index 09037ff897..7962091ca5 100644 --- a/include/deal.II/simplex/quadrature_lib.h +++ b/include/deal.II/simplex/quadrature_lib.h @@ -66,6 +66,21 @@ namespace Simplex explicit QGaussWedge(const unsigned int n_points_1D); }; + /** + * Integration rule for pyramid entities. + */ + template + class QGaussPyramid : public Quadrature + { + 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 QGaussPyramid(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 277633a373..5d3ea4df22 100644 --- a/source/simplex/quadrature_lib.cc +++ b/source/simplex/quadrature_lib.cc @@ -168,6 +168,50 @@ namespace Simplex ExcMessage("No valid quadrature points!")); } + + + template + QGaussPyramid::QGaussPyramid(const unsigned int n_points_1D) + : Quadrature() + { + AssertDimension(dim, 3); + + if (n_points_1D == 1) + { + const double Q14 = 1.0 / 4.0; + const double Q43 = 4.0 / 3.0; + + this->quadrature_points.emplace_back(0, 0, Q14); + this->weights.emplace_back(Q43); + } + else if (n_points_1D == 2) + { + // clang-format off + this->quadrature_points.emplace_back(-0.26318405556971, -0.26318405556971, 0.54415184401122); + this->quadrature_points.emplace_back(-0.50661630334979, -0.50661630334979, 0.12251482265544); + this->quadrature_points.emplace_back(-0.26318405556971, +0.26318405556971, 0.54415184401122); + this->quadrature_points.emplace_back(-0.50661630334979, +0.50661630334979, 0.12251482265544); + this->quadrature_points.emplace_back(+0.26318405556971, -0.26318405556971, 0.54415184401122); + this->quadrature_points.emplace_back(+0.50661630334979, -0.50661630334979, 0.12251482265544); + this->quadrature_points.emplace_back(+0.26318405556971, +0.26318405556971, 0.54415184401122); + this->quadrature_points.emplace_back(+0.50661630334979, +0.50661630334979, 0.12251482265544); + // clang-format on + + this->weights.emplace_back(0.10078588207983); + this->weights.emplace_back(0.23254745125351); + this->weights.emplace_back(0.10078588207983); + this->weights.emplace_back(0.23254745125351); + this->weights.emplace_back(0.10078588207983); + this->weights.emplace_back(0.23254745125351); + this->weights.emplace_back(0.10078588207983); + this->weights.emplace_back(0.23254745125351); + } + + AssertDimension(this->quadrature_points.size(), this->weights.size()); + Assert(this->quadrature_points.size() > 0, + ExcMessage("No valid quadrature points!")); + } + } // namespace Simplex @@ -177,5 +221,8 @@ template class Simplex::QGauss<3>; template class Simplex::QGaussWedge<1>; template class Simplex::QGaussWedge<2>; template class Simplex::QGaussWedge<3>; +template class Simplex::QGaussPyramid<1>; +template class Simplex::QGaussPyramid<2>; +template class Simplex::QGaussPyramid<3>; DEAL_II_NAMESPACE_CLOSE -- 2.39.5