From 4497b67c20a3a1b1fae6314b29e6c7cbbbf659da Mon Sep 17 00:00:00 2001 From: Jiaqi Zhang Date: Thu, 28 Jul 2022 10:57:28 -0400 Subject: [PATCH] no normal flux on levels --- .../numerics/vector_tools_constraints.h | 29 ++ .../vector_tools_constraints.templates.h | 427 ++++++++++++++++++ .../numerics/vector_tools_constraints.inst.in | 10 + tests/numerics/no_flux_16.cc | 105 +++++ tests/numerics/no_flux_16.output | 141 ++++++ 5 files changed, 712 insertions(+) create mode 100644 tests/numerics/no_flux_16.cc create mode 100644 tests/numerics/no_flux_16.output diff --git a/include/deal.II/numerics/vector_tools_constraints.h b/include/deal.II/numerics/vector_tools_constraints.h index 5b3aa88d08..aa36668448 100644 --- a/include/deal.II/numerics/vector_tools_constraints.h +++ b/include/deal.II/numerics/vector_tools_constraints.h @@ -21,6 +21,8 @@ #include +#include + #include #include @@ -316,6 +318,33 @@ namespace VectorTools #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 + )); + /** * Compute the constraints that correspond to boundary conditions of the * form $\vec u \times \vec n=\vec u_\Gamma \times \vec n$, i.e., tangential diff --git a/include/deal.II/numerics/vector_tools_constraints.templates.h b/include/deal.II/numerics/vector_tools_constraints.templates.h index a556af03fe..08c2a44a70 100644 --- a/include/deal.II/numerics/vector_tools_constraints.templates.h +++ b/include/deal.II/numerics/vector_tools_constraints.templates.h @@ -1293,6 +1293,433 @@ namespace VectorTools 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); + } + } + } + } + + + template void compute_normal_flux_constraints( diff --git a/source/numerics/vector_tools_constraints.inst.in b/source/numerics/vector_tools_constraints.inst.in index d23919f0dd..4c112dd8c0 100644 --- a/source/numerics/vector_tools_constraints.inst.in +++ b/source/numerics/vector_tools_constraints.inst.in @@ -34,6 +34,16 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) AffineConstraints & constraints, const Mapping &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); + template void compute_nonzero_tangential_flux_constraints( const DoFHandler &dof_handler, diff --git a/tests/numerics/no_flux_16.cc b/tests/numerics/no_flux_16.cc new file mode 100644 index 0000000000..f9acc1fe28 --- /dev/null +++ b/tests/numerics/no_flux_16.cc @@ -0,0 +1,105 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2012 - 2022 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// add tests for compute_no_normal_flux_constraints_on_level() + +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include "../tests.h" + +template +void +run(const Triangulation & triangulation, + const std::set &no_flux_boundary) +{ + FESystem fe(FE_Q(1), dim); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe); + dof_handler.distribute_mg_dofs(); + + MGConstrainedDoFs mg_constrained_dofs; + mg_constrained_dofs.initialize(dof_handler); + + const unsigned int n_levels = triangulation.n_global_levels(); + + MappingQ mapping(4); + + for (unsigned int level = 0; level < n_levels; ++level) + { + 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); + + user_level_constraints.print(deallog.get_file_stream()); + + deallog.get_file_stream() << std::flush; + user_level_constraints.close(); + + deallog << "Level " << level << " OK" << std::endl; + } +} + + +int +main() +{ + initlog(); + deallog.get_file_stream().precision(7); + deallog.get_file_stream().setf(std::ios::fixed); + + { + const unsigned int dim = 2; + 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); + } + { + 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); + } +} diff --git a/tests/numerics/no_flux_16.output b/tests/numerics/no_flux_16.output new file mode 100644 index 0000000000..5dee504c58 --- /dev/null +++ b/tests/numerics/no_flux_16.output @@ -0,0 +1,141 @@ + + 0 = 0 + 1 = 0 + 3 = 0 + 4 = 0 + 8 = 0 + 9 = 0 + 11 10: -1.0000000 + 12 = 0 + 13 = 0 +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::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