]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add support for hanging nodes to the NedelecSZ class.
authorSebastian Kinnewig <kinnewig@ifam.uni-hannover.de>
Fri, 19 May 2023 09:19:44 +0000 (11:19 +0200)
committerSebastian Kinnewig <kinnewig@ifam.uni-hannover.de>
Tue, 16 Jan 2024 18:19:11 +0000 (19:19 +0100)
doc/news/changes/major/20231128SebastianKinnewig [new file with mode: 0644]
include/deal.II/fe/fe_nedelec_sz.h
source/dofs/dof_tools_constraints.cc
source/fe/fe_nedelec_sz.cc

diff --git a/doc/news/changes/major/20231128SebastianKinnewig b/doc/news/changes/major/20231128SebastianKinnewig
new file mode 100644 (file)
index 0000000..c936511
--- /dev/null
@@ -0,0 +1,3 @@
+New: Support for hanging nodes in FE_NedelecSZ.
+<br>
+(Sebastian Kinnewig, 2023/11/28)
index de396e6bb7d34e507dffff3f5bc4b1c13dc6ecbf..234fbe9672be878bf860c83f9785ba8d93cf353e 100644 (file)
@@ -63,13 +63,22 @@ DEAL_II_NAMESPACE_OPEN
  * vertices are decided such that the second has a higher global vertex
  * numbering than the fourth.
  *
