From e73794256eb2f39d03282e97e9fb5c667ccaa4d2 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 20 Dec 2020 10:57:04 +0100 Subject: [PATCH] Implement FE_WedgeDGP and FE_PyramidDGP --- include/deal.II/fe/fe_base.h | 12 ++ include/deal.II/simplex/fe_lib.h | 102 ++++++++++++++++- source/fe/fe_data.cc | 6 +- source/simplex/fe_lib.cc | 181 ++++++++++++++++++++++++------- source/simplex/fe_lib.inst.in | 8 ++ tests/simplex/fe_lib_01.cc | 28 ++++- tests/simplex/fe_lib_01.output | 86 ++++++++++++--- 7 files changed, 360 insertions(+), 63 deletions(-) diff --git a/include/deal.II/fe/fe_base.h b/include/deal.II/fe/fe_base.h index ac3deeac3e..803776665e 100644 --- a/include/deal.II/fe/fe_base.h +++ b/include/deal.II/fe/fe_base.h @@ -673,6 +673,18 @@ public: get_first_face_quad_index(const unsigned int face_no = 0) const; }; +namespace internal +{ + /** + * Utility function to convert "dofs per object" information + * of a @p dim dimensional reference cell @p cell_type. + */ + internal::GenericDoFsPerObject + expand(const unsigned int dim, + const std::vector &dofs_per_object, + const ReferenceCell::Type cell_type); +} // namespace internal + // --------- inline and template functions --------------- diff --git a/include/deal.II/simplex/fe_lib.h b/include/deal.II/simplex/fe_lib.h index aa833316b9..0d5fe5cc8e 100644 --- a/include/deal.II/simplex/fe_lib.h +++ b/include/deal.II/simplex/fe_lib.h @@ -169,6 +169,24 @@ namespace Simplex const FiniteElement &fe_other) const override; }; + /** + * Base class of FE_WedgeP and FE_WedgeDGP. + * + * @note Only implemented for 3D. + * + * @ingroup simplex + */ + template + class FE_Wedge : public dealii::FE_Poly + { + public: + /** + * Constructor. + */ + FE_Wedge(const unsigned int degree, + const internal::GenericDoFsPerObject & dpos, + const typename FiniteElementData::Conformity conformity); + }; /** * Implementation of a scalar Lagrange finite element on a wedge that yields @@ -178,7 +196,7 @@ namespace Simplex * @ingroup simplex */ template - class FE_WedgeP : public dealii::FE_Poly + class FE_WedgeP : public FE_Wedge { public: /** @@ -229,6 +247,55 @@ namespace Simplex const unsigned int face_no = 0) const override; }; + /** + * Implementation of a scalar Lagrange finite element on a wedge that yields + * the finite element space of discontinuous, piecewise polynomials of + * degree p. + * + * @ingroup simplex + */ + template + class FE_WedgeDGP : public FE_Wedge + { + public: + /** + * Constructor. + */ + FE_WedgeDGP(const unsigned int degree); + + /** + * @copydoc dealii::FiniteElement::clone() + */ + std::unique_ptr> + clone() const override; + + /** + * Return a string that uniquely identifies a finite element. This class + * returns Simplex::FE_WedgeDGP(degree), with @p dim and @p degree + * replaced by appropriate values. + */ + std::string + get_name() const override; + }; + + /** + * Base class of FE_PyramidP and FE_PyramidDGP. + * + * @note Only implemented for 3D. + * + * @ingroup simplex + */ + template + class FE_Pyramid : public dealii::FE_Poly + { + public: + /** + * Constructor. + */ + FE_Pyramid(const unsigned int degree, + const internal::GenericDoFsPerObject & dpos, + const typename FiniteElementData::Conformity conformity); + }; /** * Implementation of a scalar Lagrange finite element on a pyramid that yields @@ -238,7 +305,7 @@ namespace Simplex * @ingroup simplex */ template - class FE_PyramidP : public dealii::FE_Poly + class FE_PyramidP : public FE_Pyramid { public: /** @@ -289,6 +356,37 @@ namespace Simplex const unsigned int face_no = 0) const override; }; + /** + * Implementation of a scalar Lagrange finite element on a pyramid that yields + * the finite element space of discontinuous, piecewise polynomials of + * degree p. + * + * @ingroup simplex + */ + template + class FE_PyramidDGP : public FE_Pyramid + { + public: + /** + * Constructor. + */ + FE_PyramidDGP(const unsigned int degree); + + /** + * @copydoc dealii::FiniteElement::clone() + */ + std::unique_ptr> + clone() const override; + + /** + * Return a string that uniquely identifies a finite element. This class + * returns Simplex::FE_PyramidDGP(degree), with @p dim and @p degree + * replaced by appropriate values. + */ + std::string + get_name() const override; + }; + } // namespace Simplex diff --git a/source/fe/fe_data.cc b/source/fe/fe_data.cc index 497f398444..333baa7f4d 100644 --- a/source/fe/fe_data.cc +++ b/source/fe/fe_data.cc @@ -19,7 +19,7 @@ DEAL_II_NAMESPACE_OPEN -namespace +namespace internal { internal::GenericDoFsPerObject expand(const unsigned int dim, @@ -118,7 +118,7 @@ namespace return result; } -} // namespace +} // namespace internal template FiniteElementData::FiniteElementData( @@ -147,7 +147,7 @@ FiniteElementData::FiniteElementData( const unsigned int degree, const Conformity conformity, const BlockIndices & block_indices) - : FiniteElementData(expand(dim, dofs_per_object, cell_type), + : FiniteElementData(internal::expand(dim, dofs_per_object, cell_type), cell_type, n_components, degree, diff --git a/source/simplex/fe_lib.cc b/source/simplex/fe_lib.cc index 0b1831d2bd..91dbf2e260 100644 --- a/source/simplex/fe_lib.cc +++ b/source/simplex/fe_lib.cc @@ -85,11 +85,10 @@ namespace Simplex * Helper function to set up the dpo vector of FE_WedgeP for a given @p degree. */ internal::GenericDoFsPerObject - get_dpo_vector_fe_wedge(const unsigned int degree) + get_dpo_vector_fe_wedge_p(const unsigned int degree) { internal::GenericDoFsPerObject dpo; - // all dofs are internal if (degree == 1) { dpo.dofs_per_object_exclusive = {{1}, {0}, {0, 0, 0, 0, 0}, {0}}; @@ -117,14 +116,33 @@ namespace Simplex } /** - * Helper function to set up the dpo vector of FE_WedgeP for a given @p degree. + * Helper function to set up the dpo vector of FE_WedgeDGP for a given @p degree. + */ + internal::GenericDoFsPerObject + get_dpo_vector_fe_wedge_dgp(const unsigned int degree) + { + unsigned int n_dofs = 0; + + if (degree == 1) + n_dofs = 6; + else if (degree == 2) + n_dofs = 18; + else + Assert(false, ExcNotImplemented()); + + return internal::expand(3, + {{0, 0, 0, n_dofs}}, + ReferenceCell::Type::Wedge); + } + + /** + * Helper function to set up the dpo vector of FE_PyramidP for a given @p degree. */ internal::GenericDoFsPerObject - get_dpo_vector_fe_pyramid(const unsigned int degree) + get_dpo_vector_fe_pyramid_p(const unsigned int degree) { internal::GenericDoFsPerObject dpo; - // all dofs are internal if (degree == 1) { dpo.dofs_per_object_exclusive = {{1}, {0}, {0, 0, 0, 0, 0}, {0}}; @@ -141,6 +159,24 @@ namespace Simplex return dpo; } + + /** + * Helper function to set up the dpo vector of FE_PyramidDGP for a given @p degree. + */ + internal::GenericDoFsPerObject + get_dpo_vector_fe_pyramid_dgp(const unsigned int degree) + { + unsigned int n_dofs = 0; + + if (degree == 1) + n_dofs = 5; + else + Assert(false, ExcNotImplemented()); + + return internal::expand(3, + {{0, 0, 0, n_dofs}}, + ReferenceCell::Type::Pyramid); + } } // namespace @@ -451,32 +487,28 @@ namespace Simplex template - FE_WedgeP::FE_WedgeP(const unsigned int degree) + FE_Wedge::FE_Wedge( + const unsigned int degree, + const internal::GenericDoFsPerObject & dpos, + const typename FiniteElementData::Conformity conformity) : dealii::FE_Poly( Simplex::ScalarWedgePolynomial(degree), - FiniteElementData(get_dpo_vector_fe_wedge(degree), + FiniteElementData(dpos, ReferenceCell::Type::Wedge, 1, degree, - FiniteElementData::H1), - std::vector(FiniteElementData(get_dpo_vector_fe_wedge( - degree), - ReferenceCell::Type::Wedge, - 1, - degree) - .dofs_per_cell, - true), + conformity), + std::vector( + FiniteElementData(dpos, ReferenceCell::Type::Wedge, 1, degree) + .dofs_per_cell, + true), std::vector( - FiniteElementData(get_dpo_vector_fe_wedge(degree), - ReferenceCell::Type::Wedge, - 1, - degree) + FiniteElementData(dpos, ReferenceCell::Type::Wedge, 1, degree) .dofs_per_cell, std::vector(1, true))) { AssertDimension(dim, 3); - if (degree == 1) { this->unit_support_points.emplace_back(1.0, 0.0, 0.0); @@ -485,16 +517,20 @@ namespace Simplex this->unit_support_points.emplace_back(1.0, 0.0, 1.0); this->unit_support_points.emplace_back(0.0, 1.0, 1.0); this->unit_support_points.emplace_back(0.0, 0.0, 1.0); - - // TODO - this->unit_face_support_points[0].emplace_back(1.0, 0.0); - this->unit_face_support_points[0].emplace_back(0.0, 1.0); - this->unit_face_support_points[0].emplace_back(0.0, 0.0); } } + template + FE_WedgeP::FE_WedgeP(const unsigned int degree) + : FE_Wedge(degree, + get_dpo_vector_fe_wedge_p(degree), + FiniteElementData::H1) + {} + + + template std::unique_ptr> FE_WedgeP::clone() const @@ -579,7 +615,6 @@ namespace Simplex { (void)fe_other; - AssertIndexRange(face_no, 5); if (face_no < 2) @@ -604,26 +639,53 @@ namespace Simplex template - FE_PyramidP::FE_PyramidP(const unsigned int degree) + FE_WedgeDGP::FE_WedgeDGP(const unsigned int degree) + : FE_Wedge(degree, + get_dpo_vector_fe_wedge_dgp(degree), + FiniteElementData::L2) + {} + + + + template + std::unique_ptr> + FE_WedgeDGP::clone() const + { + return std::make_unique>(*this); + } + + + + template + std::string + FE_WedgeDGP::get_name() const + { + std::ostringstream namebuf; + namebuf << "FE_WedgeDGP<" << dim << ">(" << this->degree << ")"; + + return namebuf.str(); + } + + + + template + FE_Pyramid::FE_Pyramid( + const unsigned int degree, + const internal::GenericDoFsPerObject & dpos, + const typename FiniteElementData::Conformity conformity) : dealii::FE_Poly( Simplex::ScalarPyramidPolynomial(degree), - FiniteElementData(get_dpo_vector_fe_pyramid(degree), + FiniteElementData(dpos, ReferenceCell::Type::Pyramid, 1, degree, - FiniteElementData::H1), - std::vector(FiniteElementData(get_dpo_vector_fe_pyramid( - degree), - ReferenceCell::Type::Pyramid, - 1, - degree) - .dofs_per_cell, - true), + conformity), + std::vector( + FiniteElementData(dpos, ReferenceCell::Type::Pyramid, 1, degree) + .dofs_per_cell, + true), std::vector( - FiniteElementData(get_dpo_vector_fe_pyramid(degree), - ReferenceCell::Type::Pyramid, - 1, - degree) + FiniteElementData(dpos, ReferenceCell::Type::Pyramid, 1, degree) .dofs_per_cell, std::vector(1, true))) { @@ -642,6 +704,15 @@ namespace Simplex + template + FE_PyramidP::FE_PyramidP(const unsigned int degree) + : FE_Pyramid(degree, + get_dpo_vector_fe_pyramid_p(degree), + FiniteElementData::H1) + {} + + + template std::unique_ptr> FE_PyramidP::clone() const @@ -750,6 +821,36 @@ namespace Simplex + template + FE_PyramidDGP::FE_PyramidDGP(const unsigned int degree) + : FE_Pyramid(degree, + get_dpo_vector_fe_pyramid_dgp(degree), + FiniteElementData::L2) + {} + + + + template + std::unique_ptr> + FE_PyramidDGP::clone() const + { + return std::make_unique>(*this); + } + + + + template + std::string + FE_PyramidDGP::get_name() const + { + std::ostringstream namebuf; + namebuf << "FE_PyramidDGP<" << dim << ">(" << this->degree << ")"; + + return namebuf.str(); + } + + + } // namespace Simplex // explicit instantiations diff --git a/source/simplex/fe_lib.inst.in b/source/simplex/fe_lib.inst.in index 36b7279cd3..be8135b313 100644 --- a/source/simplex/fe_lib.inst.in +++ b/source/simplex/fe_lib.inst.in @@ -21,9 +21,17 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS) template class Simplex::FE_Poly; template class Simplex::FE_P; template class Simplex::FE_DGP; + template class Simplex::FE_Wedge; template class Simplex::FE_WedgeP; + template class Simplex::FE_WedgeDGP; + template class Simplex::FE_Pyramid; template class Simplex::FE_PyramidP; + template class Simplex::FE_PyramidDGP; #endif } diff --git a/tests/simplex/fe_lib_01.cc b/tests/simplex/fe_lib_01.cc index e8f01a75bc..104d4d60dd 100644 --- a/tests/simplex/fe_lib_01.cc +++ b/tests/simplex/fe_lib_01.cc @@ -29,11 +29,24 @@ test(const FiniteElement &fe) { deallog << fe.get_name() << ": " << std::endl; + const auto &reference_cell = + ReferenceCell::internal::Info::get_cell(fe.reference_cell_type()); + deallog << " n_dofs_per_vertex(): " << fe.n_dofs_per_vertex() << std::endl; deallog << " n_dofs_per_line(): " << fe.n_dofs_per_line() << std::endl; - deallog << " n_dofs_per_quad(): " << fe.n_dofs_per_quad() << std::endl; + + deallog << " n_dofs_per_quad(): "; + for (unsigned int i = 0; i < (dim == 2 ? 1 : reference_cell.n_faces()); ++i) + deallog << fe.n_dofs_per_quad(i) << " "; + deallog << std::endl; + deallog << " n_dofs_per_hex(): " << fe.n_dofs_per_hex() << std::endl; - deallog << " n_dofs_per_face(): " << fe.n_dofs_per_face() << std::endl; + + deallog << " n_dofs_per_face(): "; + for (unsigned int i = 0; i < reference_cell.n_faces(); ++i) + deallog << fe.n_dofs_per_face(i) << " "; + deallog << std::endl; + deallog << " n_dofs_per_cell(): " << fe.n_dofs_per_cell() << std::endl; deallog << " tensor_degree(): " << fe.tensor_degree() << std::endl; @@ -49,8 +62,19 @@ main() test(Simplex::FE_P<2>(2)); test(Simplex::FE_P<3>(1)); test(Simplex::FE_P<3>(2)); + test(Simplex::FE_DGP<2>(1)); test(Simplex::FE_DGP<2>(2)); test(Simplex::FE_DGP<3>(1)); test(Simplex::FE_DGP<3>(2)); + + test(Simplex::FE_WedgeP<3>(1)); + test(Simplex::FE_WedgeP<3>(2)); + + test(Simplex::FE_WedgeDGP<3>(1)); + test(Simplex::FE_WedgeDGP<3>(2)); + + test(Simplex::FE_PyramidP<3>(1)); + + test(Simplex::FE_PyramidDGP<3>(1)); } diff --git a/tests/simplex/fe_lib_01.output b/tests/simplex/fe_lib_01.output index 4410d62f69..ec3b0b86fd 100644 --- a/tests/simplex/fe_lib_01.output +++ b/tests/simplex/fe_lib_01.output @@ -2,72 +2,126 @@ DEAL::FE_P<2>(1): DEAL:: n_dofs_per_vertex(): 1 DEAL:: n_dofs_per_line(): 0 -DEAL:: n_dofs_per_quad(): 0 +DEAL:: n_dofs_per_quad(): 0 DEAL:: n_dofs_per_hex(): 0 -DEAL:: n_dofs_per_face(): 2 +DEAL:: n_dofs_per_face(): 2 2 2 DEAL:: n_dofs_per_cell(): 3 DEAL:: tensor_degree(): 1 DEAL:: DEAL::FE_P<2>(2): DEAL:: n_dofs_per_vertex(): 1 DEAL:: n_dofs_per_line(): 1 -DEAL:: n_dofs_per_quad(): 0 +DEAL:: n_dofs_per_quad(): 0 DEAL:: n_dofs_per_hex(): 0 -DEAL:: n_dofs_per_face(): 3 +DEAL:: n_dofs_per_face(): 3 3 3 DEAL:: n_dofs_per_cell(): 6 DEAL:: tensor_degree(): 2 DEAL:: DEAL::FE_P<3>(1): DEAL:: n_dofs_per_vertex(): 1 DEAL:: n_dofs_per_line(): 0 -DEAL:: n_dofs_per_quad(): 0 +DEAL:: n_dofs_per_quad(): 0 0 0 0 DEAL:: n_dofs_per_hex(): 0 -DEAL:: n_dofs_per_face(): 3 +DEAL:: n_dofs_per_face(): 3 3 3 3 DEAL:: n_dofs_per_cell(): 4 DEAL:: tensor_degree(): 1 DEAL:: DEAL::FE_P<3>(2): DEAL:: n_dofs_per_vertex(): 1 DEAL:: n_dofs_per_line(): 1 -DEAL:: n_dofs_per_quad(): 0 +DEAL:: n_dofs_per_quad(): 0 0 0 0 DEAL:: n_dofs_per_hex(): 0 -DEAL:: n_dofs_per_face(): 6 +DEAL:: n_dofs_per_face(): 6 6 6 6 DEAL:: n_dofs_per_cell(): 10 DEAL:: tensor_degree(): 2 DEAL:: DEAL::FE_DGP<2>(1): DEAL:: n_dofs_per_vertex(): 0 DEAL:: n_dofs_per_line(): 0 -DEAL:: n_dofs_per_quad(): 3 +DEAL:: n_dofs_per_quad(): 3 DEAL:: n_dofs_per_hex(): 0 -DEAL:: n_dofs_per_face(): 0 +DEAL:: n_dofs_per_face(): 0 0 0 DEAL:: n_dofs_per_cell(): 3 DEAL:: tensor_degree(): 1 DEAL:: DEAL::FE_DGP<2>(2): DEAL:: n_dofs_per_vertex(): 0 DEAL:: n_dofs_per_line(): 0 -DEAL:: n_dofs_per_quad(): 6 +DEAL:: n_dofs_per_quad(): 6 DEAL:: n_dofs_per_hex(): 0 -DEAL:: n_dofs_per_face(): 0 +DEAL:: n_dofs_per_face(): 0 0 0 DEAL:: n_dofs_per_cell(): 6 DEAL:: tensor_degree(): 2 DEAL:: DEAL::FE_DGP<3>(1): DEAL:: n_dofs_per_vertex(): 0 DEAL:: n_dofs_per_line(): 0 -DEAL:: n_dofs_per_quad(): 0 +DEAL:: n_dofs_per_quad(): 0 0 0 0 DEAL:: n_dofs_per_hex(): 4 -DEAL:: n_dofs_per_face(): 0 +DEAL:: n_dofs_per_face(): 0 0 0 0 DEAL:: n_dofs_per_cell(): 4 DEAL:: tensor_degree(): 1 DEAL:: DEAL::FE_DGP<3>(2): DEAL:: n_dofs_per_vertex(): 0 DEAL:: n_dofs_per_line(): 0 -DEAL:: n_dofs_per_quad(): 0 +DEAL:: n_dofs_per_quad(): 0 0 0 0 DEAL:: n_dofs_per_hex(): 10 -DEAL:: n_dofs_per_face(): 0 +DEAL:: n_dofs_per_face(): 0 0 0 0 DEAL:: n_dofs_per_cell(): 10 DEAL:: tensor_degree(): 2 DEAL:: +DEAL::FE_WedgeP<3>(1): +DEAL:: n_dofs_per_vertex(): 1 +DEAL:: n_dofs_per_line(): 0 +DEAL:: n_dofs_per_quad(): 0 0 0 0 0 +DEAL:: n_dofs_per_hex(): 0 +DEAL:: n_dofs_per_face(): 3 3 4 4 4 +DEAL:: n_dofs_per_cell(): 6 +DEAL:: tensor_degree(): 1 +DEAL:: +DEAL::FE_WedgeP<3>(2): +DEAL:: n_dofs_per_vertex(): 1 +DEAL:: n_dofs_per_line(): 1 +DEAL:: n_dofs_per_quad(): 0 0 1 1 1 +DEAL:: n_dofs_per_hex(): 0 +DEAL:: n_dofs_per_face(): 6 6 9 9 9 +DEAL:: n_dofs_per_cell(): 18 +DEAL:: tensor_degree(): 2 +DEAL:: +DEAL::FE_WedgeDGP<3>(1): +DEAL:: n_dofs_per_vertex(): 0 +DEAL:: n_dofs_per_line(): 0 +DEAL:: n_dofs_per_quad(): 0 0 0 0 0 +DEAL:: n_dofs_per_hex(): 6 +DEAL:: n_dofs_per_face(): 0 0 0 0 0 +DEAL:: n_dofs_per_cell(): 6 +DEAL:: tensor_degree(): 1 +DEAL:: +DEAL::FE_WedgeDGP<3>(2): +DEAL:: n_dofs_per_vertex(): 0 +DEAL:: n_dofs_per_line(): 0 +DEAL:: n_dofs_per_quad(): 0 0 0 0 0 +DEAL:: n_dofs_per_hex(): 18 +DEAL:: n_dofs_per_face(): 0 0 0 0 0 +DEAL:: n_dofs_per_cell(): 18 +DEAL:: tensor_degree(): 2 +DEAL:: +DEAL::FE_PyramidP<3>(1): +DEAL:: n_dofs_per_vertex(): 1 +DEAL:: n_dofs_per_line(): 0 +DEAL:: n_dofs_per_quad(): 0 0 0 0 0 +DEAL:: n_dofs_per_hex(): 0 +DEAL:: n_dofs_per_face(): 4 3 3 3 3 +DEAL:: n_dofs_per_cell(): 5 +DEAL:: tensor_degree(): 1 +DEAL:: +DEAL::FE_PyramidDGP<3>(1): +DEAL:: n_dofs_per_vertex(): 0 +DEAL:: n_dofs_per_line(): 0 +DEAL:: n_dofs_per_quad(): 0 0 0 0 0 +DEAL:: n_dofs_per_hex(): 5 +DEAL:: n_dofs_per_face(): 0 0 0 0 0 +DEAL:: n_dofs_per_cell(): 5 +DEAL:: tensor_degree(): 1 +DEAL:: -- 2.39.5