From 55e8ca95a7e2ef0bd588598a9a430adfd8a5c057 Mon Sep 17 00:00:00 2001 From: David Wells Date: Wed, 13 Jan 2021 13:04:31 -0500 Subject: [PATCH] Add a utility function for Simplex::FE_P's support points. While we are here, fix a bug where unit face support points were assigned to a discontinous element, which doesn't make sense (see FiniteElement::get_unit_face_support_points()). --- source/simplex/fe_lib.cc | 183 +++++++++++++++++++++++---------------- 1 file changed, 106 insertions(+), 77 deletions(-) diff --git a/source/simplex/fe_lib.cc b/source/simplex/fe_lib.cc index ca3141b1cd..78bfcdd0f0 100644 --- a/source/simplex/fe_lib.cc +++ b/source/simplex/fe_lib.cc @@ -56,6 +56,107 @@ namespace Simplex return dpo; } + /** + * Set up a vector that contains the unit (reference) cell support points + * for FE_Poly and sufficiently similar elements. + */ + template + std::vector> + unit_support_points_fe_poly(const unsigned int degree) + { + std::vector> unit_points; + if (dim == 1) + { + // We don't really have dim = 1 support for simplex elements yet, but + // its convenient for populating the face array + if (degree >= 1) + { + unit_points.emplace_back(0.0); + unit_points.emplace_back(1.0); + } + if (degree == 2) + { + unit_points.emplace_back(0.5); + } + if (degree > 2) + { + Assert(false, ExcNotImplemented()); + } + } + else if (dim == 2) + { + if (degree >= 1) + { + unit_points.emplace_back(0.0, 0.0); + unit_points.emplace_back(1.0, 0.0); + unit_points.emplace_back(0.0, 1.0); + } + if (degree == 2) + { + unit_points.emplace_back(0.5, 0.0); + unit_points.emplace_back(0.5, 0.5); + unit_points.emplace_back(0.0, 0.5); + } + if (degree > 3) + { + Assert(false, ExcNotImplemented()); + } + } + else if (dim == 3) + { + if (degree >= 1) + { + unit_points.emplace_back(0.0, 0.0, 0.0); + unit_points.emplace_back(1.0, 0.0, 0.0); + unit_points.emplace_back(0.0, 1.0, 0.0); + unit_points.emplace_back(0.0, 0.0, 1.0); + } + if (degree == 2) + { + unit_points.emplace_back(0.5, 0.0, 0.0); + unit_points.emplace_back(0.5, 0.5, 0.0); + unit_points.emplace_back(0.0, 0.5, 0.0); + unit_points.emplace_back(0.0, 0.0, 0.5); + unit_points.emplace_back(0.5, 0.0, 0.5); + unit_points.emplace_back(0.0, 0.5, 0.5); + } + if (degree == 3) + { + Assert(false, ExcNotImplemented()); + } + } + else + { + Assert(false, ExcNotImplemented()); + } + + return unit_points; + } + + /** + * Set up a vector that contains the unit (reference) cell's faces support + * points for FE_Poly and sufficiently similar elements. + */ + template + std::vector>> + unit_face_support_points_fe_poly(const unsigned int degree) + { + Assert(dim == 2 || dim == 3, ExcNotImplemented()); + const auto &info = ReferenceCell::internal::Info::get_cell( + dim == 2 ? ReferenceCell::Type::Tri : ReferenceCell::Type::Tet); + std::vector>> unit_face_points; + + // all faces have the same support points + for (auto face_n : info.face_indices()) + { + (void)face_n; + unit_face_points.emplace_back( + unit_support_points_fe_poly(degree)); + } + + return unit_face_points; + } + /** * Helper function to set up the dpo vector of FE_DGP for a given @p dim and * @p degree. @@ -212,83 +313,11 @@ namespace Simplex .dofs_per_cell, std::vector(1, true))) { - this->unit_support_points.clear(); - - if (dim == 2) - { - if (degree == 1) - { - this->unit_support_points.emplace_back(0.0, 0.0); - this->unit_support_points.emplace_back(1.0, 0.0); - this->unit_support_points.emplace_back(0.0, 1.0); - - // TODO - this->unit_face_support_points[0].emplace_back(0.0); - this->unit_face_support_points[0].emplace_back(1.0); - } - else if (degree == 2) - { - this->unit_support_points.emplace_back(0.0, 0.0); - this->unit_support_points.emplace_back(1.0, 0.0); - this->unit_support_points.emplace_back(0.0, 1.0); - this->unit_support_points.emplace_back(0.5, 0.0); - this->unit_support_points.emplace_back(0.5, 0.5); - this->unit_support_points.emplace_back(0.0, 0.5); - - // TODO - this->unit_face_support_points[0].emplace_back(0.0); - this->unit_face_support_points[0].emplace_back(1.0); - this->unit_face_support_points[0].emplace_back(0.5); - } - else - { - Assert(false, ExcNotImplemented()); - } - } - else if (dim == 3) - { - if (degree == 1) - { - this->unit_support_points.emplace_back(0.0, 0.0, 0.0); - this->unit_support_points.emplace_back(1.0, 0.0, 0.0); - this->unit_support_points.emplace_back(0.0, 1.0, 0.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); - } - else if (degree == 2) - { - this->unit_support_points.emplace_back(0.0, 0.0, 0.0); - this->unit_support_points.emplace_back(1.0, 0.0, 0.0); - this->unit_support_points.emplace_back(0.0, 1.0, 0.0); - this->unit_support_points.emplace_back(0.0, 0.0, 1.0); - this->unit_support_points.emplace_back(0.5, 0.0, 0.0); - this->unit_support_points.emplace_back(0.5, 0.5, 0.0); - this->unit_support_points.emplace_back(0.0, 0.5, 0.0); - this->unit_support_points.emplace_back(0.0, 0.0, 0.5); - this->unit_support_points.emplace_back(0.5, 0.0, 0.5); - this->unit_support_points.emplace_back(0.0, 0.5, 0.5); - - // 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); - this->unit_face_support_points[0].emplace_back(0.5, 0.5); - this->unit_face_support_points[0].emplace_back(0.0, 0.5); - this->unit_face_support_points[0].emplace_back(0.5, 0.0); - } - else - { - Assert(false, ExcNotImplemented()); - } - } - else - { - Assert(false, ExcNotImplemented()); - } + this->unit_support_points = unit_support_points_fe_poly(degree); + // Discontinuous elements don't have face support points + if (conformity == FiniteElementData::Conformity::H1) + this->unit_face_support_points = + unit_face_support_points_fe_poly(degree); } -- 2.39.5