- * Note that this element does not support non-conforming meshes at this time.
+ * To overcome the sign conflict in the case of non-conforming meshes, hanging
+ * edges (and faces in 3D) must be addressed. This element chooses the
+ * orientation of hanging edges (and faces) based on the parent cell. Moreover,
+ * the associated constraint matrix that restricts the additional DoFs
+ * originating from the hanging edges and faces must be adapted to account for
+ * the orientation of the edges and faces.
  *
  * Further details on this element, including some benchmarking, can be
  * in the paper R. Kynch, P. Ledger: <b>Resolving the sign conflict
  * problem for hp-hexahedral N&eacute;d&eacute;lec elements with application to
  * eddy current problems</b>, Computers & Structures 181, 41-54, 2017 (see
  * https://doi.org/10.1016/j.compstruc.2016.05.021).
+ * For details on the implementation of the hanging node constraints, see the
+ * paper S. Kinnewig, T. Wick, S. Beuchler: <b>Algorithmic realization of the
+ * solution to the sign conflict problem for hanging nodes on hp-hexahedral
+ * Nédélec elements</b>, preprint, 2023 (https://arxiv.org/abs/2306.01416).
  */
 template <int dim, int spacedim = dim>
 class FE_NedelecSZ : public FiniteElement<dim, dim>
index 9b0066cfc405c75c8497ec3736cb078691c9740d..bb9907c83ac770be05c77f197fa83a21bf9285cb 100644 (file)
@@ -602,9 +602,8 @@ namespace DoFTools
 
     template <typename number>
     void
-    make_oldstyle_hanging_node_constraints(const DoFHandler<1> &,
-                                           AffineConstraints<number> &,
-                                           std::integral_constant<int, 1>)
+    make_hp_hanging_node_constraints(const DoFHandler<1, 2> &,
+                                     AffineConstraints<number> &)
     {
       // nothing to do for regular dof handlers in 1d
     }
@@ -612,8 +611,9 @@ namespace DoFTools
 
     template <typename number>
     void
-    make_hp_hanging_node_constraints(const DoFHandler<1, 2> &,
-                                     AffineConstraints<number> &)
+    make_hanging_node_constraints_nedelec(const dealii::DoFHandler<1, 2> &,
+                                          AffineConstraints<number> &,
+                                          std::integral_constant<int, 1>)
     {
       // nothing to do for regular dof handlers in 1d
     }
@@ -639,6 +639,17 @@ namespace DoFTools
     }
 
 
+    template <typename number, int spacedim>
+    void
+    make_hanging_node_constraints_nedelec(
+      const dealii::DoFHandler<1, spacedim> & /*dof_handler*/,
+      AffineConstraints<number> & /*constraints*/,
+      std::integral_constant<int, 1>)
+    {
+      // nothing to do for dof handlers in 1d
+    }
+
+
     template <typename number, int spacedim>
     void
     make_oldstyle_hanging_node_constraints(
@@ -1070,6 +1081,1269 @@ namespace DoFTools
 
 
 
+    template <int dim_, int spacedim, typename number>
+    void
+    make_hanging_node_constraints_nedelec(
+      const DoFHandler<dim_, spacedim> &dof_handler,
+      AffineConstraints<number>        &constraints,
+      std::integral_constant<int, 2>)
+    {
+      // Parts of this function are very similar to
+      // make_oldstyle_hanging_node_constraints.
+      // Therefore, only the parts that differ from the
+      // make_oldstyle_hanging_node_constraints are commented on.
+
+      const unsigned int dim = 2;
+
+      std::vector<types::global_dof_index> face_dof_indices;
+      std::map<types::global_dof_index, std::set<types::global_dof_index>>
+        depends_on;
+
+      // loop over all lines
+      for (const auto &cell : dof_handler.active_cell_iterators())
+        {
+          // skip artificial cells
+          if (cell->is_artificial())
+            continue;
+
+          // loop over all faces:
+          for (const unsigned int f : cell->face_indices())
+            {
+              // check if the neighbor is refined; if so, we need to
+              // treat the constraints on this interface
+              if (!cell->face(f)->has_children())
+                continue;
+
+              Assert(cell->face(f)->n_active_fe_indices() == 1,
+                     ExcInternalError());
+              Assert(cell->face(f)->fe_index_is_active(
+                       cell->active_fe_index()) == true,
+                     ExcInternalError());
+
+#ifdef DEBUG
+              for (unsigned int c = 0; c < cell->face(f)->n_children(); ++c)
+                {
+                  if (cell->neighbor_child_on_subface(f, c)->is_artificial())
+                    continue;
+
+                  Assert(cell->face(f)->child(c)->n_active_fe_indices() == 1,
+                         ExcInternalError());
+
+                  Assert(cell->face(f)->child(c)->fe_index_is_active(
+                           cell->active_fe_index()) == true,
+                         ExcNotImplemented());
+                }
+#endif // DEBUG
+
+              // Ok, start up the work:
+              const FiniteElement<dim, spacedim> &fe = cell->get_fe();
+
+              const unsigned int n_dofs = fe.n_dofs_per_line();
+              face_dof_indices.resize(n_dofs);
+
+              cell->face(f)->get_dof_indices(face_dof_indices);
+              const std::vector<types::global_dof_index> dof_on_mother_face =
+                face_dof_indices;
+
+              cell->face(f)->child(0)->get_dof_indices(face_dof_indices);
+              const std::vector<types::global_dof_index> dof_on_child_face_0 =
+                face_dof_indices;
+
+              cell->face(f)->child(1)->get_dof_indices(face_dof_indices);
+              const std::vector<types::global_dof_index> dof_on_child_face_1 =
+                face_dof_indices;
+
+              // As the Nedelec elements are oriented, we need to take care of
+              // the orientation of the lines.
+              // Remark: "false" indicates the line is not flipped.
+              //         "true" indicates the line is flipped.
+
+              // get the orientation of the faces
+              const bool direction_mother = (cell->face(f)->vertex_index(0) >
+                                             cell->face(f)->vertex_index(1)) ?
+                                              false :
+                                              true;
+              const bool direction_child_0 =
+                (cell->face(f)->child(0)->vertex_index(0) >
+                 cell->face(f)->child(0)->vertex_index(1)) ?
+                  false :
+                  true;
+              const bool direction_child_1 =
+                (cell->face(f)->child(1)->vertex_index(0) >
+                 cell->face(f)->child(1)->vertex_index(1)) ?
+                  false :
+                  true;
+
+              for (unsigned int row = 0; row < n_dofs; ++row)
+                {
+                  constraints.add_line(dof_on_child_face_0[row]);
+                  constraints.add_line(dof_on_child_face_1[row]);
+                }
+
+              for (unsigned int row = 0; row < n_dofs; ++row)
+                {
+                  for (unsigned int dof_i_on_mother = 0;
+                       dof_i_on_mother < n_dofs;
+                       ++dof_i_on_mother)
+                    {
+                      // We need to keep in mind that, if we use a FE_System
+                      // with multiple FE_NedelecSZ blocks inside, we need
+                      // to consider, that n_dofs depends on the number
+                      // of FE_NedelecSZ blocks used.
+                      unsigned int shift_0 =
+                        (direction_mother == direction_child_0) ? 0 : n_dofs;
+                      constraints.add_entry(dof_on_child_face_0[row],
+                                            dof_on_mother_face[dof_i_on_mother],
+                                            fe.constraints()(row + shift_0,
+                                                             dof_i_on_mother));
+
+                      unsigned int shift_1 =
+                        (direction_mother == direction_child_1) ? 0 : n_dofs;
+                      constraints.add_entry(dof_on_child_face_1[row],
+                                            dof_on_mother_face[dof_i_on_mother],
+                                            fe.constraints()(row + shift_1,
+                                                             dof_i_on_mother));
+                    }
+                }
+            }
+        }
+    }
+
+
+    template <int dim_, int spacedim, typename number>
+    void
+    make_hanging_node_constraints_nedelec(
+      const DoFHandler<dim_, spacedim> &dof_handler,
+      AffineConstraints<number>        &constraints,
+      std::integral_constant<int, 3>)
+    {
+      // Parts of this function are very similar to
+      // make_oldstyle_hanging_node_constraints.
+      // Therefore, only the parts that differ from the
+      // make_oldstyle_hanging_node_constraints are commented on.
+
+      const unsigned int dim = 3;
+
+      std::vector<types::global_dof_index> dofs_on_mother;
+      std::vector<types::global_dof_index> dofs_on_children;
+
+      // loop over all quads
+      for (const auto &cell : dof_handler.active_cell_iterators())
+        {
+          // skip artificial cells
+          if (cell->is_artificial())
+            continue;
+
+          // loop over all faces
+          for (const unsigned int face : cell->face_indices())
+            {
+              // skip cells without children
+              if (cell->face(face)->has_children() == false)
+                continue;
+
+              if (cell->get_fe().n_dofs_per_face(face) == 0)
+                continue;
+
+              Assert(cell->face(face)->refinement_case() ==
+                       RefinementCase<dim - 1>::isotropic_refinement,
+                     ExcNotImplemented());
+
+              AssertDimension(cell->face(face)->n_active_fe_indices(), 1);
+
+              Assert(cell->face(face)->fe_index_is_active(
+                       cell->active_fe_index()) == true,
+                     ExcInternalError());
+
+#ifdef DEBUG
+
+              for (unsigned int c = 0; c < cell->face(face)->n_children(); ++c)
+                {
+                  if (cell->neighbor_child_on_subface(face, c)->is_artificial())
+                    continue;
+
+                  AssertDimension(
+                    cell->face(face)->child(c)->n_active_fe_indices(), 1);
+
+                  Assert(cell->face(face)->child(c)->fe_index_is_active(
+                           cell->active_fe_index()) == true,
+                         ExcNotImplemented());
+
+                  for (unsigned int e = 0;
+                       e < GeometryInfo<dim>::vertices_per_face;
+                       ++e)
+                    {
+                      Assert(cell->face(face)
+                                 ->child(c)
+                                 ->line(e)
+                                 ->n_active_fe_indices() == 1,
+                             ExcNotImplemented());
+
+                      Assert(
+                        cell->face(face)->child(c)->line(e)->fe_index_is_active(
+                          cell->active_fe_index()) == true,
+                        ExcNotImplemented());
+                    }
+                }
+
+              for (unsigned int e = 0; e < GeometryInfo<dim>::vertices_per_face;
+                   ++e)
+                {
+                  Assert(cell->face(face)->line(e)->n_active_fe_indices() == 1,
+                         ExcNotImplemented());
+
+                  Assert(cell->face(face)->line(e)->fe_index_is_active(
+                           cell->active_fe_index()) == true,
+                         ExcNotImplemented());
+                }
+#endif // DEBUG
+
+              // Ok, start up the work
+              const FiniteElement<dim, spacedim> &fe = cell->get_fe();
+              const unsigned int fe_index            = cell->active_fe_index();
+
+              // get the polynomial degree
+              unsigned int degree(fe.degree);
+
+              // get the number of DoFs on mother and children;
+              // number of DoFs on the mother
+              const unsigned int n_dofs_on_mother = fe.n_dofs_per_face(face);
+              dofs_on_mother.resize(n_dofs_on_mother);
+
+              const unsigned int n_lines_on_mother =
+                GeometryInfo<dim>::lines_per_face;
+
+              // number of internal lines of the children;
+              // for more details see description of the
+              // GeometryInfo<dim> class
+              // .................
+              // .       |       .
+              // .  c2   1   c3  .
+              // .       |       .
+              // .---2---+---3---.
+              // .       |       .
+              // .  c0   0   c1  .
+              // .       |       .
+              // .................
+              const unsigned int n_internal_lines_on_children = 4;
+
+              // number of external lines of the children
+              // +---6--------7--+
+              // |       .       |
+              // 1  c2   .   c3  3
+              // |       .       |
+              // |...............|
+              // |       .       |
+              // 0  c0   .   c1  2
+              // |       .       |
+              // +---4---+---5---+
+              const unsigned int n_external_lines_on_children = 8;
+
+              const unsigned int n_lines_on_children =
+                n_internal_lines_on_children + n_external_lines_on_children;
+
+              // we only consider the isotropic case here
+              const unsigned int n_children_per_face =
+                GeometryInfo<dim>::max_children_per_face;
+              const unsigned int n_children_per_line =
+                GeometryInfo<dim - 1>::max_children_per_face;
+
+              // number of DoFs on the children
+              // Remark: Nedelec elements have no DoFs on the vertices,
+              // therefore we skip the vertices
+              const unsigned int n_dofs_on_children =
+                (n_lines_on_children * fe.n_dofs_per_line() +
+                 n_children_per_face * fe.n_dofs_per_quad(face));
+
+              dofs_on_children.clear();
+              dofs_on_children.reserve(n_dofs_on_children);
+
+              AssertDimension(n_dofs_on_mother, fe.constraints().n());
+              AssertDimension(n_dofs_on_children, fe.constraints().m());
+
+              // get the current face
+              const typename DoFHandler<dim, dim>::face_iterator this_face =
+                cell->face(face);
+
+              // fill the DoFs on the mother:
+              unsigned int next_index = 0;
+
+              // DoFs on vertices:
+              // Nedelec elements have no DoFs on the vertices
+
+              // DoFs on lines:
+              for (unsigned int line = 0;
+                   line < GeometryInfo<dim>::lines_per_face;
+                   ++line)
+                for (unsigned int dof = 0; dof != fe.n_dofs_per_line(); ++dof)
+                  dofs_on_mother[next_index++] =
+                    this_face->line(line)->dof_index(dof, fe_index);
+
+              // DoFs on the face:
+              for (unsigned int dof = 0; dof != fe.n_dofs_per_quad(face); ++dof)
+                dofs_on_mother[next_index++] =
+                  this_face->dof_index(dof, fe_index);
+
+              // check that we have added all DoFs
+              AssertDimension(next_index, dofs_on_mother.size());
+
+              // the implementation does not support anisotropic refinement
+              Assert(!dof_handler.get_triangulation()
+                        .get_anisotropic_refinement_flag(),
+                     ExcInternalError());
+
+              // fill the DoF on the children:
+              // DoFs on vertices:
+              // Nedelec elements have no DoFs on the vertices
+
+              // DoFs on lines:
+              // the DoFs on the interior lines to the children; the order
+              // of these lines is shown above (see
+              // n_internal_lines_on_children)
+              for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof)
+                dofs_on_children.push_back(
+                  this_face->child(0)->line(1)->dof_index(dof, fe_index));
+
+              for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof)
+                dofs_on_children.push_back(
+                  this_face->child(2)->line(1)->dof_index(dof, fe_index));
+
+              for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof)
+                dofs_on_children.push_back(
+                  this_face->child(0)->line(3)->dof_index(dof, fe_index));
+
+              for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof)
+                dofs_on_children.push_back(
+                  this_face->child(1)->line(3)->dof_index(dof, fe_index));
+
+              // DoFs on the bordering lines:
+              // DoFs on the exterior lines to the children; the order of
+              // these lines is shown above (see n_external_lines_on_children)
+              for (unsigned int line = 0;
+                   line < GeometryInfo<dim>::lines_per_face;
+                   ++line)
+                for (unsigned int child = 0; child < n_children_per_line;
+                     ++child)
+                  for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof)
+                    dofs_on_children.push_back(
+                      this_face->line(line)->child(child)->dof_index(dof,
+                                                                     fe_index));
+
+              // DoFs on the faces of the four children:
+              for (unsigned int child = 0; child < n_children_per_face; ++child)
+                {
+                  // skip artificial cells
+                  if (cell->neighbor_child_on_subface(face, child)
+                        ->is_artificial())
+                    continue;
+
+                  for (unsigned int dof = 0; dof < fe.n_dofs_per_quad(face);
+                       ++dof)
+                    dofs_on_children.push_back(
+                      this_face->child(child)->dof_index(dof, fe_index));
+                } // rof: child
+
+              // consistency check:
+              // note: we can get fewer DoFs when we have artificial cells
+              Assert(dofs_on_children.size() <= n_dofs_on_children,
+                     ExcInternalError());
+
+              // As the Nedelec elements are oriented, we need to take care of
+              // the orientation of the lines.
+              // Remark: "false" indicates the line is not flipped.
+              //         "true" indicates the line is flipped.
+
+              // Orientation - Lines:
+              // get the orientation from the edges from the mother cell
+              std::vector<bool> direction_mother(
+                GeometryInfo<dim>::lines_per_face, false);
+              for (unsigned int line = 0;
+                   line < GeometryInfo<dim>::lines_per_face;
+                   ++line)
+                if (this_face->line(line)->vertex_index(0) >
+                    this_face->line(line)->vertex_index(1))
+                  direction_mother[line] = true;
+
+              // get the orientation from the intern edges of the children
+              std::vector<bool> direction_child_intern(
+                n_internal_lines_on_children, false);
+
+              // get the global vertex index of vertex in the center;
+              // we need this vertex index, to compute the direction
+              // of the internal edges
+              unsigned int center = this_face->child(0)->vertex_index(3);
+
+              // compute the direction of the internal edges
+              for (unsigned int line = 0; line < n_internal_lines_on_children;
+                   ++line)
+                if (line % 2 == 0)
+                  {
+                    direction_child_intern[line] =
+                      this_face->line(line)->child(0)->vertex_index(1) <
+                          center ?
+                        false :
+                        true;
+                  }
+                else
+                  {
+                    direction_child_intern[line] =
+                      this_face->line(line)->child(0)->vertex_index(1) >
+                          center ?
+                        false :
+                        true;
+                  }
+
+              // compute the direction of the outer edges
+              std::vector<bool> direction_child(n_external_lines_on_children,
+                                                false);
+              for (unsigned int line = 0;
+                   line < GeometryInfo<dim>::lines_per_face;
+                   ++line)
+                {
+                  if (this_face->line(line)->child(0)->vertex_index(0) >
+                      this_face->line(line)->child(0)->vertex_index(1))
+                    direction_child[2 * line] = true;
+                  if (this_face->line(line)->child(1)->vertex_index(0) >
+                      this_face->line(line)->child(1)->vertex_index(1))
+                    direction_child[2 * line + 1] = true;
+                }
+
+
+              // Orientation - Faces:
+              bool              mother_flip_x  = false;
+              bool              mother_flip_y  = false;
+              bool              mother_flip_xy = false;
+              std::vector<bool> child_flip_x(n_children_per_face, false);
+              std::vector<bool> child_flip_y(n_children_per_face, false);
+              std::vector<bool> child_flip_xy(n_children_per_face, false);
+              const unsigned int
+                vertices_adjacent_on_face[GeometryInfo<dim>::vertices_per_face]
+                                         [2] = {{1, 2}, {0, 3}, {3, 0}, {2, 1}};
+
+              {
+                // Mother
+                // get the position of the vertex with the highest number
+                unsigned int current_glob = cell->face(face)->vertex_index(0);
+                unsigned int current_max  = 0;
+                for (unsigned int v = 1;
+                     v < GeometryInfo<dim>::vertices_per_face;
+                     ++v)
+                  if (current_glob < this_face->vertex_index(v))
+                    {
+                      current_max  = v;
+                      current_glob = this_face->vertex_index(v);
+                    }
+
+                // if the vertex with the highest DoF index is in the lower row
+                // of the face, the face is flipped in y direction
+                if (current_max < 2)
+                  mother_flip_y = true;
+
+                // if the vertex with the highest DoF index is on the left side
+                // of the face is flipped in x direction
+                if (current_max % 2 == 0)
+                  mother_flip_x = true;
+
+                // get the minor direction of the face of the mother
+                if (this_face->vertex_index(
+                      vertices_adjacent_on_face[current_max][0]) <
+                    this_face->vertex_index(
+                      vertices_adjacent_on_face[current_max][1]))
+                  mother_flip_xy = true;
+              }
+
+              // Children:
+              // get the orientation of the faces of the children
+              for (unsigned int child = 0; child < n_children_per_face; ++child)
+                {
+                  unsigned int current_max = 0;
+                  unsigned int current_glob =
+                    this_face->child(child)->vertex_index(0);
+
+                  for (unsigned int v = 1;
+                       v < GeometryInfo<dim>::vertices_per_face;
+                       ++v)
+                    if (current_glob < this_face->child(child)->vertex_index(v))
+                      {
+                        current_max  = v;
+                        current_glob = this_face->child(child)->vertex_index(v);
+                      }
+
+                  if (current_max < 2)
+                    child_flip_y[child] = true;
+
+                  if (current_max % 2 == 0)
+                    child_flip_x[child] = true;
+
+                  if (this_face->child(child)->vertex_index(
+                        vertices_adjacent_on_face[current_max][0]) <
+                      this_face->child(child)->vertex_index(
+                        vertices_adjacent_on_face[current_max][1]))
+                    child_flip_xy[child] = true;
+
+                  child_flip_xy[child] = mother_flip_xy;
+                }
+
+              // copy the constraint matrix, since we need to modify that matrix
+              std::vector<std::vector<double>> constraints_matrix(
+                n_lines_on_children * fe.n_dofs_per_line() +
+                  n_children_per_face * fe.n_dofs_per_quad(),
+                std::vector<double>(dofs_on_mother.size(), 0));
+
+              {
+                // copy the constraint matrix
+                // internal lines
+                for (unsigned int line = 0; line < n_internal_lines_on_children;
+                     ++line)
+                {
+                  unsigned int row_start   = line * fe.n_dofs_per_line();
+                  unsigned int line_mother = line / 2;
+                  unsigned int row_mother =
+                    (line_mother * 2) * fe.n_dofs_per_line();
+                  for (unsigned int row = 0; row < fe.n_dofs_per_line();
+                       ++row)
+                    for (unsigned int i = 0;
+                         i < n_lines_on_mother * fe.n_dofs_per_line();
+                         ++i)
+                      constraints_matrix[row + row_start][i] =
+                        fe.constraints()(row + row_mother, i);
+                }
+
+                for (unsigned int line = 0; line < n_internal_lines_on_children;
+                     ++line)
+                  {
+                    unsigned int row_start   = line * fe.n_dofs_per_line();
+                    unsigned int line_mother = line / 2;
+                    unsigned int row_mother =
+                      (line_mother * 2) * fe.n_dofs_per_line();
+                    for (unsigned int row = 0; row < fe.n_dofs_per_line();
+                         ++row)
+                      for (unsigned int i = n_lines_on_mother * fe.n_dofs_per_line();
+                           i < dofs_on_mother.size();
+                           ++i)
+                        constraints_matrix[row + row_start][i] =
+                          fe.constraints()(row + row_mother, i);
+                  }
+
+                // external lines
+                unsigned int row_offset =
+                  n_internal_lines_on_children * fe.n_dofs_per_line();
+                for (unsigned int line = 0; line < n_external_lines_on_children;
+                     line++)
+                  {
+                    unsigned int row_start   = line * fe.n_dofs_per_line();
+                    unsigned int line_mother = line / 2;
+                    unsigned int row_mother =
+                      (line_mother * 2) * fe.n_dofs_per_line();
+                    for (unsigned int row = row_offset;
+                         row < row_offset + fe.n_dofs_per_line();
+                         ++row)
+                      for (unsigned int i = 0; i < dofs_on_mother.size(); ++i)
+                        constraints_matrix[row + row_start][i] =
+                          fe.constraints()(row + row_mother, i);
+                  }
+
+                // copy the weights for the faces
+                row_offset = n_lines_on_children * fe.n_dofs_per_line();
+                for (unsigned int face = 0; face < n_children_per_face; ++face)
+                  {
+                    unsigned int row_start = face * fe.n_dofs_per_quad();
+                    for (unsigned int row = row_offset;
+                         row < row_offset + fe.n_dofs_per_quad();
+                         row++)
+                      for (unsigned int i = 0; i < dofs_on_mother.size(); ++i)
+                        constraints_matrix[row + row_start][i] =
+                          fe.constraints()(row, i);
+                  }
+              }
+
+              // Modify the matrix
+              // Edge - Edge:
+              // Interior edges: the interior edges have support on the
+              // corresponding edges and faces loop over all 4 intern edges
+              for (unsigned int i = 0;
+                   i < n_internal_lines_on_children * fe.n_dofs_per_line();
+                   ++i)
+                {
+                  unsigned int line_i = i / fe.n_dofs_per_line();
+                  unsigned int tmp_i  = i % degree;
+
+                  // loop over the edges of the mother cell
+                  for (unsigned int j = 0;
+                       j < n_lines_on_mother * fe.n_dofs_per_line();
+                       ++j)
+                    {
+                      unsigned int line_j = j / fe.n_dofs_per_line();
+                      unsigned int tmp_j  = j % degree;
+
+                      if ((line_i < 2 && line_j < 2) ||
+                          (line_i >= 2 && line_j >= 2))
+                        {
+                          if (direction_child_intern[line_i] !=
+                              direction_mother[line_j])
+                            {
+                              if ((tmp_i + tmp_j) % 2 == 1)
+                                { // anti-symmetric
+                                  constraints_matrix[i][j] *= -1.0;
+                                }
+                            }
+                        }
+                      else
+                        {
+                          if (direction_mother[line_i])
+                            {
+                              if ((tmp_i + tmp_j) % 2 == 1)
+                                { // anti-symmetric
+                                  constraints_matrix[i][j] *= -1.0;
+                                }
+                            }
+                        }
+                    }
+                }
+
+              // Exterior edges:
+              for (unsigned int i =
+                     n_internal_lines_on_children * fe.n_dofs_per_line();
+                   i < n_lines_on_children * fe.n_dofs_per_line();
+                   ++i)
+                {
+                  unsigned int line_i = (i / fe.n_dofs_per_line()) - 4;
+                  unsigned int tmp_i  = i % degree;
+
+                  // loop over the edges of the mother cell
+                  for (unsigned int j = 0;
+                       j < n_lines_on_mother * fe.n_dofs_per_line();
+                       ++j)
+                    {
+                      unsigned int line_j = j / fe.n_dofs_per_line();
+                      unsigned int tmp_j  = j % degree;
+
+                      if (direction_child[line_i] != direction_mother[line_j])
+                        {
+                          if ((tmp_i + tmp_j) % 2 == 1)
+                            { // anti-symmetric
+                              constraints_matrix[i][j] *= -1.0;
+                            }
+                        }
+                    }
+                }
+
+              // Note:
+              // We need to keep in mind that, if we use a FE_System
+              // with multiple FE_NedelecSZ blocks inside, we need
+              // to consider, that fe.n_dofs_per_line() depends on the number
+              // of FE_NedelecSZ blocks used.
+              const unsigned int n_blocks = fe.n_dofs_per_line() / degree;
+
+              // Edge - Face
+              // Interior edges: for x-direction
+              for (unsigned int i = 0; i < 2 * fe.n_dofs_per_line(); ++i)
+                {
+                  unsigned int line_i = i / fe.n_dofs_per_line();
+                  unsigned int tmp_i  = i % degree;
+
+                  unsigned int start_j =
+                    n_lines_on_mother * fe.n_dofs_per_line();
+
+                  for (unsigned int block = 0; block < n_blocks; ++block)
+                    {
+                      // Type 1:
+                      for (unsigned int jy = 0; jy < degree - 1; ++jy)
+                        for (unsigned int jx = 0; jx < degree - 1; ++jx)
+                          {
+                            unsigned int j = start_j + jx + (jy * (degree - 1));
+                            if (direction_child_intern[line_i] != mother_flip_y)
+                              {
+                                if ((jy + tmp_i) % 2 == 0)
+                                  { // anti-symmetric case
+                                    constraints_matrix[i][j] *= -1.0;
+                                  }
+                              }
+                          }
+
+                      start_j += (degree - 1) * (degree - 1);
+
+                      // Type 2:
+                      for (unsigned int jy = 0; jy < degree - 1; ++jy)
+                        for (unsigned int jx = 0; jx < degree - 1; ++jx)
+                          {
+                            unsigned int j = start_j + jx + (jy * (degree - 1));
+
+                            if (direction_child_intern[line_i] != mother_flip_y)
+                              {
+                                if ((jy + tmp_i) % 2 == 0)
+                                  { // anti-symmetric case
+                                    constraints_matrix[i][j] *= -1.0;
+                                  }
+                              }
+                          }
+                      start_j += (degree - 1) * (degree - 1);
+
+                      // Type 3.1:
+                      // nothing to do
+                      start_j += degree - 1;
+
+                      // Type 3.2:
+                      // nothing to do
+                      start_j += degree - 1;
+                    }
+                }
+
+              // Interior edges: for y-direction
+              for (unsigned int i = 2 * fe.n_dofs_per_line();
+                   i < 4 * fe.n_dofs_per_line();
+                   i++)
+                {
+                  unsigned int line_i = i / fe.n_dofs_per_line();
+                  unsigned int tmp_i  = i % degree;
+
+                  unsigned int start_j =
+                    n_lines_on_mother * fe.n_dofs_per_line();
+
+                  for (unsigned int block = 0; block < n_blocks; block++)
+                    {
+                      // Type 1:
+                      for (unsigned int jy = 0; jy < degree - 1; ++jy)
+                        for (unsigned int jx = 0; jx < degree - 1; ++jx)
+                          {
+                            unsigned int j = start_j + jx + (jy * (degree - 1));
+                            if (direction_child_intern[line_i] != mother_flip_x)
+                              {
+                                if ((jx + tmp_i) % 2 == 0)
+                                  { // anti-symmetric case
+                                    constraints_matrix[i][j] *= -1.0;
+                                  }
+                              }
+                          }
+
+                      start_j += (degree - 1) * (degree - 1);
+
+                      // Type 2:
+                      for (unsigned int jy = 0; jy < degree - 1; ++jy)
+                        for (unsigned int jx = 0; jx < degree - 1; ++jx)
+                          {
+                            unsigned int j = start_j + jx + (jy * (degree - 1));
+                            if (direction_child_intern[line_i] != mother_flip_x)
+                              {
+                                if ((jx + tmp_i) % 2 == 0)
+                                  { // anti-symmetric case
+                                    constraints_matrix[i][j] *= -1.0;
+                                  }
+                              }
+                          }
+                      start_j += (degree - 1) * (degree - 1);
+
+                      // Type 3.1:
+                      // nothing to do
+                      start_j += degree - 1;
+
+                      // Type 3.2:
+                      // nothing to do
+                      start_j += degree - 1;
+                    }
+                }
+
+              // Face - Face
+              unsigned int degree_square = (degree - 1) * (degree - 1);
+              {
+                // Face
+                unsigned int i = n_lines_on_children * fe.n_dofs_per_line();
+                for (unsigned int child_face = 0;
+                     child_face < n_children_per_face;
+                     ++child_face)
+                  for (unsigned int block = 0; block < n_blocks; ++block)
+                    {
+                      unsigned int block_size = fe.n_dofs_per_quad() / n_blocks;
+
+                      // check if the counting of the DoFs is correct:
+                      Assert((block == 0 &&
+                              i != n_lines_on_children * fe.n_dofs_per_line() +
+                                     child_face * fe.n_dofs_per_quad()) ==
+                               false,
+                             ExcInternalError());
+
+                      // Type 1:
+                      for (unsigned int iy = 0; iy < degree - 1; ++iy)
+                        for (unsigned int ix = 0; ix < degree - 1; ++ix)
+                          {
+                            // Type 1 on mother:
+                            unsigned int j =
+                              n_lines_on_mother * fe.n_dofs_per_line() +
+                              block * block_size;
+                            for (unsigned int jy = 0; jy < degree - 1; ++jy)
+                              for (unsigned int jx = 0; jx < degree - 1; ++jx)
+                                {
+                                  if (child_flip_x[child_face] !=
+                                      mother_flip_x) //  x - direction (x-flip)
+                                    {
+                                      if ((ix + jx) % 2 == 1)
+                                        { // anti-symmetric in x
+                                          constraints_matrix[i][j] *= -1.0;
+                                        }
+                                    }
+
+                                  if (child_flip_y[child_face] !=
+                                      mother_flip_y) // y - direction (y-flip)
+                                    {
+                                      if ((iy + jy) % 2 == 1)
+                                        { // anti-symmetric in y
+                                          constraints_matrix[i][j] *= -1.0;
+                                        }
+                                    }
+
+                                  j++;
+                                }
+                            i++;
+                          }
+
+                      // Type 2:
+                      for (unsigned int iy = 0; iy < degree - 1; ++iy)
+                        for (unsigned int ix = 0; ix < degree - 1; ++ix)
+                          {
+                            // Type 2 on mother:
+                            unsigned int j =
+                              n_lines_on_mother * fe.n_dofs_per_line() +
+                              degree_square + block * block_size;
+                            for (unsigned int jy = 0; jy < degree - 1; ++jy)
+                              for (unsigned int jx = 0; jx < degree - 1; ++jx)
+                                {
+                                  if (child_flip_x[child_face] !=
+                                      mother_flip_x) //  x - direction (x-flip)
+                                    {
+                                      if ((ix + jx) % 2 == 1)
+                                        { // anti-symmetric in x
+                                          constraints_matrix[i][j] *= -1.0;
+                                        }
+                                    }
+
+                                  if (child_flip_y[child_face] !=
+                                      mother_flip_y) // y - direction (y-flip)
+                                    {
+                                      if ((iy + jy) % 2 == 1)
+                                        { // anti-symmetric in y
+                                          constraints_matrix[i][j] *= -1.0;
+                                        }
+                                    }
+
+                                  j++;
+                                }
+
+                            i++;
+                          }
+
+
+                      // Type 3 (y):
+                      for (unsigned int iy = 0; iy < degree - 1; ++iy)
+                        {
+                          // Type 2 on mother:
+                          unsigned int j =
+                            n_lines_on_mother * fe.n_dofs_per_line() +
+                            degree_square + block * block_size;
+                          for (unsigned int jy = 0; jy < degree - 1; ++jy)
+                            for (unsigned int jx = 0; jx < degree - 1; ++jx)
+                              {
+                                if (child_flip_x[child_face] !=
+                                    mother_flip_x) //  x - direction (x-flip)
+                                  {
+                                    if ((jx) % 2 == 0)
+                                      { // anti-symmetric in x
+                                        constraints_matrix[i][j] *= -1.0;
+                                      }
+                                  }
+
+                                if (child_flip_y[child_face] !=
+                                    mother_flip_y) // y - direction (y-flip)
+                                  {
+                                    if ((iy + jy) % 2 == 1)
+                                      { // anti-symmetric in y
+                                        constraints_matrix[i][j] *= -1.0;
+                                      }
+                                  }
+
+                                j++;
+                              }
+
+                          // Type 3 on mother:
+                          j = n_lines_on_mother * fe.n_dofs_per_line() +
+                              2 * degree_square + block * block_size;
+                          for (unsigned int jy = 0; jy < degree - 1; ++jy)
+                            {
+                              if (child_flip_y[child_face] !=
+                                  mother_flip_y) // y - direction (y-flip)
+                                {
+                                  if ((iy + jy) % 2 == 1)
+                                    { // anti-symmetric in y
+                                      constraints_matrix[i][j] *= -1.0;
+                                    }
+                                }
+
+                              j++;
+                            }
+                          i++;
+                        }
+
+                      // Type 3 (x):
+                      for (unsigned int ix = 0; ix < degree - 1; ++ix)
+                        {
+                          // Type 2 on mother:
+                          unsigned int j =
+                            n_lines_on_mother * fe.n_dofs_per_line() +
+                            degree_square + block * block_size;
+                          for (unsigned int jy = 0; jy < degree - 1; ++jy)
+                            for (unsigned int jx = 0; jx < degree - 1; ++jx)
+                              {
+                                if (child_flip_x[child_face] !=
+                                    mother_flip_x) //  x - direction (x-flip)
+                                  {
+                                    if ((ix + jx) % 2 == 1)
+                                      { // anti-symmetric in x
+                                        constraints_matrix[i][j] *= -1.0;
+                                      }
+                                  }
+
+                                if (child_flip_y[child_face] !=
+                                    mother_flip_y) // y - direction (y-flip)
+                                  {
+                                    if ((jy) % 2 == 0)
+                                      { // anti-symmetric in y
+                                        constraints_matrix[i][j] *= -1.0;
+                                      }
+                                  }
+
+                                j++;
+                              } // rof: Dof j
+
+                          // Type 3 on mother:
+                          j = n_lines_on_mother * fe.n_dofs_per_line() +
+                              2 * degree_square + (degree - 1) +
+                              block * block_size;
+                          for (unsigned int jx = 0; jx < degree - 1; ++jx)
+                            {
+                              if (child_flip_x[child_face] !=
+                                  mother_flip_x) //  x - direction (x-flip)
+                                {
+                                  if ((ix + jx) % 2 == 1)
+                                    { // anti-symmetric in x
+                                      constraints_matrix[i][j] *= -1.0;
+                                    }
+                                }
+
+                              j++;
+                            }
+                          i++;
+                        }
+                    }
+              }
+
+              // Next, after we have adapted the signs in the constraint matrix,
+              // based on the directions of the edges, we need to modify the
+              // constraint matrix based on the orientation of the faces (i.e.
+              // if x and y direction are exchanged on the face)
+
+              // interior edges:
+              for (unsigned int i = 0;
+                   i < n_internal_lines_on_children * fe.n_dofs_per_line();
+                   ++i)
+                {
+                  // check if x and y are permuted on the parent's face
+                  if (mother_flip_xy)
+                    {
+                      // copy the constraints:
+                      std::vector<double> constraints_matrix_old(
+                        dofs_on_mother.size(), 0);
+                      for (unsigned int j = 0; j < dofs_on_mother.size(); ++j)
+                        {
+                          constraints_matrix_old[j] = constraints_matrix[i][j];
+                        }
+
+                      unsigned int j_start =
+                        n_lines_on_mother * fe.n_dofs_per_line();
+                      for (unsigned block = 0; block < n_blocks; block++)
+                        {
+                          // Type 1
+                          for (unsigned int jy = 0; jy < degree - 1; ++jy)
+                            for (unsigned int jx = 0; jx < degree - 1; ++jx)
+                              {
+                                unsigned int j_old =
+                                  j_start + jx + (jy * (degree - 1));
+                                unsigned int j_new =
+                                  j_start + jy + (jx * (degree - 1));
+                                constraints_matrix[i][j_new] =
+                                  constraints_matrix_old[j_old];
+                              }
+                          j_start += degree_square;
+
+                          // Type 2
+                          for (unsigned int jy = 0; jy < degree - 1; ++jy)
+                            for (unsigned int jx = 0; jx < degree - 1; ++jx)
+                              {
+                                unsigned int j_old =
+                                  j_start + jx + (jy * (degree - 1));
+                                unsigned int j_new =
+                                  j_start + jy + (jx * (degree - 1));
+                                constraints_matrix[i][j_new] =
+                                  -constraints_matrix_old[j_old];
+                              }
+                          j_start += degree_square;
+
+                          // Type 3
+                          for (unsigned int j = j_start;
+                               j < j_start + (degree - 1);
+                               j++)
+                            {
+                              constraints_matrix[i][j] =
+                                constraints_matrix_old[j + (degree - 1)];
+                              constraints_matrix[i][j + (degree - 1)] =
+                                constraints_matrix_old[j];
+                            }
+                          j_start += 2 * (degree - 1);
+                        }
+                    }
+                }
+
+              {
+                // faces:
+                const unsigned int deg = degree - 1;
+
+                // copy the constraints
+                std::vector<std::vector<double>> constraints_matrix_old(
+                  4 * fe.n_dofs_per_quad(),
+                  std::vector<double>(fe.n_dofs_per_quad(), 0));
+                for (unsigned int i = 0;
+                     i < n_children_per_face * fe.n_dofs_per_quad();
+                     ++i)
+                  for (unsigned int j = 0; j < fe.n_dofs_per_quad(); ++j)
+                    constraints_matrix_old[i][j] = constraints_matrix
+                      [i + (n_lines_on_children * fe.n_dofs_per_line())]
+                      [j + (n_lines_on_mother * fe.n_dofs_per_line())];
+
+                // permute rows (on child)
+                for (unsigned int child = 0; child < n_children_per_face;
+                     ++child)
+                  {
+                    if (!child_flip_xy[child])
+                      continue;
+
+                    unsigned int i_start_new =
+                      n_lines_on_children * fe.n_dofs_per_line() +
+                      (child * fe.n_dofs_per_quad());
+                    unsigned int i_start_old = child * fe.n_dofs_per_quad();
+
+                    unsigned int j_start =
+                      n_lines_on_mother * fe.n_dofs_per_line();
+
+                    for (unsigned int block = 0; block < n_blocks; block++)
+                      {
+                        // Type 1:
+                        for (unsigned int ix = 0; ix < deg; ++ix)
+                          {
+                            for (unsigned int iy = 0; iy < deg; ++iy)
+                              {
+                                for (unsigned int j = 0;
+                                     j < fe.n_dofs_per_quad();
+                                     ++j)
+                                  constraints_matrix[i_start_new + iy +
+                                                     (ix * deg)][j + j_start] =
+                                    constraints_matrix_old[i_start_old + ix +
+                                                           (iy * deg)][j];
+                              }
+                          }
+                        i_start_new += deg * deg;
+                        i_start_old += deg * deg;
+
+                        // Type 2:
+                        for (unsigned int ix = 0; ix < deg; ++ix)
+                          {
+                            for (unsigned int iy = 0; iy < deg; ++iy)
+                              {
+                                for (unsigned int j = 0;
+                                     j < fe.n_dofs_per_quad();
+                                     j++)
+                                  constraints_matrix[i_start_new + iy +
+                                                     (ix * deg)][j + j_start] =
+                                    -constraints_matrix_old[i_start_old + ix +
+                                                            (iy * deg)][j];
+                              }
+                          }
+                        i_start_new += deg * deg;
+                        i_start_old += deg * deg;
+
+                        // Type 3:
+                        for (unsigned int ix = 0; ix < deg; ++ix)
+                          {
+                            for (unsigned int j = 0; j < fe.n_dofs_per_quad();
+                                 ++j)
+                              constraints_matrix[i_start_new + ix][j +
+                                                                   j_start] =
+                                constraints_matrix_old[i_start_old + ix + deg]
+                                                      [j];
+                            for (unsigned int j = 0; j < fe.n_dofs_per_quad();
+                                 ++j)
+                              constraints_matrix[i_start_new + ix +
+                                                 deg][j + j_start] =
+                                constraints_matrix_old[i_start_old + ix][j];
+                          } // rof: ix
+
+                        i_start_new += 2 * deg;
+                        i_start_old += 2 * deg;
+                      }
+                  }
+
+                // update the constraints_old
+                for (unsigned int i = 0;
+                     i < n_children_per_face * fe.n_dofs_per_quad();
+                     i++)
+                  for (unsigned int j = 0; j < fe.n_dofs_per_quad(); j++)
+                    constraints_matrix_old[i][j] = constraints_matrix
+                      [i + (n_lines_on_children * fe.n_dofs_per_line())]
+                      [j + (n_lines_on_mother * fe.n_dofs_per_line())];
+
+                // Mother
+                if (mother_flip_xy)
+                  {
+                    unsigned int i_start =
+                      n_lines_on_children * fe.n_dofs_per_line();
+
+                    unsigned int j_start_new =
+                      n_lines_on_mother * fe.n_dofs_per_line();
+                    unsigned int j_start_old = 0;
+
+                    for (unsigned int block = 0; block < n_blocks; ++block)
+                      {
+                        // Type 1:
+                        for (unsigned int jx = 0; jx < deg; ++jx)
+                          {
+                            for (unsigned int jy = 0; jy < deg; ++jy)
+                              {
+                                for (unsigned int i = 0;
+                                     i <
+                                     n_children_per_face * fe.n_dofs_per_quad();
+                                     ++i)
+                                  constraints_matrix[i + i_start][j_start_new +
+                                                                  jy +
+                                                                  (jx * deg)] =
+                                    constraints_matrix_old[i][j_start_old + jx +
+                                                              (jy * deg)];
+                              }
+                          }
+                        j_start_new += deg * deg;
+                        j_start_old += deg * deg;
+
+                        // Type 2:
+                        for (unsigned int jx = 0; jx < deg; ++jx)
+                          {
+                            for (unsigned int jy = 0; jy < deg; ++jy)
+                              {
+                                for (unsigned int i = 0;
+                                     i <
+                                     n_children_per_face * fe.n_dofs_per_quad();
+                                     ++i)
+                                  constraints_matrix[i + i_start][j_start_new +
+                                                                  jy +
+                                                                  (jx * deg)] =
+                                    -constraints_matrix_old[i][j_start_old +
+                                                               jx + (jy * deg)];
+                              }
+                          }
+                        j_start_new += deg * deg;
+                        j_start_old += deg * deg;
+
+                        // Type 3:
+                        for (unsigned int jx = 0; jx < deg; ++jx)
+                          {
+                            for (unsigned int i = 0;
+                                 i < n_children_per_face * fe.n_dofs_per_quad();
+                                 ++i)
+                              {
+                                constraints_matrix[i + i_start][j_start_new +
+                                                                jx] =
+                                  constraints_matrix_old[i][j_start_old + jx +
+                                                            deg];
+                                constraints_matrix[i + i_start][j_start_new +
+                                                                jx + deg] =
+                                  constraints_matrix_old[i][j_start_old + jx];
+                              }
+                          }
+                        j_start_new += 2 * deg;
+                        j_start_old += 2 * deg;
+                      }
+                  }
+              }
+
+              // For each row in the AffineConstraints object for
+              // this line, add the constraint. We split this into the different
+              // cases.
+
+              // internal edges:
+              for (unsigned int line = 0; line < n_internal_lines_on_children;
+                   ++line)
+                {
+                  unsigned int row_start = line * fe.n_dofs_per_line();
+
+                  for (unsigned int row = 0; row < fe.n_dofs_per_line(); ++row)
+                    {
+                      constraints.add_line(dofs_on_children[row_start + row]);
+                      for (unsigned int i = 0; i < dofs_on_mother.size(); ++i)
+                        {
+                          constraints.add_entry(
+                            dofs_on_children[row_start + row],
+                            dofs_on_mother[i],
+                            constraints_matrix[row_start + row][i]);
+                        }
+                      constraints.set_inhomogeneity(
+                        dofs_on_children[row_start + row], 0.);
+                    }
+                }
+
+              // Exterior edges
+              for (unsigned int line = 0; line < n_external_lines_on_children;
+                   ++line)
+                {
+                  unsigned int row_start =
+                    (4 * fe.n_dofs_per_line()) + (line * fe.n_dofs_per_line());
+
+                  for (unsigned int row = 0; row < fe.n_dofs_per_line(); ++row)
+                    {
+                      constraints.add_line(dofs_on_children[row_start + row]);
+                      for (unsigned int i = 0; i < dofs_on_mother.size(); ++i)
+                        {
+                          constraints.add_entry(
+                            dofs_on_children[row_start + row],
+                            dofs_on_mother[i],
+                            constraints_matrix[row_start + row][i]);
+                        }
+                      constraints.set_inhomogeneity(
+                        dofs_on_children[row_start + row], 0.);
+                    }
+                }
+
+              // Faces:
+              for (unsigned int f = 0; f < n_children_per_face; ++f)
+                {
+                  unsigned int row_start =
+                    (n_lines_on_children * fe.n_dofs_per_line()) +
+                    (f * fe.n_dofs_per_quad());
+
+                  for (unsigned int row = 0; row < fe.n_dofs_per_quad(); ++row)
+                    {
+                      constraints.add_line(dofs_on_children[row_start + row]);
+
+                      for (unsigned int i = 0; i < dofs_on_mother.size(); ++i)
+                        {
+                          constraints.add_entry(
+                            dofs_on_children[row_start + row],
+                            dofs_on_mother[i],
+                            constraints_matrix[row_start + row][i]);
+                        }
+
+                      constraints.set_inhomogeneity(
+                        dofs_on_children[row_start + row], 0.);
+                    }
+                }
+            }
+        }
+    }
+
+
     template <int dim, int spacedim, typename number>
     void
     make_hp_hanging_node_constraints(
@@ -1837,11 +3111,16 @@ namespace DoFTools
              "The given DoFHandler does not have any DoFs. Did you forget to "
              "call dof_handler.distribute_dofs()?"));
 
