From 1e8b060ea7e36d61a3afeb65dd85802e58428b36 Mon Sep 17 00:00:00 2001 From: Jiaqi Zhang Date: Thu, 4 Aug 2022 00:59:28 -0400 Subject: [PATCH] combine into one function --- .../numerics/vector_tools_constraints.h | 41 +- .../vector_tools_constraints.templates.h | 683 ++++-------------- .../numerics/vector_tools_constraints.inst.in | 17 +- tests/numerics/no_flux_16.cc | 123 +++- tests/numerics/no_flux_16.output | 138 +--- 5 files changed, 266 insertions(+), 736 deletions(-) diff --git a/include/deal.II/numerics/vector_tools_constraints.h b/include/deal.II/numerics/vector_tools_constraints.h index aa36668448..9b2380d688 100644 --- a/include/deal.II/numerics/vector_tools_constraints.h +++ b/include/deal.II/numerics/vector_tools_constraints.h @@ -51,6 +51,9 @@ namespace VectorTools * i.e., normal flux constraints where $\vec u$ is a vector-valued solution * variable and $\vec u_\Gamma$ is a prescribed vector field whose normal * component we want to be equal to the normal component of the solution. + * This function can also be used on level meshes in the multigrid method + * if @p refinement_edge_indices and @p level are provided, and the former + * can be obtained by MGConstrainedDoFs::get_refinement_edge_indices(). * These conditions have exactly the form handled by the * AffineConstraints class, in that they relate a linear * combination of boundary degrees of freedom to a corresponding @@ -288,7 +291,9 @@ namespace VectorTools #else .ReferenceCell::get_default_linear_mapping() #endif - )); + ), + const IndexSet & refinement_edge_indices = IndexSet(), + const unsigned int level = numbers::invalid_unsigned_int); /** * This function does the same as the @@ -296,6 +301,9 @@ namespace VectorTools * information), but for the simpler case of homogeneous normal-flux * constraints, i.e., for imposing the condition * $\vec u \cdot \vec n= 0$. This function is used in step-31 and step-32. + * This function can also be used on level meshes in the multigrid method + * if @p refinement_edge_indices and @p level are provided, and the former + * can be obtained by MGConstrainedDoFs::get_refinement_edge_indices(). * * @ingroup constraints * @@ -316,34 +324,9 @@ namespace VectorTools #else .ReferenceCell::get_default_linear_mapping() #endif - )); - - /** - * This function does the same as the - * compute_no_normal_flux_constraints(), but for the case of level meshes - * in the multigrid method. - * @ingroup constraints - * - * @see - * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" - */ - template - void - compute_no_normal_flux_constraints_on_level( - const DoFHandler & dof_handler, - const MGConstrainedDoFs & mg_constrained_dofs, - const unsigned int level, - const unsigned int first_vector_component, - const std::set &boundary_ids, - AffineConstraints & constraints, - const Mapping & mapping = - (ReferenceCells::get_hypercube() -#ifndef _MSC_VER - .template get_default_linear_mapping() -#else - .ReferenceCell::get_default_linear_mapping() -#endif - )); + ), + const IndexSet & refinement_edge_indices = IndexSet(), + const unsigned int level = numbers::invalid_unsigned_int); /** * Compute the constraints that correspond to boundary conditions of the diff --git a/include/deal.II/numerics/vector_tools_constraints.templates.h b/include/deal.II/numerics/vector_tools_constraints.templates.h index 08c2a44a70..17c05a85ee 100644 --- a/include/deal.II/numerics/vector_tools_constraints.templates.h +++ b/include/deal.II/numerics/vector_tools_constraints.templates.h @@ -473,13 +473,20 @@ namespace VectorTools const std::map *> & function_map, AffineConstraints & constraints, - const Mapping &mapping) + 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 @@ -533,30 +540,34 @@ namespace VectorTools using DoFToNormalsMap = std::multimap< internal::VectorDoFTuple, std::pair, - typename DoFHandler::active_cell_iterator>>; + typename DoFHandler::cell_iterator>>; std::map, Vector> dof_vector_to_b_values; DoFToNormalsMap dof_to_normals_map; // now loop over all cells and all faces - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); std::set::iterator b_id; - for (; cell != endc; ++cell) - if (!cell->is_artificial()) + 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::face_iterator face = + 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)); - face->get_dof_indices(face_dofs, cell->active_fe_index()); + + if (level != numbers::invalid_unsigned_int) + face->get_mg_dof_indices(mesh_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 FEFaceValues &fe_values = @@ -567,115 +578,122 @@ namespace VectorTools 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; - vector_dofs.dof_indices[0] = face_dofs[i]; + // 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(), + 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(); - // 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)); + 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)); #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; + 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 - } + } } // Now do something with the collected information. To this end, loop @@ -716,7 +734,7 @@ namespace VectorTools // vectors. the values of the map are pairs of normal vectors and // number of cells that have contributed using CellToNormalsMap = - std::map::active_cell_iterator, + std::map::cell_iterator, std::pair, unsigned int>>; CellToNormalsMap cell_to_normals_map; @@ -892,9 +910,9 @@ namespace VectorTools // use a std::list instead of a std::set (which would be more // natural) because std::set requires that the stored elements // are comparable with operator< - using CellContributions = std::map< - typename DoFHandler::active_cell_iterator, - std::list>>; + using CellContributions = + std::map::cell_iterator, + std::list>>; CellContributions cell_contributions; for (typename DoFToNormalsMap::const_iterator q = @@ -1279,7 +1297,9 @@ namespace VectorTools const unsigned int first_vector_component, const std::set &boundary_ids, AffineConstraints & constraints, - const Mapping & mapping) + const Mapping & mapping, + const IndexSet & refinement_edge_indices, + const unsigned int level) { Functions::ZeroFunction zero_function(dim); std::map *> function_map; @@ -1290,432 +1310,9 @@ namespace VectorTools boundary_ids, function_map, constraints, - mapping); - } - - template - void - compute_no_normal_flux_constraints_on_level( - const DoFHandler & dof_handler, - const MGConstrainedDoFs & mg_constrained_dofs, - const unsigned int level, - const unsigned int first_vector_component, - const std::set &boundary_ids, - AffineConstraints & constraints, - const Mapping & mapping) - { - // Copied from compute_nonzero_normal_flux_constraints() - // Only changed active_cell_iterator to cell_iterator in DoFToNormalsMap, - // CellToNormalsMap, CellContributions, and set inhomogeneity to 0. - const IndexSet &refinement_edge_indices = - mg_constrained_dofs.get_refinement_edge_indices(level); - - const auto & fe = dof_handler.get_fe(); - const std::vector> &unit_support_points = - fe.get_unit_face_support_points(); - const Quadrature quadrature(unit_support_points); - const unsigned int dofs_per_face = fe.dofs_per_face; - std::vector face_dofs(dofs_per_face); - - - FEFaceValues fe_face_values(mapping, - fe, - quadrature, - update_quadrature_points | - update_normal_vectors); - - std::set::iterator b_id; - using DoFToNormalsMap = std::multimap< - internal::VectorDoFTuple, - std::pair, - typename DoFHandler::cell_iterator>>; - DoFToNormalsMap dof_to_normals_map; - 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()) - { - typename DoFHandler::level_face_iterator face = - cell->face(face_no); - face->get_mg_dof_indices(level, face_dofs); - fe_face_values.reinit(cell, face_no); - - for (unsigned int i = 0; i < face_dofs.size(); ++i) - if (fe.face_system_to_component_index(i).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])) - { - const Point position = - fe_face_values.quadrature_point(i); - internal::VectorDoFTuple vector_dofs; - vector_dofs.dof_indices[0] = face_dofs[i]; - for (unsigned int k = 0; k < dofs_per_face; ++k) - if ((k != i) && - (quadrature.point(k) == quadrature.point(i)) && - (fe.face_system_to_component_index(k).first >= - first_vector_component) && - (fe.face_system_to_component_index(k).first < - first_vector_component + dim)) - vector_dofs.dof_indices - [fe.face_system_to_component_index(k).first - - first_vector_component] = face_dofs[k]; - - Tensor<1, dim> normal_vector = - cell->face(face_no)->get_manifold().normal_vector( - cell->face(face_no), position); - - // make sure the normal vector is pointing to the right - // direction, more dietails can be found in - // compute_nonzero_normal_flux_constraints() - if (normal_vector * fe_face_values.normal_vector(i) < 0) - normal_vector *= -1; - Assert(std::fabs(normal_vector.norm() - 1) < 1e-14, - ExcInternalError()); - - // remove small entries: - 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(); - - dof_to_normals_map.insert( - std::make_pair(vector_dofs, - std::make_pair(normal_vector, cell))); -#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 - } - } - // Now do something with the collected information. To this end, loop - // through all sets of pairs (dofs,normal_vector) and identify which - // entries belong to the same set of dofs and then do as described in the - // documentation, i.e. either average the normal vector or don't for this - // particular set of dofs - typename DoFToNormalsMap::const_iterator p = dof_to_normals_map.begin(); - - while (p != dof_to_normals_map.end()) - { - // first find the range of entries in the multimap that corresponds to - // the same vector-dof tuple. as usual, we define the range - // half-open. the first entry of course is 'p' - typename DoFToNormalsMap::const_iterator same_dof_range[2] = {p}; - for (++p; p != dof_to_normals_map.end(); ++p) - if (p->first != same_dof_range[0]->first) - { - same_dof_range[1] = p; - break; - } - if (p == dof_to_normals_map.end()) - same_dof_range[1] = dof_to_normals_map.end(); - -#ifdef DEBUG_NO_NORMAL_FLUX - std::cout << "For dof indices <" << p->first - << ">, found the following normals" << std::endl; - for (typename DoFToNormalsMap::const_iterator q = same_dof_range[0]; - q != same_dof_range[1]; - ++q) - std::cout << " " << q->second.first << " from cell " - << q->second.second << std::endl; -#endif - - - // now compute the reverse mapping: for each of the cells that - // contributed to the current set of vector dofs, add up the normal - // vectors. the values of the map are pairs of normal vectors and - // number of cells that have contributed - using CellToNormalsMap = - std::map::cell_iterator, - std::pair, unsigned int>>; - - CellToNormalsMap cell_to_normals_map; - for (typename DoFToNormalsMap::const_iterator q = same_dof_range[0]; - q != same_dof_range[1]; - ++q) - if (cell_to_normals_map.find(q->second.second) == - cell_to_normals_map.end()) - cell_to_normals_map[q->second.second] = - std::make_pair(q->second.first, 1U); - else - { - const Tensor<1, dim> old_normal = - cell_to_normals_map[q->second.second].first; - const unsigned int old_count = - cell_to_normals_map[q->second.second].second; - - Assert(old_count > 0, ExcInternalError()); - - // in the same entry, store again the now averaged normal vector - // and the new count - cell_to_normals_map[q->second.second] = - std::make_pair((old_normal * old_count + q->second.first) / - (old_count + 1), - old_count + 1); - } - Assert(cell_to_normals_map.size() >= 1, ExcInternalError()); - -#ifdef DEBUG_NO_NORMAL_FLUX - std::cout << " cell_to_normals_map:" << std::endl; - for (typename CellToNormalsMap::const_iterator x = - cell_to_normals_map.begin(); - x != cell_to_normals_map.end(); - ++x) - std::cout << " " << x->first << " -> (" << x->second.first << ',' - << x->second.second << ')' << std::endl; -#endif - - // count the maximum number of contributions from each cell - unsigned int max_n_contributions_per_cell = 1; - for (typename CellToNormalsMap::const_iterator x = - cell_to_normals_map.begin(); - x != cell_to_normals_map.end(); - ++x) - max_n_contributions_per_cell = - std::max(max_n_contributions_per_cell, x->second.second); - - // verify that each cell can have only contributed at most dim times, - // since that is the maximum number of faces that come together at a - // single place - Assert(max_n_contributions_per_cell <= dim, ExcInternalError()); - - switch (max_n_contributions_per_cell) - { - case 1: - { - // compute the average normal vector from all the ones that - // have the same set of dofs. we could add them up and divide - // them by the number of additions, or simply normalize them - // right away since we want them to have unit length anyway - Tensor<1, dim> normal; - for (typename CellToNormalsMap::const_iterator x = - cell_to_normals_map.begin(); - x != cell_to_normals_map.end(); - ++x) - normal += x->second.first; - normal /= normal.norm(); - - // normalize again - for (unsigned int d = 0; d < dim; ++d) - if (std::fabs(normal[d]) < 1e-13) - normal[d] = 0; - normal /= normal.norm(); - - const auto &dof_indices = same_dof_range[0]->first; - internal::add_constraint(dof_indices, - normal, - constraints, - 0.0); - - break; - } - - case dim: - { - // assert that indeed only a single cell has contributed - Assert(cell_to_normals_map.size() == 1, ExcInternalError()); - - // check linear independence by computing the determinant of - // the matrix created from all the normal vectors. if they are - // linearly independent, then the determinant is nonzero. if - // they are orthogonal, then the matrix is in fact equal to 1 - // (since they are all unit vectors); make sure the - // determinant is larger than 1e-3 to avoid cases where cells - // are degenerate - { - Tensor<2, dim> t; - - typename DoFToNormalsMap::const_iterator x = - same_dof_range[0]; - for (unsigned int i = 0; i < dim; ++i, ++x) - for (unsigned int j = 0; j < dim; ++j) - t[i][j] = x->second.first[j]; - - Assert( - std::fabs(determinant(t)) > 1e-3, - ExcMessage( - "Found a set of normal vectors that are nearly collinear.")); - } - - // so all components of this vector dof are constrained. enter - // this into the AffineConstraints object - // - // ignore dofs already constrained - const auto &dof_indices = same_dof_range[0]->first; - - for (unsigned int i = 0; i < dim; ++i) - if (!constraints.is_constrained( - same_dof_range[0]->first.dof_indices[i]) && - constraints.can_store_line( - same_dof_range[0]->first.dof_indices[i])) - { - const types::global_dof_index line = - dof_indices.dof_indices[i]; - constraints.add_line(line); - } - - break; - } - - // this is the case of an edge contribution in 3d, i.e. the vector - // is constrained in two directions but not the third. - default: - { - Assert(dim >= 3, ExcNotImplemented()); - Assert(max_n_contributions_per_cell == 2, ExcInternalError()); - - // as described in the documentation, let us first collect - // what each of the cells contributed at the current point. we - // use a std::list instead of a std::set (which would be more - // natural) because std::set requires that the stored elements - // are comparable with operator< - using CellContributions = - std::map::cell_iterator, - std::list>>; - CellContributions cell_contributions; - - for (typename DoFToNormalsMap::const_iterator q = - same_dof_range[0]; - q != same_dof_range[1]; - ++q) - cell_contributions[q->second.second].push_back( - q->second.first); - Assert(cell_contributions.size() >= 1, ExcInternalError()); - - // now for each cell that has contributed determine the number - // of normal vectors it has contributed. we currently only - // implement if this is dim-1 for all cells (if a single cell - // has contributed dim, or if all adjacent cells have - // contributed 1 normal vector, this is already handled - // above). - // - // we only implement the case that all cells contribute - // dim-1 because we assume that we are following an edge - // of the domain (think: we are looking at a vertex - // located on one of the edges of a refined cube where the - // boundary indicators of the two adjacent faces of the - // cube are both listed in the set of boundary indicators - // passed to this function). in that case, all cells along - // that edge of the domain are assumed to have contributed - // dim-1 normal vectors. however, there are cases where - // this assumption is not justified (see the lengthy - // explanation in test no_flux_12.cc) and in those cases - // we simply ignore the cell that contributes only - // once. this is also discussed at length in the - // documentation of this function. - // - // for each contributing cell compute the tangential vector - // that remains unconstrained - std::list> tangential_vectors; - for (typename CellContributions::const_iterator contribution = - cell_contributions.begin(); - contribution != cell_contributions.end(); - ++contribution) - { -#ifdef DEBUG_NO_NORMAL_FLUX - std::cout - << " Treating edge case with dim-1 contributions." - << std::endl - << " Looking at cell " << contribution->first - << " which has contributed these normal vectors:" - << std::endl; - for (typename std::list>::const_iterator t = - contribution->second.begin(); - t != contribution->second.end(); - ++t) - std::cout << " " << *t << std::endl; -#endif - - // as mentioned above, simply ignore cells that only - // contribute once - if (contribution->second.size() < dim - 1) - continue; - - Tensor<1, dim> normals[dim - 1]; - { - unsigned int index = 0; - for (typename std::list>::const_iterator - t = contribution->second.begin(); - t != contribution->second.end(); - ++t, ++index) - normals[index] = *t; - Assert(index == dim - 1, ExcInternalError()); - } - - // calculate the tangent as the outer product of the - // normal vectors. since these vectors do not need to be - // orthogonal (think, for example, the case of the - // deal.II/no_flux_07 test: a sheared cube in 3d, with Q2 - // elements, where we have constraints from the two normal - // vectors of two faces of the sheared cube that are not - // perpendicular to each other), we have to normalize the - // outer product - Tensor<1, dim> tangent; - switch (dim) - { - case 3: - // take cross product between normals[0] and - // normals[1]. write it in the current form (with - // [dim-2]) to make sure that compilers don't warn - // about out-of-bounds accesses -- the warnings are - // bogus since we get here only for dim==3, but at - // least one isn't quite smart enough to notice this - // and warns when compiling the function in 2d - tangent = - cross_product_3d(normals[0], normals[dim - 2]); - break; - default: - Assert(false, ExcNotImplemented()); - } - - Assert( - std::fabs(tangent.norm()) > 1e-12, - ExcMessage( - "Two normal vectors from adjacent faces are almost " - "parallel.")); - tangent /= tangent.norm(); - - tangential_vectors.push_back(tangent); - } - - // go through the list of tangents and make sure that they all - // roughly point in the same direction as the first one (i.e. - // have an angle less than 90 degrees); if they don't then - // flip their sign - { - const Tensor<1, dim> first_tangent = - tangential_vectors.front(); - typename std::list>::iterator t = - tangential_vectors.begin(); - ++t; - for (; t != tangential_vectors.end(); ++t) - if (*t * first_tangent < 0) - *t *= -1; - } - - // now compute the average tangent and normalize it - Tensor<1, dim> average_tangent; - for (typename std::list>::const_iterator t = - tangential_vectors.begin(); - t != tangential_vectors.end(); - ++t) - average_tangent += *t; - average_tangent /= average_tangent.norm(); - - const auto &dof_indices = same_dof_range[0]->first; - internal::add_tangentiality_constraints(dof_indices, - average_tangent, - constraints); - } - } - } + mapping, + refinement_edge_indices, + level); } diff --git a/source/numerics/vector_tools_constraints.inst.in b/source/numerics/vector_tools_constraints.inst.in index 4c112dd8c0..9e47b54f84 100644 --- a/source/numerics/vector_tools_constraints.inst.in +++ b/source/numerics/vector_tools_constraints.inst.in @@ -32,17 +32,10 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const std::map *> & function_map, AffineConstraints & constraints, - const Mapping &mapping); + const Mapping &mapping, + const IndexSet & refinement_edge_indices, + const unsigned int level); - template void - compute_no_normal_flux_constraints_on_level( - const DoFHandler &dof_handler, - const MGConstrainedDoFs & mg_constrained_dofs, - const unsigned int level, - const unsigned int first_vector_component, - const std::set & boundary_ids, - AffineConstraints & constraints, - const Mapping & mapping); template void compute_nonzero_tangential_flux_constraints( @@ -60,7 +53,9 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const unsigned int first_vector_component, const std::set & boundary_ids, AffineConstraints & constraints, - const Mapping & mapping); + const Mapping & mapping, + const IndexSet & refinement_edge_indices, + const unsigned int level); template void compute_normal_flux_constraints( diff --git a/tests/numerics/no_flux_16.cc b/tests/numerics/no_flux_16.cc index f9acc1fe28..b979b5380d 100644 --- a/tests/numerics/no_flux_16.cc +++ b/tests/numerics/no_flux_16.cc @@ -37,10 +37,64 @@ #include "../tests.h" +bool +compare_constraints(const AffineConstraints &constraints1, + const AffineConstraints &constraints2) +{ + if (constraints1.n_constraints() != constraints2.n_constraints()) + return false; + + for (auto line : constraints1.get_lines()) + { + const unsigned int line_index = line.index; + const std::vector> + *constraint_entries_1 = constraints1.get_constraint_entries(line_index); + const std::vector> + *constraint_entries_2 = constraints2.get_constraint_entries(line_index); + if (constraint_entries_1->size() != constraint_entries_2->size() || + (constraints1.get_inhomogeneity(line_index) != + constraints2.get_inhomogeneity(line_index))) + return false; + for (unsigned int i = 0; i < constraint_entries_1->size(); ++i) + { + if ((*constraint_entries_1)[i].first != + (*constraint_entries_2)[i].first) + return false; + if ((*constraint_entries_1)[i].second != + (*constraint_entries_2)[i].second) + return false; + } + } + + return true; +} + template void -run(const Triangulation & triangulation, - const std::set &no_flux_boundary) +constraints_on_active_cells( + const Triangulation & tria, + const std::set &no_normal_flux_boundaries, + AffineConstraints & constraints) +{ + SphericalManifold spherical; + + MappingQ mapping(4); + + FESystem fe(FE_Q(1), dim); + DoFHandler dofh(tria); + + dofh.distribute_dofs(fe); + + VectorTools::compute_no_normal_flux_constraints( + dofh, 0, no_normal_flux_boundaries, constraints, mapping); +} + +template +void +constraints_on_levels( + const Triangulation & triangulation, + const std::set & no_flux_boundary, + MGLevelObject> &level_constraints) { FESystem fe(FE_Q(1), dim); DoFHandler dof_handler(triangulation); @@ -58,25 +112,62 @@ run(const Triangulation & triangulation, { AffineConstraints user_level_constraints; - VectorTools::compute_no_normal_flux_constraints_on_level( - dof_handler, - mg_constrained_dofs, - level, - 0, - no_flux_boundary, - user_level_constraints, - mapping); + const IndexSet &refinement_edge_indices = + mg_constrained_dofs.get_refinement_edge_indices(level); - user_level_constraints.print(deallog.get_file_stream()); + VectorTools::compute_no_normal_flux_constraints(dof_handler, + 0, + no_flux_boundary, + user_level_constraints, + mapping, + refinement_edge_indices, + level); - deallog.get_file_stream() << std::flush; + // user_level_constraints.print(deallog.get_file_stream()); + + // deallog.get_file_stream() << std::flush; user_level_constraints.close(); + level_constraints[level].copy_from(user_level_constraints); + } +} - deallog << "Level " << level << " OK" << std::endl; +template +void +run(Triangulation & triangulation, + std::set no_flux_boundary) +{ + const unsigned int n_levels = 2; + std::vector> ref_constraints(n_levels); + + constraints_on_active_cells(triangulation, + no_flux_boundary, + ref_constraints[0]); + for (unsigned int level = 1; level < n_levels; ++level) + { + triangulation.refine_global(1); + constraints_on_active_cells(triangulation, + no_flux_boundary, + ref_constraints[level]); + } + + MGLevelObject> mg_level_constraints(0, + n_levels - 1); + constraints_on_levels(triangulation, + no_flux_boundary, + mg_level_constraints); + + deallog << " dim " << dim << std::endl; + for (unsigned int i = 0; i < n_levels; ++i) + { + if (compare_constraints(ref_constraints[i], mg_level_constraints[i])) + deallog << "Level " << i << " OK" << std::endl; + else + deallog << "Level " << i << " failed" << std::endl; } } + int main() { @@ -89,17 +180,15 @@ main() Triangulation triangulation( Triangulation::limit_level_difference_at_vertices); GridGenerator::quarter_hyper_ball(triangulation); - triangulation.refine_global(1); std::set no_flux_boundary{0, 1}; - run(triangulation, no_flux_boundary); + run(triangulation, no_flux_boundary); } { const unsigned int dim = 3; Triangulation triangulation( Triangulation::limit_level_difference_at_vertices); GridGenerator::quarter_hyper_ball(triangulation); - triangulation.refine_global(1); std::set no_flux_boundary{0, 1}; - run(triangulation, no_flux_boundary); + run(triangulation, no_flux_boundary); } } diff --git a/tests/numerics/no_flux_16.output b/tests/numerics/no_flux_16.output index 5dee504c58..8fc6891ed7 100644 --- a/tests/numerics/no_flux_16.output +++ b/tests/numerics/no_flux_16.output @@ -1,141 +1,7 @@ - 0 = 0 - 1 = 0 - 3 = 0 - 4 = 0 - 8 = 0 - 9 = 0 - 11 10: -1.0000000 - 12 = 0 - 13 = 0 +DEAL:: dim 2 DEAL::Level 0 OK - 0 = 0 - 1 = 0 - 3 = 0 - 4 = 0 - 9 = 0 - 12 = 0 - 18 = 0 - 19 = 0 - 20 21: -0.4142136 - 23 = 0 - 27 26: -1.0000000 - 30 = 0 - 31 = 0 - 32 = 0 - 35 34: -0.4142136 DEAL::Level 1 OK - 0 = 0 - 1 = 0 - 2 = 0 - 4 = 0 - 5 = 0 - 6 = 0 - 8 = 0 - 11 = 0 - 12 = 0 - 13 = 0 - 16 = 0 - 18 = 0 - 24 = 0 - 25 = 0 - 26 = 0 - 28 27: -1.0000000 - 29 = 0 - 31 = 0 - 32 30: -1.0000000 - 35 33: -1.0000000 - 35 34: -1.0000000 - 36 = 0 - 37 = 0 - 38 = 0 - 39 = 0 - 41 40: -1.0000000 - 42 = 0 - 43 = 0 - 44 = 0 +DEAL:: dim 3 DEAL::Level 0 OK - 0 = 0 - 1 = 0 - 2 = 0 - 4 = 0 - 5 = 0 - 6 = 0 - 8 = 0 - 11 = 0 - 12 = 0 - 13 = 0 - 16 = 0 - 18 = 0 - 25 = 0 - 26 = 0 - 29 = 0 - 31 = 0 - 36 = 0 - 38 = 0 - 41 = 0 - 42 = 0 - 50 = 0 - 54 = 0 - 55 = 0 - 58 = 0 - 60 = 0 - 67 = 0 - 72 = 0 - 81 = 0 - 82 = 0 - 83 = 0 - 84 85: -0.4142136 - 86 = 0 - 88 = 0 - 89 = 0 - 92 = 0 - 93 95: -0.4142136 - 94 = 0 - 96 97: -0.4237876 - 96 98: -0.4237876 - 100 = 0 - 106 105: -1.0000000 - 107 = 0 - 110 = 0 - 112 111: -1.0000000 - 112 113: -0.4494897 - 118 = 0 - 119 117: -1.0000000 - 122 120: -1.0000000 - 122 121: -0.4494897 - 124 = 0 - 131 129: -1.0000000 - 131 130: -1.0000000 - 135 = 0 - 136 = 0 - 137 = 0 - 138 = 0 - 140 = 0 - 142 141: -0.4142136 - 143 = 0 - 146 = 0 - 147 = 0 - 148 149: -0.4142136 - 150 = 0 - 154 153: -0.4237876 - 154 155: -0.4237876 - 159 = 0 - 161 160: -1.0000000 - 162 = 0 - 167 165: -0.4494897 - 167 166: -1.0000000 - 171 = 0 - 172 = 0 - 175 = 0 - 177 = 0 - 183 = 0 - 184 = 0 - 185 = 0 - 187 = 0 - 188 186: -0.4142136 - 189 = 0 - 191 190: -0.4142136 - 194 192: -0.4237876 - 194 193: -0.4237876 DEAL::Level 1 OK -- 2.39.5