From ac0114858c2b67136df78e772f52bafd18135ca9 Mon Sep 17 00:00:00 2001 From: Jiaqi Zhang Date: Mon, 15 Aug 2022 01:03:47 -0400 Subject: [PATCH] fix a bug --- .../vector_tools_constraints.templates.h | 433 +++++++++++++----- tests/numerics/no_flux_17.cc | 46 +- tests/numerics/no_flux_17.output | 45 ++ 3 files changed, 414 insertions(+), 110 deletions(-) diff --git a/include/deal.II/numerics/vector_tools_constraints.templates.h b/include/deal.II/numerics/vector_tools_constraints.templates.h index 17c05a85ee..14c7c14987 100644 --- a/include/deal.II/numerics/vector_tools_constraints.templates.h +++ b/include/deal.II/numerics/vector_tools_constraints.templates.h @@ -461,127 +461,237 @@ namespace VectorTools Assert(false, ExcNotImplemented()); } } - } // namespace internal - - - template - void - compute_nonzero_normal_flux_constraints( - const DoFHandler & dof_handler, - const unsigned int first_vector_component, - const std::set &boundary_ids, - const std::map *> - & function_map, - AffineConstraints & constraints, - const Mapping &mapping, - const IndexSet & refinement_edge_indices, - const unsigned int level) - { - Assert(dim > 1, - ExcMessage("This function is not useful in 1d because it amounts " - "to imposing Dirichlet values on the vector-valued " - "quantity.")); - - const unsigned int mesh_level = - (level == numbers::invalid_unsigned_int) ? - dof_handler.get_triangulation().n_global_levels() - 1 : - level; - std::vector face_dofs; - - // create FE and mapping collections for all elements in use by this - // DoFHandler - const hp::FECollection &fe_collection = - dof_handler.get_fe_collection(); - hp::MappingCollection mapping_collection; - for (unsigned int i = 0; i < fe_collection.size(); ++i) - mapping_collection.push_back(mapping); - // TODO: the implementation makes the assumption that all faces have the - // same number of dofs - AssertDimension(dof_handler.get_fe().n_unique_faces(), 1); - const unsigned int face_no = 0; - // now also create a quadrature collection for the faces of a cell. fill - // it with a quadrature formula with the support points on faces for each - // FE - hp::QCollection face_quadrature_collection; - for (unsigned int i = 0; i < fe_collection.size(); ++i) - { - const std::vector> &unit_support_points = - fe_collection[i].get_unit_face_support_points(face_no); + template + void + map_dof_to_normals_on_level( + const DoFHandler & dof_handler, + const unsigned int first_vector_component, + const std::set &boundary_ids, + const std::map *> + & function_map, + hp::FEFaceValues &x_fe_face_values, + const IndexSet & refinement_edge_indices, + const unsigned int level, + std::multimap< + internal::VectorDoFTuple, + std::pair, + typename DoFHandler::cell_iterator>> + &dof_to_normals_map, + std::map, Vector> + &dof_vector_to_b_values) + { + Assert(level < dof_handler.get_triangulation().n_levels(), + ExcInternalError()); + + std::vector face_dofs; + + const auto &face_quadrature_collection = + x_fe_face_values.get_quadrature_collection(); + + // now loop over all cells and all faces + std::set::iterator b_id; + for (const auto &cell : dof_handler.cell_iterators_on_level(level)) + if (cell->level_subdomain_id() != numbers::artificial_subdomain_id && + cell->level_subdomain_id() != numbers::invalid_subdomain_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); - Assert(unit_support_points.size() == - fe_collection[i].n_dofs_per_face(face_no), - ExcInternalError()); + // get the indices of the dofs on this cell... + face_dofs.resize(fe.n_dofs_per_face(face_no)); - face_quadrature_collection.push_back( - Quadrature(unit_support_points)); - } + face->get_mg_dof_indices(level, + face_dofs, + cell->active_fe_index()); - // now create the object with which we will generate the normal vectors - hp::FEFaceValues x_fe_face_values(mapping_collection, - fe_collection, - face_quadrature_collection, - update_quadrature_points | - update_normal_vectors); + x_fe_face_values.reinit(cell, face_no); + const FEFaceValues &fe_values = + x_fe_face_values.get_present_fe_values(); + + // then identify which of them correspond to the selected set of + // vector components + for (unsigned int i = 0; i < face_dofs.size(); ++i) + if (fe.face_system_to_component_index(i, face_no).first == + first_vector_component) + // Refinement edge indices are going to be constrained to 0 + // during a multigrid cycle and do not need no-normal-flux + // constraints, so skip them: + if (!refinement_edge_indices.is_element(face_dofs[i])) + { + // find corresponding other components of vector + internal::VectorDoFTuple vector_dofs; + vector_dofs.dof_indices[0] = face_dofs[i]; + + 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 = 0; d < dim; ++d) + Assert(vector_dofs.dof_indices[d] < + dof_handler.n_dofs(), + ExcInternalError()); + + // we need the normal vector on this face. we know that + // it is a vector of length 1 but at least with higher + // order mappings it isn't always possible to guarantee + // that each component is exact up to zero tolerance. in + // particular, as shown in the deal.II/no_flux_06 test, + // if we just take the normal vector as given by the + // fe_values object, we can get entries in the normal + // vectors of the unit cube that have entries up to + // several times 1e-14. + // + // the problem with this is that this later yields + // constraints that are circular (e.g., in the testcase, + // we get constraints of the form + // + // x22 = 2.93099e-14*x21 + 2.93099e-14*x23 + // x21 = -2.93099e-14*x22 + 2.93099e-14*x21 + // + // in both of these constraints, the small numbers + // should be zero and the constraints should simply be + // x22 = x21 = 0 + // + // to achieve this, we utilize that we know that the + // normal vector has (or should have) length 1 and that + // we can simply set small elements to zero (without + // having to check that they are small *relative to + // something else*). we do this and then normalize the + // length of the vector back to one, just to be on the + // safe side + // + // one more point: we would like to use the "real" + // normal vector here, as provided by the boundary + // description and as opposed to what we get from the + // FEValues object. we do this in the immediately next + // line, but as is obvious, the boundary only has a + // vague idea which side of a cell it is on -- indicated + // by the face number. in other words, it may provide + // the inner or outer normal. by and large, there is no + // harm from this, since the tangential vector we + // compute is still the same. however, we do average + // over normal vectors from adjacent cells and if they + // have recorded normal vectors from the inside once and + // from the outside the other time, then this averaging + // is going to run into trouble. as a consequence we ask + // the mapping after all for its normal vector, but we + // only ask it so that we can possibly correct the sign + // of the normal vector provided by the boundary if they + // should point in different directions. this is the + // case in tests/deal.II/no_flux_11. + Tensor<1, dim> normal_vector = + (cell->face(face_no)->get_manifold().normal_vector( + cell->face(face_no), + fe_values.quadrature_point(i))); + if (normal_vector * fe_values.normal_vector(i) < 0) + normal_vector *= -1; + Assert(std::fabs(normal_vector.norm() - 1) < 1e-14, + ExcInternalError()); + for (unsigned int d = 0; d < dim; ++d) + if (std::fabs(normal_vector[d]) < 1e-13) + normal_vector[d] = 0; + normal_vector /= normal_vector.norm(); + + const Point &point = fe_values.quadrature_point(i); + Vector b_values(dim); + function_map.at(*b_id)->vector_value(point, b_values); + + // now enter the (dofs,(normal_vector,cell)) entry into + // the map + dof_to_normals_map.insert( + std::make_pair(vector_dofs, + std::make_pair(normal_vector, cell))); + dof_vector_to_b_values.insert( + std::make_pair(vector_dofs, b_values)); - // have a map that stores normal vectors for each vector-dof tuple we want - // to constrain. since we can get at the same vector dof tuple more than - // once (for example if it is located at a vertex that we visit from all - // adjacent cells), we will want to average later on the normal vectors - // computed on different cells as described in the documentation of this - // function. however, we can only average if the contributions came from - // different cells, whereas we want to constrain twice or more in case the - // contributions came from different faces of the same cell - // (i.e. constrain not just the *average normal direction* but *all normal - // directions* we find). consequently, we also have to store which cell a - // normal vector was computed on - using DoFToNormalsMap = std::multimap< - internal::VectorDoFTuple, - std::pair, - typename DoFHandler::cell_iterator>>; - std::map, Vector> - dof_vector_to_b_values; +#ifdef DEBUG_NO_NORMAL_FLUX + std::cout << "Adding normal vector:" << std::endl + << " dofs=" << vector_dofs << std::endl + << " cell=" << cell << " at " + << cell->center() << std::endl + << " normal=" << normal_vector << std::endl; +#endif + } + } + } - DoFToNormalsMap dof_to_normals_map; - // now loop over all cells and all faces - std::set::iterator b_id; - for (const auto &cell : dof_handler.cell_iterators_on_level(mesh_level)) - if (cell->level_subdomain_id() != numbers::artificial_subdomain_id && - cell->level_subdomain_id() != numbers::invalid_subdomain_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); - // get the indices of the dofs on this cell... - face_dofs.resize(fe.n_dofs_per_face(face_no)); + template + void + map_dof_to_normals( + const DoFHandler & dof_handler, + const unsigned int first_vector_component, + const std::set &boundary_ids, + const std::map *> + & function_map, + hp::FEFaceValues &x_fe_face_values, + std::multimap< + internal::VectorDoFTuple, + std::pair, + typename DoFHandler::cell_iterator>> + &dof_to_normals_map, + std::map, Vector> + &dof_vector_to_b_values) + { + std::vector face_dofs; + + const auto &face_quadrature_collection = + x_fe_face_values.get_quadrature_collection(); + + // now loop over all cells and all faces + std::set::iterator b_id; + for (const auto &cell : dof_handler.active_cell_iterators()) + if (!cell->is_artificial()) + 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::face_iterator face = + cell->face(face_no); - if (level != numbers::invalid_unsigned_int) - face->get_mg_dof_indices(mesh_level, - face_dofs, - cell->active_fe_index()); - else + // get the indices of the dofs on this cell... + face_dofs.resize(fe.n_dofs_per_face(face_no)); face->get_dof_indices(face_dofs, cell->active_fe_index()); - x_fe_face_values.reinit(cell, face_no); - const FEFaceValues &fe_values = - x_fe_face_values.get_present_fe_values(); + x_fe_face_values.reinit(cell, face_no); + const FEFaceValues &fe_values = + x_fe_face_values.get_present_fe_values(); - // then identify which of them correspond to the selected set of - // vector components - for (unsigned int i = 0; i < face_dofs.size(); ++i) - if (fe.face_system_to_component_index(i, face_no).first == - first_vector_component) - // Refinement edge indices are going to be constrained to 0 - // during a multigrid cycle and do not need no-normal-flux - // constraints, so skip them: - if (!refinement_edge_indices.is_element(face_dofs[i])) + // then identify which of them correspond to the selected set of + // vector components + for (unsigned int i = 0; i < face_dofs.size(); ++i) + if (fe.face_system_to_component_index(i, face_no).first == + first_vector_component) { // find corresponding other components of vector internal::VectorDoFTuple vector_dofs; @@ -694,7 +804,114 @@ namespace VectorTools << " normal=" << normal_vector << std::endl; #endif } - } + } + } + + } // namespace internal + + + + template + void + compute_nonzero_normal_flux_constraints( + const DoFHandler & dof_handler, + const unsigned int first_vector_component, + const std::set &boundary_ids, + const std::map *> + & function_map, + AffineConstraints & constraints, + const Mapping &mapping, + const IndexSet & refinement_edge_indices, + const unsigned int level) + { + Assert(dim > 1, + ExcMessage("This function is not useful in 1d because it amounts " + "to imposing Dirichlet values on the vector-valued " + "quantity.")); + + // create FE and mapping collections for all elements in use by this + // DoFHandler + const hp::FECollection &fe_collection = + dof_handler.get_fe_collection(); + hp::MappingCollection mapping_collection; + for (unsigned int i = 0; i < fe_collection.size(); ++i) + mapping_collection.push_back(mapping); + + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(dof_handler.get_fe().n_unique_faces(), 1); + const unsigned int face_no = 0; + + // now also create a quadrature collection for the faces of a cell. fill + // it with a quadrature formula with the support points on faces for each + // FE + hp::QCollection face_quadrature_collection; + for (unsigned int i = 0; i < fe_collection.size(); ++i) + { + const std::vector> &unit_support_points = + fe_collection[i].get_unit_face_support_points(face_no); + + Assert(unit_support_points.size() == + fe_collection[i].n_dofs_per_face(face_no), + ExcInternalError()); + + face_quadrature_collection.push_back( + Quadrature(unit_support_points)); + } + + // now create the object with which we will generate the normal vectors + hp::FEFaceValues x_fe_face_values(mapping_collection, + fe_collection, + face_quadrature_collection, + update_quadrature_points | + update_normal_vectors); + + // have a map that stores normal vectors for each vector-dof tuple we want + // to constrain. since we can get at the same vector dof tuple more than + // once (for example if it is located at a vertex that we visit from all + // adjacent cells), we will want to average later on the normal vectors + // computed on different cells as described in the documentation of this + // function. however, we can only average if the contributions came from + // different cells, whereas we want to constrain twice or more in case the + // contributions came from different faces of the same cell + // (i.e. constrain not just the *average normal direction* but *all normal + // directions* we find). consequently, we also have to store which cell a + // normal vector was computed on + using DoFToNormalsMap = std::multimap< + internal::VectorDoFTuple, + std::pair, + typename DoFHandler::cell_iterator>>; + std::map, Vector> + dof_vector_to_b_values; + + DoFToNormalsMap dof_to_normals_map; + + if (level == numbers::invalid_unsigned_int) + { + // active cells + internal::map_dof_to_normals(dof_handler, + first_vector_component, + boundary_ids, + function_map, + x_fe_face_values, + dof_to_normals_map, + dof_vector_to_b_values); + } + else + { // level cells + internal::map_dof_to_normals_on_level( + dof_handler, + first_vector_component, + boundary_ids, + function_map, + x_fe_face_values, + refinement_edge_indices, + level, + dof_to_normals_map, + dof_vector_to_b_values); + } + + // Now do something with the collected information. To this end, loop // through all sets of pairs (dofs,normal_vector) and identify which diff --git a/tests/numerics/no_flux_17.cc b/tests/numerics/no_flux_17.cc index 9acf3076fa..edd36d3a76 100644 --- a/tests/numerics/no_flux_17.cc +++ b/tests/numerics/no_flux_17.cc @@ -43,7 +43,12 @@ compare_constraints(const AffineConstraints &constraints1, { const double tol = 1.e-15; if (constraints1.n_constraints() != constraints2.n_constraints()) - return false; + { + deallog << "n_constraints_1 = " << constraints1.n_constraints() + << "n_constraints_2 = " << constraints2.n_constraints() + << std::endl; + return false; + } auto line1 = (constraints1.get_lines()).begin(); auto line2 = (constraints2.get_lines()).begin(); @@ -54,7 +59,11 @@ compare_constraints(const AffineConstraints &constraints1, const unsigned int line_index = line1->index; const unsigned int line_index2 = line2->index; if (line_index != line_index2) - return false; + { + deallog << " line_index_1 = " << line_index + << " line_index_2 = " << line_index2 << std::endl; + return false; + } const std::vector> *constraint_entries_1 = constraints1.get_constraint_entries(line_index); @@ -132,6 +141,9 @@ get_constraints_on_active_cells( VectorTools::compute_no_normal_flux_constraints( dofh, 0, no_normal_flux_boundaries, constraints, mapping); constraints.close(); + // deallog<<" ***active cell constraints: + // "< @@ -169,6 +181,9 @@ get_constraints_on_levels( level); user_level_constraints.close(); level_constraints[level].copy_from(user_level_constraints); + // deallog<<" ***level cell constraints: + // "< fe(FE_Q(2), dim); + + const MappingQ mapping(4); + Triangulation triangulation( + Triangulation::limit_level_difference_at_vertices); + GridGenerator::hyper_shell(triangulation, Point(), 0.5, 1.); + std::set no_flux_boundary{0}; + triangulation.begin_active()->set_refine_flag(); + triangulation.execute_coarsening_and_refinement(); + + DoFHandler dofh(triangulation); + dofh.distribute_dofs(fe); + + AffineConstraints constraints; + + VectorTools::compute_no_normal_flux_constraints( + dofh, 0, no_flux_boundary, constraints, mapping); + + constraints.print(deallog.get_file_stream()); + } } diff --git a/tests/numerics/no_flux_17.output b/tests/numerics/no_flux_17.output index 6196d0fe2d..7d7b0df189 100644 --- a/tests/numerics/no_flux_17.output +++ b/tests/numerics/no_flux_17.output @@ -58,3 +58,48 @@ DEAL::Level 1 OK DEAL:: dim 3 DEAL::Level 0 OK DEAL::Level 1 OK +DEAL:: adaptive hyper_shell: + 0 1: -0.7265425 + 3 2: -0.3249197 + 4 5: -0.7265425 + 7 6: -0.3249197 + 13 12: -0.7265425 + 15 14: -0.7265425 + 19 18: 0.3249197 + 21 20: 0.3249197 + 25 = 0 + 27 = 0 + 30 31: 0.7265425 + 32 33: 0.7265425 + 37 36: 0.7265425 + 39 38: 0.7265425 + 42 = 0 + 44 = 0 + 48 49: 0.3249197 + 50 51: 0.3249197 + 54 55: -0.7265425 + 56 57: -0.7265425 + 60 61: -0.3249197 + 62 63: -0.3249197 + 67 66: -0.3249197 + 69 68: -0.3249197 + 73 72: -0.7265425 + 75 74: -0.7265425 + 79 78: 0.3249197 + 81 80: 0.3249197 + 85 = 0 + 87 = 0 + 90 91: 0.7265425 + 92 93: 0.7265425 + 97 96: 0.7265425 + 99 98: 0.7265425 + 102 = 0 + 104 = 0 + 108 109: 0.3249197 + 110 111: 0.3249197 + 114 115: -0.3249197 + 124 125: -0.1583844 + 134 135: -0.5095254 + 140 141: -0.3249197 + 146 147: -0.1583844 + 152 153: -0.5095254 -- 2.39.5