Add get_embedding_dofs() to FE_Nedelec (2-D)

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.
+<br>
+(Jake Harmon, 2022/01/25)

include/deal.II/fe/fe_nedelec.h

@@ -324,6 +324,18 @@ class FE_Nedelec : public FE_PolyTensor<dim>
   virtual std::unique_ptr<FiniteElement<dim, dim>>
   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<unsigned int>
+  get_embedding_dofs(const unsigned int sub_degree) const;
+
 private:
   /**
    * Only for internal use. Its full name is @p get_dofs_per_object_vector

source/fe/fe_nedelec.cc

@@ -4121,7 +4121,63 @@ FE_Nedelec<dim>::memory_consumption() const
   return 0;
 }
 
+template <int dim>
+std::vector<unsigned int>
+FE_Nedelec<dim>::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<dim>::lines_per_cell * sub_degree;
+          const unsigned int n_cell_dofs_sub =
+            2 * (sub_degree - 1) * sub_degree;
+
+          std::vector<unsigned int> 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<dim>::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<unsigned int>();
+    }
+}
 //----------------------------------------------------------------------//
 
 

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 <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_nedelec.h>
+
+#include "interpolate_common.h"
+
+template <int dim>
+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<dim> fe_sub(sub_degree - 1);
+  const FE_Nedelec<dim> 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<dim> quadrature(fe_sup.degree + 1);
+
+  std::vector<double> dofs_sub(fe_sub.dofs_per_cell, 0.);
+  std::vector<double> 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<unsigned int> 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);
+  }
+}

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