From fb1286c0815c1bc72547c1ab5f70ac5553ba8efe Mon Sep 17 00:00:00 2001 From: Jake Harmon Date: Thu, 30 Dec 2021 15:07:12 -0700 Subject: [PATCH] Add get_embedding_dofs() to FE_Nedelec (2-D) --- doc/news/changes/minor/20220125Harmon | 5 ++ include/deal.II/fe/fe_nedelec.h | 12 ++++ source/fe/fe_nedelec.cc | 56 +++++++++++++++++ tests/fe/nedelec_embed_dofs.cc | 91 +++++++++++++++++++++++++++ tests/fe/nedelec_embed_dofs.output | 41 ++++++++++++ 5 files changed, 205 insertions(+) create mode 100644 doc/news/changes/minor/20220125Harmon create mode 100644 tests/fe/nedelec_embed_dofs.cc create mode 100644 tests/fe/nedelec_embed_dofs.output diff --git a/doc/news/changes/minor/20220125Harmon b/doc/news/changes/minor/20220125Harmon new file mode 100644 index 0000000000..c72da63e28 --- /dev/null +++ b/doc/news/changes/minor/20220125Harmon @@ -0,0 +1,5 @@ +New: FE_Nedelec<2> now includes the function get_embedding_dofs(), +which permits direct identification of DoF indices between two +realizations of the FE_Nedelec cell. +
+(Jake Harmon, 2022/01/25) diff --git a/include/deal.II/fe/fe_nedelec.h b/include/deal.II/fe/fe_nedelec.h index 5857a056aa..f0e28fd73a 100644 --- a/include/deal.II/fe/fe_nedelec.h +++ b/include/deal.II/fe/fe_nedelec.h @@ -324,6 +324,18 @@ public: virtual std::unique_ptr> clone() const override; + /** + * For a finite element of degree larger than @p sub_degree, we return a + * vector which maps the numbering on an FE of degree @p sub_degree into the + * numbering on this element. + * + * Note that for the Nedelec element, by @p sub_degree, + * we refer to the maximal polynomial degree (in any coordinate direction) as + * opposed to the Nedelec degree. + */ + std::vector + get_embedding_dofs(const unsigned int sub_degree) const; + private: /** * Only for internal use. Its full name is @p get_dofs_per_object_vector diff --git a/source/fe/fe_nedelec.cc b/source/fe/fe_nedelec.cc index 677fa68973..9b8a7690f7 100644 --- a/source/fe/fe_nedelec.cc +++ b/source/fe/fe_nedelec.cc @@ -4121,7 +4121,63 @@ FE_Nedelec::memory_consumption() const return 0; } +template +std::vector +FE_Nedelec::get_embedding_dofs(const unsigned int sub_degree) const +{ + Assert((sub_degree > 0) && (sub_degree <= this->degree), + ExcIndexRange(sub_degree, 1, this->degree)); + + switch (dim) + { + case 2: + { + // The Nedelec cell has only Face (Line) and Cell DoFs... + const unsigned int n_face_dofs_sub = + GeometryInfo::lines_per_cell * sub_degree; + const unsigned int n_cell_dofs_sub = + 2 * (sub_degree - 1) * sub_degree; + + std::vector embedding_dofs(n_face_dofs_sub + + n_cell_dofs_sub); + + unsigned int i = 0; + + // Identify the Face/Line DoFs + while (i < n_face_dofs_sub) + { + const unsigned int face_index = i / sub_degree; + embedding_dofs[i] = i % sub_degree + face_index * this->degree; + ++i; + } + + // Identify the Cell DoFs + if (sub_degree >= 2) + { + const unsigned int n_face_dofs = + GeometryInfo::lines_per_cell * this->degree; + // For the first component + for (unsigned ku = 0; ku < sub_degree; ++ku) + for (unsigned kv = 2; kv <= sub_degree; ++kv) + embedding_dofs[i++] = + n_face_dofs + ku * (this->degree - 1) + (kv - 2); + + // For the second component + for (unsigned ku = 2; ku <= sub_degree; ++ku) + for (unsigned kv = 0; kv < sub_degree; ++kv) + embedding_dofs[i++] = n_face_dofs + + this->degree * (this->degree - 1) + + (ku - 2) * (this->degree) + kv; + } + Assert(i == (n_face_dofs_sub + n_cell_dofs_sub), ExcInternalError()); + return embedding_dofs; + } + default: + Assert(false, ExcNotImplemented()); + return std::vector(); + } +} //----------------------------------------------------------------------// diff --git a/tests/fe/nedelec_embed_dofs.cc b/tests/fe/nedelec_embed_dofs.cc new file mode 100644 index 0000000000..451c09c3c2 --- /dev/null +++ b/tests/fe/nedelec_embed_dofs.cc @@ -0,0 +1,91 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 - 2021 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. +// +// --------------------------------------------------------------------- + +// Checks if get_embedding_dofs() correctly maps the local DoF indices +// of an FE_Nedelec element of sub_degree to an element of sup_degree, +// where sub_degree <= sup_degree + +#include + +#include + +#include "interpolate_common.h" + +template +void +check_nedelec_embed(const unsigned int sub_degree, + const unsigned int sup_degree) +{ + Assert((sub_degree > 0) && (sub_degree <= sup_degree), + ExcIndexRange(sub_degree, 1, sup_degree)); + + const FE_Nedelec fe_sub(sub_degree - 1); + const FE_Nedelec fe_sup(sup_degree - 1); + + deallog << fe_sub.get_name() << ' ' << fe_sup.get_name() << ' ' + << fe_sub.dofs_per_cell << ' ' << fe_sup.dofs_per_cell; + + // Set the quadrature sufficiently high for the enriched finite element space + const QGauss quadrature(fe_sup.degree + 1); + + std::vector dofs_sub(fe_sub.dofs_per_cell, 0.); + std::vector dofs_sup(fe_sup.dofs_per_cell, 0.); + + // Assign arbitrary values to dofs_sub + for (unsigned int i = 0; i < dofs_sub.size(); ++i) + dofs_sub[i] = std::sin(i + 0.5); + + + // Map the DoFs from fe_sub to fe_sup + const std::vector embedding_dofs = + fe_sup.get_embedding_dofs(fe_sub.degree); + for (unsigned int i = 0; i < dofs_sub.size(); ++i) + dofs_sup[embedding_dofs[i]] = dofs_sub[i]; + + // Compare the values at each quadrature point + double result = 0.; + for (unsigned int k = 0; k < quadrature.size(); ++k) + for (unsigned int comp = 0; comp < fe_sup.n_components(); ++comp) + { + double eval_sub = 0., eval_sup = 0.; + + for (unsigned int i = 0; i < dofs_sub.size(); ++i) + eval_sub += + dofs_sub[i] * + fe_sub.shape_value_component(i, quadrature.point(k), comp); + + for (unsigned int i = 0; i < dofs_sup.size(); ++i) + eval_sup += + dofs_sup[i] * + fe_sup.shape_value_component(i, quadrature.point(k), comp); + + const double diff = std::abs(eval_sub - eval_sup); + result = std::max(result, diff); + } + deallog << " max diff " << result << std::endl; +} + +int +main() +{ + initlog(); + + // For 2-D + { + for (unsigned int low_deg = 1; low_deg <= 5; ++low_deg) + for (unsigned int high_deg = low_deg; high_deg <= 10; ++high_deg) + check_nedelec_embed<2>(low_deg, high_deg); + } +} diff --git a/tests/fe/nedelec_embed_dofs.output b/tests/fe/nedelec_embed_dofs.output new file mode 100644 index 0000000000..1162ee76e4 --- /dev/null +++ b/tests/fe/nedelec_embed_dofs.output @@ -0,0 +1,41 @@ + +DEAL::FE_Nedelec<2>(0) FE_Nedelec<2>(0) 4 4 max diff 0.00000 +DEAL::FE_Nedelec<2>(0) FE_Nedelec<2>(1) 4 12 max diff 0.00000 +DEAL::FE_Nedelec<2>(0) FE_Nedelec<2>(2) 4 24 max diff 0.00000 +DEAL::FE_Nedelec<2>(0) FE_Nedelec<2>(3) 4 40 max diff 0.00000 +DEAL::FE_Nedelec<2>(0) FE_Nedelec<2>(4) 4 60 max diff 0.00000 +DEAL::FE_Nedelec<2>(0) FE_Nedelec<2>(5) 4 84 max diff 0.00000 +DEAL::FE_Nedelec<2>(0) FE_Nedelec<2>(6) 4 112 max diff 0.00000 +DEAL::FE_Nedelec<2>(0) FE_Nedelec<2>(7) 4 144 max diff 0.00000 +DEAL::FE_Nedelec<2>(0) FE_Nedelec<2>(8) 4 180 max diff 0.00000 +DEAL::FE_Nedelec<2>(0) FE_Nedelec<2>(9) 4 220 max diff 0.00000 +DEAL::FE_Nedelec<2>(1) FE_Nedelec<2>(1) 12 12 max diff 0.00000 +DEAL::FE_Nedelec<2>(1) FE_Nedelec<2>(2) 12 24 max diff 0.00000 +DEAL::FE_Nedelec<2>(1) FE_Nedelec<2>(3) 12 40 max diff 0.00000 +DEAL::FE_Nedelec<2>(1) FE_Nedelec<2>(4) 12 60 max diff 0.00000 +DEAL::FE_Nedelec<2>(1) FE_Nedelec<2>(5) 12 84 max diff 0.00000 +DEAL::FE_Nedelec<2>(1) FE_Nedelec<2>(6) 12 112 max diff 0.00000 +DEAL::FE_Nedelec<2>(1) FE_Nedelec<2>(7) 12 144 max diff 0.00000 +DEAL::FE_Nedelec<2>(1) FE_Nedelec<2>(8) 12 180 max diff 0.00000 +DEAL::FE_Nedelec<2>(1) FE_Nedelec<2>(9) 12 220 max diff 0.00000 +DEAL::FE_Nedelec<2>(2) FE_Nedelec<2>(2) 24 24 max diff 0.00000 +DEAL::FE_Nedelec<2>(2) FE_Nedelec<2>(3) 24 40 max diff 0.00000 +DEAL::FE_Nedelec<2>(2) FE_Nedelec<2>(4) 24 60 max diff 0.00000 +DEAL::FE_Nedelec<2>(2) FE_Nedelec<2>(5) 24 84 max diff 0.00000 +DEAL::FE_Nedelec<2>(2) FE_Nedelec<2>(6) 24 112 max diff 0.00000 +DEAL::FE_Nedelec<2>(2) FE_Nedelec<2>(7) 24 144 max diff 0.00000 +DEAL::FE_Nedelec<2>(2) FE_Nedelec<2>(8) 24 180 max diff 0.00000 +DEAL::FE_Nedelec<2>(2) FE_Nedelec<2>(9) 24 220 max diff 0.00000 +DEAL::FE_Nedelec<2>(3) FE_Nedelec<2>(3) 40 40 max diff 0.00000 +DEAL::FE_Nedelec<2>(3) FE_Nedelec<2>(4) 40 60 max diff 0.00000 +DEAL::FE_Nedelec<2>(3) FE_Nedelec<2>(5) 40 84 max diff 0.00000 +DEAL::FE_Nedelec<2>(3) FE_Nedelec<2>(6) 40 112 max diff 0.00000 +DEAL::FE_Nedelec<2>(3) FE_Nedelec<2>(7) 40 144 max diff 0.00000 +DEAL::FE_Nedelec<2>(3) FE_Nedelec<2>(8) 40 180 max diff 0.00000 +DEAL::FE_Nedelec<2>(3) FE_Nedelec<2>(9) 40 220 max diff 0.00000 +DEAL::FE_Nedelec<2>(4) FE_Nedelec<2>(4) 60 60 max diff 0.00000 +DEAL::FE_Nedelec<2>(4) FE_Nedelec<2>(5) 60 84 max diff 0.00000 +DEAL::FE_Nedelec<2>(4) FE_Nedelec<2>(6) 60 112 max diff 0.00000 +DEAL::FE_Nedelec<2>(4) FE_Nedelec<2>(7) 60 144 max diff 0.00000 +DEAL::FE_Nedelec<2>(4) FE_Nedelec<2>(8) 60 180 max diff 0.00000 +DEAL::FE_Nedelec<2>(4) FE_Nedelec<2>(9) 60 220 max diff 0.00000 -- 2.39.5