-    // Decide whether to use the new or old make_hanging_node_constraints
+    // Decide whether to use make_hanging_node_constraints_nedelec,
+    // the new or old make_hanging_node_constraints
     // function. If all the FiniteElement or all elements in a FECollection
     // support the new face constraint matrix, the new code will be used.
     // Otherwise, the old implementation is used for the moment.
-    if (dof_handler.get_fe_collection().hp_constraints_are_implemented())
+    if (dof_handler.get_fe().get_name().find("FE_NedelecSZ") !=
+        std::string::npos)
+      internal::make_hanging_node_constraints_nedelec(
+        dof_handler, constraints, std::integral_constant<int, dim>());
+    else if (dof_handler.get_fe_collection().hp_constraints_are_implemented())
       internal::make_hp_hanging_node_constraints(dof_handler, constraints);
     else
       internal::make_oldstyle_hanging_node_constraints(
index dc6e6fa5f87947f0464d3e659e5f565139e558ca..134107d4d0b6914052528d60b050adf466251bd7 100644 (file)
@@ -15,6 +15,7 @@
 
 
 #include <deal.II/fe/fe_nedelec_sz.h>
+#include <deal.II/fe/fe_tools.h>
 
 #include <memory>
 
@@ -45,6 +46,104 @@ FE_NedelecSZ<dim, spacedim>::FE_NedelecSZ(const unsigned int order)
 
   // Generate the 1-D polynomial basis.
   create_polynomials(order);
+
+  // Compute the face embedding.
+  FullMatrix<double> face_embeddings[GeometryInfo<dim>::max_children_per_face];
+
+  // The implementation assumes that all faces have the same
+  // number of DoFs.
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 1;
+  for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
+    {
+      face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
+                                this->n_dofs_per_face(face_no));
+    }
+
+  FETools::compute_face_embedding_matrices<dim, double>(
+    *this, face_embeddings, 0, 0, 1.e-15 * std::exp(std::pow(order, 1.075)));
+
+  switch (dim)
+    {
+      case 1:
+        {
+          this->interface_constraints.reinit(0, 0);
+          break;
+        }
+
+      case 2:
+        {
+          this->interface_constraints.reinit(2 * this->n_dofs_per_face(face_no),
+                                             this->n_dofs_per_face(face_no));
+          for (unsigned int i = 0; i < GeometryInfo<2>::max_children_per_face;
+               ++i)
+            {
+              for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
+                {
+                  for (unsigned int k = 0; k < this->n_dofs_per_face(face_no);
+                       ++k)
+                    {
+                      this->interface_constraints(
+                        i * this->n_dofs_per_face(face_no) + j, k) =
+                        face_embeddings[i](j, k);
+                    }
+                }
+            }
+          break;
+        }
+
+      case 3:
+        {
+          this->interface_constraints.reinit(
+            4 * (this->n_dofs_per_face(face_no) - this->degree),
+            this->n_dofs_per_face(face_no));
+          unsigned int target_row = 0;
+          for (unsigned int i = 0; i < 2; ++i)
+            for (unsigned int j = this->degree; j < 2 * this->degree;
+                 ++j, ++target_row)
+              for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
+                this->interface_constraints(target_row, k) =
+                  face_embeddings[2 * i](j, k);
+          for (unsigned int i = 0; i < 2; ++i)
+            for (unsigned int j = 3 * this->degree;
+                 j < GeometryInfo<3>::lines_per_face * this->degree;
+                 ++j, ++target_row)
+              for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
+                this->interface_constraints(target_row, k) =
+                  face_embeddings[i](j, k);
+          for (unsigned int i = 0; i < 2; ++i)
+            for (unsigned int j = 0; j < 2; ++j)
+              for (unsigned int k = i * this->degree;
+                   k < (i + 1) * this->degree;
+                   ++k, ++target_row)
+                for (unsigned int l = 0; l < this->n_dofs_per_face(face_no);
+                     ++l)
+                  this->interface_constraints(target_row, l) =
+                    face_embeddings[i + 2 * j](k, l);
+          for (unsigned int i = 0; i < 2; ++i)
+            for (unsigned int j = 0; j < 2; ++j)
+              for (unsigned int k = (i + 2) * this->degree;
+                   k < (i + 3) * this->degree;
+                   ++k, ++target_row)
+                for (unsigned int l = 0; l < this->n_dofs_per_face(face_no);
+                     ++l)
+                  this->interface_constraints(target_row, l) =
+                    face_embeddings[2 * i + j](k, l);
+          for (unsigned int i = 0; i < GeometryInfo<3>::max_children_per_face;
+               ++i)
+            for (unsigned int j =
+                   GeometryInfo<3>::lines_per_face * this->degree;
+                 j < this->n_dofs_per_face(face_no);
+                 ++j, ++target_row)
+              for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
+                this->interface_constraints(target_row, k) =
+                  face_embeddings[i](j, k);
+          break;
+        }
+
+      default:
+        Assert(false, ExcNotImplemented());
+    }
 }
 
 
