From 104ecad246dbb84f87b2db616493d84660ea49b6 Mon Sep 17 00:00:00 2001 From: buerg Date: Fri, 10 Jan 2014 22:35:44 +0000 Subject: [PATCH] Compute hessians for FE_Nedelec correctly. git-svn-id: https://svn.dealii.org/trunk@32191 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/fe/fe.h | 62 +---------------- deal.II/source/fe/fe.cc | 104 +++------------------------- deal.II/source/fe/fe_poly_tensor.cc | 62 +++++++++++++++-- 3 files changed, 67 insertions(+), 161 deletions(-) diff --git a/deal.II/include/deal.II/fe/fe.h b/deal.II/include/deal.II/fe/fe.h index 40b4c742c5..55898960f7 100644 --- a/deal.II/include/deal.II/fe/fe.h +++ b/deal.II/include/deal.II/fe/fe.h @@ -931,65 +931,6 @@ public: std::pair 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::InternalDataBase &mapping_internal, InternalDataBase &fe_internal, - FEValuesData &data) const; + FEValuesData &data, + const MappingType mapping_type = mapping_covariant) const; /** * Given the pattern of nonzero components for each shape function, compute diff --git a/deal.II/source/fe/fe.cc b/deal.II/source/fe/fe.cc index df61dbc73a..a9f9e0aeca 100644 --- a/deal.II/source/fe/fe.cc +++ b/deal.II/source/fe/fe.cc @@ -544,96 +544,6 @@ block_mask (const ComponentMask &component_mask) const -template -unsigned int -FiniteElement:: -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::faces_per_cell, - ExcIndexRange(face, 0, GeometryInfo::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::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::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 unsigned int FiniteElement::adjust_quad_dof_index_for_face_orientation (const unsigned int index, @@ -1195,7 +1105,8 @@ FiniteElement<1,2>::compute_2nd ( const unsigned int, Mapping<1,2>::InternalDataBase &, InternalDataBase &, - FEValuesData<1,2> &) const + FEValuesData<1,2> &, + MappingType) const { Assert(false, ExcNotImplemented()); @@ -1210,7 +1121,8 @@ FiniteElement<1,3>::compute_2nd ( const unsigned int, Mapping<1,3>::InternalDataBase &, InternalDataBase &, - FEValuesData<1,3> &) const + FEValuesData<1,3> &, + MappingType) const { Assert(false, ExcNotImplemented()); @@ -1226,7 +1138,8 @@ FiniteElement<2,3>::compute_2nd ( const unsigned int, Mapping<2,3>::InternalDataBase &, InternalDataBase &, - FEValuesData<2,3> &) const + FEValuesData<2,3> &, + MappingType) const { Assert(false, ExcNotImplemented()); @@ -1242,7 +1155,8 @@ FiniteElement::compute_2nd ( const unsigned int offset, typename Mapping::InternalDataBase &mapping_internal, InternalDataBase &fe_internal, - FEValuesData &data) const + FEValuesData &data, + MappingType mapping_type) const { Assert ((fe_internal.update_each | fe_internal.update_once) & update_hessians, @@ -1368,7 +1282,7 @@ FiniteElement::compute_2nd ( // cell for (unsigned int d=0; d::fill_fe_values ( } if (flags & update_hessians && cell_similarity != CellSimilarity::translation) - this->compute_2nd (mapping, cell, - typename QProjector::DataSetDescriptor().cell(), - mapping_data, fe_data, data); + if (mapping_type == mapping_nedelec) + { + this->compute_2nd (mapping, cell, + typename QProjector::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::DataSetDescriptor().cell(), + mapping_data, fe_data, data); } @@ -791,7 +809,23 @@ FE_PolyTensor::fill_fe_face_values ( } 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); } @@ -989,7 +1023,7 @@ FE_PolyTensor::fill_fe_subface_values ( 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; } @@ -1001,7 +1035,23 @@ FE_PolyTensor::fill_fe_subface_values ( } 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); } -- 2.39.5