// at quadrature points.
boundary_function.vector_value_list(quadrature_points, values);
- // Find the group of vector components (dim of them,
- // starting at first_vector_component) are within an FESystem.
- //
- // If not using FESystem then must be using FE_Nedelec,
- // which has one base element and one copy of it (with 3 components).
+ // Find the group of vector components we want to project onto
+ // (dim of them, starting at first_vector_component) within the
+ // overall finite element (which may be an FESystem).
std::pair<unsigned int, unsigned int> base_indices(0, 0);
if (dynamic_cast<const FESystem<dim> *>(&cell->get_fe()) != nullptr)
{
base_indices.second = (first_vector_component - fe_index_old) /
fe.base_element(i).n_components();
}
- // Store degree as fe.degree-1
- // For nedelec elements FE_Nedelec<dim> (0) returns fe.degree = 1.
- // For FESystem get the degree from the base_element
- // indicated by the first_vector_component
+ else
+ // The only other element we know how to deal with (so far) is
+ // FE_Nedelec, which has one base element and one copy of it
+ // (with 3 components). In that case, the values of
+ // 'base_indices' as initialized above are correct.
+ Assert((dynamic_cast<const FE_Nedelec<dim> *>(&cell->get_fe()) !=
+ nullptr) ||
+ (dynamic_cast<const FE_NedelecSZ<dim> *>(&cell->get_fe()) !=
+ nullptr),
+ ExcNotImplemented());
+
+
+ // Store the 'degree' of the Nedelec element as fe.degree-1. For
+ // Nedelec elements, FE_Nedelec<dim>(0) returns fe.degree = 1
+ // because fe.degree stores the *polynomial* degree, not the
+ // degree of the element (which is typically defined based on
+ // the largest polynomial space that is *complete* within the
+ // finite element).
const unsigned int degree =
fe.base_element(base_indices.first).degree - 1;
- // Find DoFs we want to constrain:
- // There are fe.dofs_per_line DoFs associated with the
- // given line on the given face on the given cell.
+ // Find DoFs we want to constrain: There are
+ // fe.base_element(base_indices.first).dofs_per_line DoFs
+ // associated with the given line on the given face on the given
+ // cell.
//
- // Want to know which of these DoFs (there are degree+1 of interest)
+ // We need to know which of these DoFs (there are degree+1 of interest)
// are associated with the components given by first_vector_component.
// Then we can make a map from the associated line DoFs to the face DoFs.
//
// element and the index within this ordering.
//
// We call the map associated_edge_dof_to_face_dof
- std::vector<unsigned int> associated_edge_dof_to_face_dof(degree + 1);
+ std::vector<unsigned int> associated_edge_dof_to_face_dof(
+ degree + 1, numbers::invalid_dof_index);
// Lowest DoF in the base element allowed for this edge:
const unsigned int lower_bound =
.face_to_cell_index((line + 1) * (degree + 1) - 1, face);
unsigned int associated_edge_dof_index = 0;
- // for (unsigned int face_idx = 0; face_idx < fe.dofs_per_face;
- // ++face_idx)
- for (unsigned int line_idx = 0; line_idx < fe.dofs_per_line; ++line_idx)
+ for (unsigned int line_dof_idx = 0; line_dof_idx < fe.dofs_per_line;
+ ++line_dof_idx)
{
- // Assuming DoFs on a face are numbered in order by lines then faces.
+ // For each DoF associated with the (interior of) the line, we need
+ // to figure out which base element it belongs to and then if
+ // that's the correct base element. This is complicated by the
+ // fact that the FiniteElement class has functions that translate
+ // from face to cell, but not from edge to cell index systems. So
+ // we have to do that step by step.
+ //
+ // DoFs on a face in 3d are numbered in order by vertices then lines
+ // then faces.
// i.e. line 0 has degree+1 dofs numbered 0,..,degree
// line 1 has degree+1 dofs numbered (degree+1),..,2*(degree+1)
// and so on.
- const unsigned int face_idx = line * fe.dofs_per_line + line_idx;
+
+ const unsigned int face_dof_idx =
+ GeometryInfo<dim>::vertices_per_face * fe.dofs_per_vertex +
+ line * fe.dofs_per_line + line_dof_idx;
+
// Note, assuming that the edge orientations are "standard"
// i.e. cell->line_orientation(line) = true.
- const unsigned int cell_idx = fe.face_to_cell_index(face_idx, face);
-
- // Check this cell_idx belongs to the correct base_element, component
- // and line:
- if (((dynamic_cast<const FESystem<dim> *>(&fe) != nullptr) &&
- (fe.system_to_base_index(cell_idx).first == base_indices) &&
- (lower_bound <= fe.system_to_base_index(cell_idx).second) &&
- (fe.system_to_base_index(cell_idx).second <= upper_bound)) ||
- (((dynamic_cast<const FE_Nedelec<dim> *>(&fe) != nullptr) ||
- (dynamic_cast<const FE_NedelecSZ<dim> *>(&fe) != nullptr)) &&
- (line * (degree + 1) <= face_idx) &&
- (face_idx <= (line + 1) * (degree + 1) - 1)))
+ Assert(cell->line_orientation(line),
+ ExcMessage("Edge orientation doesnot meet expectation."));
+ // Next, translate from face to cell. Note, this might be assuming
+ // that the edge orientations are "standard" (not sure any more at
+ // this time), i.e.
+ // cell->line_orientation(line) = true.
+ const unsigned int cell_dof_idx =
+ fe.face_to_cell_index(face_dof_idx, face);
+
+ // Check that this cell_idx belongs to the correct base_element,
+ // component and line. We do this for each of the supported elements
+ // separately
+ bool dof_is_of_interest = false;
+ if (dynamic_cast<const FESystem<dim> *>(&fe) != nullptr)
+ {
+ dof_is_of_interest =
+ (fe.system_to_base_index(cell_dof_idx).first == base_indices) &&
+ (lower_bound <= fe.system_to_base_index(cell_dof_idx).second) &&
+ (fe.system_to_base_index(cell_dof_idx).second <= upper_bound);
+ }
+ else if ((dynamic_cast<const FE_Nedelec<dim> *>(&fe) != nullptr) ||
+ (dynamic_cast<const FE_NedelecSZ<dim> *>(&fe) != nullptr))
+ {
+ Assert((line * (degree + 1) <= face_dof_idx) &&
+ (face_dof_idx < (line + 1) * (degree + 1)),
+ ExcInternalError());
+ dof_is_of_interest = true;
+ }
+ else
+ Assert(false, ExcNotImplemented());
+
+ if (dof_is_of_interest)
{
associated_edge_dof_to_face_dof[associated_edge_dof_index] =
- face_idx;
+ face_dof_idx;
++associated_edge_dof_index;
}
}
// Sanity check:
- const unsigned int associated_edge_dofs = associated_edge_dof_index;
- Assert(associated_edge_dofs == degree + 1,
+ const unsigned int n_associated_edge_dofs = associated_edge_dof_index;
+ Assert(n_associated_edge_dofs == degree + 1,
ExcMessage("Error: Unexpected number of 3D edge DoFs"));
// Matrix and RHS vectors to store linear system:
// The RHS entries are:
// \int_{edge}
// (tangential*boundary_value)*(tangential*edge_shape_function_i) dS.
- for (unsigned int j = 0; j < associated_edge_dofs; ++j)
+ for (unsigned int j = 0; j < n_associated_edge_dofs; ++j)
{
const unsigned int j_face_idx =
associated_edge_dof_to_face_dof[j];
const unsigned int j_cell_idx =
fe.face_to_cell_index(j_face_idx, face);
- for (unsigned int i = 0; i < associated_edge_dofs; ++i)
+ for (unsigned int i = 0; i < n_associated_edge_dofs; ++i)
{
const unsigned int i_face_idx =
associated_edge_dof_to_face_dof[i];
edge_matrix_inv.vmult(edge_solution, edge_rhs);
// Store computed DoFs
- for (unsigned int i = 0; i < associated_edge_dofs; ++i)
+ for (unsigned int i = 0; i < n_associated_edge_dofs; ++i)
{
dof_values[associated_edge_dof_to_face_dof[i]] = edge_solution(i);
dofs_processed[associated_edge_dof_to_face_dof[i]] = true;
base_indices.second = (first_vector_component - fe_index_old) /
fe.base_element(i).n_components();
}
+ else
+ {
+ // Assert that the FE is in fact an FE_Nedelec, so that the default
+ // base_indices == (0,0) is correct.
+ Assert((dynamic_cast<const FE_Nedelec<dim> *>(&cell->get_fe()) !=
+ nullptr) ||
+ (dynamic_cast<const FE_NedelecSZ<dim> *>(&cell->get_fe()) !=
+ nullptr),
+ ExcNotImplemented());
+ }
const unsigned int degree =
fe.base_element(base_indices.first).degree - 1;
fe.base_element(base_indices.first)
.face_to_cell_index((line + 1) * (degree + 1) - 1, face);
unsigned int associated_edge_dof_index = 0;
- for (unsigned int line_idx = 0; line_idx < fe.dofs_per_line;
- ++line_idx)
+
+ for (unsigned int line_dof_idx = 0;
+ line_dof_idx < fe.dofs_per_line;
+ ++line_dof_idx)
{
- const unsigned int face_idx =
- line * fe.dofs_per_line + line_idx;
- const unsigned int cell_idx =
- fe.face_to_cell_index(face_idx, face);
- // Check this cell_idx belongs to the correct
- // base_element, component and line:
- if (((dynamic_cast<const FESystem<dim> *>(&fe) !=
- nullptr) &&
- (fe.system_to_base_index(cell_idx).first ==
- base_indices) &&
- (lower_bound <=
- fe.system_to_base_index(cell_idx).second) &&
- (fe.system_to_base_index(cell_idx).second <=
- upper_bound)) ||
- (((dynamic_cast<const FE_Nedelec<dim> *>(&fe) !=
- nullptr) ||
- (dynamic_cast<const FE_NedelecSZ<dim> *>(&fe) !=
- nullptr)) &&
- (line * (degree + 1) <= face_idx) &&
- (face_idx <= (line + 1) * (degree + 1) - 1)))
+ // For each DoF associated with the (interior of) the
+ // line, we need to figure out which base element it
+ // belongs to and then if that's the correct base element.
+ // This is complicated by the fact that the FiniteElement
+ // class has functions that translate from face to cell,
+ // but not from edge to cell index systems. So we have to
+ // do that step by step.
+ //
+ // DoFs on a face in 3d are numbered in order by vertices
+ // then lines then faces. i.e. line 0 has degree+1 dofs
+ // numbered 0,..,degree
+ // line 1 has degree+1 dofs numbered
+ // (degree+1),..,2*(degree+1) and so on.
+ const unsigned int face_dof_idx =
+ GeometryInfo<dim>::vertices_per_face *
+ fe.dofs_per_vertex +
+ line * fe.dofs_per_line + line_dof_idx;
+
+ // Next, translate from face to cell. Note, this might be
+ // assuming that the edge orientations are "standard" (not
+ // sure any more at this time), i.e.
+ // cell->line_orientation(line) = true.
+ const unsigned int cell_dof_idx =
+ fe.face_to_cell_index(face_dof_idx, face);
+
+ // Check that this cell_idx belongs to the correct
+ // base_element, component and line. We do this for each
+ // of the supported elements separately
+ bool dof_is_of_interest = false;
+ if (dynamic_cast<const FESystem<dim> *>(&fe) != nullptr)
+ {
+ dof_is_of_interest =
+ (fe.system_to_base_index(cell_dof_idx).first ==
+ base_indices) &&
+ (lower_bound <=
+ fe.system_to_base_index(cell_dof_idx).second) &&
+ (fe.system_to_base_index(cell_dof_idx).second <=
+ upper_bound);
+ }
+ else if ((dynamic_cast<const FE_Nedelec<dim> *>(&fe) !=
+ nullptr) ||
+ (dynamic_cast<const FE_NedelecSZ<dim> *>(&fe) !=
+ nullptr))
+ {
+ Assert((line * (degree + 1) <= face_dof_idx) &&
+ (face_dof_idx < (line + 1) * (degree + 1)),
+ ExcInternalError());
+ dof_is_of_interest = true;
+ }
+ else
+ Assert(false, ExcNotImplemented());
+
+ if (dof_is_of_interest)
{
associated_edge_dof_to_face_dof
- [line][associated_edge_dof_index] = face_idx;
+ [line][associated_edge_dof_index] = face_dof_idx;
++associated_edge_dof_index;
}
}
// Sanity check:
associated_edge_dofs[line] = associated_edge_dof_index;
Assert(associated_edge_dofs[line] == degree + 1,
- ExcMessage(
- "Error: Unexpected number of 3D edge DoFs"));
+ ExcInternalError());
}
// Next find the face DoFs associated with the vector components
// we're interested in. There are 2*degree*(degree+1) DoFs
- // associated with each face (not including edges!).
+ // associated with the interior of each face (not including
+ // edges!).
//
// Create a map mapping from the consecutively numbered
// associated_dofs to the face DoF (which can be transferred to a
std::vector<unsigned int> associated_face_dof_to_face_dof(
2 * degree * (degree + 1));
- // Skip the edge DoFs, so we start at
- // lines_per_face*(fe.dofs_per_line).
+ // Loop over these quad-interior dofs.
unsigned int associated_face_dof_index = 0;
- for (unsigned int face_idx = lines_per_face * (fe.dofs_per_line);
- face_idx < fe.dofs_per_face;
- ++face_idx)
+ for (unsigned int quad_dof_idx = 0;
+ quad_dof_idx < fe.dofs_per_quad;
+ ++quad_dof_idx)
{
+ const unsigned int face_idx =
+ GeometryInfo<dim>::vertices_per_face * fe.dofs_per_vertex +
+ lines_per_face * fe.dofs_per_line + quad_dof_idx;
const unsigned int cell_idx =
fe.face_to_cell_index(face_idx, face);
if (((dynamic_cast<const FESystem<dim> *>(&fe) != nullptr) &&
(dynamic_cast<const FE_Nedelec<dim> *>(&fe) != nullptr) ||
(dynamic_cast<const FE_NedelecSZ<dim> *>(&fe) != nullptr))
{
+ AssertIndexRange(associated_face_dof_index,
+ associated_face_dof_to_face_dof.size());
associated_face_dof_to_face_dof
- [associated_face_dof_index] = face_idx;
+ [associated_face_dof_index] = quad_dof_idx;
++associated_face_dof_index;
}
}
projects / dealii.git / commitdiff