#include <deal.II/base/config.h>
+#include <deal.II/base/ndarray.h>
#include <deal.II/base/polynomials_barycentric.h>
#include <deal.II/base/scalar_polynomials_base.h>
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<unsigned int, 6, 2> 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<unsigned int, 18, 2> 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.
// ---------------------------------------------------------------------
-#include <deal.II/base/ndarray.h>
#include <deal.II/base/polynomials_barycentric.h>
#include <deal.II/base/polynomials_wedge.h>
{}
-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<unsigned int, 6, 2> 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<unsigned int, 18, 2> 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 <int dim>
double
ScalarLagrangePolynomialWedge<dim>::compute_value(const unsigned int i,
const Point<dim> & 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);
ScalarLagrangePolynomialWedge<dim>::compute_grad(const unsigned int i,
const Point<dim> & 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);
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<dim>(i));
+ this->reference_cell().template vertex<dim>(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<dim - 1>(i));
+ }
}
else
Assert(false, ExcNotImplemented());
{
AssertDimension(dim, 3);
- if (degree == 1)
- {
- for (const unsigned int i : ReferenceCells::Wedge.vertex_indices())
- this->unit_support_points.emplace_back(
- ReferenceCells::Wedge.vertex<dim>(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());
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_wedge_p.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim, int spacedim>
+void
+test(const FiniteElement<dim, spacedim> &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<dim - 1> fe_q(fe.degree);
+ FE_SimplexP<dim - 1> 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<FiniteElement<dim - 1> &>(fe_q) :
+ static_cast<FiniteElement<dim - 1> &>(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));
+}
--- /dev/null
+
+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!
projects / dealii.git / commitdiff