From 70407c1fb32b0f48baa05e52512b41a718210784 Mon Sep 17 00:00:00 2001 From: Sebastian Kinnewig Date: Fri, 19 May 2023 11:19:44 +0200 Subject: [PATCH] Add support for hanging nodes to the NedelecSZ class. --- .../changes/major/20231128SebastianKinnewig | 3 + include/deal.II/fe/fe_nedelec_sz.h | 11 +- source/dofs/dof_tools_constraints.cc | 1293 ++++++++++++++++- source/fe/fe_nedelec_sz.cc | 370 ++++- 4 files changed, 1640 insertions(+), 37 deletions(-) create mode 100644 doc/news/changes/major/20231128SebastianKinnewig diff --git a/doc/news/changes/major/20231128SebastianKinnewig b/doc/news/changes/major/20231128SebastianKinnewig new file mode 100644 index 0000000000..c936511c0c --- /dev/null +++ b/doc/news/changes/major/20231128SebastianKinnewig @@ -0,0 +1,3 @@ +New: Support for hanging nodes in FE_NedelecSZ. +
+(Sebastian Kinnewig, 2023/11/28) diff --git a/include/deal.II/fe/fe_nedelec_sz.h b/include/deal.II/fe/fe_nedelec_sz.h index de396e6bb7..234fbe9672 100644 --- a/include/deal.II/fe/fe_nedelec_sz.h +++ b/include/deal.II/fe/fe_nedelec_sz.h @@ -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: Resolving the sign conflict * problem for hp-hexahedral Nédélec elements with application to * eddy current problems, 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: Algorithmic realization of the + * solution to the sign conflict problem for hanging nodes on hp-hexahedral + * Nédélec elements, preprint, 2023 (https://arxiv.org/abs/2306.01416). */ template class FE_NedelecSZ : public FiniteElement diff --git a/source/dofs/dof_tools_constraints.cc b/source/dofs/dof_tools_constraints.cc index 9b0066cfc4..bb9907c83a 100644 --- a/source/dofs/dof_tools_constraints.cc +++ b/source/dofs/dof_tools_constraints.cc @@ -602,9 +602,8 @@ namespace DoFTools template void - make_oldstyle_hanging_node_constraints(const DoFHandler<1> &, - AffineConstraints &, - std::integral_constant) + make_hp_hanging_node_constraints(const DoFHandler<1, 2> &, + AffineConstraints &) { // nothing to do for regular dof handlers in 1d } @@ -612,8 +611,9 @@ namespace DoFTools template void - make_hp_hanging_node_constraints(const DoFHandler<1, 2> &, - AffineConstraints &) + make_hanging_node_constraints_nedelec(const dealii::DoFHandler<1, 2> &, + AffineConstraints &, + std::integral_constant) { // nothing to do for regular dof handlers in 1d } @@ -639,6 +639,17 @@ namespace DoFTools } + template + void + make_hanging_node_constraints_nedelec( + const dealii::DoFHandler<1, spacedim> & /*dof_handler*/, + AffineConstraints & /*constraints*/, + std::integral_constant) + { + // nothing to do for dof handlers in 1d + } + + template void make_oldstyle_hanging_node_constraints( @@ -1070,6 +1081,1269 @@ namespace DoFTools + template + void + make_hanging_node_constraints_nedelec( + const DoFHandler &dof_handler, + AffineConstraints &constraints, + std::integral_constant) + { + // 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 face_dof_indices; + std::map> + 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 &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 dof_on_mother_face = + face_dof_indices; + + cell->face(f)->child(0)->get_dof_indices(face_dof_indices); + const std::vector dof_on_child_face_0 = + face_dof_indices; + + cell->face(f)->child(1)->get_dof_indices(face_dof_indices); + const std::vector 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 + void + make_hanging_node_constraints_nedelec( + const DoFHandler &dof_handler, + AffineConstraints &constraints, + std::integral_constant) + { + // 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 dofs_on_mother; + std::vector 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::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::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::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 &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::lines_per_face; + + // number of internal lines of the children; + // for more details see description of the + // GeometryInfo 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::max_children_per_face; + const unsigned int n_children_per_line = + GeometryInfo::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::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::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::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 direction_mother( + GeometryInfo::lines_per_face, false); + for (unsigned int line = 0; + line < GeometryInfo::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 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 direction_child(n_external_lines_on_children, + false); + for (unsigned int line = 0; + line < GeometryInfo::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 child_flip_x(n_children_per_face, false); + std::vector child_flip_y(n_children_per_face, false); + std::vector child_flip_xy(n_children_per_face, false); + const unsigned int + vertices_adjacent_on_face[GeometryInfo::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::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::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> constraints_matrix( + n_lines_on_children * fe.n_dofs_per_line() + + n_children_per_face * fe.n_dofs_per_quad(), + std::vector(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 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> constraints_matrix_old( + 4 * fe.n_dofs_per_quad(), + std::vector(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 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()); + 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( diff --git a/source/fe/fe_nedelec_sz.cc b/source/fe/fe_nedelec_sz.cc index dc6e6fa5f8..134107d4d0 100644 --- a/source/fe/fe_nedelec_sz.cc +++ b/source/fe/fe_nedelec_sz.cc @@ -15,6 +15,7 @@ #include +#include #include @@ -45,6 +46,104 @@ FE_NedelecSZ::FE_NedelecSZ(const unsigned int order) // Generate the 1-D polynomial basis. create_polynomials(order); + + // Compute the face embedding. + FullMatrix face_embeddings[GeometryInfo::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::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( + *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::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::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::face_to_cell_lines(f, m); + + unsigned int v0_loc = + GeometryInfo::line_to_cell_vertices(parent_m, 0); + unsigned int v1_loc = + GeometryInfo::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 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::face_to_cell_lines(f, 0); + unsigned int v0_loc = + GeometryInfo::line_to_cell_vertices(parent_m, + 0); + unsigned int v1_loc = + GeometryInfo::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::face_to_cell_lines(f, 1); + unsigned int v0_loc = + GeometryInfo::line_to_cell_vertices(parent_m, + 0); + unsigned int v1_loc = + GeometryInfo::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::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::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::face_to_cell_vertices(m, v))) + unsigned int current_glob = cell->parent()->vertex_index( + GeometryInfo::face_to_cell_vertices(m, 0)); + for (unsigned int v = 1; v < vertices_per_face; ++v) + { + if (current_glob < + cell->parent()->vertex_index( + GeometryInfo::face_to_cell_vertices(m, v))) + { + current_max = v; + current_glob = cell->parent()->vertex_index( + GeometryInfo::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::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::face_to_cell_vertices(m, v)); + if (current_glob < + cell->vertex_index( + GeometryInfo::face_to_cell_vertices(m, v))) + { + current_max = v; + current_glob = cell->vertex_index( + GeometryInfo::face_to_cell_vertices(m, v)); + } } } + face_orientation[m][0] = GeometryInfo::face_to_cell_vertices(m, current_max); @@ -2111,29 +2388,63 @@ FE_NedelecSZ::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::face_to_cell_vertices( - m, vertices_adjacent_on_face[current_max][0])) > - cell->vertex_index(GeometryInfo::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::face_to_cell_vertices( - m, vertices_adjacent_on_face[current_max][0]); - face_orientation[m][3] = - GeometryInfo::face_to_cell_vertices( - m, vertices_adjacent_on_face[current_max][1]); + if (cell->parent()->vertex_index( + GeometryInfo::face_to_cell_vertices( + m, vertices_adjacent_on_face[current_max][0])) > + cell->parent()->vertex_index( + GeometryInfo::face_to_cell_vertices( + m, vertices_adjacent_on_face[current_max][1]))) + { + face_orientation[m][1] = + GeometryInfo::face_to_cell_vertices( + m, vertices_adjacent_on_face[current_max][0]); + face_orientation[m][3] = + GeometryInfo::face_to_cell_vertices( + m, vertices_adjacent_on_face[current_max][1]); + } + else + { + face_orientation[m][1] = + GeometryInfo::face_to_cell_vertices( + m, vertices_adjacent_on_face[current_max][1]); + face_orientation[m][3] = + GeometryInfo::face_to_cell_vertices( + m, vertices_adjacent_on_face[current_max][0]); + } } else { - face_orientation[m][1] = - GeometryInfo::face_to_cell_vertices( - m, vertices_adjacent_on_face[current_max][1]); - face_orientation[m][3] = - GeometryInfo::face_to_cell_vertices( - m, vertices_adjacent_on_face[current_max][0]); + if (cell->vertex_index( + GeometryInfo::face_to_cell_vertices( + m, vertices_adjacent_on_face[current_max][0])) > + cell->vertex_index( + GeometryInfo::face_to_cell_vertices( + m, vertices_adjacent_on_face[current_max][1]))) + { + face_orientation[m][1] = + GeometryInfo::face_to_cell_vertices( + m, vertices_adjacent_on_face[current_max][0]); + face_orientation[m][3] = + GeometryInfo::face_to_cell_vertices( + m, vertices_adjacent_on_face[current_max][1]); + } + else + { + face_orientation[m][1] = + GeometryInfo::face_to_cell_vertices( + m, vertices_adjacent_on_face[current_max][1]); + face_orientation[m][3] = + GeometryInfo::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> face_xi_values( @@ -2211,7 +2522,8 @@ FE_NedelecSZ::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. -- 2.39.5