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
}
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
}
}
+ 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(
+ 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(
"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(
#include <deal.II/fe/fe_nedelec_sz.h>
+#include <deal.II/fe/fe_tools.h>
#include <memory>
// 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());
+ }
}
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
}
}
+ // 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}}
//
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);
// 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(
// 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.
projects / dealii.git / commitdiff