void
FE_Nedelec<dim>::initialize_quad_dof_index_permutation_and_sign_change()
{
- // for 1d and 2d, do nothing
- if (dim < 3)
+ // The order of the Nedelec elements equals the tensor degree minus one,
+ // k = n - 1. In the three-dimensional space the Nedelec elements of the
+ // lowermost order, k = 0, have only 12 line (edge) dofs. The Nedelec
+ // elements of the higher orders, k > 0, have 3*(k+1)*(k+2)^2 dofs in
+ // total if dim=3. The dofs in a cell are distributed between lines
+ // (edges), quads (faces), and the hex (the interior of the cell) as the
+ // following:
+ //
+ // 12*(k+1) line dofs; (k+1) dofs per line.
+ // 2*6*k*(k+1) quad dofs; 2*k*(k+1) dofs per quad.
+ // 3*(k+2)^2*(k+1) hex dofs;
+ //
+ // The dofs are indexed in the following order: first all line dofs,
+ // then all quad dofs, and then all hex dofs.
+ //
+ // The line dofs need only sign adjustments. No permutation of line
+ // dofs is needed. The line dofs are treated by
+ // internal::FE_PolyTensor::get_dof_sign_change_nedelec(...)
+ // in fe_poly_tensor.cc.
+ //
+ // The hex dofs need no adjustments: they are not shared between
+ // neighbouring mesh cells.
+ //
+ // The two-dimensional Nedelec finite elements share no quad dofs between
+ // neighbouring mesh cells. The zero-order three-dimensional Nedelec
+ // finite elements have no quad dofs. Consequently, here we treat only
+ // quad dofs of the three-dimensional Nedelec finite elements of the
+ // higher orders, k>0. The questions how the curl looks like in the
+ // higher-dimensional spaces and what does it mean to be curl-conforming
+ // if dim>3 we leave unanswered.
+ //
+ // In this function we need to change some entries in the following two
+ // vectors of tables:
+ // adjust_quad_dof_index_for_face_orientation_table
+ // and
+ // adjust_quad_dof_sign_for_face_orientation_table.
+ // These tables specify the permutations and sign adjustments of the quad
+ // dofs only. The tables are already filled with zeros meaning no
+ // permutations or sign change are required. We need to change some
+ // entries of the tables such that the shape functions that correspond to
+ // the quad dofs and are shared between neighbouring cells have consistent
+ // orientations.
+ //
+ // The swap tables below describe the dof permutations and sign changes
+ // that need to be done. The function
+ // FE_Nedelec<dim>::initialize_quad_dof_index_permutation_and_sign_change()
+ // simply reads the information in the swap tables below and puts it into
+ // tables
+ // adjust_quad_dof_index_for_face_orientation_table
+ // and
+ // adjust_quad_dof_sign_for_face_orientation_table.
+ // A good question is: why don't we put the information into the tables of
+ // deal.II right away? The answer is the following. The information on the
+ // necessary dof permutations and sign changes is derived by plotting the
+ // shape functions and observing them on faces of different orientations.
+ // It is convenient to put the observations first in the format of the
+ // swap tables below and then convert the swap tables into the format used
+ // by deal.II.
+ //
+ // The dofs on a quad are indexed as the following:
+ //
+ // | x0, x1, x2, x3, ..., xk | y0, y1, y2, y3 ..., yk |
+ // | | |
+ // |<------ k*(k+1) --------->|<------ k*(k+1) -------->|
+ // | |
+ // |<------------------- 2*k*(k+1) -------------------->|
+ //
+ // Only one type of dof permutation is needed: swap between two dofs; one
+ // dof being xi, another yj. That is, if x4 is replaced with y7,
+ // then y7 must be replaced with x4. Such swaps can be ordered as
+ // illustrated by the following example:
+ //
+ // *
+ // y0, y9, y1, y2, ..., yk
+ // --------------------------- (swap)
+ // x0, x1, x2, x3, ..., xk
+ // *
+ //
+ // An x-dof below the line is swapped with the corresponding y-dof above
+ // the line. A dof marked by the asterisk must change its sign before the
+ // swap.
+ //
+ // The x-dofs are assumed to have the normal order. There is no need to
+ // encode it. Therefore, the swap tables need to encode the following
+ // information: indices of the y-dofs, the sign change of the x-dofs, and
+ // sign change of the y-dofs. The swap above is encoded as the following:
+ //
+ // swap = { 0, 9, 1, 2, ...., yk, // indices of the y-dofs
+ // 1, 0, 0, 0, ...., 0, // sign change of the x-dofs,
+ // 0, 1, 0, 0, ...., 0}; // sign change of the y-dofs.
+ //
+ // If no swap is needed, -1 is placed instead of the y-dof index.
+ //
+ // Such swaps are assembled into the swap table:
+ //
+ // swap_table = {swap_0, swap_1, ... swap_7};
+ //
+ // Each swap table contains eight swaps - one swap for each possible quad
+ // orientation. The deal.II encodes the orientation of a quad using
+ // three boolean parameters:
+ // face_orientation - true if face is in standard orientation
+ // and false otherwise;
+ // face_rotation - rotation by 90 deg counterclockwise if true;
+ // face_flip - rotation by 180 deg counterclockwise if true.
+ // See the documentation of GeometryInfo<dim>.
+ //
+ // The combined face orientation is computes as
+ // orientation_no = face_flip*4 + face_rotation*2 + face_orientation*1;
+ // See tria_orientation.h.
+ //
+ // The parameter orientation_no (0...7) indexes the swaps in a swap table.
+ //
+ // Nedelec elements of order k have their own swap table, swap_table_k.
+ // Recall, the swap_table_0 is empty as the Nedelec finite elements of the
+ // lowermost order have no quad dofs.
+
+ static const int c_swap_table_0 = 0;
+
+ static const int c_swap_table_1[8][3][2] = { // 0 1
+ {{0, 1}, // 0
+ {0, 0},
+ {0, 0}},
+ {{-1, -1}, // 1
+ {0, 0},
+ {0, 0}},
+ {{-1, -1}, // 2
+ {0, 0},
+ {1, 0}},
+ {{0, 1}, // 3
+ {1, 0},
+ {0, 0}},
+ {{0, 1}, // 4
+ {1, 0},
+ {1, 0}},
+ {{-1, -1}, // 5
+ {1, 0},
+ {1, 0}},
+ {{-1, -1}, // 6
+ {1, 0},
+ {0, 0}},
+ {{0, 1}, // 7
+ {0, 0},
+ {1, 0}}};
+
+ static const int c_swap_table_2[8][3][6] = {// 0 1 2 3 4 5
+ {{0, 3, 1, 4, 2, 5}, // 0
+ {0, 0, 0, 0, 0, 0},
+ {0, 0, 0, 0, 0, 0}},
+ {{-1, -1, -1, -1, -1, -1}, // 1
+ {0, 0, 0, 0, 0, 0},
+ {0, 0, 0, 0, 0, 0}},
+ {{-1, -1, -1, -1, -1, -1}, // 2
+ {0, 1, 0, 1, 0, 1},
+ {1, 0, 1, 1, 0, 1}},
+ {{0, 3, 1, 4, 2, 5}, // 3
+ {1, 1, 0, 0, 1, 1},
+ {0, 0, 0, 1, 1, 1}},
+ {{0, 3, 1, 4, 2, 5}, // 4
+ {1, 0, 0, 1, 1, 0},
+ {1, 0, 1, 0, 1, 0}},
+ {{-1, -1, -1, -1, -1, -1}, // 5
+ {1, 0, 0, 1, 1, 0},
+ {1, 0, 1, 0, 1, 0}},
+ {{-1, -1, -1, -1, -1, -1}, // 6
+ {1, 1, 0, 0, 1, 1},
+ {0, 0, 0, 1, 1, 1}},
+ {{0, 3, 1, 4, 2, 5}, // 7
+ {0, 1, 0, 1, 0, 1},
+ {1, 0, 1, 1, 0, 1}}};
+
+ static const int c_swap_table_3[8][3][12] = {
+ // 0 1 2 3 4 5 6 7 8 9 10 11
+ {{0, 4, 8, 1, 5, 9, 2, 6, 10, 3, 7, 11}, // 0
+ {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+ {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}},
+ {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // 1
+ {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+ {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}},
+ {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // 2
+ {0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0},
+ {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0}},
+ {{0, 4, 8, 1, 5, 9, 2, 6, 10, 3, 7, 11}, // 3
+ {1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0},
+ {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0}},
+ {{0, 4, 8, 1, 5, 9, 2, 6, 10, 3, 7, 11}, // 4
+ {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0},
+ {1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0}},
+ {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // 5
+ {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0},
+ {1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0}},
+ {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // 6
+ {1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0},
+ {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0}},
+ {{0, 4, 8, 1, 5, 9, 2, 6, 10, 3, 7, 11}, // 7
+ {0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0},
+ {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0}}};
+
+ static const int c_swap_table_4[8][3][20] = {
+ // Swap sign_X and sign_Y rows if k=4. Why?...
+ // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
+ // 19
+ {{0, 5, 10, 15, 1, 6, 11, 16, 2, 7,
+ 12, 17, 3, 8, 13, 18, 4, 9, 14, 19}, // 0
+ {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+ {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}},
+ {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
+ -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // 1
+ {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+ {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}},
+ {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
+ -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // 2
+ {1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1},
+ {0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1}},
+ {{0, 5, 10, 15, 1, 6, 11, 16, 2, 7,
+ 12, 17, 3, 8, 13, 18, 4, 9, 14, 19}, // 3
+ {0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1},
+ {1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1}},
+ {{0, 5, 10, 15, 1, 6, 11, 16, 2, 7,
+ 12, 17, 3, 8, 13, 18, 4, 9, 14, 19}, // 4
+ {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0},
+ {1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0}},
+ {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
+ -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // 5
+ {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0},
+ {1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0}},
+ {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
+ -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // 6
+ {0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1},
+ {1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1}},
+ {{0, 5, 10, 15, 1, 6, 11, 16, 2, 7,
+ 12, 17, 3, 8, 13, 18, 4, 9, 14, 19}, // 7
+ {1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1},
+ {0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1}}};
+
+ static const int *swap_table_array[5] = {&c_swap_table_0,
+ &c_swap_table_1[0][0][0],
+ &c_swap_table_2[0][0][0],
+ &c_swap_table_3[0][0][0],
+ &c_swap_table_4[0][0][0]};
+
+ static const int row_length[5] = {0, 2, 6, 12, 20};
+ static const int table_size[5] = {
+ 0, 8 * 3 * 2, 8 * 3 * 6, 8 * 3 * 12, 8 * 3 * 20};
+
+ // Only three-dimensional Nedelec finite elements are treated. The
+ // two-dimensional Nedelec finite elements only need sign adjustments of the
+ // line dofs. These adjustments are done by
+ // internal::FE_PolyTensor::get_dof_sign_change_nedelec(...)
+ // in fe_poly_tensor.cc. The notions of curl and curl-conforming finite
+ // elements in higher-dimensional spaces, dim >3, are somewhat unclear as
+ // curl, strictly peaking, exists only in the three-dimensional space.
+ if (dim != 3)
return;
- // TODO: Implement this for this class
+ const unsigned int k = this->tensor_degree() - 1;
+
+ // The Nedelec finite elements of the lowermost order have no quad dofs.
+ if (k == 0)
+ return;
+
+ // The finite element orders > 4 are not implemented.
+ AssertThrow(k < 5, ExcNotImplemented());
+
+ // TODO: the implementation makes the assumption that all quads have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
+ Assert(
+ this->adjust_quad_dof_index_for_face_orientation_table[0].n_elements() ==
+ this->reference_cell().n_face_orientations(face_no) *
+ this->n_dofs_per_quad(face_no),
+ ExcInternalError());
+
+ Assert(
+ this->adjust_quad_dof_sign_for_face_orientation_table[0].n_elements() ==
+ this->reference_cell().n_face_orientations(face_no) *
+ this->n_dofs_per_quad(face_no),
+ ExcInternalError());
+
+ // The 3D Nedelec finite elements have 2*k*(k+1) dofs per each quad.
+ Assert(2 * k * (k + 1) == this->n_dofs_per_quad(face_no), ExcInternalError());
+
+ const int *swap_table = swap_table_array[k];
+
+ const unsigned int half_dofs = k * (k + 1); // see below;
+
+ const int rl = row_length[k];
+
+ int offset = table_size[0];
+ // The assignment above is to prevent the compiler from complaining about the
+ // unused table_size.
+
+ int value = 0;
+
+ for (const bool face_orientation : {false, true})
+ for (const bool face_rotation : {false, true})
+ for (const bool face_flip : {false, true})
+ {
+ const auto case_no =
+ internal::combined_face_orientation(face_orientation,
+ face_rotation,
+ face_flip);
+
+ // The dofs on a quad are indexed as the following:
+ //
+ // | x0, x1, x2, x3, ..., xk | y0, y1, y2, y3 ..., yk |
+ // | | |
+ // |-- half_ dofs = k*(k+1) --|-- half_dofs = k*(k+1) --|
+ // | |
+ // |-------------------- 2*k*(k+1) ---------------------|
+
+ for (unsigned int indx_x = 0; indx_x < half_dofs; indx_x++)
+ {
+ offset = 3 * rl * case_no + 0 * rl + indx_x;
+ Assert(offset < table_size[k], ExcInternalError());
+ value = *(swap_table + offset);
+ Assert(value < table_size[k], ExcInternalError());
+ Assert(value > -2, ExcInternalError());
+
+ if (value != -1)
+ {
+ const unsigned int indx_y =
+ half_dofs + static_cast<unsigned int>(value);
+
+ // dofs swap
+ this
+ ->adjust_quad_dof_index_for_face_orientation_table[face_no](
+ indx_x, case_no) = indx_y - indx_x;
+
+ this
+ ->adjust_quad_dof_index_for_face_orientation_table[face_no](
+ indx_y, case_no) = indx_x - indx_y;
+ }
+
+ // dof sign change
+ offset = 3 * rl * case_no + 1 * rl + indx_x;
+ Assert(offset < table_size[k], ExcInternalError());
+ value = *(swap_table + offset);
+ Assert((value == 0) || (value == 1), ExcInternalError());
+
+ this->adjust_quad_dof_sign_for_face_orientation_table[face_no](
+ indx_x, case_no) = static_cast<bool>(value);
+
+
+ offset = 3 * rl * case_no + 2 * rl + indx_x;
+ Assert(offset < table_size[k], ExcInternalError());
+ value = *(swap_table + offset);
+ Assert((value == 0) || (value == 1), ExcInternalError());
+
+ this->adjust_quad_dof_sign_for_face_orientation_table[face_no](
+ indx_x + half_dofs, case_no) = static_cast<bool>(value);
+ }
+ }
+
return;
}
// nothing to do in 1d
}
-
-
template <int spacedim>
void
get_dof_sign_change_nedelec(
const typename dealii::Triangulation<2, spacedim>::cell_iterator &cell,
- const FiniteElement<2, spacedim> & /*fe*/,
+ const FiniteElement<2, spacedim> &fe,
const std::vector<MappingKind> &mapping_kind,
std::vector<double> &line_dof_sign)
{
- const unsigned int dim = 2;
- // TODO: This fixes only lowest order
- for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
- if ((cell->line_orientation(l) ==
- numbers::reverse_line_orientation) &&
+ // The Nedelec finite elements in two spacial dimensions have two types
+ // of dofs: the line dofs and the quad dofs. The line dofs are
+ // associated with the edges. They are shared between the neighbouring
+ // cells and, consequently, need to be adjusted to compensate for a
+ // possible mismatch in the edge orientation. The quad dofs, they are
+ // associated with the interiors of the cells, are not shared between
+ // the cells in two spacial dimensions and need no adjustments.
+ //
+ // The Nedelec finite elements in two spacial dimensions have
+ // 2*(k+1)*(k+2) dofs per cell. All dofs are distributed between the
+ // line and quad dofs as the following:
+ //
+ // 4*(k+1) line dofs; (k+1) dofs per line.
+ // 2*k*(k+1) quad dofs.
+ //
+ // The dofs are indexed in the following order: first all line dofs,
+ // then all quad dofs.
+ //
+ // Here we adjust the line dofs. The sign of a line dof needs to be
+ // changed if the edge on which the dof resides points in the opposite
+ // direction.
+ const unsigned int k = fe.tensor_degree() - 1;
+
+ for (unsigned int l = 0; l < GeometryInfo<2>::lines_per_cell; ++l)
+ if (!(cell->line_orientation(l)) &&
mapping_kind[0] == mapping_nedelec)
- line_dof_sign[l] = -1.0;
+ {
+ if (k == 0)
+ {
+ // The lowest order element (k=0) is straightforward, because
+ // there is a single dof per edge, which needs to be flipped:
+ line_dof_sign[l] = -1.0;
+ }
+ else
+ {
+ // The case k > 0 is a bit more complicated. As we adjust
+ // only edge dofs in this function, we need to concern
+ // ourselves with the first 4*(k+1) entries in line_dof_sign
+ // vector ignoring the rest. There are (k+1) dofs per edge.
+ // Let us consider the local dof indices on one edge,
+ // local_line_dof = 0...k. The shape functions with even
+ // indices are asymmetric. The corresponding dofs need sign
+ // adjustment if the edge points in the opposite direction.
+ // The shape functions with odd indices are symmetric. The
+ // corresponding dofs need no sign adjustment even if the edge
+ // points in the opposite direction. In the current context
+ // the notion of symmetry of a shape function means that a
+ // shape function looks exactly the same if it is looked upon
+ // from the centers of the two neighbouring cells that share
+ // it.
+ for (unsigned int local_line_dof = 0;
+ local_line_dof < (k + 1);
+ local_line_dof++)
+ if (local_line_dof % 2 == 0)
+ line_dof_sign[local_line_dof + l * (k + 1)] = -1.0;
+ }
+ }
}
-
template <int spacedim>
void
get_dof_sign_change_nedelec(
const typename dealii::Triangulation<3, spacedim>::cell_iterator &cell,
- const FiniteElement<3, spacedim> & /*fe*/,
+ const FiniteElement<3, spacedim> &fe,
const std::vector<MappingKind> &mapping_kind,
std::vector<double> &line_dof_sign)
{
- const unsigned int dim = 3;
- // TODO: This is probably only going to work for those elements for
- // which all dofs are face dofs
- for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
- if ((cell->line_orientation(l) ==
- numbers::reverse_line_orientation) &&
+ // This function does half of the job - it adjusts the sign of the
+ // line (edge) dofs. In the three-dimensional space the quad (face) dofs
+ // need to be adjusted as well. The quad dofs are treated by
+ // FE_Nedelec<dim>::initialize_quad_dof_index_permutation_and_sign_change()
+ // in fe_nedelec.cc. The dofs associated with the interior of the cells,
+ // the hex dofs, need no adjustments as they are not shared between the
+ // neighboring cells.
+
+ const unsigned int k = fe.tensor_degree() - 1;
+ // The order of the Nedelec elements equals the tensor degree minus one,
+ // k = n - 1. In the three-dimensional space the Nedelec elements of the
+ // lowermost order, k = 0, have only 12 line (edge) dofs. The Nedelec
+ // elements of the higher orders, k > 0, have 3*(k+1)*(k+2)^2 dofs in
+ // total if dim=3. The dofs in a cell are distributed between lines
+ // (edges), quads (faces), and the hex (the interior of the cell) as the
+ // following:
+ //
+ // 12*(k+1) line dofs; (k+1) dofs per line.
+ // 2*6*k*(k+1) quad dofs; 2*k*(k+1) dofs per quad.
+ // 3*(k+2)^2*(k+1) hex dofs.
+ //
+ // The dofs are indexed in the following order: first all line dofs,
+ // then all quad dofs, and then all hex dofs.
+ //
+ // Here we adjust only the line (edge) dofs. The line dofs need only
+ // sign adjustment. That is, no permutation of the line dofs is needed.
+ for (unsigned int l = 0; l < GeometryInfo<3>::lines_per_cell; ++l)
+ if (!(cell->line_orientation(l)) &&
mapping_kind[0] == mapping_nedelec)
- line_dof_sign[l] = -1.0;
+ {
+ if (k == 0)
+ {
+ // The lowest order element (k=0) is straightforward, because
+ // there is a single dof per edge, which needs to be flipped:
+ line_dof_sign[l] = -1.0;
+ }
+ else
+ {
+ // The case k > 0 is a bit more complicated. As we adjust
+ // only edge dofs in this function, we need to concern
+ // ourselves with the first 12*(k+1) entries in line_dof_sign
+ // vector ignoring the rest. There are (k+1) dofs per edge.
+ // Let us consider the local dof indices on one edge,
+ // local_line_dof = 0...k. The shape functions with even
+ // indices are asymmetric. The corresponding dofs need sign
+ // adjustment if the edge points in the opposite direction.
+ // The shape functions with odd indices are symmetric. The
+ // corresponding dofs need no sign adjustment even if the edge
+ // points in the opposite direction. In the current context
+ // the notion of symmetry of a shape function means that a
+ // shape function looks exactly the same if it is looked upon
+ // from the centers of the two neighbouring cells that share
+ // it.
+ for (unsigned int local_line_dof = 0;
+ local_line_dof < (k + 1);
+ local_line_dof++)
+ if (local_line_dof % 2 == 0)
+ line_dof_sign[local_line_dof + l * (k + 1)] = -1.0;
+ }
+ }
}
} // namespace
} // namespace FE_PolyTensor
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2013 - 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+// In general, vertices of a face can be incorporated into two neighbouring
+// cells such that the face appears to have different orientations in the two
+// cells. The code of deal.II compensates for such mismatch in face orientation.
+// This code tests if deal.II does the compensation well for the Nedelec finite
+// elements in two- and three-dimensional spaces on globally refined meshes.
+//
+// In the two-dimensional space the Nedelec finite elements have dofs associated
+// with lines (edges) and the quad (the interior of the cell). Only line
+// dofs are shared between the cells in 2D. In the three-dimensional space the
+// Nedelec finite elements have dofs that are associated with lines (edges),
+// quads (faces), and the hex (the interior of the cell). The dofs associated
+// with lines and quads are shared between the cells in 3D. Therefore, the dofs
+// affected by the mismatch in face orientation are:
+// (i) the line dofs in 2D,
+// (ii) the line dofs in 3D,
+// (iii) the quad dofs in 3D.
+//
+// The shape functions that correspond to these three kinds of dofs share the
+// same feature: they exhibit reflection symmetry in the shared face if the
+// two cells that share the face are cubes (squares) of equal volume (area).
+// The code below exploits this feature.
+//
+// First, a mesh with a particular face mismatch is generated. Second, all shape
+// functions are sampled at sampling points. Indices of the sampling points in
+// one cell is the mirrored image of the indices in the other cell. Third, the
+// sampled shape functions are sorted in two categories: (i) shared between two
+// cells and (ii) locked inside a cell. The locked shaped functions are ignored.
+// The reflection symmetry of the shared shape functions is tested by computing
+// a score coefficient which is a cumulative l2_norm of a difference between a
+// sampled vector and its reflection. The score coefficient should be zero. The
+// program logs the score coefficients and two versions of the total amount of
+// shared dofs - one is theoretically expected, another is actually counted
+// during the execution of the code.
+//
+// In the three-dimensional space deal.II encodes the mismatch in face
+// orientation by three boolean parameters:
+//
+// (i) face_orientation - true if face is in standard orientation.
+// (ii) face_rotation - rotation by 90 deg counterclockwise if true;
+// (iii) face_flip - rotation by 180 deg counterclockwise if true.
+//
+// See the documentation of GeometryInfo<dim>.
+//
+// The combined face orientation is computes as
+//
+// orientation_no = face_flip*4 + face_rotation*2 + face_orientation*1;
+//
+// See tria_orientation.h.
+//
+// In the two-dimensional space the face orientation is encoded by
+// one boolean parameter line_orientation. The parameter equals true if the
+// faces have the same orientation.
+//
+// In this program the one parameter is used to encode the face orientation in
+// two- and three- dimensions. This parameter is
+//
+// unsigned char combined_face_orientation;
+//
+// Interpretation of this parameter in three dimensions is straightforward:
+// it is the combined face orientation as it is known in deal.II. In the
+// two-dimensional space this parameter is interpreted as the following.
+// The odd parameter (LSB = 1) implies that the shared faces (lines) are
+// aligned. The even parameter (LSB = 0) implies that the shared faces point in
+// opposite directions. That is, the least significant bit (LSB) of the
+// combined_face_orientation parameter is interpreted as the
+// line_orientation parameter which is used in deal.II on two-dimensional
+// meshes.
+//
+// Test meshes in two and three dimensions are generated by functions
+// GridGenerator::non_standard_orientation_mesh().
+// The three-dimensional mesh is used as it is. Only the leftmost shared face
+// is used in two dimensions. That is, shape functions are sampled only on
+// the leftmost and the middle cells of the two-dimensional mesh.
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_orientation.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/fe_field_function.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <cmath>
+#include <string>
+
+#include "../tests.h"
+
+template <int dim>
+void
+create_triangulation(Triangulation<dim> &triangulation,
+ const unsigned char combined_face_orientation);
+
+template <int dim>
+void
+create_sampling_points(std::vector<Point<dim>> &sampling_points_left,
+ std::vector<Point<dim>> &sampling_points_right,
+ unsigned int N);
+
+template <int dim>
+void
+save_sampling_points(std::vector<Point<dim>> &sampling_points_left,
+ std::vector<Point<dim>> &sampling_points_right);
+
+template <int dim>
+void
+print_shared_face(Triangulation<dim> &triangulation);
+
+template <>
+void
+create_triangulation<2>(Triangulation<2> &triangulation,
+ const unsigned char combined_face_orientation)
+
+{
+ // Only the leftmost shared face is considered. Even values of
+ // combined_face_orientation produce misaligned faces (the middle square is
+ // rotated by 2*[pi/2]). Odd values of combined_face_orientation produce
+ // aligned faces (the middle square is not rotated).
+ const unsigned int n_rotate_middle_square =
+ (combined_face_orientation % 2 == 0) ? 2 : 0;
+
+ GridGenerator::non_standard_orientation_mesh(triangulation,
+ n_rotate_middle_square);
+}
+
+template <>
+void
+create_triangulation<3>(Triangulation<3> &triangulation,
+ const unsigned char combined_face_orientation)
+
+{
+ bool face_orientation;
+ bool face_rotation;
+ bool face_flip;
+
+ std::tie(face_orientation, face_rotation, face_flip) =
+ dealii::internal::split_face_orientation(combined_face_orientation);
+
+ GridGenerator::non_standard_orientation_mesh(
+ triangulation, face_orientation, face_flip, face_rotation, false);
+}
+
+template <>
+void
+create_sampling_points<2>(std::vector<Point<2>> &sampling_points_left,
+ std::vector<Point<2>> &sampling_points_right,
+ unsigned int N)
+{
+ const double offset = 0.1;
+ const double step = (1.0 - 2 * offset) / (N - 1);
+
+ unsigned int indx = 0;
+
+ for (unsigned int j = 0; j < N; j++)
+ for (unsigned int i = 0; i < N; i++)
+ {
+ indx = i + N * j;
+
+ sampling_points_left.at(indx)[0] = offset + step * i;
+ sampling_points_left.at(indx)[1] = offset + step * j;
+
+ sampling_points_right.at(indx)[0] =
+ 2.0 - sampling_points_left.at(indx)[0];
+ sampling_points_right.at(indx)[1] = sampling_points_left.at(indx)[1];
+ }
+}
+
+template <>
+void
+create_sampling_points<3>(std::vector<Point<3>> &sampling_points_left,
+ std::vector<Point<3>> &sampling_points_right,
+ unsigned int N)
+{
+ const double offset = 0.1;
+ const double step = (1.0 - 2 * offset) / (N - 1);
+
+ unsigned int indx = 0;
+
+ for (unsigned int k = 0; k < N; k++)
+ for (unsigned int j = 0; j < N; j++)
+ for (unsigned int i = 0; i < N; i++)
+ {
+ indx = i + N * j + N * N * k;
+
+ sampling_points_left.at(indx)[0] = offset + step * i;
+ sampling_points_left.at(indx)[1] = offset + step * j;
+ sampling_points_left.at(indx)[2] = offset + step * k;
+
+ sampling_points_right.at(indx)[0] =
+ 2.0 - sampling_points_left.at(indx)[0];
+ sampling_points_right.at(indx)[1] = sampling_points_left.at(indx)[1];
+ sampling_points_right.at(indx)[2] = sampling_points_left.at(indx)[2];
+ }
+}
+
+template <int dim>
+void
+save_mesh(Triangulation<dim> &triangulation)
+{
+ GridOut gridout;
+ GridOutFlags::Msh msh_flags(true, true);
+ gridout.set_flags(msh_flags);
+
+ std::ofstream ofs("mesh_" + std::to_string(dim) + "D.msh");
+ gridout.write_msh(triangulation, ofs);
+}
+
+template <>
+void
+save_sampling_points<2>(std::vector<Point<2>> &sampling_points_left,
+ std::vector<Point<2>> &sampling_points_right)
+{
+ std::ofstream ofs_left("sampling_points_left_2D.csv");
+ std::ofstream ofs_right("sampling_points_right_2D.csv");
+
+ for (unsigned int i = 0; i < sampling_points_left.size(); i++)
+ {
+ ofs_left << sampling_points_left.at(i)[0] << " , "
+ << sampling_points_left.at(i)[1] << std::endl;
+
+ ofs_right << sampling_points_right.at(i)[0] << " , "
+ << sampling_points_right.at(i)[1] << std::endl;
+ }
+
+ ofs_left.close();
+ ofs_right.close();
+}
+
+template <>
+void
+save_sampling_points<3>(std::vector<Point<3>> &sampling_points_left,
+ std::vector<Point<3>> &sampling_points_right)
+{
+ std::ofstream ofs_left("sampling_points_left_3D.csv");
+ std::ofstream ofs_right("sampling_points_right_3D.csv");
+
+ for (unsigned int i = 0; i < sampling_points_left.size(); i++)
+ {
+ ofs_left << sampling_points_left.at(i)[0] << " , "
+ << sampling_points_left.at(i)[1] << " , "
+ << sampling_points_left.at(i)[2] << std::endl;
+
+ ofs_right << sampling_points_right.at(i)[0] << " , "
+ << sampling_points_right.at(i)[1] << " , "
+ << sampling_points_right.at(i)[2] << std::endl;
+ }
+
+ ofs_left.close();
+ ofs_right.close();
+}
+
+template <>
+void
+print_shared_face(Triangulation<2> &triangulation)
+{
+ const double eps = 1e-12;
+
+ for (const auto &cell : triangulation.active_cell_iterators())
+ for (unsigned int f = 0; f < 4; f++)
+ if (cell->neighbor_index(f) != -1)
+ if (((cell->face(f)->center()(0) - 1.0) < eps) &&
+ ((cell->face(f)->center()(1) - 0.5) < eps))
+ {
+ std::cout << "Cell center = (" << cell->center() << ")"
+ << std::endl;
+ std::cout << "The local index of the shared face: " << f
+ << std::endl;
+ std::cout << "Shared face - orientation: "
+ << cell->line_orientation(f) << std::endl;
+ }
+}
+
+template <>
+void
+print_shared_face<3>(Triangulation<3> &triangulation)
+{
+ for (const auto &cell : triangulation.active_cell_iterators())
+ {
+ unsigned int face = -1;
+
+ for (unsigned int f = 0; f < 6; f++)
+ if (cell->neighbor_index(f) != -1)
+ face = f;
+
+ std::cout << "Cell id = " << cell->id() << ": " << std::endl;
+
+ std::cout << "The local index of the shared face: " << face << std::endl;
+
+ std::cout << "Shared face - flip: " << cell->face_flip(face) << std::endl;
+ std::cout << "Shared face - rotation: " << cell->face_rotation(face)
+ << std::endl;
+ std::cout << "Shared face - orientation: " << cell->face_orientation(face)
+ << std::endl;
+ std::cout
+ << "Shared face - combined orientation: "
+ << static_cast<unsigned int>(
+ internal::combined_face_orientation(cell->face_orientation(face),
+ cell->face_rotation(face),
+ cell->face_flip(face)))
+ << std::endl;
+ }
+}
+
+struct SFData
+{
+ bool shared;
+ bool error;
+ double cum_norm_left;
+ double cum_norm_right;
+ double cum_norm_diff;
+};
+
+void
+classify(const std::vector<Vector<double>> &shape_function_left,
+ const std::vector<Vector<double>> &shape_function_right,
+ SFData &result)
+{
+ result.shared = false;
+ result.error = false;
+ result.cum_norm_left = 0.0;
+ result.cum_norm_right = 0.0;
+ result.cum_norm_diff = 0.0;
+
+ double th = 1e-3;
+ double eps = 1e-12;
+
+ unsigned int v_size = shape_function_left.size();
+ for (unsigned int i = 0; i < v_size; i++)
+ {
+ result.cum_norm_left += shape_function_left.at(i).l2_norm();
+ result.cum_norm_right += shape_function_right.at(i).l2_norm();
+
+ Vector<double> diff(shape_function_left.at(i));
+
+ diff.add(-1.0, shape_function_right.at(i));
+
+ result.cum_norm_diff += diff.l2_norm();
+ }
+
+ double s_v_size = std::sqrt(v_size);
+
+ result.cum_norm_left /= s_v_size;
+ result.cum_norm_right /= s_v_size;
+ result.cum_norm_diff /= s_v_size;
+
+ if ((result.cum_norm_left > th) && (result.cum_norm_right > th))
+ result.shared = true;
+
+ if (result.shared)
+ if (result.cum_norm_diff > eps)
+ result.error = true;
+}
+
+template <int dim, unsigned int order>
+void
+run(unsigned char combined_face_orientation)
+{
+ Triangulation<dim> triangulation;
+
+ create_triangulation<dim>(triangulation, combined_face_orientation);
+
+ // save_mesh<dim>(triangulation);
+
+ // print_shared_face<dim>(triangulation);
+
+ FE_Nedelec<dim> fe(order);
+
+ DoFHandler<dim> dof_handler(triangulation);
+
+ dof_handler.distribute_dofs(fe);
+
+ unsigned int N = order + 2;
+ unsigned int M = static_cast<unsigned int>(std::pow(N, dim));
+
+ std::vector<Point<dim>> sampling_points_left(M, Point<dim>());
+ std::vector<Point<dim>> sampling_points_right(M, Point<dim>());
+
+ std::vector<Vector<double>> shape_function_left(M, Vector<double>(dim));
+ std::vector<Vector<double>> shape_function_right(M, Vector<double>(dim));
+
+ create_sampling_points<dim>(sampling_points_left, sampling_points_right, N);
+
+ // save_sampling_points<dim>(sampling_points_left, sampling_points_right);
+
+ Vector<double> solution(dof_handler.n_dofs());
+
+ dealii::Functions::FEFieldFunction field_function(dof_handler, solution);
+
+ unsigned int n_shared_dofs_theory;
+
+ if (dim == 2)
+ n_shared_dofs_theory = order + 1;
+
+ if (dim == 3)
+ n_shared_dofs_theory = (2 * order + 1) * (order + 1);
+
+ unsigned int Q = solution.size();
+ unsigned int n_shared_dofs = 0;
+ for (unsigned int i = 0; i < Q; i++)
+ {
+ solution.reinit(Q);
+ solution(i) = 1.0;
+
+ field_function.vector_value_list(sampling_points_left,
+ shape_function_left);
+ field_function.vector_value_list(sampling_points_right,
+ shape_function_right);
+
+ SFData result;
+ classify(shape_function_left, shape_function_right, result);
+
+ if (result.error)
+ deallog << "error: " << result.cum_norm_left << " "
+ << result.cum_norm_right << " " << result.cum_norm_diff
+ << std::endl;
+
+ if ((result.shared) && (!result.error))
+ {
+ n_shared_dofs++;
+ deallog << result.cum_norm_diff << std::endl;
+ }
+
+
+ if (n_shared_dofs == n_shared_dofs_theory)
+ break;
+ }
+
+ deallog << n_shared_dofs << " : " << n_shared_dofs_theory << std::endl;
+}