@@ -1531,6 +1630,16 @@ FE_NedelecSZ<dim, spacedim>::fill_edge_values(
                   unsigned int v0_glob = cell->vertex_index(v0_loc);
                   unsigned int v1_glob = cell->vertex_index(v1_loc);
 
+                  // Check for hanging edges on the current face. If we
+                  // encounter a hanging edge, we use the vertex indices
+                  // from the parent.
+                  if (cell->face(m)->at_boundary() == false)
+                    if (cell->neighbor_is_coarser(m))
+                      {
+                        v0_glob = cell->parent()->vertex_index(v0_loc);
+                        v1_glob = cell->parent()->vertex_index(v1_loc);
+                      }
+
                   if (v0_glob > v1_glob)
                     {
                       // Opposite to global numbering on our reference element
@@ -1779,6 +1888,142 @@ FE_NedelecSZ<dim, spacedim>::fill_edge_values(
                     }
                 }
 
+              // Check for hanging faces. If we encounter hanging faces,
+              // we use the vertex indices from the parent.
+              for (unsigned int f = 0; f < 6; ++f)
+                if (cell->face(f)->at_boundary() == false)
+                  if (cell->neighbor_is_coarser(f))
+                    for (unsigned int m = 0; m < 4; ++m)
+                      {
+                        unsigned int parent_m =
+                          GeometryInfo<dim>::face_to_cell_lines(f, m);
+
+                        unsigned int v0_loc =
+                          GeometryInfo<dim>::line_to_cell_vertices(parent_m, 0);
+                        unsigned int v1_loc =
+                          GeometryInfo<dim>::line_to_cell_vertices(parent_m, 1);
+
+                        unsigned int v0_glob =
+                          cell->parent()->vertex_index(v0_loc);
+                        unsigned int v1_glob =
+                          cell->parent()->vertex_index(v1_loc);
+
+                        if (v0_glob > v1_glob)
+                          {
+                            // Opposite to global numbering on our reference
+                            // element
+                            edge_sign[parent_m] = -1.0;
+                          }
+                        else
+                          {
+                            // Aligns with global numbering on our reference
+                            // element.
+                            edge_sign[parent_m] = 1.0;
+                          }
+                      }
+
+              // Next, we cover the case where we encounter a hanging edge but
+              // not a hanging face. This is always the case if the cell
+              // currently considered shares all faces with cells of the same
+              // refinement level but shares one edge with a coarser cell,
+              // e.g.:
+              //            *----*----*---------*
+              //          /    /    /         / |
+              //        *----*----*         /   |
+              //      /    /    /         /     |
+              //    *----*----*---------*       *
+              //    |    |    |         |     / |
+              //    *----*----*         |   *   *
+              //    |    |    |         | / | / |
+              //    *----*----*----*----*   *   *
+              //    |    |    |    |    | / | /
+              //    *----*----*----*----*   *
+              //    |    |    |    |    | /
+              //    *----*----*----*----*
+              // where the cell at the left bottom is the currently
+              // considered cell.
+              //
+              // In that case, we determine the direction of the hanging edge
+              // based on the vertex indices from the parent cell.
+              //
+              // Note: We assume here that at most eight cells are adjacent to
+              // a single vertex.
+              std::vector<unsigned int> adjacent_faces = {2, 2, 4, 4, 0, 0};
+              for (unsigned int f = 0; f < 6; ++f)
+                {
+                  if (!cell->face(f)->at_boundary() &&
+                      !cell->face(adjacent_faces[f])->at_boundary())
+                    if (!cell->neighbor_is_coarser(f) &&
+                        !cell->neighbor_is_coarser(adjacent_faces[f]))
+                      if (!cell->neighbor(f)
+                             ->face(adjacent_faces[f])
+                             ->at_boundary())
+                        if (cell->neighbor(f)->neighbor_is_coarser(
+                              adjacent_faces[f]))
+                          {
+                            unsigned int parent_m =
+                              GeometryInfo<dim>::face_to_cell_lines(f, 0);
+                            unsigned int v0_loc =
+                              GeometryInfo<dim>::line_to_cell_vertices(parent_m,
+                                                                       0);
+                            unsigned int v1_loc =
+                              GeometryInfo<dim>::line_to_cell_vertices(parent_m,
+                                                                       1);
+                            unsigned int v0_glob =
+                              cell->parent()->vertex_index(v0_loc);
+                            unsigned int v1_glob =
+                              cell->parent()->vertex_index(v1_loc);
+                            if (v0_glob > v1_glob)
+                              {
+                                // Opposite to global numbering on our reference
+                                // element
+                                edge_sign[parent_m] = -1.0;
+                              }
+                            else
+                              {
+                                // Aligns with global numbering on our reference
+                                // element.
+                                edge_sign[parent_m] = 1.0;
+                              }
+                          }
+
+                  if (!cell->face(f)->at_boundary() &&
+                      !cell->face(adjacent_faces[f] + 1)->at_boundary())
+                    if (!cell->neighbor_is_coarser(f) &&
+                        !cell->neighbor_is_coarser(adjacent_faces[f] + 1))
+                      if (!cell->neighbor(f)
+                             ->face(adjacent_faces[f] + 1)
+                             ->at_boundary())
+                        if (cell->neighbor(f)->neighbor_is_coarser(
+                              adjacent_faces[f] + 1))
+                          {
+                            unsigned int parent_m =
+                              GeometryInfo<dim>::face_to_cell_lines(f, 1);
+                            unsigned int v0_loc =
+                              GeometryInfo<dim>::line_to_cell_vertices(parent_m,
+                                                                       0);
+                            unsigned int v1_loc =
+                              GeometryInfo<dim>::line_to_cell_vertices(parent_m,
+                                                                       1);
+                            unsigned int v0_glob =
+                              cell->parent()->vertex_index(v0_loc);
+                            unsigned int v1_glob =
+                              cell->parent()->vertex_index(v1_loc);
+                            if (v0_glob > v1_glob)
+                              {
+                                // Opposite to global numbering on our reference
+                                // element
+                                edge_sign[parent_m] = -1.0;
+                              }
+                            else
+                              {
+                                // Aligns with global numbering on our reference
+                                // element.
+                                edge_sign[parent_m] = 1.0;
+                              }
+                          }
+                }
+
               // Define \sigma_{m} = sigma_{e^{m}_{1}} - sigma_{e^{m}_{2}}
               //        \lambda_{m} = \lambda_{e^{m}_{1}} + \lambda_{e^{m}_{2}}
               //
@@ -2082,26 +2327,58 @@ FE_NedelecSZ<dim, spacedim>::fill_face_values(
 
           const unsigned int
             vertices_adjacent_on_face[GeometryInfo<3>::vertices_per_face][2] = {
-              {1, 2}, {0, 3}, {0, 3}, {1, 2}};
+              {1, 2}, {0, 3}, {3, 0}, {2, 1}};
 
           for (unsigned int m = 0; m < faces_per_cell; ++m)
             {
+              // Check, if we are on a hanging face.
+              bool cell_has_coarser_neighbor = false;
+              if (cell->face(m)->at_boundary() == false)
+                if (cell->neighbor_is_coarser(m))
+                  cell_has_coarser_neighbor = true;
+
               // Find the local vertex on this face with the highest global
               // numbering. This is f^m_0.
-              unsigned int current_max  = 0;
-              unsigned int current_glob = cell->vertex_index(
-                GeometryInfo<dim>::face_to_cell_vertices(m, 0));
-              for (unsigned int v = 1; v < vertices_per_face; ++v)
+              unsigned int current_max = 0;
+
+              // We start with the hanging face case, where the face
+              // orientation is determined based on the vertex indices
+              // of the parent cell.
+              if (cell_has_coarser_neighbor)
                 {
-                  if (current_glob <
-                      cell->vertex_index(
-                        GeometryInfo<dim>::face_to_cell_vertices(m, v)))
+                  unsigned int current_glob = cell->parent()->vertex_index(
+                    GeometryInfo<dim>::face_to_cell_vertices(m, 0));
+                  for (unsigned int v = 1; v < vertices_per_face; ++v)
+                    {
+                      if (current_glob <
+                          cell->parent()->vertex_index(
+                            GeometryInfo<dim>::face_to_cell_vertices(m, v)))
+                        {
+                          current_max  = v;
+                          current_glob = cell->parent()->vertex_index(
+                            GeometryInfo<dim>::face_to_cell_vertices(m, v));
+                        }
+                    }
+                }
+              // Otherwise, the face orientation is based on its own
+              // vertex indices.
+              else
+                {
+                  unsigned int current_glob = cell->vertex_index(
+                    GeometryInfo<dim>::face_to_cell_vertices(m, 0));
+                  for (unsigned int v = 1; v < vertices_per_face; ++v)
                     {
-                      current_max  = v;
-                      current_glob = cell->vertex_index(
-                        GeometryInfo<dim>::face_to_cell_vertices(m, v));
+                      if (current_glob <
+                          cell->vertex_index(
+                            GeometryInfo<dim>::face_to_cell_vertices(m, v)))
+                        {
+                          current_max  = v;
+                          current_glob = cell->vertex_index(
+                            GeometryInfo<dim>::face_to_cell_vertices(m, v));
+                        }
                     }
                 }
+
               face_orientation[m][0] =
                 GeometryInfo<dim>::face_to_cell_vertices(m, current_max);
 
@@ -2111,29 +2388,63 @@ FE_NedelecSZ<dim, spacedim>::fill_face_values(
 
               // Finally, f^m_1 is the vertex with the greater global numbering
               // of the remaining two local vertices. Then, f^m_3 is the other.
-              if (cell->vertex_index(GeometryInfo<dim>::face_to_cell_vertices(
-                    m, vertices_adjacent_on_face[current_max][0])) >
-                  cell->vertex_index(GeometryInfo<dim>::face_to_cell_vertices(
-                    m, vertices_adjacent_on_face[current_max][1])))
+              // Again, we need to distinguish between the hanging face and the
+              // non-hanging face cases. In the case of hanging faces, we
+              // consider the vertex indices from the parent. Otherwise, we
+              // consider the vertex indices of the face itself.
+              if (cell_has_coarser_neighbor)
                 {
-                  face_orientation[m][1] =
-                    GeometryInfo<dim>::face_to_cell_vertices(
-                      m, vertices_adjacent_on_face[current_max][0]);
-                  face_orientation[m][3] =
-                    GeometryInfo<dim>::face_to_cell_vertices(
-                      m, vertices_adjacent_on_face[current_max][1]);
+                  if (cell->parent()->vertex_index(
+                        GeometryInfo<dim>::face_to_cell_vertices(
+                          m, vertices_adjacent_on_face[current_max][0])) >
+                      cell->parent()->vertex_index(
+                        GeometryInfo<dim>::face_to_cell_vertices(
+                          m, vertices_adjacent_on_face[current_max][1])))
+                    {
+                      face_orientation[m][1] =
+                        GeometryInfo<dim>::face_to_cell_vertices(
+                          m, vertices_adjacent_on_face[current_max][0]);
+                      face_orientation[m][3] =
+                        GeometryInfo<dim>::face_to_cell_vertices(
+                          m, vertices_adjacent_on_face[current_max][1]);
+                    }
+                  else
+                    {
+                      face_orientation[m][1] =
+                        GeometryInfo<dim>::face_to_cell_vertices(
+                          m, vertices_adjacent_on_face[current_max][1]);
+                      face_orientation[m][3] =
+                        GeometryInfo<dim>::face_to_cell_vertices(
+                          m, vertices_adjacent_on_face[current_max][0]);
+                    }
                 }
               else
                 {
-                  face_orientation[m][1] =
-                    GeometryInfo<dim>::face_to_cell_vertices(
-                      m, vertices_adjacent_on_face[current_max][1]);
-                  face_orientation[m][3] =
-                    GeometryInfo<dim>::face_to_cell_vertices(
-                      m, vertices_adjacent_on_face[current_max][0]);
+                  if (cell->vertex_index(
+                        GeometryInfo<dim>::face_to_cell_vertices(
+                          m, vertices_adjacent_on_face[current_max][0])) >
+                      cell->vertex_index(
+                        GeometryInfo<dim>::face_to_cell_vertices(
+                          m, vertices_adjacent_on_face[current_max][1])))
+                    {
+                      face_orientation[m][1] =
+                        GeometryInfo<dim>::face_to_cell_vertices(
+                          m, vertices_adjacent_on_face[current_max][0]);
+                      face_orientation[m][3] =
+                        GeometryInfo<dim>::face_to_cell_vertices(
+                          m, vertices_adjacent_on_face[current_max][1]);
+                    }
+                  else
+                    {
+                      face_orientation[m][1] =
+                        GeometryInfo<dim>::face_to_cell_vertices(
+                          m, vertices_adjacent_on_face[current_max][1]);
+                      face_orientation[m][3] =
+                        GeometryInfo<dim>::face_to_cell_vertices(
+                          m, vertices_adjacent_on_face[current_max][0]);
+                    }
                 }
             }
-
           // Now we know the face orientation on the current cell, we can
           // generate the parameterisation:
           std::vector<std::vector<double>> face_xi_values(
@@ -2211,7 +2522,8 @@ FE_NedelecSZ<dim, spacedim>::fill_face_values(
               // Type-3:
               //
               // \phi^{F_m,3}_{i} = L_{i+2}(\eta_{F_{m}}) \lambda_{F_{m}}
-              // \nabla\xi_{F_{m}} \phi^{F_m,3}_{i+p} = L_{i+2}(\xi_{F_{m}})
+              // \nabla\xi_{F_{m}}
+              // \phi^{F_m,3}_{i+p} = L_{i+2}(\xi_{F_{m}})
               // \lambda_{F_{m}} \nabla\eta_{F_{m}}
               //
               // 0 <= i < degree.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.