From: Peter Munch Date: Sat, 28 Nov 2020 08:17:39 +0000 (+0100) Subject: Introduce FE_Wedge X-Git-Tag: v9.3.0-rc1~819^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=943eb56101a55bdbbf878a02c8e552fd4833c02e;p=dealii.git Introduce FE_Wedge --- diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h index 32d1c49b72..5402c73267 100644 --- a/include/deal.II/grid/reference_cell.h +++ b/include/deal.II/grid/reference_cell.h @@ -822,8 +822,6 @@ namespace ReferenceCell standard_line_to_face_and_line_index( const unsigned int line) const override { - Assert(false, ExcNotImplemented()); - static const std::array table[6] = { {{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 1}}, {{1, 2}}, {{2, 1}}}; @@ -986,10 +984,15 @@ namespace ReferenceCell standard_line_to_face_and_line_index( const unsigned int line) const override { - Assert(false, ExcNotImplemented()); - - static const std::array table[6] = { - {{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 1}}, {{1, 2}}, {{2, 1}}}; + static const std::array table[9] = {{{0, 0}}, + {{0, 2}}, + {{0, 1}}, + {{1, 0}}, + {{1, 1}}, + {{1, 2}}, + {{2, 0}}, + {{2, 1}}, + {{3, 1}}}; return table[line]; } @@ -1000,19 +1003,26 @@ namespace ReferenceCell const unsigned int face, const unsigned char face_orientation) const override { - Assert(false, ExcNotImplemented()); - - (void)face; - - static const std::array, 6> table = { - {{{2, 1, 0}}, - {{0, 1, 2}}, - {{1, 2, 0}}, - {{0, 2, 1}}, - {{1, 0, 2}}, - {{2, 0, 1}}}}; + if (face > 1) // QUAD + { + return GeometryInfo<3>::standard_to_real_face_line( + line, + get_bit(face_orientation, 0), + get_bit(face_orientation, 2), + get_bit(face_orientation, 1)); + } + else // TRI + { + static const std::array, 6> table = { + {{{2, 1, 0}}, + {{0, 1, 2}}, + {{1, 2, 0}}, + {{0, 2, 1}}, + {{1, 0, 2}}, + {{2, 0, 1}}}}; - return table[face_orientation][line]; + return table[face_orientation][line]; + } } bool diff --git a/include/deal.II/simplex/fe_lib.h b/include/deal.II/simplex/fe_lib.h index f618ff8e5b..e7c45a8ed9 100644 --- a/include/deal.II/simplex/fe_lib.h +++ b/include/deal.II/simplex/fe_lib.h @@ -161,6 +161,67 @@ namespace Simplex const FiniteElement &fe_other) const override; }; + + /** + * Implementation of a scalar Lagrange finite element on a wedge that yields + * the finite element space of continuous, piecewise polynomials of + * degree p. + * + * @ingroup simplex + */ + template + class FE_WedgeP : public dealii::FE_Poly + { + public: + /** + * Constructor. + */ + FE_WedgeP(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_WedgeP(degree), with @p dim and @p degree + * replaced by appropriate values. + */ + std::string + get_name() const override; + + /** + * @copydoc dealii::FiniteElement::compare_for_domination() + */ + FiniteElementDomination::Domination + compare_for_domination(const FiniteElement &fe_other, + const unsigned int codim) const override; + + /** + * @copydoc dealii::FiniteElement::hp_vertex_dof_identities() + */ + std::vector> + hp_vertex_dof_identities( + const FiniteElement &fe_other) const override; + + /** + * @copydoc dealii::FiniteElement::hp_line_dof_identities() + */ + std::vector> + hp_line_dof_identities( + const FiniteElement &fe_other) const override; + + /** + * @copydoc dealii::FiniteElement::hp_quad_dof_identities() + */ + std::vector> + hp_quad_dof_identities(const FiniteElement &fe_other, + const unsigned int face_no = 0) const override; + }; + + } // namespace Simplex DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/simplex/polynomials.h b/include/deal.II/simplex/polynomials.h index 7561d28a3a..a1c484aaa5 100644 --- a/include/deal.II/simplex/polynomials.h +++ b/include/deal.II/simplex/polynomials.h @@ -145,6 +145,127 @@ namespace Simplex clone() const override; }; + + + /** + * Polynomials defined on wedge entities. This class is basis of + * Simplex::FE_WedgeP. + * + * The polynomials are created via a tensor product of a + * Simplex::ScalarPolynomial<2>(degree) and a + * Simplex::ScalarPolynomial<1>(degree), however, are re-numerated to better + * match the definition of FiniteElement. + */ + template + class ScalarWedgePolynomial : public ScalarPolynomialsBase + { + public: + /** + * Make the dimension available to the outside. + */ + static const unsigned int dimension = dim; + + /* + * Constructor taking the polynomial @p degree as input. + * + * @note Currently, only linear (degree=1) and quadratic polynomials + * (degree=2) are implemented. + */ + ScalarWedgePolynomial(const unsigned int degree); + + /** + * @copydoc ScalarPolynomialsBase::evaluate() + * + * @note Currently, only the vectors @p values and @p grads are filled. + */ + void + evaluate(const Point & unit_point, + std::vector & values, + std::vector> &grads, + std::vector> &grad_grads, + std::vector> &third_derivatives, + std::vector> &fourth_derivatives) const override; + + /** + * @copydoc ScalarPolynomialsBase::compute_value() + */ + double + compute_value(const unsigned int i, const Point &p) const override; + + /** + * @copydoc ScalarPolynomialsBase::compute_derivative() + * + * @note Currently, only implemented for first derivative. + */ + template + Tensor + compute_derivative(const unsigned int i, const Point &p) const; + + /** + * @copydoc ScalarPolynomialsBase::compute_1st_derivative() + */ + Tensor<1, dim> + compute_1st_derivative(const unsigned int i, + const Point & p) const override; + + /** + * @copydoc ScalarPolynomialsBase::compute_2nd_derivative() + * + * @note Not implemented yet. + */ + Tensor<2, dim> + compute_2nd_derivative(const unsigned int i, + const Point & p) const override; + + /** + * @copydoc ScalarPolynomialsBase::compute_3rd_derivative() + * + * @note Not implemented yet. + */ + Tensor<3, dim> + compute_3rd_derivative(const unsigned int i, + const Point & p) const override; + + /** + * @copydoc ScalarPolynomialsBase::compute_4th_derivative() + * + * @note Not implemented yet. + */ + Tensor<4, dim> + compute_4th_derivative(const unsigned int i, + const Point & p) const override; + + /** + * @copydoc ScalarPolynomialsBase::compute_grad() + * + * @note Not implemented yet. + */ + Tensor<1, dim> + compute_grad(const unsigned int i, const Point &p) const override; + + /** + * @copydoc ScalarPolynomialsBase::compute_grad_grad() + * + * @note Not implemented yet. + */ + Tensor<2, dim> + compute_grad_grad(const unsigned int i, const Point &p) const override; + + /** + * @copydoc ScalarPolynomialsBase::name() + */ + std::string + name() const override; + + /** + * @copydoc ScalarPolynomialsBase::clone() + */ + virtual std::unique_ptr> + clone() const override; + }; + + + // template functions template template @@ -162,6 +283,26 @@ namespace Simplex return der; } + + + + template + template + Tensor + ScalarWedgePolynomial::compute_derivative(const unsigned int i, + const Point & p) const + { + Tensor der; + + AssertDimension(order, 1); + const auto grad = compute_grad(i, p); + + for (unsigned int i = 0; i < dim; i++) + der[i] = grad[i]; + + return der; + } + } // namespace Simplex DEAL_II_NAMESPACE_CLOSE diff --git a/source/simplex/fe_lib.cc b/source/simplex/fe_lib.cc index a034abd234..abd6e59042 100644 --- a/source/simplex/fe_lib.cc +++ b/source/simplex/fe_lib.cc @@ -80,6 +80,41 @@ namespace Simplex return dpo; } + + /** + * 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) + { + internal::GenericDoFsPerObject dpo; + + // all dofs are internal + if (degree == 1) + { + dpo.dofs_per_object_exclusive = {{1}, {0}, {0, 0, 0, 0, 0}, {0}}; + dpo.dofs_per_object_inclusive = {{1}, {2}, {3, 3, 4, 4, 4}, {6}}; + dpo.object_index = {{}, {6}, {6}, {6}}; + dpo.first_object_index_on_face = {{}, + {3, 3, 4, 4, 4}, + {3, 3, 4, 4, 4}}; + } + else if (degree == 2) + { + dpo.dofs_per_object_exclusive = {{1}, {1}, {0, 0, 1, 1, 1}, {0}}; + dpo.dofs_per_object_inclusive = {{1}, {3}, {6, 6, 9, 9, 9}, {18}}; + dpo.object_index = {{}, {6}, {15, 15, 15, 16, 17}, {18}}; + dpo.first_object_index_on_face = {{}, + {3, 3, 4, 4, 4}, + {6, 6, 8, 8, 8}}; + } + else + { + Assert(false, ExcNotImplemented()); + } + + return dpo; + } } // namespace @@ -371,6 +406,159 @@ namespace Simplex + template + FE_WedgeP::FE_WedgeP(const unsigned int degree) + : dealii::FE_Poly( + Simplex::ScalarWedgePolynomial(degree), + FiniteElementData(get_dpo_vector_fe_wedge(degree), + ReferenceCell::Type::Wedge, + 1, + degree, + FiniteElementData::L2), + std::vector(FiniteElementData(get_dpo_vector_fe_wedge( + degree), + ReferenceCell::Type::Wedge, + 1, + degree) + .dofs_per_cell, + true), + std::vector( + FiniteElementData(get_dpo_vector_fe_wedge(degree), + 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); + this->unit_support_points.emplace_back(0.0, 1.0, 0.0); + this->unit_support_points.emplace_back(0.0, 0.0, 0.0); + 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 + std::unique_ptr> + FE_WedgeP::clone() const + { + return std::make_unique>(*this); + } + + + + template + std::string + FE_WedgeP::get_name() const + { + std::ostringstream namebuf; + namebuf << "FE_WedgeP<" << dim << ">(" << this->degree << ")"; + + return namebuf.str(); + } + + + + template + FiniteElementDomination::Domination + FE_WedgeP::compare_for_domination( + const FiniteElement &fe_other, + const unsigned int codim) const + { + (void)fe_other; + (void)codim; + + Assert((dynamic_cast *>(&fe_other)) || + (dynamic_cast *>(&fe_other)), + ExcNotImplemented()); + + return FiniteElementDomination::this_element_dominates; + } + + + + template + std::vector> + FE_WedgeP::hp_vertex_dof_identities( + const FiniteElement &fe_other) const + { + (void)fe_other; + + Assert((dynamic_cast *>(&fe_other)) || + (dynamic_cast *>(&fe_other)), + ExcNotImplemented()); + + return {{0, 0}}; + } + + + + template + std::vector> + FE_WedgeP::hp_line_dof_identities( + const FiniteElement &fe_other) const + { + (void)fe_other; + + Assert((dynamic_cast *>(&fe_other)) || + (dynamic_cast *>(&fe_other)), + ExcNotImplemented()); + + std::vector> result; + + for (unsigned int i = 0; i < this->degree - 1; ++i) + result.emplace_back(i, i); + + return result; + } + + + + template + std::vector> + FE_WedgeP::hp_quad_dof_identities( + const FiniteElement &fe_other, + const unsigned int face_no) const + { + (void)fe_other; + + + AssertIndexRange(face_no, 5); + + if (face_no < 2) + { + Assert((dynamic_cast *>(&fe_other)), + ExcNotImplemented()); + } + else + { + Assert((dynamic_cast *>(&fe_other)), + ExcNotImplemented()); + } + + std::vector> result; + + for (unsigned int i = 0; i < this->n_dofs_per_quad(face_no); ++i) + result.emplace_back(i, i); + + return result; + } + + + } // namespace Simplex // explicit instantiations diff --git a/source/simplex/fe_lib.inst.in b/source/simplex/fe_lib.inst.in index 0f29f676d1..160d876f95 100644 --- a/source/simplex/fe_lib.inst.in +++ b/source/simplex/fe_lib.inst.in @@ -21,5 +21,7 @@ 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_WedgeP; #endif } diff --git a/source/simplex/polynomials.cc b/source/simplex/polynomials.cc index 5c33c826dc..bae98e5f95 100644 --- a/source/simplex/polynomials.cc +++ b/source/simplex/polynomials.cc @@ -51,6 +51,23 @@ namespace Simplex return 0; } + + unsigned int + compute_n_polynomials_wedge(const unsigned int dim, + const unsigned int degree) + { + if (dim == 3) + { + if (degree == 1) + return 6; + if (degree == 2) + return 18; + } + + Assert(false, ExcNotImplemented()); + + return 0; + } } // namespace @@ -485,9 +502,221 @@ namespace Simplex return std::make_unique>(*this); } + + + template + ScalarWedgePolynomial::ScalarWedgePolynomial(const unsigned int degree) + : ScalarPolynomialsBase(degree, + compute_n_polynomials_wedge(dim, degree)) + {} + + + namespace + { + /** + * TODO + */ + static const constexpr std::array, 6> + wedge_table_1{ + {{{0, 0}}, {{1, 0}}, {{2, 0}}, {{0, 1}}, {{1, 1}}, {{2, 1}}}}; + + /** + * TODO + */ + static const constexpr std::array, 18> + 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 + ScalarWedgePolynomial::compute_value(const unsigned int i, + const Point & p) const + { + AssertDimension(dim, 3); + AssertIndexRange(this->degree(), 3); + + const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i]; + + const ScalarPolynomial<2> poly_tri(this->degree()); + const Point<2> p_tri(p[0], p[1]); + const auto v_tri = poly_tri.compute_value(pair[0], p_tri); + + const ScalarPolynomial<1> poly_line(this->degree()); + const Point<1> p_line(p[2]); + const auto v_line = poly_line.compute_value(pair[1], p_line); + + return v_tri * v_line; + } + + + + template + Tensor<1, dim> + ScalarWedgePolynomial::compute_grad(const unsigned int i, + const Point & p) const + { + AssertDimension(dim, 3); + AssertIndexRange(this->degree(), 3); + + const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i]; + + const ScalarPolynomial<2> poly_tri(this->degree()); + const Point<2> p_tri(p[0], p[1]); + const auto v_tri = poly_tri.compute_value(pair[0], p_tri); + const auto g_tri = poly_tri.compute_grad(pair[0], p_tri); + + const ScalarPolynomial<1> poly_line(this->degree()); + const Point<1> p_line(p[2]); + const auto v_line = poly_line.compute_value(pair[1], p_line); + const auto g_line = poly_line.compute_grad(pair[1], p_line); + + Tensor<1, dim> grad; + grad[0] = g_tri[0] * v_line; + grad[1] = g_tri[1] * v_line; + grad[2] = v_tri * g_line[0]; + + return grad; + } + + + + template + Tensor<2, dim> + ScalarWedgePolynomial::compute_grad_grad(const unsigned int i, + const Point & p) const + { + (void)i; + (void)p; + + Assert(false, ExcNotImplemented()); + return Tensor<2, dim>(); + } + + + + template + void + ScalarWedgePolynomial::evaluate( + const Point & unit_point, + std::vector & values, + std::vector> &grads, + std::vector> &grad_grads, + std::vector> &third_derivatives, + std::vector> &fourth_derivatives) const + { + (void)grads; + (void)grad_grads; + (void)third_derivatives; + (void)fourth_derivatives; + + if (values.size() == this->n()) + for (unsigned int i = 0; i < this->n(); i++) + values[i] = compute_value(i, unit_point); + + if (grads.size() == this->n()) + for (unsigned int i = 0; i < this->n(); i++) + grads[i] = compute_grad(i, unit_point); + } + + + + template + Tensor<1, dim> + ScalarWedgePolynomial::compute_1st_derivative(const unsigned int i, + const Point & p) const + { + return compute_grad(i, p); + } + + + + template + Tensor<2, dim> + ScalarWedgePolynomial::compute_2nd_derivative(const unsigned int i, + const Point & p) const + { + (void)i; + (void)p; + + Assert(false, ExcNotImplemented()); + + return {}; + } + + + + template + Tensor<3, dim> + ScalarWedgePolynomial::compute_3rd_derivative(const unsigned int i, + const Point & p) const + { + (void)i; + (void)p; + + Assert(false, ExcNotImplemented()); + + return {}; + } + + + + template + Tensor<4, dim> + ScalarWedgePolynomial::compute_4th_derivative(const unsigned int i, + const Point & p) const + { + (void)i; + (void)p; + + Assert(false, ExcNotImplemented()); + + return {}; + } + + + + template + std::string + ScalarWedgePolynomial::name() const + { + return "ScalarWedgePolynomial"; + } + + + + template + std::unique_ptr> + ScalarWedgePolynomial::clone() const + { + return std::make_unique>(*this); + } + + + template class ScalarPolynomial<1>; template class ScalarPolynomial<2>; template class ScalarPolynomial<3>; + template class ScalarWedgePolynomial<1>; + template class ScalarWedgePolynomial<2>; + template class ScalarWedgePolynomial<3>; } // namespace Simplex diff --git a/tests/simplex/wedge_01.cc b/tests/simplex/wedge_01.cc new file mode 100644 index 0000000000..5629248bee --- /dev/null +++ b/tests/simplex/wedge_01.cc @@ -0,0 +1,113 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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. +// +// --------------------------------------------------------------------- + + +// Distribute Simplex::FE_Wedge on a DoFHandler. + +#include + +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +void +test_3() +{ + const int dim = 3; + const int spacedim = 3; + + std::vector> vertices; + std::vector> cells; + +#if true + vertices.emplace_back(0.0, 0.0, 0.0); + vertices.emplace_back(1.0, 0.0, 0.0); + vertices.emplace_back(0.0, 1.0, 0.0); + vertices.emplace_back(0.0, 0.0, 1.0); + vertices.emplace_back(1.0, 0.0, 1.0); + vertices.emplace_back(0.0, 1.0, 1.0); + + { + CellData cell; + cell.vertices = {0, 1, 2, 3, 4, 5}; + cells.push_back(cell); + } +#else + vertices.emplace_back(0.0, 0.0, 0.0); + vertices.emplace_back(1.0, 0.0, 0.0); + vertices.emplace_back(0.0, 1.0, 0.0); + vertices.emplace_back(1.0, 1.0, 0.0); + vertices.emplace_back(0.0, 0.0, 1.0); + vertices.emplace_back(1.0, 0.0, 1.0); + vertices.emplace_back(0.0, 1.0, 1.0); + vertices.emplace_back(1.0, 1.0, 1.0); + + { + CellData cell; + cell.vertices = {0, 1, 2, 4, 5, 6}; + cells.push_back(cell); + } + { + CellData cell; + cell.vertices = {2, 1, 3, 6, 5, 7}; + cells.push_back(cell); + } +#endif + + Triangulation tria; + tria.create_triangulation(vertices, cells, SubCellData()); + + GridOut grid_out; + std::ofstream out("mesh.vtk"); + grid_out.write_vtk(tria, out); + + Simplex::FE_WedgeP fe(2); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + deallog << dof_handler.n_dofs() << std::endl; + + std::vector dof_indices; + for (const auto &cell : dof_handler.active_cell_iterators()) + { + for (const auto face : cell->face_indices()) + { +#if false + dof_indices.resize(face <= 1 ? 3 : 4); +#else + dof_indices.resize(fe.n_dofs_per_face(face)); +#endif + cell->face(face)->get_dof_indices(dof_indices); + + for (const auto i : dof_indices) + deallog << i << " "; + deallog << std::endl; + } + } +} + +int +main() +{ + initlog(); + + test_3(); +} diff --git a/tests/simplex/wedge_01.with_simplex_support=on.output b/tests/simplex/wedge_01.with_simplex_support=on.output new file mode 100644 index 0000000000..ad3705a248 --- /dev/null +++ b/tests/simplex/wedge_01.with_simplex_support=on.output @@ -0,0 +1,7 @@ + +DEAL::18 +DEAL::1 0 2 6 8 7 +DEAL::3 4 5 9 10 11 +DEAL::0 1 3 4 12 13 6 9 15 +DEAL::1 2 4 5 13 14 7 10 16 +DEAL::2 0 5 3 14 12 11 8 17 \ No newline at end of file