From 87410017dfddb5b234d71b85f8d41922ede0e6fe Mon Sep 17 00:00:00 2001 From: Jiaqi Zhang Date: Thu, 1 Sep 2022 11:07:29 -0400 Subject: [PATCH] face dof to vector dof --- .../vector_tools_constraints.templates.h | 180 ++++++++++++++++-- 1 file changed, 162 insertions(+), 18 deletions(-) diff --git a/include/deal.II/numerics/vector_tools_constraints.templates.h b/include/deal.II/numerics/vector_tools_constraints.templates.h index d15bd41817..e94598f9f5 100644 --- a/include/deal.II/numerics/vector_tools_constraints.templates.h +++ b/include/deal.II/numerics/vector_tools_constraints.templates.h @@ -85,6 +85,52 @@ namespace VectorTools + /** + * Copied from VectorDoFTuple but only stores + * active fe index, face number, DoF index. + */ + struct FaceDoFInfo + { + static const unsigned int size = 3; + // store active fe indeX, face number, DoF index + unsigned int indices[size]; + + FaceDoFInfo() + { + for (unsigned int i = 0; i < size; ++i) + indices[i] = numbers::invalid_unsigned_int; + } + + bool + operator<(const FaceDoFInfo &other) const + { + for (unsigned int i = 0; i < size; ++i) + if (indices[i] < other.indices[i]) + return true; + else if (indices[i] > other.indices[i]) + return false; + return false; + } + + bool + operator==(const FaceDoFInfo &other) const + { + for (unsigned int i = 0; i < size; ++i) + if (indices[i] != other.indices[i]) + return false; + + return true; + } + + bool + operator!=(const FaceDoFInfo &other) const + { + return !(*this == other); + } + }; + + + /** * Add the constraint $\vec n \cdot \vec u = inhom$ to the list of * constraints. @@ -463,6 +509,89 @@ namespace VectorTools } + + /** + * Map active fe index, face number and + * local degrees of freedom index to the dim degrees of freedom indices that + * correspond to a vector-valued quantity at the same support point. + */ + template + void + map_face_dof_to_vector_dof( + const typename DoFHandler::cell_iterator &cell, + const unsigned int first_vector_component, + const std::set &boundary_ids, + hp::FEFaceValues & x_fe_face_values, + const IndexSet & refinement_edge_indices, + const unsigned int level, + std::map> + &dof_to_vector_dof) + { + std::set::iterator b_id; + for (const unsigned int face_no : cell->face_indices()) + if ((b_id = boundary_ids.find(cell->face(face_no)->boundary_id())) != + boundary_ids.end()) + { + const FiniteElement &fe = cell->get_fe(); + typename DoFHandler::level_face_iterator face = + cell->face(face_no); + + std::vector face_dofs; + // get the indices of the dofs on this cell... + face_dofs.resize(fe.n_dofs_per_face(face_no)); + + if (level != numbers::invalid_unsigned_int) + face->get_mg_dof_indices(level, + face_dofs, + cell->active_fe_index()); + else + face->get_dof_indices(face_dofs, cell->active_fe_index()); + + x_fe_face_values.reinit(cell, face_no); + + const auto &face_quadrature_collection = + x_fe_face_values.get_quadrature_collection(); + + std::array local_vector_indices; + for (unsigned int i = 0; i < face_dofs.size(); ++i) + { + if (fe.face_system_to_component_index(i, face_no).first == + first_vector_component) + if (level == numbers::invalid_unsigned_int || + !refinement_edge_indices.is_element(face_dofs[i])) + { + local_vector_indices[0] = i; + internal::FaceDoFInfo local_face_dof_info; + local_face_dof_info.indices[0] = cell->active_fe_index(); + local_face_dof_info.indices[1] = face_no; + local_face_dof_info.indices[2] = i; + for (unsigned int k = 0; k < fe.n_dofs_per_face(face_no); + ++k) + if ((k != i) && + (face_quadrature_collection[cell->active_fe_index()] + .point(k) == + face_quadrature_collection[cell->active_fe_index()] + .point(i)) && + (fe.face_system_to_component_index(k, face_no) + .first >= first_vector_component) && + (fe.face_system_to_component_index(k, face_no) + .first < first_vector_component + dim)) + { + local_vector_indices + [fe.face_system_to_component_index(k, face_no) + .first - + first_vector_component] = k; + } + dof_to_vector_dof.insert( + std::make_pair(local_face_dof_info, + local_vector_indices)); + } + } + } + } + + + /** * Compute the mappings from vector degrees of freedom to normal vectors @p dof_to_normals_map * and vector degrees of freedom to prescribed normal fluxes @p dof_vector_to_b_values. @@ -487,6 +616,21 @@ namespace VectorTools std::map, Vector> &dof_vector_to_b_values) { + // mapping from (active_fe_index, face_no and local + // dof index) to dim vector indices of the same + // supporting point. + std::map> + dof_to_vector_dof; + internal::map_face_dof_to_vector_dof( + cell, + first_vector_component, + boundary_ids, + x_fe_face_values, + refinement_edge_indices, + level, + dof_to_vector_dof); + + std::set::iterator b_id; for (const unsigned int face_no : cell->face_indices()) if ((b_id = boundary_ids.find(cell->face(face_no)->boundary_id())) != @@ -511,9 +655,6 @@ namespace VectorTools const FEFaceValues &fe_values = x_fe_face_values.get_present_fe_values(); - const auto &face_quadrature_collection = - x_fe_face_values.get_quadrature_collection(); - // then identify which of them correspond to the selected set of // vector components for (unsigned int i = 0; i < face_dofs.size(); ++i) @@ -529,26 +670,29 @@ namespace VectorTools internal::VectorDoFTuple vector_dofs; vector_dofs.dof_indices[0] = face_dofs[i]; + internal::FaceDoFInfo face_dof_info; + face_dof_info.indices[0] = cell->active_fe_index(); + face_dof_info.indices[1] = face_no; + face_dof_info.indices[2] = i; + + std::array local_vector_indices; + for (unsigned int d = 0; d < dim; ++d) + local_vector_indices[d] = numbers::invalid_unsigned_int; + + const auto it = dof_to_vector_dof.find(face_dof_info); + if (it != dof_to_vector_dof.end()) + local_vector_indices = it->second; + Assert( first_vector_component + dim <= fe.n_components(), ExcMessage( "Error: the finite element does not have enough components " "to define a normal direction.")); - - for (unsigned int k = 0; k < fe.n_dofs_per_face(face_no); - ++k) - if ((k != i) && - (face_quadrature_collection[cell->active_fe_index()] - .point(k) == - face_quadrature_collection[cell->active_fe_index()] - .point(i)) && - (fe.face_system_to_component_index(k, face_no) - .first >= first_vector_component) && - (fe.face_system_to_component_index(k, face_no).first < - first_vector_component + dim)) - vector_dofs.dof_indices - [fe.face_system_to_component_index(k, face_no).first - - first_vector_component] = face_dofs[k]; + for (unsigned int d = 1; d < dim; ++d) + { + vector_dofs.dof_indices[d] = + face_dofs[local_vector_indices[d]]; + } for (unsigned int d = 0; d < dim; ++d) Assert(vector_dofs.dof_indices[d] < n_dofs, -- 2.39.5