From 729fd212b26cfc37edc7dc462b826471dec0ace4 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 8 May 2022 17:13:39 +0200 Subject: [PATCH] unit_support_points and unit_face_support_points for wedges and pyramids --- include/deal.II/base/polynomials_wedge.h | 39 ++++++++ source/base/polynomials_wedge.cc | 43 +-------- source/fe/fe_pyramid_p.cc | 16 +++- source/fe/fe_wedge_p.cc | 39 +++++--- ...rt_points.cc => unit_support_points_01.cc} | 0 ...s.output => unit_support_points_01.output} | 0 tests/bits/unit_support_points_03.cc | 96 +++++++++++++++++++ tests/bits/unit_support_points_03.output | 23 +++++ 8 files changed, 204 insertions(+), 52 deletions(-) rename tests/bits/{unit_support_points.cc => unit_support_points_01.cc} (100%) rename tests/bits/{unit_support_points.output => unit_support_points_01.output} (100%) create mode 100644 tests/bits/unit_support_points_03.cc create mode 100644 tests/bits/unit_support_points_03.output diff --git a/include/deal.II/base/polynomials_wedge.h b/include/deal.II/base/polynomials_wedge.h index c29410fc0a..dfd86334f7 100644 --- a/include/deal.II/base/polynomials_wedge.h +++ b/include/deal.II/base/polynomials_wedge.h @@ -19,11 +19,50 @@ #include +#include #include #include DEAL_II_NAMESPACE_OPEN + +namespace internal +{ + /** + * Decompose the shape-function index of a linear wedge into an index + * to access the right shape function within the triangle and and within + * the line. + */ + static const constexpr dealii::ndarray wedge_table_1{ + {{{0, 0}}, {{1, 0}}, {{2, 0}}, {{0, 1}}, {{1, 1}}, {{2, 1}}}}; + + /** + * Decompose the shape-function index of a quadratic wedge into an index + * to access the right shape function within the triangle and and within + * the line. + */ + static const constexpr dealii::ndarray wedge_table_2{ + {{{0, 0}}, + {{1, 0}}, + {{2, 0}}, + {{0, 1}}, + {{1, 1}}, + {{2, 1}}, + {{3, 0}}, + {{4, 0}}, + {{5, 0}}, + {{3, 1}}, + {{4, 1}}, + {{5, 1}}, + {{0, 2}}, + {{1, 2}}, + {{2, 2}}, + {{3, 2}}, + {{4, 2}}, + {{5, 2}}}}; +} // namespace internal + + /** * Polynomials defined on wedge entities. This class is basis of * FE_WedgeP. diff --git a/source/base/polynomials_wedge.cc b/source/base/polynomials_wedge.cc index 402fbf954e..186b228f15 100644 --- a/source/base/polynomials_wedge.cc +++ b/source/base/polynomials_wedge.cc @@ -14,7 +14,6 @@ // --------------------------------------------------------------------- -#include #include #include @@ -50,49 +49,14 @@ ScalarLagrangePolynomialWedge::ScalarLagrangePolynomialWedge( {} -namespace -{ - /** - * Decompose the shape-function index of a linear wedge into an index - * to access the right shape function within the triangle and and within - * the line. - */ - static const constexpr ndarray wedge_table_1{ - {{{0, 0}}, {{1, 0}}, {{2, 0}}, {{0, 1}}, {{1, 1}}, {{2, 1}}}}; - - /** - * Decompose the shape-function index of a quadratic wedge into an index - * to access the right shape function within the triangle and and within - * the line. - */ - static const constexpr ndarray wedge_table_2{{{{0, 0}}, - {{1, 0}}, - {{2, 0}}, - {{0, 1}}, - {{1, 1}}, - {{2, 1}}, - {{3, 0}}, - {{4, 0}}, - {{5, 0}}, - {{3, 1}}, - {{4, 1}}, - {{5, 1}}, - {{0, 2}}, - {{1, 2}}, - {{2, 2}}, - {{3, 2}}, - {{4, 2}}, - {{5, 2}}}}; -} // namespace - - template double ScalarLagrangePolynomialWedge::compute_value(const unsigned int i, const Point & p) const { - const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i]; + const auto pair = this->degree() == 1 ? internal::wedge_table_1[i] : + internal::wedge_table_2[i]; const Point<2> p_tri(p[0], p[1]); const auto v_tri = poly_tri.compute_value(pair[0], p_tri); @@ -110,7 +74,8 @@ Tensor<1, dim> ScalarLagrangePolynomialWedge::compute_grad(const unsigned int i, const Point & p) const { - const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i]; + const auto pair = this->degree() == 1 ? internal::wedge_table_1[i] : + internal::wedge_table_2[i]; const Point<2> p_tri(p[0], p[1]); const auto v_tri = poly_tri.compute_value(pair[0], p_tri); diff --git a/source/fe/fe_pyramid_p.cc b/source/fe/fe_pyramid_p.cc index 60bd0d0a03..2ef5822bfb 100644 --- a/source/fe/fe_pyramid_p.cc +++ b/source/fe/fe_pyramid_p.cc @@ -97,9 +97,21 @@ FE_PyramidPoly::FE_PyramidPoly( if (degree == 1) { - for (const unsigned int i : ReferenceCells::Pyramid.vertex_indices()) + for (const auto i : this->reference_cell().vertex_indices()) this->unit_support_points.emplace_back( - ReferenceCells::Pyramid.vertex(i)); + this->reference_cell().template vertex(i)); + + this->unit_face_support_points.resize(this->reference_cell().n_faces()); + + for (const auto f : this->reference_cell().face_indices()) + { + const auto face_reference_cell = + this->reference_cell().face_reference_cell(f); + + for (const auto i : face_reference_cell.vertex_indices()) + this->unit_face_support_points[f].emplace_back( + face_reference_cell.template vertex(i)); + } } else Assert(false, ExcNotImplemented()); diff --git a/source/fe/fe_wedge_p.cc b/source/fe/fe_wedge_p.cc index 42e6c1dcec..263699c606 100644 --- a/source/fe/fe_wedge_p.cc +++ b/source/fe/fe_wedge_p.cc @@ -105,19 +105,36 @@ FE_WedgePoly::FE_WedgePoly( { AssertDimension(dim, 3); - if (degree == 1) - { - for (const unsigned int i : ReferenceCells::Wedge.vertex_indices()) - this->unit_support_points.emplace_back( - ReferenceCells::Wedge.vertex(i)); - } - else + Assert(1 <= degree && degree <= 2, ExcNotImplemented()); + + FE_SimplexP<2> fe_triangle(degree); + FE_Q<1> fe_line(degree); + FE_Q<2> fe_quad(degree); + + for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i) { - // TODO: Wedge elements work for higher degrees, but we don't currently - // fill their support points. Leaving the array empty is valid, however, - // and will simply result in an error when someone tries to access the - // array. + const auto pair = this->degree == 1 ? internal::wedge_table_1[i] : + internal::wedge_table_2[i]; + + this->unit_support_points.emplace_back( + fe_triangle.get_unit_support_points()[pair[0]][0], + fe_triangle.get_unit_support_points()[pair[0]][1], + fe_line.get_unit_support_points()[pair[1]][0]); } + + this->unit_face_support_points.resize(this->reference_cell().n_faces()); + + for (const auto f : this->reference_cell().face_indices()) + if (this->reference_cell().face_reference_cell(f) == + ReferenceCells::Triangle) + for (const auto &p : fe_triangle.get_unit_support_points()) + this->unit_face_support_points[f].emplace_back(p[0], p[1]); + else if (this->reference_cell().face_reference_cell(f) == + ReferenceCells::Quadrilateral) + for (const auto &p : fe_quad.get_unit_support_points()) + this->unit_face_support_points[f].emplace_back(p[0], p[1]); + else + Assert(false, ExcInternalError()); } diff --git a/tests/bits/unit_support_points.cc b/tests/bits/unit_support_points_01.cc similarity index 100% rename from tests/bits/unit_support_points.cc rename to tests/bits/unit_support_points_01.cc diff --git a/tests/bits/unit_support_points.output b/tests/bits/unit_support_points_01.output similarity index 100% rename from tests/bits/unit_support_points.output rename to tests/bits/unit_support_points_01.output diff --git a/tests/bits/unit_support_points_03.cc b/tests/bits/unit_support_points_03.cc new file mode 100644 index 0000000000..b630f5fbfc --- /dev/null +++ b/tests/bits/unit_support_points_03.cc @@ -0,0 +1,96 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Test FiniteElement::unit_support_points() and ::unit_face_support_points() +// for different shapes. + +#include +#include +#include +#include + +#include "../tests.h" + +using namespace dealii; + +template +void +test(const FiniteElement &fe) +{ + deallog << fe.get_name() << std::endl; + + const auto &unit_support_points = fe.get_unit_support_points(); + + AssertDimension(unit_support_points.size(), fe.n_dofs_per_cell()); + + for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) + for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j) + { + const auto value = fe.shape_value(i, unit_support_points[j]); + Assert((i == j ? (std::abs(value - 1.0) < 1e-8) : + (std::abs(value) < 1e-8)), + ExcInternalError()); + } + + FE_Q fe_q(fe.degree); + FE_SimplexP fe_p(fe.degree); + + for (const auto f : fe.reference_cell().face_indices()) + { + const auto &unit_face_support_points = fe.get_unit_face_support_points(f); + + AssertDimension(unit_face_support_points.size(), fe.n_dofs_per_face(f)); + + const auto &fe_face = + fe.reference_cell().face_reference_cell(f).is_hyper_cube() ? + static_cast &>(fe_q) : + static_cast &>(fe_p); + + for (unsigned int i = 0; i < fe.n_dofs_per_face(f); ++i) + for (unsigned int j = 0; j < fe.n_dofs_per_face(f); ++j) + { + const auto value = + fe_face.shape_value(i, unit_face_support_points[j]); + Assert((i == j ? (std::abs(value - 1.0) < 1e-8) : + (std::abs(value) < 1e-8)), + ExcInternalError()); + } + } + + + + deallog << "OK!" << std::endl; +} + +int +main() +{ + initlog(); + + test(FE_Q<2>(1)); + test(FE_Q<2>(2)); + test(FE_Q<3>(1)); + test(FE_Q<3>(2)); + + test(FE_SimplexP<2>(1)); + test(FE_SimplexP<2>(2)); + test(FE_SimplexP<3>(1)); + test(FE_SimplexP<3>(2)); + + test(FE_PyramidP<3>(1)); + + test(FE_WedgeP<3>(1)); + test(FE_WedgeP<3>(2)); +} diff --git a/tests/bits/unit_support_points_03.output b/tests/bits/unit_support_points_03.output new file mode 100644 index 0000000000..8c41550fd3 --- /dev/null +++ b/tests/bits/unit_support_points_03.output @@ -0,0 +1,23 @@ + +DEAL::FE_Q<2>(1) +DEAL::OK! +DEAL::FE_Q<2>(2) +DEAL::OK! +DEAL::FE_Q<3>(1) +DEAL::OK! +DEAL::FE_Q<3>(2) +DEAL::OK! +DEAL::FE_SimplexP<2>(1) +DEAL::OK! +DEAL::FE_SimplexP<2>(2) +DEAL::OK! +DEAL::FE_SimplexP<3>(1) +DEAL::OK! +DEAL::FE_SimplexP<3>(2) +DEAL::OK! +DEAL::FE_PyramidP<3>(1) +DEAL::OK! +DEAL::FE_WedgeP<3>(1) +DEAL::OK! +DEAL::FE_WedgeP<3>(2) +DEAL::OK! -- 2.39.5