std::pair<unsigned int, unsigned int>
face_system_to_component_index (const unsigned int index) const;
- /**
- * Given an index in the natural ordering of indices on a face, return the
- * index of the same degree of freedom on the cell.
- *
- * To explain the concept, consider the case where we would like to know
- * whether a degree of freedom on a face, for example as part of an FESystem
- * element, is primitive. Unfortunately, the
- * is_primitive() function in the FiniteElement class takes a cell index, so
- * we would need to find the cell index of the shape function that
- * corresponds to the present face index. This function does that.
- *
- * Code implementing this would then look like this:
- * @code
- * for (i=0; i<dofs_per_face; ++i)
- * if (fe.is_primitive(fe.face_to_equivalent_cell_index(i, some_face_no)))
- * ... do whatever
- * @endcode
- * The function takes additional arguments that account for the fact that
- * actual faces can be in their standard ordering with respect to the cell
- * under consideration, or can be flipped, oriented, etc.
- *
- * @param face_dof_index The index of the degree of freedom on a face.
- * This index must be between zero and dofs_per_face.
- * @param face The number of the face this degree of freedom lives on.
- * This number must be between zero and GeometryInfo::faces_per_cell.
- * @param face_orientation One part of the description of the orientation
- * of the face. See @ref GlossFaceOrientation .
- * @param face_flip One part of the description of the orientation
- * of the face. See @ref GlossFaceOrientation .
- * @param face_rotation One part of the description of the orientation
- * of the face. See @ref GlossFaceOrientation .
- * @return The index of this degree of freedom within the set
- * of degrees of freedom on the entire cell. The returned value
- * will be between zero and dofs_per_cell.
- *
- * @note This function exists in this class because that is where it
- * was first implemented. However, it can't really work in the most
- * general case without knowing what element we have. The reason is that
- * when a face is flipped or rotated, we also need to know whether we
- * need to swap the degrees of freedom on this face, or whether they
- * are immune from this. For this, consider the situation of a $Q_3$
- * element in 2d. If face_flip is true, then we need to consider
- * the two degrees of freedom on the edge in reverse order. On the other
- * hand, if the element were a $Q_1^2$, then because the two degrees of
- * freedom on this edge belong to different vector components, they
- * should not be considered in reverse order. What all of this shows is
- * that the function can't work if there are more than one degree of
- * freedom per line or quad, and that in these cases the function will
- * throw an exception pointing out that this functionality will need
- * to be provided by a derived class that knows what degrees of freedom
- * actually represent.
- */
- virtual
- unsigned int face_to_cell_index (const unsigned int face_dof_index,
- const unsigned int face,
- const bool face_orientation = true,
- const bool face_flip = false,
- const bool face_rotation = false) const;
-
/**
* For faces with non-standard face_orientation in 3D, the dofs on faces
* (quads) have to be permuted in order to be combined with the correct
const unsigned int offset,
typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
InternalDataBase &fe_internal,
- FEValuesData<dim,spacedim> &data) const;
+ FEValuesData<dim,spacedim> &data,
+ const MappingType mapping_type = mapping_covariant) const;
/**
* Given the pattern of nonzero components for each shape function, compute
-template <int dim, int spacedim>
-unsigned int
-FiniteElement<dim,spacedim>::
-face_to_cell_index (const unsigned int face_index,
- const unsigned int face,
- const bool face_orientation,
- const bool face_flip,
- const bool face_rotation) const
-{
- Assert (face_index < this->dofs_per_face,
- ExcIndexRange(face_index, 0, this->dofs_per_face));
- Assert (face < GeometryInfo<dim>::faces_per_cell,
- ExcIndexRange(face, 0, GeometryInfo<dim>::faces_per_cell));
-
-//TODO: we could presumably solve the 3d case below using the
-// adjust_quad_dof_index_for_face_orientation_table field. for the
-// 2d case, we can't use adjust_line_dof_index_for_line_orientation_table
-// since that array is empty (presumably because we thought that
-// there are no flipped edges in 2d, but these can happen in
-// DoFTools::make_periodicity_constraints, for example). so we
-// would need to either fill this field, or rely on derived classes
-// implementing this function, as we currently do
-
- // see the function's documentation for an explanation of this
- // assertion -- in essence, derived classes have to implement
- // an overloaded version of this function if we are to use any
- // other than standard orientation
- if ((face_orientation != true) || (face_flip != false) || (face_rotation != false))
- Assert ((this->dofs_per_line <= 1) && (this->dofs_per_quad <= 1),
- ExcMessage ("The function in this base class can not handle this case. "
- "Rather, the derived class you are using must provide "
- "an overloaded version but apparently hasn't done so. See "
- "the documentation of this function for more information."));
-
- // we need to distinguish between DoFs on vertices, lines and in 3d quads.
- // do so in a sequence of if-else statements
- if (face_index < this->first_face_line_index)
- // DoF is on a vertex
- {
- // get the number of the vertex on the face that corresponds to this DoF,
- // along with the number of the DoF on this vertex
- const unsigned int face_vertex = face_index / this->dofs_per_vertex;
- const unsigned int dof_index_on_vertex = face_index % this->dofs_per_vertex;
-
- // then get the number of this vertex on the cell and translate
- // this to a DoF number on the cell
- return (GeometryInfo<dim>::face_to_cell_vertices(face, face_vertex,
- face_orientation,
- face_flip,
- face_rotation)
- * this->dofs_per_vertex
- +
- dof_index_on_vertex);
- }
- else if (face_index < this->first_face_quad_index)
- // DoF is on a face
- {
- // do the same kind of translation as before. we need to only consider
- // DoFs on the lines, i.e., ignoring those on the vertices
- const unsigned int index = face_index - this->first_face_line_index;
-
- const unsigned int face_line = index / this->dofs_per_line;
- const unsigned int dof_index_on_line = index % this->dofs_per_line;
-
- return (this->first_line_index
- + GeometryInfo<dim>::face_to_cell_lines(face, face_line,
- face_orientation,
- face_flip,
- face_rotation)
- * this->dofs_per_line
- +
- dof_index_on_line);
- }
- else
- // DoF is on a quad
- {
- Assert (dim >= 3, ExcInternalError());
-
- // ignore vertex and line dofs
- const unsigned int index = face_index - this->first_face_quad_index;
-
- return (this->first_quad_index
- + face * this->dofs_per_quad
- + index);
- }
-}
-
-
-
-
template <int dim, int spacedim>
unsigned int
FiniteElement<dim,spacedim>::adjust_quad_dof_index_for_face_orientation (const unsigned int index,
const unsigned int,
Mapping<1,2>::InternalDataBase &,
InternalDataBase &,
- FEValuesData<1,2> &) const
+ FEValuesData<1,2> &,
+ MappingType) const
{
Assert(false, ExcNotImplemented());
const unsigned int,
Mapping<1,3>::InternalDataBase &,
InternalDataBase &,
- FEValuesData<1,3> &) const
+ FEValuesData<1,3> &,
+ MappingType) const
{
Assert(false, ExcNotImplemented());
const unsigned int,
Mapping<2,3>::InternalDataBase &,
InternalDataBase &,
- FEValuesData<2,3> &) const
+ FEValuesData<2,3> &,
+ MappingType) const
{
Assert(false, ExcNotImplemented());
const unsigned int offset,
typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
InternalDataBase &fe_internal,
- FEValuesData<dim,spacedim> &data) const
+ FEValuesData<dim,spacedim> &data,
+ MappingType mapping_type) const
{
Assert ((fe_internal.update_each | fe_internal.update_once)
& update_hessians,
// cell
for (unsigned int d=0; d<spacedim; ++d)
{
- mapping.transform (diff_quot[d], diff_quot2, mapping_internal, mapping_covariant);
+ mapping.transform (diff_quot[d], diff_quot2, mapping_internal, mapping_type);
for (unsigned int q=0; q<n_q_points; ++q)
for (unsigned int d1=0; d1<spacedim; ++d1)
}
if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
- this->compute_2nd (mapping, cell,
- typename QProjector<dim>::DataSetDescriptor().cell(),
- mapping_data, fe_data, data);
+ if (mapping_type == mapping_nedelec)
+ {
+ this->compute_2nd (mapping, cell,
+ typename QProjector<dim>::DataSetDescriptor().cell(),
+ mapping_data, fe_data, data, mapping_covariant_gradient);
+
+ for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
+ {
+ const unsigned int first = data.shape_function_to_row_table[i * this->n_components() +
+ this->get_nonzero_components(i).first_selected_component()];
+
+ for (unsigned int k = 0; k < n_q_points; ++k)
+ for (unsigned int d = 0; d < dim; ++d)
+ data.shape_hessians[first + d][k] *= sign_change[i];
+ }
+ }
+
+ else
+ this->compute_2nd (mapping, cell,
+ typename QProjector<dim>::DataSetDescriptor().cell(),
+ mapping_data, fe_data, data);
}
}
if (flags & update_hessians)
- this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
+ if (mapping_type == mapping_nedelec)
+ {
+ this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data, mapping_covariant_gradient);
+
+ for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
+ {
+ const unsigned int first = data.shape_function_to_row_table[i * this->n_components() +
+ this->get_nonzero_components(i).first_selected_component()];
+
+ for (unsigned int k = 0; k < n_q_points; ++k)
+ for (unsigned int d = 0; d < dim; ++d)
+ data.shape_hessians[first + d][k] *= sign_change[i];
+ }
+ }
+
+ else
+ this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
}
for (unsigned int k = 0; k < n_q_points; ++k)
for (unsigned int d = 0; d < dim; ++d)
- data.shape_gradients[first + d][k] = shape_grads1[k][d];
+ data.shape_gradients[first + d][k] = sign_change[i] * shape_grads1[k][d];
break;
}
}
if (flags & update_hessians)
- this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
+ if (mapping_type == mapping_nedelec)
+ {
+ this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data, mapping_covariant_gradient);
+
+ for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
+ {
+ const unsigned int first = data.shape_function_to_row_table[i * this->n_components() +
+ this->get_nonzero_components(i).first_selected_component()];
+
+ for (unsigned int k = 0; k < n_q_points; ++k)
+ for (unsigned int d = 0; d < dim; ++d)
+ data.shape_hessians[first + d][k] *= sign_change[i];
+ }
+ }
+
+ else
+ this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
}