From 17da853c7ae48dccac4f014ab28bc45ac9c2d9ef Mon Sep 17 00:00:00 2001 From: David Wells Date: Thu, 19 Oct 2023 10:05:57 -0400 Subject: [PATCH] Move topological GridTools functions to a new file. --- include/deal.II/grid/grid_tools.h | 417 +---- include/deal.II/grid/grid_tools_topology.h | 460 +++++ source/grid/CMakeLists.txt | 2 + source/grid/grid_tools.cc | 1875 ------------------- source/grid/grid_tools.inst.in | 68 - source/grid/grid_tools_topology.cc | 1920 ++++++++++++++++++++ source/grid/grid_tools_topology.inst.in | 122 ++ 7 files changed, 2506 insertions(+), 2358 deletions(-) create mode 100644 include/deal.II/grid/grid_tools_topology.h create mode 100644 source/grid/grid_tools_topology.cc create mode 100644 source/grid/grid_tools_topology.inst.in diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index d8fe636706..e51dc09fb5 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -33,6 +33,7 @@ #include #include +#include #include #include #include @@ -51,8 +52,6 @@ #include #include -#include - #ifdef DEAL_II_WITH_ZLIB # include # include @@ -61,7 +60,7 @@ #endif #include -#include +#include #include #ifdef DEAL_II_HAVE_CXX20 @@ -133,150 +132,10 @@ namespace internal */ namespace GridTools { - /** - * @name Functions supporting the creation of meshes - */ - /** @{ */ - - /** - * Remove vertices that are not referenced by any of the cells. This - * function is called by all GridIn::read_* functions to eliminate - * vertices that are listed in the input files but are not used by the cells - * in the input file. While these vertices should not be in the input from - * the beginning, they sometimes are, most often when some cells have been - * removed by hand without wanting to update the vertex lists, as they might - * be lengthy. - * - * This function is called by all GridIn::read_* functions as the - * triangulation class requires them to be called with used vertices only. - * This is so, since the vertices are copied verbatim by that class, so we - * have to eliminate unused vertices beforehand. - * - * Not implemented for the codimension one case. - */ - template - void - delete_unused_vertices(std::vector> &vertices, - std::vector> &cells, - SubCellData &subcelldata); - - /** - * Remove vertices that are duplicated, due to the input of a structured - * grid, for example. If these vertices are not removed, the faces bounded - * by these vertices become part of the boundary, even if they are in the - * interior of the mesh. - * - * This function is called by some GridIn::read_* functions. Only - * the vertices with indices in @p considered_vertices are tested for - * equality. This speeds up the algorithm, which is, for worst-case hyper - * cube geometries $O(N^{3/2})$ in 2d and $O(N^{5/3})$ in 3d: quite slow. - * However, if you wish to consider all vertices, simply pass an empty - * vector. In that case, the function fills @p considered_vertices with all - * vertices. - * - * Two vertices are considered equal if their difference in each coordinate - * direction is less than @p tol. This implies that nothing happens if - * the tolerance is set to zero. - */ - template - void - delete_duplicated_vertices(std::vector> &all_vertices, - std::vector> &cells, - SubCellData &subcelldata, - std::vector &considered_vertices, - const double tol = 1e-12); - - /** - * Remove vertices that are duplicated. - * - * Two vertices are considered equal if their difference in each coordinate - * direction is less than @p tol. This implies that nothing happens if - * the tolerance is set to zero. - */ - template - void - delete_duplicated_vertices(std::vector> &vertices, - const double tol = 1e-12); - - /** - * Grids generated by grid generators may have an orientation of cells which - * is the inverse of the orientation required by deal.II. - * - * In 2d and 3d this function checks whether all cells have negative or - * positive measure/volume. In the former case, all cells are inverted. It - * does nothing in 1d. - * - * The inversion of cells might also work when only a subset of all cells - * have negative volume. However, grids consisting of a mixture of negative - * and positively oriented cells are very likely to be broken. Therefore, an - * exception is thrown, in case cells are not uniformly oriented. - * - * @note This function should be called before GridTools::consistently_order_cells(). - * - * @param all_vertices The vertices of the mesh. - * @param cells The array of CellData objects that describe the mesh's topology. - */ - template - void - invert_all_negative_measure_cells( - const std::vector> &all_vertices, - std::vector> &cells); - - /** - * Check the given cells and inverts any cell that is considered to have - * negative measure/volume in the orientation required by deal.II. - * - * This function is identical to invert_all_negative_measure_cells() except it - * does not throw an error if only some of the cells are inverted. Instead, - * this function returns how many cells were inverted. Additionally, it will - * always throw an exception outside of codimension 0. - */ - template - std::size_t - invert_cells_with_negative_measure( - const std::vector> &all_vertices, - std::vector> &cells); - - /** - * Given a vector of CellData objects describing a mesh, reorder their - * vertices so that all lines are consistently oriented. - * - * The expectations on orientation and a discussion of this function are - * available in the - * @ref reordering "reordering module". - * - * @param cells The array of CellData objects that describe the mesh's topology. - * @ingroup reordering - */ - template - void - consistently_order_cells(std::vector> &cells); - - /** @} */ /** * @name Rotating, stretching and otherwise transforming meshes */ /** @{ */ - /** - * Return the arrays that define the coarse mesh of a Triangulation. This - * function is the inverse of Triangulation::create_triangulation() in the - * sense that if one called this function on a triangulation, then that - * triangulation could be recreated by some kind of refinement from the - * results of this function. - * - * The return value is a tuple with the vector of vertices, the vector of - * cells, and a SubCellData structure. The latter contains additional - * information about faces and lines. These three objects are exactly - * the arguments to Triangulation::create_triangulation(). - * - * This function is useful in cases where one needs to deconstruct a - * Triangulation or manipulate the numbering of the vertices in some way: an - * example is GridGenerator::merge_triangulations(). - */ - template - std:: - tuple>, std::vector>, SubCellData> - get_coarse_mesh_description(const Triangulation &tria); /** * Transform the vertices of the given triangulation by applying the @@ -497,15 +356,6 @@ namespace GridTools const Function *coefficient = nullptr, const bool solve_for_absolute_positions = false); - /** - * Return a std::map with all vertices of faces located in the boundary - * - * @param[in] tria The Triangulation object. - */ - template - std::map> - get_all_vertices_at_boundary(const Triangulation &tria); - /** * Scale the entire triangulation by the given factor. To preserve the * orientation of the triangulation, the factor must be positive. @@ -545,75 +395,6 @@ namespace GridTools const bool keep_boundary = true, const unsigned int seed = boost::random::mt19937::default_seed); - /** - * Remove hanging nodes from a grid. If the @p isotropic parameter is set - * to @p false (default) this function detects cells with hanging nodes and - * refines the neighbours in the direction that removes hanging nodes. - * If the @p isotropic parameter is set - * to @p true, the neighbours refinement is made in each directions. - * In order to remove all hanging nodes this procedure has to be repeated: - * this could require a large number of iterations. - * To avoid this a max number (@p max_iterations) of iteration is provided. - * - * Consider the following grid: - * @image html remove_hanging_nodes-hanging.png - * - * @p isotropic == @p false would return: - * @image html remove_hanging_nodes-aniso.png - * - * @p isotropic == @p true would return: - * @image html remove_hanging_nodes-isotro.png - * - * @param[in,out] tria Triangulation to refine. - * - * @param[in] isotropic If true refine cells in each directions, otherwise - * (default value) refine the cell in the direction that removes hanging node. - * - * @param[in] max_iterations At each step only closest cells to hanging nodes - * are refined. The code may require a lot of iterations to remove all - * hanging nodes. @p max_iterations is the maximum number of iteration - * allowed. If @p max_iterations == numbers::invalid_unsigned_int this - * function continues refining until there are no hanging nodes. - * - * @note In the case of parallel codes, this function should be combined - * with GridGenerator::flatten_triangulation. - */ - template - void - remove_hanging_nodes(Triangulation &tria, - const bool isotropic = false, - const unsigned int max_iterations = 100); - - /** - * Refine a mesh anisotropically such that the resulting mesh is composed by - * cells with maximum ratio between dimensions less than @p max_ratio. - * This procedure requires an algorithm that may not terminate. Consequently, - * it is possible to set a maximum number of iterations through the - * @p max_iterations parameter. - * - * Starting from a cell like this: - * @image html remove_anisotropy-coarse.png - * - * This function would return: - * @image html remove_anisotropy-refined.png - * - * @param[in,out] tria Triangulation to refine. - * - * @param[in] max_ratio Maximum value allowed among the ratio between - * the dimensions of each cell. - * - * @param[in] max_iterations Maximum number of iterations allowed. - * - * @note In the case of parallel codes, this function should be combined - * with GridGenerator::flatten_triangulation and - * GridTools::remove_hanging_nodes. - */ - template - void - remove_anisotropy(Triangulation &tria, - const double max_ratio = 1.6180339887, - const unsigned int max_iterations = 5); - /** * Analyze the boundary cells of a mesh, and if one cell is found at * a corner position (with dim adjacent faces on the boundary), and its @@ -1146,55 +927,6 @@ namespace GridTools } // namespace internal - /** - * Return a map `vertex index -> Point` containing the used - * vertices of the given `container`. The key of the returned map (i.e., - * the first element of the pair above) is the global index in the - * triangulation, whereas the value of each pair is the physical - * location of the corresponding vertex. The used vertices are obtained by - * looping over all cells, - * and querying for each cell where its vertices are through the (optional) - * `mapping` argument. - * - * In serial Triangulation objects and parallel::shared::Triangulation - * objects, the size of the returned map - * equals Triangulation::n_used_vertices() (not Triangulation::n_vertices()). - * Note that in parallel::distributed::Triangulation objects, only vertices in - * locally owned cells and ghost cells are returned, as for all other vertices - * their real location might not be known (e.g. for distributed computations - * using MappingQEulerian). - * - * If you use the default `mapping`, the returned map satisfies the following - * equality: - * - * @code - * const auto used_vertices = extract_used_vertices(tria); - * auto all_vertices = tria.get_vertices(); - * - * for(const auto &id_and_v : used_vertices) - * all_vertices[id_and_v.first] == id_and_v.second; // true - * @endcode - * - * Notice that the above is not satisfied for mappings that change the - * location of vertices, like MappingQEulerian. - * - * @ref ConceptMeshType "MeshType concept". - * @param container The container to extract vertices from. - * @param mapping The mapping to use to compute the points locations. - */ - template - std::map> - extract_used_vertices( - const Triangulation &container, - const Mapping &mapping = - (ReferenceCells::get_hypercube() -#ifndef _MSC_VER - .template get_default_linear_mapping() -#else - .ReferenceCell::get_default_linear_mapping() -#endif - )); - /** * Find and return the index of the closest vertex to a given point in the * map of vertices passed as the first argument. @@ -2040,20 +1772,6 @@ namespace GridTools const RTree, unsigned int>> &covering_rtree, const std::vector> &points); - - /** - * Return the adjacent cells of all the vertices. If a vertex is also a - * hanging node, the associated coarse cell is also returned. The vertices - * are ordered by the vertex index. This is the number returned by the - * function cell-@>vertex_index(). Notice that only the indices - * marked in the array returned by - * Triangulation::get_used_vertices() are used. - */ - template - std::vector< - std::set::active_cell_iterator>> - vertex_to_cell_map(const Triangulation &triangulation); - /** * Return a vector of normalized tensors for each vertex-cell combination of * the output of GridTools::vertex_to_cell_map() (which is expected as input @@ -2118,49 +1836,6 @@ namespace GridTools */ /** @{ */ - /** - * Produce a sparsity pattern in which nonzero entries indicate that two - * cells are connected via a common face. The diagonal entries of the - * sparsity pattern are also set. - * - * The rows and columns refer to the cells as they are traversed in their - * natural order using cell iterators. - */ - template - void - get_face_connectivity_of_cells( - const Triangulation &triangulation, - DynamicSparsityPattern &connectivity); - - /** - * Produce a sparsity pattern in which nonzero entries indicate that two - * cells are connected via a common vertex. The diagonal entries of the - * sparsity pattern are also set. - * - * The rows and columns refer to the cells as they are traversed in their - * natural order using cell iterators. - */ - template - void - get_vertex_connectivity_of_cells( - const Triangulation &triangulation, - DynamicSparsityPattern &connectivity); - - /** - * Produce a sparsity pattern for a given level mesh in which nonzero entries - * indicate that two cells are connected via a common vertex. The diagonal - * entries of the sparsity pattern are also set. - * - * The rows and columns refer to the cells as they are traversed in their - * natural order using cell iterators. - */ - template - void - get_vertex_connectivity_of_cells_on_level( - const Triangulation &triangulation, - const unsigned int level, - DynamicSparsityPattern &connectivity); - /** * Use graph partitioner to partition the active cells making up the entire * domain. After calling this function, the subdomain ids of all active cells @@ -2393,85 +2068,6 @@ namespace GridTools std::vector get_locally_owned_vertices(const Triangulation &triangulation); - /** @} */ - /** - * @name Comparing different meshes - */ - /** @{ */ - - /** - * Given two meshes (i.e. objects of type Triangulation or DoFHandler) that - * are based on the same coarse mesh, this function figures out a set of cells - * that are matched between the two meshes and where at most one of the meshes - * is more refined on this cell. In other words, it finds the smallest cells - * that are common to both meshes, and that together completely cover the - * domain. - * - * This function is useful, for example, in time-dependent or nonlinear - * application, where one has to integrate a solution defined on one mesh - * (e.g., the one from the previous time step or nonlinear iteration) - * against the shape functions of another mesh (the next time step, the next - * nonlinear iteration). If, for example, the new mesh is finer, then one - * has to obtain the solution on the coarse mesh (mesh_1) and interpolate it - * to the children of the corresponding cell of mesh_2. Conversely, if the - * new mesh is coarser, one has to express the coarse cell shape function by - * a linear combination of fine cell shape functions. In either case, one - * needs to loop over the finest cells that are common to both - * triangulations. This function returns a list of pairs of matching - * iterators to cells in the two meshes that can be used to this end. - * - * Note that the list of these iterators is not necessarily ordered, and - * does also not necessarily coincide with the order in which cells are - * traversed in one, or both, of the meshes given as arguments. - * - * @tparam MeshType A type that satisfies the requirements of the - * @ref ConceptMeshType "MeshType concept". - * - * @note This function can only be used with - * parallel::distributed::Triangulation when both meshes use the same - * Triangulation since, with a distributed Triangulation, not all cells are - * stored locally, so the resulting list may not cover the entire domain. - * - * @dealiiConceptRequires{concepts::is_triangulation_or_dof_handler} - */ - template - DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler) - std::list> get_finest_common_cells(const MeshType - &mesh_1, - const MeshType - &mesh_2); - - /** - * Return true if the two triangulations are based on the same coarse mesh. - * This is determined by checking whether they have the same number of cells - * on the coarsest level, and then checking that they have the same - * vertices. - * - * The two meshes may have different refinement histories beyond the coarse - * mesh. - */ - template - bool - have_same_coarse_mesh(const Triangulation &mesh_1, - const Triangulation &mesh_2); - - /** - * The same function as above, but working on arguments of type DoFHandler. - * This function is provided to allow calling have_same_coarse_mesh for all - * types of containers representing triangulations or the classes built on - * triangulations. - * - * @tparam MeshType A type that satisfies the requirements of the - * @ref ConceptMeshType "MeshType concept". - * - * @dealiiConceptRequires{concepts::is_triangulation_or_dof_handler} - */ - template - DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler) - bool have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2); - /** @} */ /** * @name Dealing with distorted cells @@ -3527,15 +3123,6 @@ namespace GridTools << "The given vertex with index " << arg1 << " is not used in the given triangulation."); - /** - * An exception that is thrown whenever the edges of a mesh are not - * orientable. - * - * @ingroup Exceptions - */ - DeclExceptionMsg(ExcMeshNotOrientable, - "The edges of the mesh are not consistently orientable."); - /** @} */ } /*namespace GridTools*/ diff --git a/include/deal.II/grid/grid_tools_topology.h b/include/deal.II/grid/grid_tools_topology.h new file mode 100644 index 0000000000..6075e074fa --- /dev/null +++ b/include/deal.II/grid/grid_tools_topology.h @@ -0,0 +1,460 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2001 - 2023 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii_grid_tools_topology_h +#define dealii_grid_tools_topology_h + +#include + +#include +#include + +#include +#include + +#include + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace GridTools +{ + /** + * @name Querying or modifying topological information + */ + /** @{ */ + /** + * Return the arrays that define the coarse mesh of a Triangulation. This + * function is the inverse of Triangulation::create_triangulation() in the + * sense that if one called this function on a triangulation, then that + * triangulation could be recreated by some kind of refinement from the + * results of this function. + * + * The return value is a tuple with the vector of vertices, the vector of + * cells, and a SubCellData structure. The latter contains additional + * information about faces and lines. These three objects are exactly + * the arguments to Triangulation::create_triangulation(). + * + * This function is useful in cases where one needs to deconstruct a + * Triangulation or manipulate the numbering of the vertices in some way: an + * example is GridGenerator::merge_triangulations(). + */ + template + std:: + tuple>, std::vector>, SubCellData> + get_coarse_mesh_description(const Triangulation &tria); + + /** + * Remove vertices that are not referenced by any of the cells. This + * function is called by all GridIn::read_* functions to eliminate + * vertices that are listed in the input files but are not used by the cells + * in the input file. While these vertices should not be in the input from + * the beginning, they sometimes are, most often when some cells have been + * removed by hand without wanting to update the vertex lists, as they might + * be lengthy. + * + * This function is called by all GridIn::read_* functions as the + * triangulation class requires them to be called with used vertices only. + * This is so, since the vertices are copied verbatim by that class, so we + * have to eliminate unused vertices beforehand. + * + * Not implemented for the codimension one case. + */ + template + void + delete_unused_vertices(std::vector> &vertices, + std::vector> &cells, + SubCellData &subcelldata); + + /** + * Remove vertices that are duplicated, due to the input of a structured + * grid, for example. If these vertices are not removed, the faces bounded + * by these vertices become part of the boundary, even if they are in the + * interior of the mesh. + * + * This function is called by some GridIn::read_* functions. Only + * the vertices with indices in @p considered_vertices are tested for + * equality. This speeds up the algorithm, which is, for worst-case hyper + * cube geometries $O(N^{3/2})$ in 2d and $O(N^{5/3})$ in 3d: quite slow. + * However, if you wish to consider all vertices, simply pass an empty + * vector. In that case, the function fills @p considered_vertices with all + * vertices. + * + * Two vertices are considered equal if their difference in each coordinate + * direction is less than @p tol. This implies that nothing happens if + * the tolerance is set to zero. + */ + template + void + delete_duplicated_vertices(std::vector> &all_vertices, + std::vector> &cells, + SubCellData &subcelldata, + std::vector &considered_vertices, + const double tol = 1e-12); + + /** + * Remove vertices that are duplicated. + * + * Two vertices are considered equal if their difference in each coordinate + * direction is less than @p tol. This implies that nothing happens if + * the tolerance is set to zero. + */ + template + void + delete_duplicated_vertices(std::vector> &vertices, + const double tol = 1e-12); + + /** + * Grids generated by grid generators may have an orientation of cells which + * is the inverse of the orientation required by deal.II. + * + * In 2d and 3d this function checks whether all cells have negative or + * positive measure/volume. In the former case, all cells are inverted. It + * does nothing in 1d. + * + * The inversion of cells might also work when only a subset of all cells + * have negative volume. However, grids consisting of a mixture of negative + * and positively oriented cells are very likely to be broken. Therefore, an + * exception is thrown, in case cells are not uniformly oriented. + * + * @note This function should be called before GridTools::consistently_order_cells(). + * + * @param all_vertices The vertices of the mesh. + * @param cells The array of CellData objects that describe the mesh's topology. + */ + template + void + invert_all_negative_measure_cells( + const std::vector> &all_vertices, + std::vector> &cells); + + /** + * Check the given cells and inverts any cell that is considered to have + * negative measure/volume in the orientation required by deal.II. + * + * This function is identical to invert_all_negative_measure_cells() except it + * does not throw an error if only some of the cells are inverted. Instead, + * this function returns how many cells were inverted. Additionally, it will + * always throw an exception outside of codimension 0. + */ + template + std::size_t + invert_cells_with_negative_measure( + const std::vector> &all_vertices, + std::vector> &cells); + + /** + * Given a vector of CellData objects describing a mesh, reorder their + * vertices so that all lines are consistently oriented. + * + * The expectations on orientation and a discussion of this function are + * available in the + * @ref reordering "reordering module". + * + * @param cells The array of CellData objects that describe the mesh's topology. + * @ingroup reordering + */ + template + void + consistently_order_cells(std::vector> &cells); + + /** + * Return a std::map with all vertices of faces located in the boundary + * + * @param[in] tria The Triangulation object. + */ + template + std::map> + get_all_vertices_at_boundary(const Triangulation &tria); + + /** + * Remove hanging nodes from a grid. If the @p isotropic parameter is set + * to @p false (default) this function detects cells with hanging nodes and + * refines the neighbours in the direction that removes hanging nodes. + * If the @p isotropic parameter is set + * to @p true, the neighbours refinement is made in each directions. + * In order to remove all hanging nodes this procedure has to be repeated: + * this could require a large number of iterations. + * To avoid this a max number (@p max_iterations) of iteration is provided. + * + * Consider the following grid: + * @image html remove_hanging_nodes-hanging.png + * + * @p isotropic == @p false would return: + * @image html remove_hanging_nodes-aniso.png + * + * @p isotropic == @p true would return: + * @image html remove_hanging_nodes-isotro.png + * + * @param[in,out] tria Triangulation to refine. + * + * @param[in] isotropic If true refine cells in each directions, otherwise + * (default value) refine the cell in the direction that removes hanging node. + * + * @param[in] max_iterations At each step only closest cells to hanging nodes + * are refined. The code may require a lot of iterations to remove all + * hanging nodes. @p max_iterations is the maximum number of iteration + * allowed. If @p max_iterations == numbers::invalid_unsigned_int this + * function continues refining until there are no hanging nodes. + * + * @note In the case of parallel codes, this function should be combined + * with GridGenerator::flatten_triangulation. + */ + template + void + remove_hanging_nodes(Triangulation &tria, + const bool isotropic = false, + const unsigned int max_iterations = 100); + + /** + * Refine a mesh anisotropically such that the resulting mesh is composed by + * cells with maximum ratio between dimensions less than @p max_ratio. + * This procedure requires an algorithm that may not terminate. Consequently, + * it is possible to set a maximum number of iterations through the + * @p max_iterations parameter. + * + * Starting from a cell like this: + * @image html remove_anisotropy-coarse.png + * + * This function would return: + * @image html remove_anisotropy-refined.png + * + * @param[in,out] tria Triangulation to refine. + * + * @param[in] max_ratio Maximum value allowed among the ratio between + * the dimensions of each cell. + * + * @param[in] max_iterations Maximum number of iterations allowed. + * + * @note In the case of parallel codes, this function should be combined + * with GridGenerator::flatten_triangulation and + * GridTools::remove_hanging_nodes. + */ + template + void + remove_anisotropy(Triangulation &tria, + const double max_ratio = 1.6180339887, + const unsigned int max_iterations = 5); + + /** + * Return a map `vertex index -> Point` containing the used + * vertices of the given `container`. The key of the returned map (i.e., + * the first element of the pair above) is the global index in the + * triangulation, whereas the value of each pair is the physical + * location of the corresponding vertex. The used vertices are obtained by + * looping over all cells, + * and querying for each cell where its vertices are through the (optional) + * `mapping` argument. + * + * In serial Triangulation objects and parallel::shared::Triangulation + * objects, the size of the returned map + * equals Triangulation::n_used_vertices() (not Triangulation::n_vertices()). + * Note that in parallel::distributed::Triangulation objects, only vertices in + * locally owned cells and ghost cells are returned, as for all other vertices + * their real location might not be known (e.g. for distributed computations + * using MappingQEulerian). + * + * If you use the default `mapping`, the returned map satisfies the following + * equality: + * + * @code + * const auto used_vertices = extract_used_vertices(tria); + * auto all_vertices = tria.get_vertices(); + * + * for(const auto &id_and_v : used_vertices) + * all_vertices[id_and_v.first] == id_and_v.second; // true + * @endcode + * + * Notice that the above is not satisfied for mappings that change the + * location of vertices, like MappingQEulerian. + * + * @ref ConceptMeshType "MeshType concept". + * @param container The container to extract vertices from. + * @param mapping The mapping to use to compute the points locations. + */ + template + std::map> + extract_used_vertices( + const Triangulation &container, + const Mapping &mapping = + (ReferenceCells::get_hypercube() +#ifndef _MSC_VER + .template get_default_linear_mapping() +#else + .ReferenceCell::get_default_linear_mapping() +#endif + )); + + /** + * Return the adjacent cells of all the vertices. If a vertex is also a + * hanging node, the associated coarse cell is also returned. The vertices + * are ordered by the vertex index. This is the number returned by the + * function cell-@>vertex_index(). Notice that only the indices + * marked in the array returned by + * Triangulation::get_used_vertices() are used. + */ + template + std::vector< + std::set::active_cell_iterator>> + vertex_to_cell_map(const Triangulation &triangulation); + + /** + * Produce a sparsity pattern in which nonzero entries indicate that two + * cells are connected via a common face. The diagonal entries of the + * sparsity pattern are also set. + * + * The rows and columns refer to the cells as they are traversed in their + * natural order using cell iterators. + */ + template + void + get_face_connectivity_of_cells( + const Triangulation &triangulation, + DynamicSparsityPattern &connectivity); + + /** + * Produce a sparsity pattern in which nonzero entries indicate that two + * cells are connected via a common vertex. The diagonal entries of the + * sparsity pattern are also set. + * + * The rows and columns refer to the cells as they are traversed in their + * natural order using cell iterators. + */ + template + void + get_vertex_connectivity_of_cells( + const Triangulation &triangulation, + DynamicSparsityPattern &connectivity); + + /** + * Produce a sparsity pattern for a given level mesh in which nonzero entries + * indicate that two cells are connected via a common vertex. The diagonal + * entries of the sparsity pattern are also set. + * + * The rows and columns refer to the cells as they are traversed in their + * natural order using cell iterators. + */ + template + void + get_vertex_connectivity_of_cells_on_level( + const Triangulation &triangulation, + const unsigned int level, + DynamicSparsityPattern &connectivity); + + /** @} */ + + /** + * @name Comparing different meshes + */ + /** @{ */ + + /** + * Given two meshes (i.e. objects of type Triangulation or DoFHandler) that + * are based on the same coarse mesh, this function figures out a set of cells + * that are matched between the two meshes and where at most one of the meshes + * is more refined on this cell. In other words, it finds the smallest cells + * that are common to both meshes, and that together completely cover the + * domain. + * + * This function is useful, for example, in time-dependent or nonlinear + * application, where one has to integrate a solution defined on one mesh + * (e.g., the one from the previous time step or nonlinear iteration) + * against the shape functions of another mesh (the next time step, the next + * nonlinear iteration). If, for example, the new mesh is finer, then one + * has to obtain the solution on the coarse mesh (mesh_1) and interpolate it + * to the children of the corresponding cell of mesh_2. Conversely, if the + * new mesh is coarser, one has to express the coarse cell shape function by + * a linear combination of fine cell shape functions. In either case, one + * needs to loop over the finest cells that are common to both + * triangulations. This function returns a list of pairs of matching + * iterators to cells in the two meshes that can be used to this end. + * + * Note that the list of these iterators is not necessarily ordered, and + * does also not necessarily coincide with the order in which cells are + * traversed in one, or both, of the meshes given as arguments. + * + * @tparam MeshType A type that satisfies the requirements of the + * @ref ConceptMeshType "MeshType concept". + * + * @note This function can only be used with + * parallel::distributed::Triangulation when both meshes use the same + * Triangulation since, with a distributed Triangulation, not all cells are + * stored locally, so the resulting list may not cover the entire domain. + * + * @dealiiConceptRequires{concepts::is_triangulation_or_dof_handler} + */ + template + DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler) + std::list> get_finest_common_cells(const MeshType + &mesh_1, + const MeshType + &mesh_2); + + /** + * Return true if the two triangulations are based on the same coarse mesh. + * This is determined by checking whether they have the same number of cells + * on the coarsest level, and then checking that they have the same + * vertices. + * + * The two meshes may have different refinement histories beyond the coarse + * mesh. + */ + template + bool + have_same_coarse_mesh(const Triangulation &mesh_1, + const Triangulation &mesh_2); + + /** + * The same function as above, but working on arguments of type DoFHandler. + * This function is provided to allow calling have_same_coarse_mesh for all + * types of containers representing triangulations or the classes built on + * triangulations. + * + * @tparam MeshType A type that satisfies the requirements of the + * @ref ConceptMeshType "MeshType concept". + * + * @dealiiConceptRequires{concepts::is_triangulation_or_dof_handler} + */ + template + DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler) + bool have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2); + + /** @} */ + + /** + * @name Exceptions + */ + /** @{ */ + + /** + * An exception that is thrown whenever the edges of a mesh are not + * orientable. + * + * @ingroup Exceptions + */ + DeclExceptionMsg(ExcMeshNotOrientable, + "The edges of the mesh are not consistently orientable."); + + /** @} */ +} // namespace GridTools + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/grid/CMakeLists.txt b/source/grid/CMakeLists.txt index 4c31e658aa..fa7745790c 100644 --- a/source/grid/CMakeLists.txt +++ b/source/grid/CMakeLists.txt @@ -57,6 +57,7 @@ set(_separate_src grid_tools_dof_handlers.cc grid_tools_geometry.cc grid_tools_nontemplates.cc + grid_tools_topology.cc tria.cc ) @@ -82,6 +83,7 @@ set(_inst grid_tools_cache.inst.in grid_tools_dof_handlers.inst.in grid_tools_geometry.inst.in + grid_tools_topology.inst.in intergrid_map.inst.in manifold.inst.in manifold_lib.inst.in diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index e21c012af9..43b778a793 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -106,1615 +106,6 @@ DEAL_II_NAMESPACE_OPEN namespace GridTools { - // Generic functions for appending face data in 2d or 3d. TODO: we can - // remove these once we have 'if constexpr'. - namespace internal - { - inline void - append_face_data(const CellData<1> &face_data, SubCellData &subcell_data) - { - subcell_data.boundary_lines.push_back(face_data); - } - - - - inline void - append_face_data(const CellData<2> &face_data, SubCellData &subcell_data) - { - subcell_data.boundary_quads.push_back(face_data); - } - - - - // Lexical comparison for sorting CellData objects. - template - struct CellDataComparator - { - bool - operator()(const CellData &a, - const CellData &b) const - { - // Check vertices: - if (std::lexicographical_compare(std::begin(a.vertices), - std::end(a.vertices), - std::begin(b.vertices), - std::end(b.vertices))) - return true; - // it should never be necessary to check the material or manifold - // ids as a 'tiebreaker' (since they must be equal if the vertex - // indices are equal). Assert it anyway: -#ifdef DEBUG - if (std::equal(std::begin(a.vertices), - std::end(a.vertices), - std::begin(b.vertices))) - { - Assert(a.material_id == b.material_id && - a.manifold_id == b.manifold_id, - ExcMessage( - "Two CellData objects with equal vertices must " - "have the same material/boundary ids and manifold " - "ids.")); - } -#endif - return false; - } - }; - - - /** - * get_coarse_mesh_description() needs to store face data for dim>1, but - * we can not have this code in the function, as this requires either an - * instantiation of CellData, or constexpr if. We use a class with - * specialization instead for now. - * - * Data on faces is added with insert_face_data() and then retrieved with - * get(). - */ - template - class FaceDataHelper - { - public: - /** - * Store the data about the given face @p face. - */ - template - void - insert_face_data(const FaceIteratorType &face) - { - CellData face_cell_data(face->n_vertices()); - for (unsigned int vertex_n = 0; vertex_n < face->n_vertices(); - ++vertex_n) - face_cell_data.vertices[vertex_n] = face->vertex_index(vertex_n); - face_cell_data.boundary_id = face->boundary_id(); - face_cell_data.manifold_id = face->manifold_id(); - - face_data.insert(std::move(face_cell_data)); - } - - /** - * Return the @p subcell_data with the stored information. - */ - SubCellData - get() - { - SubCellData subcell_data; - - for (const CellData &face_cell_data : face_data) - internal::append_face_data(face_cell_data, subcell_data); - return subcell_data; - } - - - private: - std::set, internal::CellDataComparator> - face_data; - }; - - - // Do nothing for dim=1: - template <> - class FaceDataHelper<1> - { - public: - template - void - insert_face_data(const FaceIteratorType &) - {} - - SubCellData - get() - { - return SubCellData(); - } - }; - } // namespace internal - - - - template - std:: - tuple>, std::vector>, SubCellData> - get_coarse_mesh_description(const Triangulation &tria) - { - Assert(tria.n_levels() >= 1, - ExcMessage("The input triangulation must be non-empty.")); - - std::vector> vertices = tria.get_vertices(); - std::vector> cells; - - internal::FaceDataHelper face_data; - std::set, internal::CellDataComparator<1>> - line_data; // only used in 3d - - for (const auto &cell : tria.cell_iterators_on_level(0)) - { - // Save cell data - CellData cell_data(cell->n_vertices()); - for (const unsigned int cell_vertex_n : cell->vertex_indices()) - { - Assert(cell->vertex_index(cell_vertex_n) < vertices.size(), - ExcInternalError()); - cell_data.vertices[cell_vertex_n] = - cell->vertex_index(cell_vertex_n); - } - cell_data.material_id = cell->material_id(); - cell_data.manifold_id = cell->manifold_id(); - cells.emplace_back(std::move(cell_data)); - - // Save face data - if (dim > 1) - { - for (const unsigned int face_n : cell->face_indices()) - // We don't need to insert anything if we have default values - { - const auto face = cell->face(face_n); - if (face->boundary_id() != numbers::internal_face_boundary_id || - face->manifold_id() != numbers::flat_manifold_id) - face_data.insert_face_data(face); - } - } - // Save line data - if (dim == 3) - { - for (unsigned int line_n = 0; line_n < cell->n_lines(); ++line_n) - { - const auto line = cell->line(line_n); - // We don't need to insert anything if we have default values - if (line->boundary_id() != numbers::internal_face_boundary_id || - line->manifold_id() != numbers::flat_manifold_id) - { - CellData<1> line_cell_data(line->n_vertices()); - for (const unsigned int vertex_n : line->vertex_indices()) - line_cell_data.vertices[vertex_n] = - line->vertex_index(vertex_n); - line_cell_data.boundary_id = line->boundary_id(); - line_cell_data.manifold_id = line->manifold_id(); - line_data.insert(std::move(line_cell_data)); - } - } - } - } - - SubCellData subcell_data = face_data.get(); - - if (dim == 3) - for (const CellData<1> &face_line_data : line_data) - subcell_data.boundary_lines.push_back(face_line_data); - - // We end up with a 'vertices' array that uses some of the entries, - // but not all -- specifically, all vertices referenced by level-0 - // cells. We can compress the array: - GridTools::delete_unused_vertices(vertices, cells, subcell_data); - - return std::tuple>, - std::vector>, - SubCellData>(std::move(vertices), - std::move(cells), - std::move(subcell_data)); - } - - - - template - void - delete_unused_vertices(std::vector> &vertices, - std::vector> &cells, - SubCellData &subcelldata) - { - Assert( - subcelldata.check_consistency(dim), - ExcMessage( - "Invalid SubCellData supplied according to ::check_consistency(). " - "This is caused by data containing objects for the wrong dimension.")); - - // first check which vertices are actually used - std::vector vertex_used(vertices.size(), false); - for (unsigned int c = 0; c < cells.size(); ++c) - for (unsigned int v = 0; v < cells[c].vertices.size(); ++v) - { - Assert(cells[c].vertices[v] < vertices.size(), - ExcMessage("Invalid vertex index encountered! cells[" + - Utilities::int_to_string(c) + "].vertices[" + - Utilities::int_to_string(v) + "]=" + - Utilities::int_to_string(cells[c].vertices[v]) + - " is invalid, because only " + - Utilities::int_to_string(vertices.size()) + - " vertices were supplied.")); - vertex_used[cells[c].vertices[v]] = true; - } - - - // then renumber the vertices that are actually used in the same order as - // they were beforehand - const unsigned int invalid_vertex = numbers::invalid_unsigned_int; - std::vector new_vertex_numbers(vertices.size(), - invalid_vertex); - unsigned int next_free_number = 0; - for (unsigned int i = 0; i < vertices.size(); ++i) - if (vertex_used[i] == true) - { - new_vertex_numbers[i] = next_free_number; - ++next_free_number; - } - - // next replace old vertex numbers by the new ones - for (unsigned int c = 0; c < cells.size(); ++c) - for (auto &v : cells[c].vertices) - v = new_vertex_numbers[v]; - - // same for boundary data - for (unsigned int c = 0; c < subcelldata.boundary_lines.size(); // NOLINT - ++c) - for (unsigned int v = 0; - v < subcelldata.boundary_lines[c].vertices.size(); - ++v) - { - Assert(subcelldata.boundary_lines[c].vertices[v] < - new_vertex_numbers.size(), - ExcMessage( - "Invalid vertex index in subcelldata.boundary_lines. " - "subcelldata.boundary_lines[" + - Utilities::int_to_string(c) + "].vertices[" + - Utilities::int_to_string(v) + "]=" + - Utilities::int_to_string( - subcelldata.boundary_lines[c].vertices[v]) + - " is invalid, because only " + - Utilities::int_to_string(vertices.size()) + - " vertices were supplied.")); - subcelldata.boundary_lines[c].vertices[v] = - new_vertex_numbers[subcelldata.boundary_lines[c].vertices[v]]; - } - - for (unsigned int c = 0; c < subcelldata.boundary_quads.size(); // NOLINT - ++c) - for (unsigned int v = 0; - v < subcelldata.boundary_quads[c].vertices.size(); - ++v) - { - Assert(subcelldata.boundary_quads[c].vertices[v] < - new_vertex_numbers.size(), - ExcMessage( - "Invalid vertex index in subcelldata.boundary_quads. " - "subcelldata.boundary_quads[" + - Utilities::int_to_string(c) + "].vertices[" + - Utilities::int_to_string(v) + "]=" + - Utilities::int_to_string( - subcelldata.boundary_quads[c].vertices[v]) + - " is invalid, because only " + - Utilities::int_to_string(vertices.size()) + - " vertices were supplied.")); - - subcelldata.boundary_quads[c].vertices[v] = - new_vertex_numbers[subcelldata.boundary_quads[c].vertices[v]]; - } - - // finally copy over the vertices which we really need to a new array and - // replace the old one by the new one - std::vector> tmp; - tmp.reserve(std::count(vertex_used.begin(), vertex_used.end(), true)); - for (unsigned int v = 0; v < vertices.size(); ++v) - if (vertex_used[v] == true) - tmp.push_back(vertices[v]); - swap(vertices, tmp); - } - - - - template - void - delete_duplicated_vertices(std::vector> &vertices, - std::vector> &cells, - SubCellData &subcelldata, - std::vector &considered_vertices, - const double tol) - { - if (tol == 0.0) - return; // nothing to do per definition - - AssertIndexRange(2, vertices.size()); - std::vector new_vertex_numbers(vertices.size()); - std::iota(new_vertex_numbers.begin(), new_vertex_numbers.end(), 0); - - // if the considered_vertices vector is empty, consider all vertices - if (considered_vertices.empty()) - considered_vertices = new_vertex_numbers; - Assert(considered_vertices.size() <= vertices.size(), ExcInternalError()); - - // The algorithm below improves upon the naive O(n^2) algorithm by first - // sorting vertices by their value in one component and then only - // comparing vertices for equality which are nearly equal in that - // component. For example, if @p vertices form a cube, then we will only - // compare points that have the same x coordinate when we try to find - // duplicated vertices. - - // Start by finding the longest coordinate direction. This minimizes the - // number of points that need to be compared against each-other in a - // single set for typical geometries. - const BoundingBox bbox(vertices); - - unsigned int longest_coordinate_direction = 0; - double longest_coordinate_length = bbox.side_length(0); - for (unsigned int d = 1; d < spacedim; ++d) - { - const double coordinate_length = bbox.side_length(d); - if (longest_coordinate_length < coordinate_length) - { - longest_coordinate_length = coordinate_length; - longest_coordinate_direction = d; - } - } - - // Sort vertices (while preserving their vertex numbers) along that - // coordinate direction: - std::vector>> sorted_vertices; - sorted_vertices.reserve(vertices.size()); - for (const unsigned int vertex_n : considered_vertices) - { - AssertIndexRange(vertex_n, vertices.size()); - sorted_vertices.emplace_back(vertex_n, vertices[vertex_n]); - } - std::sort(sorted_vertices.begin(), - sorted_vertices.end(), - [&](const std::pair> &a, - const std::pair> &b) { - return a.second[longest_coordinate_direction] < - b.second[longest_coordinate_direction]; - }); - - auto within_tolerance = [=](const Point &a, - const Point &b) { - for (unsigned int d = 0; d < spacedim; ++d) - if (std::abs(a[d] - b[d]) > tol) - return false; - return true; - }; - - // Find a range of numbers that have the same component in the longest - // coordinate direction: - auto range_start = sorted_vertices.begin(); - while (range_start != sorted_vertices.end()) - { - auto range_end = range_start + 1; - while (range_end != sorted_vertices.end() && - std::abs(range_end->second[longest_coordinate_direction] - - range_start->second[longest_coordinate_direction]) < - tol) - ++range_end; - - // preserve behavior with older versions of this function by replacing - // higher vertex numbers by lower vertex numbers - std::sort(range_start, - range_end, - [](const std::pair> &a, - const std::pair> &b) { - return a.first < b.first; - }); - - // Now de-duplicate [range_start, range_end) - // - // We have identified all points that are within a strip of width 'tol' - // in one coordinate direction. Now we need to figure out which of these - // are also close in other coordinate directions. If two are close, we - // can mark the second one for deletion. - for (auto reference = range_start; reference != range_end; ++reference) - { - if (reference->first != numbers::invalid_unsigned_int) - for (auto it = reference + 1; it != range_end; ++it) - { - if (within_tolerance(reference->second, it->second)) - { - new_vertex_numbers[it->first] = reference->first; - // skip the replaced vertex in the future - it->first = numbers::invalid_unsigned_int; - } - } - } - range_start = range_end; - } - - // now we got a renumbering list. simply renumber all vertices - // (non-duplicate vertices get renumbered to themselves, so nothing bad - // happens). after that, the duplicate vertices will be unused, so call - // delete_unused_vertices() to do that part of the job. - for (auto &cell : cells) - for (auto &vertex_index : cell.vertices) - vertex_index = new_vertex_numbers[vertex_index]; - for (auto &quad : subcelldata.boundary_quads) - for (auto &vertex_index : quad.vertices) - vertex_index = new_vertex_numbers[vertex_index]; - for (auto &line : subcelldata.boundary_lines) - for (auto &vertex_index : line.vertices) - vertex_index = new_vertex_numbers[vertex_index]; - - delete_unused_vertices(vertices, cells, subcelldata); - } - - - - template - void - delete_duplicated_vertices(std::vector> &vertices, - const double tol) - { - if (vertices.empty()) - return; - - // 1) map point to local vertex index - std::map, unsigned int, FloatingPointComparator> - map_point_to_local_vertex_index{FloatingPointComparator(tol)}; - - // 2) initialize map with existing points uniquely - for (unsigned int i = 0; i < vertices.size(); ++i) - map_point_to_local_vertex_index[vertices[i]] = i; - - // no duplicate points are found - if (map_point_to_local_vertex_index.size() == vertices.size()) - return; - - // 3) remove duplicate entries from vertices - vertices.resize(map_point_to_local_vertex_index.size()); - { - unsigned int j = 0; - for (const auto &p : map_point_to_local_vertex_index) - vertices[j++] = p.first; - } - } - - - - template - std::size_t - invert_cells_with_negative_measure( - const std::vector> &all_vertices, - std::vector> &cells) - { - // This function is presently only implemented for volumetric (codimension - // 0) elements. - - if (dim == 1) - return 0; - if (dim == 2 && spacedim == 3) - Assert(false, ExcNotImplemented()); - - std::size_t n_negative_cells = 0; - std::size_t cell_no = 0; - for (auto &cell : cells) - { - const ArrayView vertices(cell.vertices); - // Some pathologically twisted cells can have exactly zero measure but - // we can still fix them - if (GridTools::cell_measure(all_vertices, vertices) <= 0) - { - ++n_negative_cells; - const auto reference_cell = - ReferenceCell::n_vertices_to_type(dim, vertices.size()); - - if (reference_cell.is_hyper_cube()) - { - if (dim == 2) - { - // flip the cell across the y = x line in 2d - std::swap(cell.vertices[1], cell.vertices[2]); - } - else if (dim == 3) - { - // swap the front and back faces in 3d - std::swap(cell.vertices[0], cell.vertices[2]); - std::swap(cell.vertices[1], cell.vertices[3]); - std::swap(cell.vertices[4], cell.vertices[6]); - std::swap(cell.vertices[5], cell.vertices[7]); - } - } - else if (reference_cell.is_simplex()) - { - // By basic rules for computing determinants we can just swap - // two vertices to fix a negative volume. Arbitrarily pick the - // last two. - std::swap(cell.vertices[cell.vertices.size() - 2], - cell.vertices[cell.vertices.size() - 1]); - } - else if (reference_cell == ReferenceCells::Wedge) - { - // swap the two triangular faces - std::swap(cell.vertices[0], cell.vertices[3]); - std::swap(cell.vertices[1], cell.vertices[4]); - std::swap(cell.vertices[2], cell.vertices[5]); - } - else if (reference_cell == ReferenceCells::Pyramid) - { - // Try swapping two vertices in the base - perhaps things were - // read in the UCD (counter-clockwise) order instead of lexical - std::swap(cell.vertices[2], cell.vertices[3]); - } - else - { - AssertThrow(false, ExcNotImplemented()); - } - // Check whether the resulting cell is now ok. - // If not, then the grid is seriously broken and - // we just give up. - AssertThrow(GridTools::cell_measure(all_vertices, vertices) > 0, - ExcGridHasInvalidCell(cell_no)); - } - ++cell_no; - } - return n_negative_cells; - } - - - template - void - invert_all_negative_measure_cells( - const std::vector> &all_vertices, - std::vector> &cells) - { - const std::size_t n_negative_cells = - invert_cells_with_negative_measure(all_vertices, cells); - - // We assume that all cells of a grid have - // either positive or negative volumes but - // not both mixed. Although above reordering - // might work also on single cells, grids - // with both kind of cells are very likely to - // be broken. Check for this here. - AssertThrow(n_negative_cells == 0 || n_negative_cells == cells.size(), - ExcMessage( - std::string( - "This function assumes that either all cells have positive " - "volume, or that all cells have been specified in an " - "inverted vertex order so that their volume is negative. " - "(In the latter case, this class automatically inverts " - "every cell.) However, the mesh you have specified " - "appears to have both cells with positive and cells with " - "negative volume. You need to check your mesh which " - "cells these are and how they got there.\n" - "As a hint, of the total ") + - std::to_string(cells.size()) + " cells in the mesh, " + - std::to_string(n_negative_cells) + - " appear to have a negative volume.")); - } - - - - // Functions and classes for consistently_order_cells - namespace - { - /** - * A simple data structure denoting an edge, i.e., the ordered pair - * of its vertex indices. This is only used in the is_consistent() - * function. - */ - struct CheapEdge - { - /** - * Construct an edge from the global indices of its two vertices. - */ - CheapEdge(const unsigned int v0, const unsigned int v1) - : v0(v0) - , v1(v1) - {} - - /** - * Comparison operator for edges. It compares based on the - * lexicographic ordering of the two vertex indices. - */ - bool - operator<(const CheapEdge &e) const - { - return ((v0 < e.v0) || ((v0 == e.v0) && (v1 < e.v1))); - } - - private: - /** - * The global indices of the vertices that define the edge. - */ - const unsigned int v0, v1; - }; - - - /** - * A function that determines whether the edges in a mesh are - * already consistently oriented. It does so by adding all edges - * of all cells into a set (which automatically eliminates - * duplicates) but before that checks whether the reverse edge is - * already in the set -- which would imply that a neighboring cell - * is inconsistently oriented. - */ - template - bool - is_consistent(const std::vector> &cells) - { - std::set edges; - - for (typename std::vector>::const_iterator c = - cells.begin(); - c != cells.end(); - ++c) - { - // construct the edges in reverse order. for each of them, - // ensure that the reverse edge is not yet in the list of - // edges (return false if the reverse edge already *is* in - // the list) and then add the actual edge to it; std::set - // eliminates duplicates automatically - for (unsigned int l = 0; l < GeometryInfo::lines_per_cell; ++l) - { - const CheapEdge reverse_edge( - c->vertices[GeometryInfo::line_to_cell_vertices(l, 1)], - c->vertices[GeometryInfo::line_to_cell_vertices(l, 0)]); - if (edges.find(reverse_edge) != edges.end()) - return false; - - - // ok, not. insert edge in correct order - const CheapEdge correct_edge( - c->vertices[GeometryInfo::line_to_cell_vertices(l, 0)], - c->vertices[GeometryInfo::line_to_cell_vertices(l, 1)]); - edges.insert(correct_edge); - } - } - - // no conflicts found, so return true - return true; - } - - - /** - * A structure that describes some properties of parallel edges - * such as what starter edges are (i.e., representative elements - * of the sets of parallel edges within a cell) and what the set - * of parallel edges to each edge is. - */ - template - struct ParallelEdges - { - /** - * An array that contains the indices of dim edges that can - * serve as (arbitrarily chosen) starting points for the - * dim sets of parallel edges within each cell. - */ - static const unsigned int starter_edges[dim]; - - /** - * Number and indices of all of those edges parallel to each of the - * edges in a cell. - */ - static const unsigned int n_other_parallel_edges = (1 << (dim - 1)) - 1; - static const unsigned int - parallel_edges[GeometryInfo::lines_per_cell] - [n_other_parallel_edges]; - }; - - template <> - const unsigned int ParallelEdges<2>::starter_edges[2] = {0, 2}; - - template <> - const unsigned int ParallelEdges<2>::parallel_edges[4][1] = {{1}, - {0}, - {3}, - {2}}; - - template <> - const unsigned int ParallelEdges<3>::starter_edges[3] = {0, 2, 8}; - - template <> - const unsigned int ParallelEdges<3>::parallel_edges[12][3] = { - {1, 4, 5}, // line 0 - {0, 4, 5}, // line 1 - {3, 6, 7}, // line 2 - {2, 6, 7}, // line 3 - {0, 1, 5}, // line 4 - {0, 1, 4}, // line 5 - {2, 3, 7}, // line 6 - {2, 3, 6}, // line 7 - {9, 10, 11}, // line 8 - {8, 10, 11}, // line 9 - {8, 9, 11}, // line 10 - {8, 9, 10} // line 11 - }; - - - /** - * A structure that store the index of a cell and, crucially, how a - * given edge relates to this cell. - */ - struct AdjacentCell - { - /** - * Default constructor. Initialize the fields with invalid values. - */ - AdjacentCell() - : cell_index(numbers::invalid_unsigned_int) - , edge_within_cell(numbers::invalid_unsigned_int) - {} - - /** - * Constructor. Initialize the fields with the given values. - */ - AdjacentCell(const unsigned int cell_index, - const unsigned int edge_within_cell) - : cell_index(cell_index) - , edge_within_cell(edge_within_cell) - {} - - - unsigned int cell_index; - unsigned int edge_within_cell; - }; - - - - template - class AdjacentCells; - - /** - * A class that represents all of the cells adjacent to a given edge. - * This class corresponds to the 2d case where each edge has at most - * two adjacent cells. - */ - template <> - class AdjacentCells<2> - { - public: - /** - * An iterator that allows iterating over all cells adjacent - * to the edge represented by the current object. - */ - using const_iterator = const AdjacentCell *; - - /** - * Add the given cell to the collection of cells adjacent to - * the edge this object corresponds to. Since we are covering - * the 2d case, the set of adjacent cells currently - * represented by this object must have either zero or - * one element already, since we can not add more than two - * adjacent cells for each edge. - */ - void - push_back(const AdjacentCell &adjacent_cell) - { - if (adjacent_cells[0].cell_index == numbers::invalid_unsigned_int) - adjacent_cells[0] = adjacent_cell; - else - { - Assert(adjacent_cells[1].cell_index == - numbers::invalid_unsigned_int, - ExcInternalError()); - adjacent_cells[1] = adjacent_cell; - } - } - - - /** - * Return an iterator to the first valid cell stored as adjacent to the - * edge represented by the current object. - */ - const_iterator - begin() const - { - return adjacent_cells; - } - - - /** - * Return an iterator to the element past the last valid cell stored - * as adjacent to the edge represented by the current object. - * @return - */ - const_iterator - end() const - { - // check whether the current object stores zero, one, or two - // adjacent cells, and use this to point to the element past the - // last valid one - if (adjacent_cells[0].cell_index == numbers::invalid_unsigned_int) - return adjacent_cells; - else if (adjacent_cells[1].cell_index == numbers::invalid_unsigned_int) - return adjacent_cells + 1; - else - return adjacent_cells + 2; - } - - private: - /** - * References to the (at most) two cells that are adjacent to - * the edge this object corresponds to. Unused elements are - * default-initialized and have invalid values; in particular, - * their cell_index field equals numbers::invalid_unsigned_int. - */ - AdjacentCell adjacent_cells[2]; - }; - - - - /** - * A class that represents all of the cells adjacent to a given edge. - * This class corresponds to the 3d case where each edge can have an - * arbitrary number of adjacent cells. We represent this as a - * std::vector, from which class the current one is - * derived and from which it inherits all of its member functions. - */ - template <> - class AdjacentCells<3> : public std::vector - {}; - - - /** - * A class that describes all of the relevant properties of an - * edge. For the purpose of what we do here, that includes the - * indices of the two vertices, and the indices of the adjacent - * cells (together with a description *where* in each of the - * adjacent cells the edge is located). It also includes the - * (global) direction of the edge: either from the first vertex to - * the second, the other way around, or so far undetermined. - */ - template - class Edge - { - public: - /** - * Constructor. Create the edge based on the information given - * in @p cell, and selecting the edge with number @p edge_number - * within this cell. Initialize the edge as unoriented. - */ - Edge(const CellData &cell, const unsigned int edge_number) - : orientation_status(not_oriented) - { - Assert(edge_number < GeometryInfo::lines_per_cell, - ExcInternalError()); - - // copy vertices for this particular line - vertex_indices[0] = - cell - .vertices[GeometryInfo::line_to_cell_vertices(edge_number, 0)]; - vertex_indices[1] = - cell - .vertices[GeometryInfo::line_to_cell_vertices(edge_number, 1)]; - - // bring them into standard orientation - if (vertex_indices[0] > vertex_indices[1]) - std::swap(vertex_indices[0], vertex_indices[1]); - } - - /** - * Comparison operator for edges. It compares based on the - * lexicographic ordering of the two vertex indices. - */ - bool - operator<(const Edge &e) const - { - return ((vertex_indices[0] < e.vertex_indices[0]) || - ((vertex_indices[0] == e.vertex_indices[0]) && - (vertex_indices[1] < e.vertex_indices[1]))); - } - - /** - * Compare two edges for equality based on their vertex indices. - */ - bool - operator==(const Edge &e) const - { - return ((vertex_indices[0] == e.vertex_indices[0]) && - (vertex_indices[1] == e.vertex_indices[1])); - } - - /** - * The global indices of the two vertices that bound this edge. These - * will be ordered so that the first index is less than the second. - */ - unsigned int vertex_indices[2]; - - /** - * An enum that indicates the direction of this edge with - * regard to the two vertices that bound it. - */ - enum OrientationStatus - { - not_oriented, - forward, - backward - }; - - OrientationStatus orientation_status; - - /** - * Store the set of cells adjacent to this edge (these cells then - * also store *where* in the cell the edge is located). - */ - AdjacentCells adjacent_cells; - }; - - - - /** - * A data structure that represents a cell with all of its vertices - * and edges. - */ - template - struct Cell - { - /** - * Construct a Cell object from a CellData object. Also take a - * (sorted) list of edges and to point the edges of the current - * object into this list of edges. - */ - Cell(const CellData &c, const std::vector> &edge_list) - { - for (const unsigned int i : GeometryInfo::vertex_indices()) - vertex_indices[i] = c.vertices[i]; - - // now for each of the edges of this cell, find the location inside the - // given edge_list array and store than index - for (unsigned int l = 0; l < GeometryInfo::lines_per_cell; ++l) - { - const Edge e(c, l); - edge_indices[l] = - (std::lower_bound(edge_list.begin(), edge_list.end(), e) - - edge_list.begin()); - Assert(edge_indices[l] < edge_list.size(), ExcInternalError()); - Assert(edge_list[edge_indices[l]] == e, ExcInternalError()); - } - } - - /** - * A list of global indices for the vertices that bound this cell. - */ - unsigned int vertex_indices[GeometryInfo::vertices_per_cell]; - - /** - * A list of indices into the 'edge_list' array passed to the constructor - * for the edges of the current cell. - */ - unsigned int edge_indices[GeometryInfo::lines_per_cell]; - }; - - - - template - class EdgeDeltaSet; - - /** - * A class that represents by how much the set of parallel edges - * grows in each step. In the graph orientation paper, this set is - * called $\Delta_k$, thus the name. - * - * In 2d, this set can only include zero, one, or two elements. - * Consequently, the appropriate data structure is one in which - * we store at most 2 elements in a fixed sized data structure. - */ - template <> - class EdgeDeltaSet<2> - { - public: - /** - * Iterator type for the elements of the set. - */ - using const_iterator = const unsigned int *; - - /** - * Default constructor. Initialize both slots as unused, corresponding - * to an empty set. - */ - EdgeDeltaSet() - { - edge_indices[0] = edge_indices[1] = numbers::invalid_unsigned_int; - } - - - /** - * Delete the elements of the set by marking both slots as unused. - */ - void - clear() - { - edge_indices[0] = edge_indices[1] = numbers::invalid_unsigned_int; - } - - /** - * Insert one element into the set. This will fail if the set already - * has two elements. - */ - void - insert(const unsigned int edge_index) - { - if (edge_indices[0] == numbers::invalid_unsigned_int) - edge_indices[0] = edge_index; - else - { - Assert(edge_indices[1] == numbers::invalid_unsigned_int, - ExcInternalError()); - edge_indices[1] = edge_index; - } - } - - - /** - * Return an iterator pointing to the first element of the set. - */ - const_iterator - begin() const - { - return edge_indices; - } - - - /** - * Return an iterator pointing to the element past the last used one. - */ - const_iterator - end() const - { - // check whether the current object stores zero, one, or two - // indices, and use this to point to the element past the - // last valid one - if (edge_indices[0] == numbers::invalid_unsigned_int) - return edge_indices; - else if (edge_indices[1] == numbers::invalid_unsigned_int) - return edge_indices + 1; - else - return edge_indices + 2; - } - - private: - /** - * Storage space to store the indices of at most two edges. - */ - unsigned int edge_indices[2]; - }; - - - - /** - * A class that represents by how much the set of parallel edges - * grows in each step. In the graph orientation paper, this set is - * called $\Delta_k$, thus the name. - * - * In 3d, this set can have arbitrarily many elements, unlike the - * 2d case specialized above. Consequently, we simply represent - * the data structure with a std::set. Class derivation ensures - * that we simply inherit all of the member functions of the - * base class. - */ - template <> - class EdgeDeltaSet<3> : public std::set - {}; - - - - /** - * From a list of cells, build a sorted vector that contains all of the - * edges that exist in the mesh. - */ - template - std::vector> - build_edges(const std::vector> &cells) - { - // build the edge list for all cells. because each cell has - // GeometryInfo::lines_per_cell edges, the total number - // of edges is this many times the number of cells. of course - // some of them will be duplicates, and we throw them out below - std::vector> edge_list; - edge_list.reserve(cells.size() * GeometryInfo::lines_per_cell); - for (unsigned int i = 0; i < cells.size(); ++i) - for (unsigned int l = 0; l < GeometryInfo::lines_per_cell; ++l) - edge_list.emplace_back(cells[i], l); - - // next sort the edge list and then remove duplicates - std::sort(edge_list.begin(), edge_list.end()); - edge_list.erase(std::unique(edge_list.begin(), edge_list.end()), - edge_list.end()); - - return edge_list; - } - - - - /** - * Build the cell list. Update the edge array to let edges know - * which cells are adjacent to them. - */ - template - std::vector> - build_cells_and_connect_edges(const std::vector> &cells, - std::vector> &edges) - { - std::vector> cell_list; - cell_list.reserve(cells.size()); - for (unsigned int i = 0; i < cells.size(); ++i) - { - // create our own data structure for the cells and let it - // connect to the edges array - cell_list.emplace_back(cells[i], edges); - - // then also inform the edges that they are adjacent - // to the current cell, and where within this cell - for (unsigned int l = 0; l < GeometryInfo::lines_per_cell; ++l) - edges[cell_list.back().edge_indices[l]].adjacent_cells.push_back( - AdjacentCell(i, l)); - } - Assert(cell_list.size() == cells.size(), ExcInternalError()); - - return cell_list; - } - - - - /** - * Return the index within 'cells' of the first cell that has at least one - * edge that is not yet oriented. - */ - template - unsigned int - get_next_unoriented_cell(const std::vector> &cells, - const std::vector> &edges, - const unsigned int current_cell) - { - for (unsigned int c = current_cell; c < cells.size(); ++c) - for (unsigned int l = 0; l < GeometryInfo::lines_per_cell; ++l) - if (edges[cells[c].edge_indices[l]].orientation_status == - Edge::not_oriented) - return c; - - return numbers::invalid_unsigned_int; - } - - - - /** - * Given a set of cells and edges, orient all edges that are - * (global) parallel to the one identified by the @p cell and - * within it the one with index @p local_edge. - */ - template - void - orient_one_set_of_parallel_edges(const std::vector> &cells, - std::vector> &edges, - const unsigned int cell, - const unsigned int local_edge) - { - // choose the direction of the first edge. we have free choice - // here and could simply choose "forward" if that's what pleases - // us. however, for backward compatibility with the previous - // implementation used till 2016, let us just choose the - // direction so that it matches what we have in the given cell. - // - // in fact, in what can only be assumed to be a bug in the - // original implementation, after orienting all edges, the code - // that rotates the cells so that they match edge orientations - // (see the rotate_cell() function below) rotated the cell two - // more times by 90 degrees. this is ok -- it simply flips all - // edge orientations, which leaves them valid. rather than do - // the same in the current implementation, we can achieve the - // same effect by modifying the rule above to choose the - // direction of the starting edge of this parallel set - // *opposite* to what it looks like in the current cell - // - // this bug only existed in the 2d implementation since there - // were different implementations for 2d and 3d. consequently, - // only replicate it for the 2d case and be "intuitive" in 3d. - if (edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] == - cells[cell].vertex_indices[GeometryInfo::line_to_cell_vertices( - local_edge, 0)]) - // orient initial edge *opposite* to the way it is in the cell - // (see above for the reason) - edges[cells[cell].edge_indices[local_edge]].orientation_status = - (dim == 2 ? Edge::backward : Edge::forward); - else - { - Assert( - edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] == - cells[cell].vertex_indices - [GeometryInfo::line_to_cell_vertices(local_edge, 1)], - ExcInternalError()); - Assert( - edges[cells[cell].edge_indices[local_edge]].vertex_indices[1] == - cells[cell].vertex_indices - [GeometryInfo::line_to_cell_vertices(local_edge, 0)], - ExcInternalError()); - - // orient initial edge *opposite* to the way it is in the cell - // (see above for the reason) - edges[cells[cell].edge_indices[local_edge]].orientation_status = - (dim == 2 ? Edge::forward : Edge::backward); - } - - // walk outward from the given edge as described in - // the algorithm in the paper that documents all of - // this - // - // note that in 2d, each of the Deltas can at most - // contain two elements, whereas in 3d it can be arbitrarily many - EdgeDeltaSet Delta_k; - EdgeDeltaSet Delta_k_minus_1; - Delta_k_minus_1.insert(cells[cell].edge_indices[local_edge]); - - while (Delta_k_minus_1.begin() != - Delta_k_minus_1.end()) // while set is not empty - { - Delta_k.clear(); - - for (typename EdgeDeltaSet::const_iterator delta = - Delta_k_minus_1.begin(); - delta != Delta_k_minus_1.end(); - ++delta) - { - Assert(edges[*delta].orientation_status != - Edge::not_oriented, - ExcInternalError()); - - // now go through the cells adjacent to this edge - for (typename AdjacentCells::const_iterator adjacent_cell = - edges[*delta].adjacent_cells.begin(); - adjacent_cell != edges[*delta].adjacent_cells.end(); - ++adjacent_cell) - { - const unsigned int K = adjacent_cell->cell_index; - const unsigned int delta_is_edge_in_K = - adjacent_cell->edge_within_cell; - - // figure out the direction of delta with respect to the cell - // K (in the orientation in which the user has given it to us) - const unsigned int first_edge_vertex = - (edges[*delta].orientation_status == Edge::forward ? - edges[*delta].vertex_indices[0] : - edges[*delta].vertex_indices[1]); - const unsigned int first_edge_vertex_in_K = - cells[K] - .vertex_indices[GeometryInfo::line_to_cell_vertices( - delta_is_edge_in_K, 0)]; - Assert( - first_edge_vertex == first_edge_vertex_in_K || - first_edge_vertex == - cells[K].vertex_indices[GeometryInfo< - dim>::line_to_cell_vertices(delta_is_edge_in_K, 1)], - ExcInternalError()); - - // now figure out which direction the each of the "opposite" - // edges needs to be oriented into. - for (unsigned int o_e = 0; - o_e < ParallelEdges::n_other_parallel_edges; - ++o_e) - { - // get the index of the opposite edge and select which its - // first vertex needs to be based on how the current edge - // is oriented in the current cell - const unsigned int opposite_edge = - cells[K].edge_indices[ParallelEdges< - dim>::parallel_edges[delta_is_edge_in_K][o_e]]; - const unsigned int first_opposite_edge_vertex = - cells[K].vertex_indices - [GeometryInfo::line_to_cell_vertices( - ParallelEdges< - dim>::parallel_edges[delta_is_edge_in_K][o_e], - (first_edge_vertex == first_edge_vertex_in_K ? 0 : - 1))]; - - // then determine the orientation of the edge based on - // whether the vertex we want to be the edge's first - // vertex is already the first vertex of the edge, or - // whether it points in the opposite direction - const typename Edge::OrientationStatus - opposite_edge_orientation = - (edges[opposite_edge].vertex_indices[0] == - first_opposite_edge_vertex ? - Edge::forward : - Edge::backward); - - // see if the opposite edge (there is only one in 2d) has - // already been oriented. - if (edges[opposite_edge].orientation_status == - Edge::not_oriented) - { - // the opposite edge is not yet oriented. do orient it - // and add it to Delta_k - edges[opposite_edge].orientation_status = - opposite_edge_orientation; - Delta_k.insert(opposite_edge); - } - else - { - // this opposite edge has already been oriented. it - // should be consistent with the current one in 2d, - // while in 3d it may in fact be mis-oriented, and in - // that case the mesh will not be orientable. indicate - // this by throwing an exception that we can catch - // further up; this has the advantage that we can - // propagate through a couple of functions without - // having to do error checking and without modifying - // the 'cells' array that the user gave us - if (dim == 2) - { - Assert(edges[opposite_edge].orientation_status == - opposite_edge_orientation, - ExcMeshNotOrientable()); - } - else if (dim == 3) - { - if (edges[opposite_edge].orientation_status != - opposite_edge_orientation) - throw ExcMeshNotOrientable(); - } - else - Assert(false, ExcNotImplemented()); - } - } - } - } - - // finally copy the new set to the previous one - // (corresponding to increasing 'k' by one in the - // algorithm) - Delta_k_minus_1 = Delta_k; - } - } - - - /** - * Given data structures @p cell_list and @p edge_list, where - * all edges are already oriented, rotate the cell with - * index @p cell_index in such a way that its local coordinate - * system matches the ones of the adjacent edges. Store the - * rotated order of vertices in raw_cells[cell_index]. - */ - template - void - rotate_cell(const std::vector> &cell_list, - const std::vector> &edge_list, - const unsigned int cell_index, - std::vector> &raw_cells) - { - // find the first vertex of the cell. this is the vertex where dim edges - // originate, so for each of the edges record which the starting vertex is - unsigned int starting_vertex_of_edge[GeometryInfo::lines_per_cell]; - for (unsigned int e = 0; e < GeometryInfo::lines_per_cell; ++e) - { - Assert(edge_list[cell_list[cell_index].edge_indices[e]] - .orientation_status != Edge::not_oriented, - ExcInternalError()); - if (edge_list[cell_list[cell_index].edge_indices[e]] - .orientation_status == Edge::forward) - starting_vertex_of_edge[e] = - edge_list[cell_list[cell_index].edge_indices[e]] - .vertex_indices[0]; - else - starting_vertex_of_edge[e] = - edge_list[cell_list[cell_index].edge_indices[e]] - .vertex_indices[1]; - } - - // find the vertex number that appears dim times. this will then be - // the vertex at which we want to locate the origin of the cell's - // coordinate system (i.e., vertex 0) - unsigned int origin_vertex_of_cell = numbers::invalid_unsigned_int; - switch (dim) - { - case 2: - { - // in 2d, we can simply enumerate the possibilities where the - // origin may be located because edges zero and one don't share - // any vertices, and the same for edges two and three - if ((starting_vertex_of_edge[0] == starting_vertex_of_edge[2]) || - (starting_vertex_of_edge[0] == starting_vertex_of_edge[3])) - origin_vertex_of_cell = starting_vertex_of_edge[0]; - else if ((starting_vertex_of_edge[1] == - starting_vertex_of_edge[2]) || - (starting_vertex_of_edge[1] == - starting_vertex_of_edge[3])) - origin_vertex_of_cell = starting_vertex_of_edge[1]; - else - Assert(false, ExcInternalError()); - - break; - } - - case 3: - { - // one could probably do something similar in 3d, but that seems - // more complicated than one wants to write down. just go - // through the list of possible starting vertices and check - for (origin_vertex_of_cell = 0; - origin_vertex_of_cell < GeometryInfo::vertices_per_cell; - ++origin_vertex_of_cell) - if (std::count(starting_vertex_of_edge, - starting_vertex_of_edge + - GeometryInfo::lines_per_cell, - cell_list[cell_index] - .vertex_indices[origin_vertex_of_cell]) == dim) - break; - Assert(origin_vertex_of_cell < - GeometryInfo::vertices_per_cell, - ExcInternalError()); - - break; - } - - default: - Assert(false, ExcNotImplemented()); - } - - // now rotate raw_cells[cell_index] in such a way that its orientation - // matches that of cell_list[cell_index] - switch (dim) - { - case 2: - { - // in 2d, we can literally rotate the cell until its origin - // matches the one that we have determined above should be - // the origin vertex - // - // when doing a rotation, take into account the ordering of - // vertices (not in clockwise or counter-clockwise sense) - while (raw_cells[cell_index].vertices[0] != origin_vertex_of_cell) - { - const unsigned int tmp = raw_cells[cell_index].vertices[0]; - raw_cells[cell_index].vertices[0] = - raw_cells[cell_index].vertices[1]; - raw_cells[cell_index].vertices[1] = - raw_cells[cell_index].vertices[3]; - raw_cells[cell_index].vertices[3] = - raw_cells[cell_index].vertices[2]; - raw_cells[cell_index].vertices[2] = tmp; - } - break; - } - - case 3: - { - // in 3d, the situation is a bit more complicated. from above, we - // now know which vertex is at the origin (because 3 edges - // originate from it), but that still leaves 3 possible rotations - // of the cube. the important realization is that we can choose - // any of them: in all 3 rotations, all edges originate from the - // one vertex, and that fixes the directions of all 12 edges in - // the cube because these 3 cover all 3 equivalence classes! - // consequently, we can select an arbitrary one among the - // permutations -- for example the following ones: - static const unsigned int cube_permutations[8][8] = { - {0, 1, 2, 3, 4, 5, 6, 7}, - {1, 5, 3, 7, 0, 4, 2, 6}, - {2, 6, 0, 4, 3, 7, 1, 5}, - {3, 2, 1, 0, 7, 6, 5, 4}, - {4, 0, 6, 2, 5, 1, 7, 3}, - {5, 4, 7, 6, 1, 0, 3, 2}, - {6, 7, 4, 5, 2, 3, 0, 1}, - {7, 3, 5, 1, 6, 2, 4, 0}}; - - unsigned int - temp_vertex_indices[GeometryInfo::vertices_per_cell]; - for (const unsigned int v : GeometryInfo::vertex_indices()) - temp_vertex_indices[v] = - raw_cells[cell_index] - .vertices[cube_permutations[origin_vertex_of_cell][v]]; - for (const unsigned int v : GeometryInfo::vertex_indices()) - raw_cells[cell_index].vertices[v] = temp_vertex_indices[v]; - - break; - } - - default: - { - Assert(false, ExcNotImplemented()); - } - } - } - - - /** - * Given a set of cells, find globally unique edge orientations - * and then rotate cells so that the coordinate system of the cell - * coincides with the coordinate systems of the adjacent edges. - */ - template - void - reorient(std::vector> &cells) - { - // first build the arrays that connect cells to edges and the other - // way around - std::vector> edge_list = build_edges(cells); - std::vector> cell_list = - build_cells_and_connect_edges(cells, edge_list); - - // then loop over all cells and start orienting parallel edge sets - // of cells that still have non-oriented edges - unsigned int next_cell_with_unoriented_edge = 0; - while ((next_cell_with_unoriented_edge = get_next_unoriented_cell( - cell_list, edge_list, next_cell_with_unoriented_edge)) != - numbers::invalid_unsigned_int) - { - // see which edge sets are still not oriented - // - // we do not need to look at each edge because if we orient edge - // 0, we will end up with edge 1 also oriented (in 2d; in 3d, there - // will be 3 other edges that are also oriented). there are only - // dim independent sets of edges, so loop over these. - // - // we need to check whether each one of these starter edges may - // already be oriented because the line (sheet) that connects - // globally parallel edges may be self-intersecting in the - // current cell - for (unsigned int l = 0; l < dim; ++l) - if (edge_list[cell_list[next_cell_with_unoriented_edge] - .edge_indices[ParallelEdges::starter_edges[l]]] - .orientation_status == Edge::not_oriented) - orient_one_set_of_parallel_edges( - cell_list, - edge_list, - next_cell_with_unoriented_edge, - ParallelEdges::starter_edges[l]); - - // ensure that we have really oriented all edges now, not just - // the starter edges - for (unsigned int l = 0; l < GeometryInfo::lines_per_cell; ++l) - Assert(edge_list[cell_list[next_cell_with_unoriented_edge] - .edge_indices[l]] - .orientation_status != Edge::not_oriented, - ExcInternalError()); - } - - // now that we have oriented all edges, we need to rotate cells - // so that the edges point in the right direction with the now - // rotated coordinate system - for (unsigned int c = 0; c < cells.size(); ++c) - rotate_cell(cell_list, edge_list, c, cells); - } - - - // overload of the function above for 1d -- there is nothing - // to orient in that case - void - reorient(std::vector> &) - {} - } // namespace - - template - void - consistently_order_cells(std::vector> &cells) - { - Assert(cells.size() != 0, - ExcMessage( - "List of elements to orient must have at least one cell")); - - // there is nothing for us to do in 1d - if (dim == 1) - return; - - // check if grids are already consistent. if so, do - // nothing. if not, then do the reordering - if (!is_consistent(cells)) - try - { - reorient(cells); - } - catch (const ExcMeshNotOrientable &) - { - // the mesh is not orientable. this is acceptable if we are in 3d, - // as class Triangulation knows how to handle this, but it is - // not in 2d; in that case, re-throw the exception - if (dim < 3) - throw; - } - } - - // define some transformations namespace internal { @@ -2015,34 +406,6 @@ namespace GridTools } } - template - std::map> - get_all_vertices_at_boundary(const Triangulation &tria) - { - std::map> vertex_map; - typename Triangulation::active_cell_iterator - cell = tria.begin_active(), - endc = tria.end(); - for (; cell != endc; ++cell) - { - for (const unsigned int i : cell->face_indices()) - { - const typename Triangulation::face_iterator &face = - cell->face(i); - if (face->at_boundary()) - { - for (unsigned j = 0; j < face->n_vertices(); ++j) - { - const Point &vertex = face->vertex(j); - const unsigned int vertex_index = face->vertex_index(j); - vertex_map[vertex_index] = vertex; - } - } - } - } - return vertex_map; - } - /** * Distort a triangulation in * some random way. @@ -3044,64 +1407,6 @@ namespace GridTools } - template - std::vector< - std::set::active_cell_iterator>> - vertex_to_cell_map(const Triangulation &triangulation) - { - std::vector< - std::set::active_cell_iterator>> - vertex_to_cell_map(triangulation.n_vertices()); - typename Triangulation::active_cell_iterator - cell = triangulation.begin_active(), - endc = triangulation.end(); - for (; cell != endc; ++cell) - for (const unsigned int i : cell->vertex_indices()) - vertex_to_cell_map[cell->vertex_index(i)].insert(cell); - - // Check if mesh has hanging nodes. Do this only locally to - // prevent communication and possible deadlock. - if (triangulation.Triangulation::has_hanging_nodes()) - { - Assert(triangulation.all_reference_cells_are_hyper_cube(), - ExcNotImplemented()); - - // Take care of hanging nodes - cell = triangulation.begin_active(); - for (; cell != endc; ++cell) - { - for (const unsigned int i : cell->face_indices()) - { - if ((cell->at_boundary(i) == false) && - (cell->neighbor(i)->is_active())) - { - typename Triangulation::active_cell_iterator - adjacent_cell = cell->neighbor(i); - for (unsigned int j = 0; j < cell->face(i)->n_vertices(); - ++j) - vertex_to_cell_map[cell->face(i)->vertex_index(j)].insert( - adjacent_cell); - } - } - - // in 3d also loop over the edges - if (dim == 3) - { - for (unsigned int i = 0; i < cell->n_lines(); ++i) - if (cell->line(i)->has_children()) - // the only place where this vertex could have been - // hiding is on the mid-edge point of the edge we - // are looking at - vertex_to_cell_map[cell->line(i)->child(0)->vertex_index(1)] - .insert(cell); - } - } - } - - return vertex_to_cell_map; - } - - template std::map @@ -3458,105 +1763,6 @@ namespace GridTools - template - void - get_face_connectivity_of_cells( - const Triangulation &triangulation, - DynamicSparsityPattern &cell_connectivity) - { - cell_connectivity.reinit(triangulation.n_active_cells(), - triangulation.n_active_cells()); - - // loop over all cells and their neighbors to build the sparsity - // pattern. note that it's a bit hard to enter all the connections when a - // neighbor has children since we would need to find out which of its - // children is adjacent to the current cell. this problem can be omitted - // if we only do something if the neighbor has no children -- in that case - // it is either on the same or a coarser level than we are. in return, we - // have to add entries in both directions for both cells - for (const auto &cell : triangulation.active_cell_iterators()) - { - const unsigned int index = cell->active_cell_index(); - cell_connectivity.add(index, index); - for (auto f : cell->face_indices()) - if ((cell->at_boundary(f) == false) && - (cell->neighbor(f)->has_children() == false)) - { - const unsigned int other_index = - cell->neighbor(f)->active_cell_index(); - cell_connectivity.add(index, other_index); - cell_connectivity.add(other_index, index); - } - } - } - - - - template - void - get_vertex_connectivity_of_cells( - const Triangulation &triangulation, - DynamicSparsityPattern &cell_connectivity) - { - std::vector> vertex_to_cell( - triangulation.n_vertices()); - for (const auto &cell : triangulation.active_cell_iterators()) - { - for (const unsigned int v : cell->vertex_indices()) - vertex_to_cell[cell->vertex_index(v)].push_back( - cell->active_cell_index()); - } - - cell_connectivity.reinit(triangulation.n_active_cells(), - triangulation.n_active_cells()); - for (const auto &cell : triangulation.active_cell_iterators()) - { - for (const unsigned int v : cell->vertex_indices()) - for (unsigned int n = 0; - n < vertex_to_cell[cell->vertex_index(v)].size(); - ++n) - cell_connectivity.add(cell->active_cell_index(), - vertex_to_cell[cell->vertex_index(v)][n]); - } - } - - - template - void - get_vertex_connectivity_of_cells_on_level( - const Triangulation &triangulation, - const unsigned int level, - DynamicSparsityPattern &cell_connectivity) - { - std::vector> vertex_to_cell( - triangulation.n_vertices()); - for (typename Triangulation::cell_iterator cell = - triangulation.begin(level); - cell != triangulation.end(level); - ++cell) - { - for (const unsigned int v : cell->vertex_indices()) - vertex_to_cell[cell->vertex_index(v)].push_back(cell->index()); - } - - cell_connectivity.reinit(triangulation.n_cells(level), - triangulation.n_cells(level)); - for (typename Triangulation::cell_iterator cell = - triangulation.begin(level); - cell != triangulation.end(level); - ++cell) - { - for (const unsigned int v : cell->vertex_indices()) - for (unsigned int n = 0; - n < vertex_to_cell[cell->vertex_index(v)].size(); - ++n) - cell_connectivity.add(cell->index(), - vertex_to_cell[cell->vertex_index(v)][n]); - } - } - - - template void partition_triangulation(const unsigned int n_partitions, @@ -4803,68 +3009,6 @@ namespace GridTools - template - void - remove_hanging_nodes(Triangulation &tria, - const bool isotropic, - const unsigned int max_iterations) - { - unsigned int iter = 0; - bool continue_refinement = true; - - while (continue_refinement && (iter < max_iterations)) - { - if (max_iterations != numbers::invalid_unsigned_int) - iter++; - continue_refinement = false; - - for (const auto &cell : tria.active_cell_iterators()) - for (const unsigned int j : cell->face_indices()) - if (cell->at_boundary(j) == false && - cell->neighbor(j)->has_children()) - { - if (isotropic) - { - cell->set_refine_flag(); - continue_refinement = true; - } - else - continue_refinement |= cell->flag_for_face_refinement(j); - } - - tria.execute_coarsening_and_refinement(); - } - } - - template - void - remove_anisotropy(Triangulation &tria, - const double max_ratio, - const unsigned int max_iterations) - { - unsigned int iter = 0; - bool continue_refinement = true; - - while (continue_refinement && (iter < max_iterations)) - { - iter++; - continue_refinement = false; - for (const auto &cell : tria.active_cell_iterators()) - { - std::pair info = - GridTools::get_longest_direction(cell); - if (info.second > max_ratio) - { - cell->set_refine_flag( - RefinementCase::cut_axis(info.first)); - continue_refinement = true; - } - } - tria.execute_coarsening_and_refinement(); - } - } - - template void regularize_corner_cells(Triangulation &tria, @@ -6422,25 +4566,6 @@ namespace GridTools - template - std::map> - extract_used_vertices(const Triangulation &container, - const Mapping &mapping) - { - std::map> result; - for (const auto &cell : container.active_cell_iterators()) - { - if (!cell->is_artificial()) - { - const auto vs = mapping.get_vertices(cell); - for (unsigned int i = 0; i < vs.size(); ++i) - result[cell->vertex_index(i)] = vs[i]; - } - } - return result; - } - - template unsigned int find_closest_vertex(const std::map> &vertices, diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index 3f4fd209e7..ea2f87577a 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -82,11 +82,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) deal_II_space_dimension> &triangulation); - template std::map> - extract_used_vertices( - const Triangulation &mesh, - const Mapping &mapping); - template std::pair< typename Triangulation::active_cell_iterator, @@ -263,39 +258,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) #if deal_II_dimension <= deal_II_space_dimension namespace GridTools \{ - template std::tuple>, - std::vector>, - SubCellData> - get_coarse_mesh_description( - const Triangulation &tria); - - template void - delete_unused_vertices(std::vector> &, - std::vector> &, - SubCellData &); - - template void - delete_duplicated_vertices(std::vector> &, - std::vector> &, - SubCellData &, - std::vector &, - double); - - template void - invert_all_negative_measure_cells( - const std::vector> &, - std::vector> &); - - template std::size_t - invert_cells_with_negative_measure( - const std::vector> &, - std::vector> &); - -# if deal_II_dimension == deal_II_space_dimension - template void - consistently_order_cells(std::vector> &); -# endif - template void shift( const Tensor<1, deal_II_space_dimension> &, @@ -334,25 +296,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const bool, const unsigned int); - template void - get_face_connectivity_of_cells( - const Triangulation - &triangulation, - DynamicSparsityPattern &cell_connectivity); - - template void - get_vertex_connectivity_of_cells( - const Triangulation - &triangulation, - DynamicSparsityPattern &cell_connectivity); - - template void - get_vertex_connectivity_of_cells_on_level( - const Triangulation - &triangulation, - const unsigned int level, - DynamicSparsityPattern &cell_connectivity); - template void partition_triangulation( const unsigned int, @@ -410,17 +353,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) get_locally_owned_vertices( const Triangulation &); - template std::map> - get_all_vertices_at_boundary( - const Triangulation &tria); - - template std::vector< - std::set::active_cell_iterator>> - vertex_to_cell_map( - const Triangulation - &triangulation); - template std::vector>> vertex_to_cell_centers_directions( const Triangulation &mesh, diff --git a/source/grid/grid_tools_topology.cc b/source/grid/grid_tools_topology.cc new file mode 100644 index 0000000000..3e0b329e58 --- /dev/null +++ b/source/grid/grid_tools_topology.cc @@ -0,0 +1,1920 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2001 - 2023 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. +// +// --------------------------------------------------------------------- + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace GridTools +{ + // Generic functions for appending face data in 2d or 3d. TODO: we can + // remove these once we have 'if constexpr'. + namespace internal + { + inline void + append_face_data(const CellData<1> &face_data, SubCellData &subcell_data) + { + subcell_data.boundary_lines.push_back(face_data); + } + + + + inline void + append_face_data(const CellData<2> &face_data, SubCellData &subcell_data) + { + subcell_data.boundary_quads.push_back(face_data); + } + + + + // Lexical comparison for sorting CellData objects. + template + struct CellDataComparator + { + bool + operator()(const CellData &a, + const CellData &b) const + { + // Check vertices: + if (std::lexicographical_compare(std::begin(a.vertices), + std::end(a.vertices), + std::begin(b.vertices), + std::end(b.vertices))) + return true; + // it should never be necessary to check the material or manifold + // ids as a 'tiebreaker' (since they must be equal if the vertex + // indices are equal). Assert it anyway: +#ifdef DEBUG + if (std::equal(std::begin(a.vertices), + std::end(a.vertices), + std::begin(b.vertices))) + { + Assert(a.material_id == b.material_id && + a.manifold_id == b.manifold_id, + ExcMessage( + "Two CellData objects with equal vertices must " + "have the same material/boundary ids and manifold " + "ids.")); + } +#endif + return false; + } + }; + + + /** + * get_coarse_mesh_description() needs to store face data for dim>1, but + * we can not have this code in the function, as this requires either an + * instantiation of CellData, or constexpr if. We use a class with + * specialization instead for now. + * + * Data on faces is added with insert_face_data() and then retrieved with + * get(). + */ + template + class FaceDataHelper + { + public: + /** + * Store the data about the given face @p face. + */ + template + void + insert_face_data(const FaceIteratorType &face) + { + CellData face_cell_data(face->n_vertices()); + for (unsigned int vertex_n = 0; vertex_n < face->n_vertices(); + ++vertex_n) + face_cell_data.vertices[vertex_n] = face->vertex_index(vertex_n); + face_cell_data.boundary_id = face->boundary_id(); + face_cell_data.manifold_id = face->manifold_id(); + + face_data.insert(std::move(face_cell_data)); + } + + /** + * Return the @p subcell_data with the stored information. + */ + SubCellData + get() + { + SubCellData subcell_data; + + for (const CellData &face_cell_data : face_data) + internal::append_face_data(face_cell_data, subcell_data); + return subcell_data; + } + + + private: + std::set, internal::CellDataComparator> + face_data; + }; + + + // Do nothing for dim=1: + template <> + class FaceDataHelper<1> + { + public: + template + void + insert_face_data(const FaceIteratorType &) + {} + + SubCellData + get() + { + return SubCellData(); + } + }; + } // namespace internal + + + + template + std:: + tuple>, std::vector>, SubCellData> + get_coarse_mesh_description(const Triangulation &tria) + { + Assert(tria.n_levels() >= 1, + ExcMessage("The input triangulation must be non-empty.")); + + std::vector> vertices = tria.get_vertices(); + std::vector> cells; + + internal::FaceDataHelper face_data; + std::set, internal::CellDataComparator<1>> + line_data; // only used in 3d + + for (const auto &cell : tria.cell_iterators_on_level(0)) + { + // Save cell data + CellData cell_data(cell->n_vertices()); + for (const unsigned int cell_vertex_n : cell->vertex_indices()) + { + Assert(cell->vertex_index(cell_vertex_n) < vertices.size(), + ExcInternalError()); + cell_data.vertices[cell_vertex_n] = + cell->vertex_index(cell_vertex_n); + } + cell_data.material_id = cell->material_id(); + cell_data.manifold_id = cell->manifold_id(); + cells.emplace_back(std::move(cell_data)); + + // Save face data + if (dim > 1) + { + for (const unsigned int face_n : cell->face_indices()) + // We don't need to insert anything if we have default values + { + const auto face = cell->face(face_n); + if (face->boundary_id() != numbers::internal_face_boundary_id || + face->manifold_id() != numbers::flat_manifold_id) + face_data.insert_face_data(face); + } + } + // Save line data + if (dim == 3) + { + for (unsigned int line_n = 0; line_n < cell->n_lines(); ++line_n) + { + const auto line = cell->line(line_n); + // We don't need to insert anything if we have default values + if (line->boundary_id() != numbers::internal_face_boundary_id || + line->manifold_id() != numbers::flat_manifold_id) + { + CellData<1> line_cell_data(line->n_vertices()); + for (const unsigned int vertex_n : line->vertex_indices()) + line_cell_data.vertices[vertex_n] = + line->vertex_index(vertex_n); + line_cell_data.boundary_id = line->boundary_id(); + line_cell_data.manifold_id = line->manifold_id(); + line_data.insert(std::move(line_cell_data)); + } + } + } + } + + SubCellData subcell_data = face_data.get(); + + if (dim == 3) + for (const CellData<1> &face_line_data : line_data) + subcell_data.boundary_lines.push_back(face_line_data); + + // We end up with a 'vertices' array that uses some of the entries, + // but not all -- specifically, all vertices referenced by level-0 + // cells. We can compress the array: + GridTools::delete_unused_vertices(vertices, cells, subcell_data); + + return std::tuple>, + std::vector>, + SubCellData>(std::move(vertices), + std::move(cells), + std::move(subcell_data)); + } + + + + template + void + delete_unused_vertices(std::vector> &vertices, + std::vector> &cells, + SubCellData &subcelldata) + { + Assert( + subcelldata.check_consistency(dim), + ExcMessage( + "Invalid SubCellData supplied according to ::check_consistency(). " + "This is caused by data containing objects for the wrong dimension.")); + + // first check which vertices are actually used + std::vector vertex_used(vertices.size(), false); + for (unsigned int c = 0; c < cells.size(); ++c) + for (unsigned int v = 0; v < cells[c].vertices.size(); ++v) + { + Assert(cells[c].vertices[v] < vertices.size(), + ExcMessage("Invalid vertex index encountered! cells[" + + Utilities::int_to_string(c) + "].vertices[" + + Utilities::int_to_string(v) + "]=" + + Utilities::int_to_string(cells[c].vertices[v]) + + " is invalid, because only " + + Utilities::int_to_string(vertices.size()) + + " vertices were supplied.")); + vertex_used[cells[c].vertices[v]] = true; + } + + + // then renumber the vertices that are actually used in the same order as + // they were beforehand + const unsigned int invalid_vertex = numbers::invalid_unsigned_int; + std::vector new_vertex_numbers(vertices.size(), + invalid_vertex); + unsigned int next_free_number = 0; + for (unsigned int i = 0; i < vertices.size(); ++i) + if (vertex_used[i] == true) + { + new_vertex_numbers[i] = next_free_number; + ++next_free_number; + } + + // next replace old vertex numbers by the new ones + for (unsigned int c = 0; c < cells.size(); ++c) + for (auto &v : cells[c].vertices) + v = new_vertex_numbers[v]; + + // same for boundary data + for (unsigned int c = 0; c < subcelldata.boundary_lines.size(); // NOLINT + ++c) + for (unsigned int v = 0; + v < subcelldata.boundary_lines[c].vertices.size(); + ++v) + { + Assert(subcelldata.boundary_lines[c].vertices[v] < + new_vertex_numbers.size(), + ExcMessage( + "Invalid vertex index in subcelldata.boundary_lines. " + "subcelldata.boundary_lines[" + + Utilities::int_to_string(c) + "].vertices[" + + Utilities::int_to_string(v) + "]=" + + Utilities::int_to_string( + subcelldata.boundary_lines[c].vertices[v]) + + " is invalid, because only " + + Utilities::int_to_string(vertices.size()) + + " vertices were supplied.")); + subcelldata.boundary_lines[c].vertices[v] = + new_vertex_numbers[subcelldata.boundary_lines[c].vertices[v]]; + } + + for (unsigned int c = 0; c < subcelldata.boundary_quads.size(); // NOLINT + ++c) + for (unsigned int v = 0; + v < subcelldata.boundary_quads[c].vertices.size(); + ++v) + { + Assert(subcelldata.boundary_quads[c].vertices[v] < + new_vertex_numbers.size(), + ExcMessage( + "Invalid vertex index in subcelldata.boundary_quads. " + "subcelldata.boundary_quads[" + + Utilities::int_to_string(c) + "].vertices[" + + Utilities::int_to_string(v) + "]=" + + Utilities::int_to_string( + subcelldata.boundary_quads[c].vertices[v]) + + " is invalid, because only " + + Utilities::int_to_string(vertices.size()) + + " vertices were supplied.")); + + subcelldata.boundary_quads[c].vertices[v] = + new_vertex_numbers[subcelldata.boundary_quads[c].vertices[v]]; + } + + // finally copy over the vertices which we really need to a new array and + // replace the old one by the new one + std::vector> tmp; + tmp.reserve(std::count(vertex_used.begin(), vertex_used.end(), true)); + for (unsigned int v = 0; v < vertices.size(); ++v) + if (vertex_used[v] == true) + tmp.push_back(vertices[v]); + swap(vertices, tmp); + } + + + + template + void + delete_duplicated_vertices(std::vector> &vertices, + std::vector> &cells, + SubCellData &subcelldata, + std::vector &considered_vertices, + const double tol) + { + if (tol == 0.0) + return; // nothing to do per definition + + AssertIndexRange(2, vertices.size()); + std::vector new_vertex_numbers(vertices.size()); + std::iota(new_vertex_numbers.begin(), new_vertex_numbers.end(), 0); + + // if the considered_vertices vector is empty, consider all vertices + if (considered_vertices.empty()) + considered_vertices = new_vertex_numbers; + Assert(considered_vertices.size() <= vertices.size(), ExcInternalError()); + + // The algorithm below improves upon the naive O(n^2) algorithm by first + // sorting vertices by their value in one component and then only + // comparing vertices for equality which are nearly equal in that + // component. For example, if @p vertices form a cube, then we will only + // compare points that have the same x coordinate when we try to find + // duplicated vertices. + + // Start by finding the longest coordinate direction. This minimizes the + // number of points that need to be compared against each-other in a + // single set for typical geometries. + const BoundingBox bbox(vertices); + + unsigned int longest_coordinate_direction = 0; + double longest_coordinate_length = bbox.side_length(0); + for (unsigned int d = 1; d < spacedim; ++d) + { + const double coordinate_length = bbox.side_length(d); + if (longest_coordinate_length < coordinate_length) + { + longest_coordinate_length = coordinate_length; + longest_coordinate_direction = d; + } + } + + // Sort vertices (while preserving their vertex numbers) along that + // coordinate direction: + std::vector>> sorted_vertices; + sorted_vertices.reserve(vertices.size()); + for (const unsigned int vertex_n : considered_vertices) + { + AssertIndexRange(vertex_n, vertices.size()); + sorted_vertices.emplace_back(vertex_n, vertices[vertex_n]); + } + std::sort(sorted_vertices.begin(), + sorted_vertices.end(), + [&](const std::pair> &a, + const std::pair> &b) { + return a.second[longest_coordinate_direction] < + b.second[longest_coordinate_direction]; + }); + + auto within_tolerance = [=](const Point &a, + const Point &b) { + for (unsigned int d = 0; d < spacedim; ++d) + if (std::abs(a[d] - b[d]) > tol) + return false; + return true; + }; + + // Find a range of numbers that have the same component in the longest + // coordinate direction: + auto range_start = sorted_vertices.begin(); + while (range_start != sorted_vertices.end()) + { + auto range_end = range_start + 1; + while (range_end != sorted_vertices.end() && + std::abs(range_end->second[longest_coordinate_direction] - + range_start->second[longest_coordinate_direction]) < + tol) + ++range_end; + + // preserve behavior with older versions of this function by replacing + // higher vertex numbers by lower vertex numbers + std::sort(range_start, + range_end, + [](const std::pair> &a, + const std::pair> &b) { + return a.first < b.first; + }); + + // Now de-duplicate [range_start, range_end) + // + // We have identified all points that are within a strip of width 'tol' + // in one coordinate direction. Now we need to figure out which of these + // are also close in other coordinate directions. If two are close, we + // can mark the second one for deletion. + for (auto reference = range_start; reference != range_end; ++reference) + { + if (reference->first != numbers::invalid_unsigned_int) + for (auto it = reference + 1; it != range_end; ++it) + { + if (within_tolerance(reference->second, it->second)) + { + new_vertex_numbers[it->first] = reference->first; + // skip the replaced vertex in the future + it->first = numbers::invalid_unsigned_int; + } + } + } + range_start = range_end; + } + + // now we got a renumbering list. simply renumber all vertices + // (non-duplicate vertices get renumbered to themselves, so nothing bad + // happens). after that, the duplicate vertices will be unused, so call + // delete_unused_vertices() to do that part of the job. + for (auto &cell : cells) + for (auto &vertex_index : cell.vertices) + vertex_index = new_vertex_numbers[vertex_index]; + for (auto &quad : subcelldata.boundary_quads) + for (auto &vertex_index : quad.vertices) + vertex_index = new_vertex_numbers[vertex_index]; + for (auto &line : subcelldata.boundary_lines) + for (auto &vertex_index : line.vertices) + vertex_index = new_vertex_numbers[vertex_index]; + + delete_unused_vertices(vertices, cells, subcelldata); + } + + + + template + void + delete_duplicated_vertices(std::vector> &vertices, + const double tol) + { + if (vertices.empty()) + return; + + // 1) map point to local vertex index + std::map, unsigned int, FloatingPointComparator> + map_point_to_local_vertex_index{FloatingPointComparator(tol)}; + + // 2) initialize map with existing points uniquely + for (unsigned int i = 0; i < vertices.size(); ++i) + map_point_to_local_vertex_index[vertices[i]] = i; + + // no duplicate points are found + if (map_point_to_local_vertex_index.size() == vertices.size()) + return; + + // 3) remove duplicate entries from vertices + vertices.resize(map_point_to_local_vertex_index.size()); + { + unsigned int j = 0; + for (const auto &p : map_point_to_local_vertex_index) + vertices[j++] = p.first; + } + } + + + + template + std::size_t + invert_cells_with_negative_measure( + const std::vector> &all_vertices, + std::vector> &cells) + { + // This function is presently only implemented for volumetric (codimension + // 0) elements. + + if (dim == 1) + return 0; + if (dim == 2 && spacedim == 3) + Assert(false, ExcNotImplemented()); + + std::size_t n_negative_cells = 0; + std::size_t cell_no = 0; + for (auto &cell : cells) + { + const ArrayView vertices(cell.vertices); + // Some pathologically twisted cells can have exactly zero measure but + // we can still fix them + if (GridTools::cell_measure(all_vertices, vertices) <= 0) + { + ++n_negative_cells; + const auto reference_cell = + ReferenceCell::n_vertices_to_type(dim, vertices.size()); + + if (reference_cell.is_hyper_cube()) + { + if (dim == 2) + { + // flip the cell across the y = x line in 2d + std::swap(cell.vertices[1], cell.vertices[2]); + } + else if (dim == 3) + { + // swap the front and back faces in 3d + std::swap(cell.vertices[0], cell.vertices[2]); + std::swap(cell.vertices[1], cell.vertices[3]); + std::swap(cell.vertices[4], cell.vertices[6]); + std::swap(cell.vertices[5], cell.vertices[7]); + } + } + else if (reference_cell.is_simplex()) + { + // By basic rules for computing determinants we can just swap + // two vertices to fix a negative volume. Arbitrarily pick the + // last two. + std::swap(cell.vertices[cell.vertices.size() - 2], + cell.vertices[cell.vertices.size() - 1]); + } + else if (reference_cell == ReferenceCells::Wedge) + { + // swap the two triangular faces + std::swap(cell.vertices[0], cell.vertices[3]); + std::swap(cell.vertices[1], cell.vertices[4]); + std::swap(cell.vertices[2], cell.vertices[5]); + } + else if (reference_cell == ReferenceCells::Pyramid) + { + // Try swapping two vertices in the base - perhaps things were + // read in the UCD (counter-clockwise) order instead of lexical + std::swap(cell.vertices[2], cell.vertices[3]); + } + else + { + AssertThrow(false, ExcNotImplemented()); + } + // Check whether the resulting cell is now ok. + // If not, then the grid is seriously broken and + // we just give up. + AssertThrow(GridTools::cell_measure(all_vertices, vertices) > 0, + ExcGridHasInvalidCell(cell_no)); + } + ++cell_no; + } + return n_negative_cells; + } + + + template + void + invert_all_negative_measure_cells( + const std::vector> &all_vertices, + std::vector> &cells) + { + const std::size_t n_negative_cells = + invert_cells_with_negative_measure(all_vertices, cells); + + // We assume that all cells of a grid have + // either positive or negative volumes but + // not both mixed. Although above reordering + // might work also on single cells, grids + // with both kind of cells are very likely to + // be broken. Check for this here. + AssertThrow(n_negative_cells == 0 || n_negative_cells == cells.size(), + ExcMessage( + std::string( + "This function assumes that either all cells have positive " + "volume, or that all cells have been specified in an " + "inverted vertex order so that their volume is negative. " + "(In the latter case, this class automatically inverts " + "every cell.) However, the mesh you have specified " + "appears to have both cells with positive and cells with " + "negative volume. You need to check your mesh which " + "cells these are and how they got there.\n" + "As a hint, of the total ") + + std::to_string(cells.size()) + " cells in the mesh, " + + std::to_string(n_negative_cells) + + " appear to have a negative volume.")); + } + + + + // Functions and classes for consistently_order_cells + namespace + { + /** + * A simple data structure denoting an edge, i.e., the ordered pair + * of its vertex indices. This is only used in the is_consistent() + * function. + */ + struct CheapEdge + { + /** + * Construct an edge from the global indices of its two vertices. + */ + CheapEdge(const unsigned int v0, const unsigned int v1) + : v0(v0) + , v1(v1) + {} + + /** + * Comparison operator for edges. It compares based on the + * lexicographic ordering of the two vertex indices. + */ + bool + operator<(const CheapEdge &e) const + { + return ((v0 < e.v0) || ((v0 == e.v0) && (v1 < e.v1))); + } + + private: + /** + * The global indices of the vertices that define the edge. + */ + const unsigned int v0, v1; + }; + + + /** + * A function that determines whether the edges in a mesh are + * already consistently oriented. It does so by adding all edges + * of all cells into a set (which automatically eliminates + * duplicates) but before that checks whether the reverse edge is + * already in the set -- which would imply that a neighboring cell + * is inconsistently oriented. + */ + template + bool + is_consistent(const std::vector> &cells) + { + std::set edges; + + for (typename std::vector>::const_iterator c = + cells.begin(); + c != cells.end(); + ++c) + { + // construct the edges in reverse order. for each of them, + // ensure that the reverse edge is not yet in the list of + // edges (return false if the reverse edge already *is* in + // the list) and then add the actual edge to it; std::set + // eliminates duplicates automatically + for (unsigned int l = 0; l < GeometryInfo::lines_per_cell; ++l) + { + const CheapEdge reverse_edge( + c->vertices[GeometryInfo::line_to_cell_vertices(l, 1)], + c->vertices[GeometryInfo::line_to_cell_vertices(l, 0)]); + if (edges.find(reverse_edge) != edges.end()) + return false; + + + // ok, not. insert edge in correct order + const CheapEdge correct_edge( + c->vertices[GeometryInfo::line_to_cell_vertices(l, 0)], + c->vertices[GeometryInfo::line_to_cell_vertices(l, 1)]); + edges.insert(correct_edge); + } + } + + // no conflicts found, so return true + return true; + } + + + /** + * A structure that describes some properties of parallel edges + * such as what starter edges are (i.e., representative elements + * of the sets of parallel edges within a cell) and what the set + * of parallel edges to each edge is. + */ + template + struct ParallelEdges + { + /** + * An array that contains the indices of dim edges that can + * serve as (arbitrarily chosen) starting points for the + * dim sets of parallel edges within each cell. + */ + static const unsigned int starter_edges[dim]; + + /** + * Number and indices of all of those edges parallel to each of the + * edges in a cell. + */ + static const unsigned int n_other_parallel_edges = (1 << (dim - 1)) - 1; + static const unsigned int + parallel_edges[GeometryInfo::lines_per_cell] + [n_other_parallel_edges]; + }; + + template <> + const unsigned int ParallelEdges<2>::starter_edges[2] = {0, 2}; + + template <> + const unsigned int ParallelEdges<2>::parallel_edges[4][1] = {{1}, + {0}, + {3}, + {2}}; + + template <> + const unsigned int ParallelEdges<3>::starter_edges[3] = {0, 2, 8}; + + template <> + const unsigned int ParallelEdges<3>::parallel_edges[12][3] = { + {1, 4, 5}, // line 0 + {0, 4, 5}, // line 1 + {3, 6, 7}, // line 2 + {2, 6, 7}, // line 3 + {0, 1, 5}, // line 4 + {0, 1, 4}, // line 5 + {2, 3, 7}, // line 6 + {2, 3, 6}, // line 7 + {9, 10, 11}, // line 8 + {8, 10, 11}, // line 9 + {8, 9, 11}, // line 10 + {8, 9, 10} // line 11 + }; + + + /** + * A structure that store the index of a cell and, crucially, how a + * given edge relates to this cell. + */ + struct AdjacentCell + { + /** + * Default constructor. Initialize the fields with invalid values. + */ + AdjacentCell() + : cell_index(numbers::invalid_unsigned_int) + , edge_within_cell(numbers::invalid_unsigned_int) + {} + + /** + * Constructor. Initialize the fields with the given values. + */ + AdjacentCell(const unsigned int cell_index, + const unsigned int edge_within_cell) + : cell_index(cell_index) + , edge_within_cell(edge_within_cell) + {} + + + unsigned int cell_index; + unsigned int edge_within_cell; + }; + + + + template + class AdjacentCells; + + /** + * A class that represents all of the cells adjacent to a given edge. + * This class corresponds to the 2d case where each edge has at most + * two adjacent cells. + */ + template <> + class AdjacentCells<2> + { + public: + /** + * An iterator that allows iterating over all cells adjacent + * to the edge represented by the current object. + */ + using const_iterator = const AdjacentCell *; + + /** + * Add the given cell to the collection of cells adjacent to + * the edge this object corresponds to. Since we are covering + * the 2d case, the set of adjacent cells currently + * represented by this object must have either zero or + * one element already, since we can not add more than two + * adjacent cells for each edge. + */ + void + push_back(const AdjacentCell &adjacent_cell) + { + if (adjacent_cells[0].cell_index == numbers::invalid_unsigned_int) + adjacent_cells[0] = adjacent_cell; + else + { + Assert(adjacent_cells[1].cell_index == + numbers::invalid_unsigned_int, + ExcInternalError()); + adjacent_cells[1] = adjacent_cell; + } + } + + + /** + * Return an iterator to the first valid cell stored as adjacent to the + * edge represented by the current object. + */ + const_iterator + begin() const + { + return adjacent_cells; + } + + + /** + * Return an iterator to the element past the last valid cell stored + * as adjacent to the edge represented by the current object. + * @return + */ + const_iterator + end() const + { + // check whether the current object stores zero, one, or two + // adjacent cells, and use this to point to the element past the + // last valid one + if (adjacent_cells[0].cell_index == numbers::invalid_unsigned_int) + return adjacent_cells; + else if (adjacent_cells[1].cell_index == numbers::invalid_unsigned_int) + return adjacent_cells + 1; + else + return adjacent_cells + 2; + } + + private: + /** + * References to the (at most) two cells that are adjacent to + * the edge this object corresponds to. Unused elements are + * default-initialized and have invalid values; in particular, + * their cell_index field equals numbers::invalid_unsigned_int. + */ + AdjacentCell adjacent_cells[2]; + }; + + + + /** + * A class that represents all of the cells adjacent to a given edge. + * This class corresponds to the 3d case where each edge can have an + * arbitrary number of adjacent cells. We represent this as a + * std::vector, from which class the current one is + * derived and from which it inherits all of its member functions. + */ + template <> + class AdjacentCells<3> : public std::vector + {}; + + + /** + * A class that describes all of the relevant properties of an + * edge. For the purpose of what we do here, that includes the + * indices of the two vertices, and the indices of the adjacent + * cells (together with a description *where* in each of the + * adjacent cells the edge is located). It also includes the + * (global) direction of the edge: either from the first vertex to + * the second, the other way around, or so far undetermined. + */ + template + class Edge + { + public: + /** + * Constructor. Create the edge based on the information given + * in @p cell, and selecting the edge with number @p edge_number + * within this cell. Initialize the edge as unoriented. + */ + Edge(const CellData &cell, const unsigned int edge_number) + : orientation_status(not_oriented) + { + Assert(edge_number < GeometryInfo::lines_per_cell, + ExcInternalError()); + + // copy vertices for this particular line + vertex_indices[0] = + cell + .vertices[GeometryInfo::line_to_cell_vertices(edge_number, 0)]; + vertex_indices[1] = + cell + .vertices[GeometryInfo::line_to_cell_vertices(edge_number, 1)]; + + // bring them into standard orientation + if (vertex_indices[0] > vertex_indices[1]) + std::swap(vertex_indices[0], vertex_indices[1]); + } + + /** + * Comparison operator for edges. It compares based on the + * lexicographic ordering of the two vertex indices. + */ + bool + operator<(const Edge &e) const + { + return ((vertex_indices[0] < e.vertex_indices[0]) || + ((vertex_indices[0] == e.vertex_indices[0]) && + (vertex_indices[1] < e.vertex_indices[1]))); + } + + /** + * Compare two edges for equality based on their vertex indices. + */ + bool + operator==(const Edge &e) const + { + return ((vertex_indices[0] == e.vertex_indices[0]) && + (vertex_indices[1] == e.vertex_indices[1])); + } + + /** + * The global indices of the two vertices that bound this edge. These + * will be ordered so that the first index is less than the second. + */ + unsigned int vertex_indices[2]; + + /** + * An enum that indicates the direction of this edge with + * regard to the two vertices that bound it. + */ + enum OrientationStatus + { + not_oriented, + forward, + backward + }; + + OrientationStatus orientation_status; + + /** + * Store the set of cells adjacent to this edge (these cells then + * also store *where* in the cell the edge is located). + */ + AdjacentCells adjacent_cells; + }; + + + + /** + * A data structure that represents a cell with all of its vertices + * and edges. + */ + template + struct Cell + { + /** + * Construct a Cell object from a CellData object. Also take a + * (sorted) list of edges and to point the edges of the current + * object into this list of edges. + */ + Cell(const CellData &c, const std::vector> &edge_list) + { + for (const unsigned int i : GeometryInfo::vertex_indices()) + vertex_indices[i] = c.vertices[i]; + + // now for each of the edges of this cell, find the location inside the + // given edge_list array and store than index + for (unsigned int l = 0; l < GeometryInfo::lines_per_cell; ++l) + { + const Edge e(c, l); + edge_indices[l] = + (std::lower_bound(edge_list.begin(), edge_list.end(), e) - + edge_list.begin()); + Assert(edge_indices[l] < edge_list.size(), ExcInternalError()); + Assert(edge_list[edge_indices[l]] == e, ExcInternalError()); + } + } + + /** + * A list of global indices for the vertices that bound this cell. + */ + unsigned int vertex_indices[GeometryInfo::vertices_per_cell]; + + /** + * A list of indices into the 'edge_list' array passed to the constructor + * for the edges of the current cell. + */ + unsigned int edge_indices[GeometryInfo::lines_per_cell]; + }; + + + + template + class EdgeDeltaSet; + + /** + * A class that represents by how much the set of parallel edges + * grows in each step. In the graph orientation paper, this set is + * called $\Delta_k$, thus the name. + * + * In 2d, this set can only include zero, one, or two elements. + * Consequently, the appropriate data structure is one in which + * we store at most 2 elements in a fixed sized data structure. + */ + template <> + class EdgeDeltaSet<2> + { + public: + /** + * Iterator type for the elements of the set. + */ + using const_iterator = const unsigned int *; + + /** + * Default constructor. Initialize both slots as unused, corresponding + * to an empty set. + */ + EdgeDeltaSet() + { + edge_indices[0] = edge_indices[1] = numbers::invalid_unsigned_int; + } + + + /** + * Delete the elements of the set by marking both slots as unused. + */ + void + clear() + { + edge_indices[0] = edge_indices[1] = numbers::invalid_unsigned_int; + } + + /** + * Insert one element into the set. This will fail if the set already + * has two elements. + */ + void + insert(const unsigned int edge_index) + { + if (edge_indices[0] == numbers::invalid_unsigned_int) + edge_indices[0] = edge_index; + else + { + Assert(edge_indices[1] == numbers::invalid_unsigned_int, + ExcInternalError()); + edge_indices[1] = edge_index; + } + } + + + /** + * Return an iterator pointing to the first element of the set. + */ + const_iterator + begin() const + { + return edge_indices; + } + + + /** + * Return an iterator pointing to the element past the last used one. + */ + const_iterator + end() const + { + // check whether the current object stores zero, one, or two + // indices, and use this to point to the element past the + // last valid one + if (edge_indices[0] == numbers::invalid_unsigned_int) + return edge_indices; + else if (edge_indices[1] == numbers::invalid_unsigned_int) + return edge_indices + 1; + else + return edge_indices + 2; + } + + private: + /** + * Storage space to store the indices of at most two edges. + */ + unsigned int edge_indices[2]; + }; + + + + /** + * A class that represents by how much the set of parallel edges + * grows in each step. In the graph orientation paper, this set is + * called $\Delta_k$, thus the name. + * + * In 3d, this set can have arbitrarily many elements, unlike the + * 2d case specialized above. Consequently, we simply represent + * the data structure with a std::set. Class derivation ensures + * that we simply inherit all of the member functions of the + * base class. + */ + template <> + class EdgeDeltaSet<3> : public std::set + {}; + + + + /** + * From a list of cells, build a sorted vector that contains all of the + * edges that exist in the mesh. + */ + template + std::vector> + build_edges(const std::vector> &cells) + { + // build the edge list for all cells. because each cell has + // GeometryInfo::lines_per_cell edges, the total number + // of edges is this many times the number of cells. of course + // some of them will be duplicates, and we throw them out below + std::vector> edge_list; + edge_list.reserve(cells.size() * GeometryInfo::lines_per_cell); + for (unsigned int i = 0; i < cells.size(); ++i) + for (unsigned int l = 0; l < GeometryInfo::lines_per_cell; ++l) + edge_list.emplace_back(cells[i], l); + + // next sort the edge list and then remove duplicates + std::sort(edge_list.begin(), edge_list.end()); + edge_list.erase(std::unique(edge_list.begin(), edge_list.end()), + edge_list.end()); + + return edge_list; + } + + + + /** + * Build the cell list. Update the edge array to let edges know + * which cells are adjacent to them. + */ + template + std::vector> + build_cells_and_connect_edges(const std::vector> &cells, + std::vector> &edges) + { + std::vector> cell_list; + cell_list.reserve(cells.size()); + for (unsigned int i = 0; i < cells.size(); ++i) + { + // create our own data structure for the cells and let it + // connect to the edges array + cell_list.emplace_back(cells[i], edges); + + // then also inform the edges that they are adjacent + // to the current cell, and where within this cell + for (unsigned int l = 0; l < GeometryInfo::lines_per_cell; ++l) + edges[cell_list.back().edge_indices[l]].adjacent_cells.push_back( + AdjacentCell(i, l)); + } + Assert(cell_list.size() == cells.size(), ExcInternalError()); + + return cell_list; + } + + + + /** + * Return the index within 'cells' of the first cell that has at least one + * edge that is not yet oriented. + */ + template + unsigned int + get_next_unoriented_cell(const std::vector> &cells, + const std::vector> &edges, + const unsigned int current_cell) + { + for (unsigned int c = current_cell; c < cells.size(); ++c) + for (unsigned int l = 0; l < GeometryInfo::lines_per_cell; ++l) + if (edges[cells[c].edge_indices[l]].orientation_status == + Edge::not_oriented) + return c; + + return numbers::invalid_unsigned_int; + } + + + + /** + * Given a set of cells and edges, orient all edges that are + * (global) parallel to the one identified by the @p cell and + * within it the one with index @p local_edge. + */ + template + void + orient_one_set_of_parallel_edges(const std::vector> &cells, + std::vector> &edges, + const unsigned int cell, + const unsigned int local_edge) + { + // choose the direction of the first edge. we have free choice + // here and could simply choose "forward" if that's what pleases + // us. however, for backward compatibility with the previous + // implementation used till 2016, let us just choose the + // direction so that it matches what we have in the given cell. + // + // in fact, in what can only be assumed to be a bug in the + // original implementation, after orienting all edges, the code + // that rotates the cells so that they match edge orientations + // (see the rotate_cell() function below) rotated the cell two + // more times by 90 degrees. this is ok -- it simply flips all + // edge orientations, which leaves them valid. rather than do + // the same in the current implementation, we can achieve the + // same effect by modifying the rule above to choose the + // direction of the starting edge of this parallel set + // *opposite* to what it looks like in the current cell + // + // this bug only existed in the 2d implementation since there + // were different implementations for 2d and 3d. consequently, + // only replicate it for the 2d case and be "intuitive" in 3d. + if (edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] == + cells[cell].vertex_indices[GeometryInfo::line_to_cell_vertices( + local_edge, 0)]) + // orient initial edge *opposite* to the way it is in the cell + // (see above for the reason) + edges[cells[cell].edge_indices[local_edge]].orientation_status = + (dim == 2 ? Edge::backward : Edge::forward); + else + { + Assert( + edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] == + cells[cell].vertex_indices + [GeometryInfo::line_to_cell_vertices(local_edge, 1)], + ExcInternalError()); + Assert( + edges[cells[cell].edge_indices[local_edge]].vertex_indices[1] == + cells[cell].vertex_indices + [GeometryInfo::line_to_cell_vertices(local_edge, 0)], + ExcInternalError()); + + // orient initial edge *opposite* to the way it is in the cell + // (see above for the reason) + edges[cells[cell].edge_indices[local_edge]].orientation_status = + (dim == 2 ? Edge::forward : Edge::backward); + } + + // walk outward from the given edge as described in + // the algorithm in the paper that documents all of + // this + // + // note that in 2d, each of the Deltas can at most + // contain two elements, whereas in 3d it can be arbitrarily many + EdgeDeltaSet Delta_k; + EdgeDeltaSet Delta_k_minus_1; + Delta_k_minus_1.insert(cells[cell].edge_indices[local_edge]); + + while (Delta_k_minus_1.begin() != + Delta_k_minus_1.end()) // while set is not empty + { + Delta_k.clear(); + + for (typename EdgeDeltaSet::const_iterator delta = + Delta_k_minus_1.begin(); + delta != Delta_k_minus_1.end(); + ++delta) + { + Assert(edges[*delta].orientation_status != + Edge::not_oriented, + ExcInternalError()); + + // now go through the cells adjacent to this edge + for (typename AdjacentCells::const_iterator adjacent_cell = + edges[*delta].adjacent_cells.begin(); + adjacent_cell != edges[*delta].adjacent_cells.end(); + ++adjacent_cell) + { + const unsigned int K = adjacent_cell->cell_index; + const unsigned int delta_is_edge_in_K = + adjacent_cell->edge_within_cell; + + // figure out the direction of delta with respect to the cell + // K (in the orientation in which the user has given it to us) + const unsigned int first_edge_vertex = + (edges[*delta].orientation_status == Edge::forward ? + edges[*delta].vertex_indices[0] : + edges[*delta].vertex_indices[1]); + const unsigned int first_edge_vertex_in_K = + cells[K] + .vertex_indices[GeometryInfo::line_to_cell_vertices( + delta_is_edge_in_K, 0)]; + Assert( + first_edge_vertex == first_edge_vertex_in_K || + first_edge_vertex == + cells[K].vertex_indices[GeometryInfo< + dim>::line_to_cell_vertices(delta_is_edge_in_K, 1)], + ExcInternalError()); + + // now figure out which direction the each of the "opposite" + // edges needs to be oriented into. + for (unsigned int o_e = 0; + o_e < ParallelEdges::n_other_parallel_edges; + ++o_e) + { + // get the index of the opposite edge and select which its + // first vertex needs to be based on how the current edge + // is oriented in the current cell + const unsigned int opposite_edge = + cells[K].edge_indices[ParallelEdges< + dim>::parallel_edges[delta_is_edge_in_K][o_e]]; + const unsigned int first_opposite_edge_vertex = + cells[K].vertex_indices + [GeometryInfo::line_to_cell_vertices( + ParallelEdges< + dim>::parallel_edges[delta_is_edge_in_K][o_e], + (first_edge_vertex == first_edge_vertex_in_K ? 0 : + 1))]; + + // then determine the orientation of the edge based on + // whether the vertex we want to be the edge's first + // vertex is already the first vertex of the edge, or + // whether it points in the opposite direction + const typename Edge::OrientationStatus + opposite_edge_orientation = + (edges[opposite_edge].vertex_indices[0] == + first_opposite_edge_vertex ? + Edge::forward : + Edge::backward); + + // see if the opposite edge (there is only one in 2d) has + // already been oriented. + if (edges[opposite_edge].orientation_status == + Edge::not_oriented) + { + // the opposite edge is not yet oriented. do orient it + // and add it to Delta_k + edges[opposite_edge].orientation_status = + opposite_edge_orientation; + Delta_k.insert(opposite_edge); + } + else + { + // this opposite edge has already been oriented. it + // should be consistent with the current one in 2d, + // while in 3d it may in fact be mis-oriented, and in + // that case the mesh will not be orientable. indicate + // this by throwing an exception that we can catch + // further up; this has the advantage that we can + // propagate through a couple of functions without + // having to do error checking and without modifying + // the 'cells' array that the user gave us + if (dim == 2) + { + Assert(edges[opposite_edge].orientation_status == + opposite_edge_orientation, + ExcMeshNotOrientable()); + } + else if (dim == 3) + { + if (edges[opposite_edge].orientation_status != + opposite_edge_orientation) + throw ExcMeshNotOrientable(); + } + else + Assert(false, ExcNotImplemented()); + } + } + } + } + + // finally copy the new set to the previous one + // (corresponding to increasing 'k' by one in the + // algorithm) + Delta_k_minus_1 = Delta_k; + } + } + + + /** + * Given data structures @p cell_list and @p edge_list, where + * all edges are already oriented, rotate the cell with + * index @p cell_index in such a way that its local coordinate + * system matches the ones of the adjacent edges. Store the + * rotated order of vertices in raw_cells[cell_index]. + */ + template + void + rotate_cell(const std::vector> &cell_list, + const std::vector> &edge_list, + const unsigned int cell_index, + std::vector> &raw_cells) + { + // find the first vertex of the cell. this is the vertex where dim edges + // originate, so for each of the edges record which the starting vertex is + unsigned int starting_vertex_of_edge[GeometryInfo::lines_per_cell]; + for (unsigned int e = 0; e < GeometryInfo::lines_per_cell; ++e) + { + Assert(edge_list[cell_list[cell_index].edge_indices[e]] + .orientation_status != Edge::not_oriented, + ExcInternalError()); + if (edge_list[cell_list[cell_index].edge_indices[e]] + .orientation_status == Edge::forward) + starting_vertex_of_edge[e] = + edge_list[cell_list[cell_index].edge_indices[e]] + .vertex_indices[0]; + else + starting_vertex_of_edge[e] = + edge_list[cell_list[cell_index].edge_indices[e]] + .vertex_indices[1]; + } + + // find the vertex number that appears dim times. this will then be + // the vertex at which we want to locate the origin of the cell's + // coordinate system (i.e., vertex 0) + unsigned int origin_vertex_of_cell = numbers::invalid_unsigned_int; + switch (dim) + { + case 2: + { + // in 2d, we can simply enumerate the possibilities where the + // origin may be located because edges zero and one don't share + // any vertices, and the same for edges two and three + if ((starting_vertex_of_edge[0] == starting_vertex_of_edge[2]) || + (starting_vertex_of_edge[0] == starting_vertex_of_edge[3])) + origin_vertex_of_cell = starting_vertex_of_edge[0]; + else if ((starting_vertex_of_edge[1] == + starting_vertex_of_edge[2]) || + (starting_vertex_of_edge[1] == + starting_vertex_of_edge[3])) + origin_vertex_of_cell = starting_vertex_of_edge[1]; + else + Assert(false, ExcInternalError()); + + break; + } + + case 3: + { + // one could probably do something similar in 3d, but that seems + // more complicated than one wants to write down. just go + // through the list of possible starting vertices and check + for (origin_vertex_of_cell = 0; + origin_vertex_of_cell < GeometryInfo::vertices_per_cell; + ++origin_vertex_of_cell) + if (std::count(starting_vertex_of_edge, + starting_vertex_of_edge + + GeometryInfo::lines_per_cell, + cell_list[cell_index] + .vertex_indices[origin_vertex_of_cell]) == dim) + break; + Assert(origin_vertex_of_cell < + GeometryInfo::vertices_per_cell, + ExcInternalError()); + + break; + } + + default: + Assert(false, ExcNotImplemented()); + } + + // now rotate raw_cells[cell_index] in such a way that its orientation + // matches that of cell_list[cell_index] + switch (dim) + { + case 2: + { + // in 2d, we can literally rotate the cell until its origin + // matches the one that we have determined above should be + // the origin vertex + // + // when doing a rotation, take into account the ordering of + // vertices (not in clockwise or counter-clockwise sense) + while (raw_cells[cell_index].vertices[0] != origin_vertex_of_cell) + { + const unsigned int tmp = raw_cells[cell_index].vertices[0]; + raw_cells[cell_index].vertices[0] = + raw_cells[cell_index].vertices[1]; + raw_cells[cell_index].vertices[1] = + raw_cells[cell_index].vertices[3]; + raw_cells[cell_index].vertices[3] = + raw_cells[cell_index].vertices[2]; + raw_cells[cell_index].vertices[2] = tmp; + } + break; + } + + case 3: + { + // in 3d, the situation is a bit more complicated. from above, we + // now know which vertex is at the origin (because 3 edges + // originate from it), but that still leaves 3 possible rotations + // of the cube. the important realization is that we can choose + // any of them: in all 3 rotations, all edges originate from the + // one vertex, and that fixes the directions of all 12 edges in + // the cube because these 3 cover all 3 equivalence classes! + // consequently, we can select an arbitrary one among the + // permutations -- for example the following ones: + static const unsigned int cube_permutations[8][8] = { + {0, 1, 2, 3, 4, 5, 6, 7}, + {1, 5, 3, 7, 0, 4, 2, 6}, + {2, 6, 0, 4, 3, 7, 1, 5}, + {3, 2, 1, 0, 7, 6, 5, 4}, + {4, 0, 6, 2, 5, 1, 7, 3}, + {5, 4, 7, 6, 1, 0, 3, 2}, + {6, 7, 4, 5, 2, 3, 0, 1}, + {7, 3, 5, 1, 6, 2, 4, 0}}; + + unsigned int + temp_vertex_indices[GeometryInfo::vertices_per_cell]; + for (const unsigned int v : GeometryInfo::vertex_indices()) + temp_vertex_indices[v] = + raw_cells[cell_index] + .vertices[cube_permutations[origin_vertex_of_cell][v]]; + for (const unsigned int v : GeometryInfo::vertex_indices()) + raw_cells[cell_index].vertices[v] = temp_vertex_indices[v]; + + break; + } + + default: + { + Assert(false, ExcNotImplemented()); + } + } + } + + + /** + * Given a set of cells, find globally unique edge orientations + * and then rotate cells so that the coordinate system of the cell + * coincides with the coordinate systems of the adjacent edges. + */ + template + void + reorient(std::vector> &cells) + { + // first build the arrays that connect cells to edges and the other + // way around + std::vector> edge_list = build_edges(cells); + std::vector> cell_list = + build_cells_and_connect_edges(cells, edge_list); + + // then loop over all cells and start orienting parallel edge sets + // of cells that still have non-oriented edges + unsigned int next_cell_with_unoriented_edge = 0; + while ((next_cell_with_unoriented_edge = get_next_unoriented_cell( + cell_list, edge_list, next_cell_with_unoriented_edge)) != + numbers::invalid_unsigned_int) + { + // see which edge sets are still not oriented + // + // we do not need to look at each edge because if we orient edge + // 0, we will end up with edge 1 also oriented (in 2d; in 3d, there + // will be 3 other edges that are also oriented). there are only + // dim independent sets of edges, so loop over these. + // + // we need to check whether each one of these starter edges may + // already be oriented because the line (sheet) that connects + // globally parallel edges may be self-intersecting in the + // current cell + for (unsigned int l = 0; l < dim; ++l) + if (edge_list[cell_list[next_cell_with_unoriented_edge] + .edge_indices[ParallelEdges::starter_edges[l]]] + .orientation_status == Edge::not_oriented) + orient_one_set_of_parallel_edges( + cell_list, + edge_list, + next_cell_with_unoriented_edge, + ParallelEdges::starter_edges[l]); + + // ensure that we have really oriented all edges now, not just + // the starter edges + for (unsigned int l = 0; l < GeometryInfo::lines_per_cell; ++l) + Assert(edge_list[cell_list[next_cell_with_unoriented_edge] + .edge_indices[l]] + .orientation_status != Edge::not_oriented, + ExcInternalError()); + } + + // now that we have oriented all edges, we need to rotate cells + // so that the edges point in the right direction with the now + // rotated coordinate system + for (unsigned int c = 0; c < cells.size(); ++c) + rotate_cell(cell_list, edge_list, c, cells); + } + + + // overload of the function above for 1d -- there is nothing + // to orient in that case + void + reorient(std::vector> &) + {} + } // namespace + + + + template + void + consistently_order_cells(std::vector> &cells) + { + Assert(cells.size() != 0, + ExcMessage( + "List of elements to orient must have at least one cell")); + + // there is nothing for us to do in 1d + if (dim == 1) + return; + + // check if grids are already consistent. if so, do + // nothing. if not, then do the reordering + if (!is_consistent(cells)) + try + { + reorient(cells); + } + catch (const ExcMeshNotOrientable &) + { + // the mesh is not orientable. this is acceptable if we are in 3d, + // as class Triangulation knows how to handle this, but it is + // not in 2d; in that case, re-throw the exception + if (dim < 3) + throw; + } + } + + + + template + std::map> + get_all_vertices_at_boundary(const Triangulation &tria) + { + std::map> vertex_map; + typename Triangulation::active_cell_iterator + cell = tria.begin_active(), + endc = tria.end(); + for (; cell != endc; ++cell) + { + for (const unsigned int i : cell->face_indices()) + { + const typename Triangulation::face_iterator &face = + cell->face(i); + if (face->at_boundary()) + { + for (unsigned j = 0; j < face->n_vertices(); ++j) + { + const Point &vertex = face->vertex(j); + const unsigned int vertex_index = face->vertex_index(j); + vertex_map[vertex_index] = vertex; + } + } + } + } + return vertex_map; + } + + + + template + void + remove_hanging_nodes(Triangulation &tria, + const bool isotropic, + const unsigned int max_iterations) + { + unsigned int iter = 0; + bool continue_refinement = true; + + while (continue_refinement && (iter < max_iterations)) + { + if (max_iterations != numbers::invalid_unsigned_int) + iter++; + continue_refinement = false; + + for (const auto &cell : tria.active_cell_iterators()) + for (const unsigned int j : cell->face_indices()) + if (cell->at_boundary(j) == false && + cell->neighbor(j)->has_children()) + { + if (isotropic) + { + cell->set_refine_flag(); + continue_refinement = true; + } + else + continue_refinement |= cell->flag_for_face_refinement(j); + } + + tria.execute_coarsening_and_refinement(); + } + } + + + + template + void + remove_anisotropy(Triangulation &tria, + const double max_ratio, + const unsigned int max_iterations) + { + unsigned int iter = 0; + bool continue_refinement = true; + + while (continue_refinement && (iter < max_iterations)) + { + iter++; + continue_refinement = false; + for (const auto &cell : tria.active_cell_iterators()) + { + std::pair info = + GridTools::get_longest_direction(cell); + if (info.second > max_ratio) + { + cell->set_refine_flag( + RefinementCase::cut_axis(info.first)); + continue_refinement = true; + } + } + tria.execute_coarsening_and_refinement(); + } + } + + + + template + std::map> + extract_used_vertices(const Triangulation &container, + const Mapping &mapping) + { + std::map> result; + for (const auto &cell : container.active_cell_iterators()) + { + if (!cell->is_artificial()) + { + const auto vs = mapping.get_vertices(cell); + for (unsigned int i = 0; i < vs.size(); ++i) + result[cell->vertex_index(i)] = vs[i]; + } + } + return result; + } + + + + template + std::vector< + std::set::active_cell_iterator>> + vertex_to_cell_map(const Triangulation &triangulation) + { + std::vector< + std::set::active_cell_iterator>> + vertex_to_cell_map(triangulation.n_vertices()); + typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + for (; cell != endc; ++cell) + for (const unsigned int i : cell->vertex_indices()) + vertex_to_cell_map[cell->vertex_index(i)].insert(cell); + + // Check if mesh has hanging nodes. Do this only locally to + // prevent communication and possible deadlock. + if (triangulation.Triangulation::has_hanging_nodes()) + { + Assert(triangulation.all_reference_cells_are_hyper_cube(), + ExcNotImplemented()); + + // Take care of hanging nodes + cell = triangulation.begin_active(); + for (; cell != endc; ++cell) + { + for (const unsigned int i : cell->face_indices()) + { + if ((cell->at_boundary(i) == false) && + (cell->neighbor(i)->is_active())) + { + typename Triangulation::active_cell_iterator + adjacent_cell = cell->neighbor(i); + for (unsigned int j = 0; j < cell->face(i)->n_vertices(); + ++j) + vertex_to_cell_map[cell->face(i)->vertex_index(j)].insert( + adjacent_cell); + } + } + + // in 3d also loop over the edges + if (dim == 3) + { + for (unsigned int i = 0; i < cell->n_lines(); ++i) + if (cell->line(i)->has_children()) + // the only place where this vertex could have been + // hiding is on the mid-edge point of the edge we + // are looking at + vertex_to_cell_map[cell->line(i)->child(0)->vertex_index(1)] + .insert(cell); + } + } + } + + return vertex_to_cell_map; + } + + + + template + void + get_face_connectivity_of_cells( + const Triangulation &triangulation, + DynamicSparsityPattern &cell_connectivity) + { + cell_connectivity.reinit(triangulation.n_active_cells(), + triangulation.n_active_cells()); + + // loop over all cells and their neighbors to build the sparsity + // pattern. note that it's a bit hard to enter all the connections when a + // neighbor has children since we would need to find out which of its + // children is adjacent to the current cell. this problem can be omitted + // if we only do something if the neighbor has no children -- in that case + // it is either on the same or a coarser level than we are. in return, we + // have to add entries in both directions for both cells + for (const auto &cell : triangulation.active_cell_iterators()) + { + const unsigned int index = cell->active_cell_index(); + cell_connectivity.add(index, index); + for (auto f : cell->face_indices()) + if ((cell->at_boundary(f) == false) && + (cell->neighbor(f)->has_children() == false)) + { + const unsigned int other_index = + cell->neighbor(f)->active_cell_index(); + cell_connectivity.add(index, other_index); + cell_connectivity.add(other_index, index); + } + } + } + + + + template + void + get_vertex_connectivity_of_cells( + const Triangulation &triangulation, + DynamicSparsityPattern &cell_connectivity) + { + std::vector> vertex_to_cell( + triangulation.n_vertices()); + for (const auto &cell : triangulation.active_cell_iterators()) + { + for (const unsigned int v : cell->vertex_indices()) + vertex_to_cell[cell->vertex_index(v)].push_back( + cell->active_cell_index()); + } + + cell_connectivity.reinit(triangulation.n_active_cells(), + triangulation.n_active_cells()); + for (const auto &cell : triangulation.active_cell_iterators()) + { + for (const unsigned int v : cell->vertex_indices()) + for (unsigned int n = 0; + n < vertex_to_cell[cell->vertex_index(v)].size(); + ++n) + cell_connectivity.add(cell->active_cell_index(), + vertex_to_cell[cell->vertex_index(v)][n]); + } + } + + + template + void + get_vertex_connectivity_of_cells_on_level( + const Triangulation &triangulation, + const unsigned int level, + DynamicSparsityPattern &cell_connectivity) + { + std::vector> vertex_to_cell( + triangulation.n_vertices()); + for (typename Triangulation::cell_iterator cell = + triangulation.begin(level); + cell != triangulation.end(level); + ++cell) + { + for (const unsigned int v : cell->vertex_indices()) + vertex_to_cell[cell->vertex_index(v)].push_back(cell->index()); + } + + cell_connectivity.reinit(triangulation.n_cells(level), + triangulation.n_cells(level)); + for (typename Triangulation::cell_iterator cell = + triangulation.begin(level); + cell != triangulation.end(level); + ++cell) + { + for (const unsigned int v : cell->vertex_indices()) + for (unsigned int n = 0; + n < vertex_to_cell[cell->vertex_index(v)].size(); + ++n) + cell_connectivity.add(cell->index(), + vertex_to_cell[cell->vertex_index(v)][n]); + } + } +} /* namespace GridTools */ + +// explicit instantiations +#include "grid_tools_topology.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/grid/grid_tools_topology.inst.in b/source/grid/grid_tools_topology.inst.in new file mode 100644 index 0000000000..ec66d3370f --- /dev/null +++ b/source/grid/grid_tools_topology.inst.in @@ -0,0 +1,122 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2001 - 2023 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. +// +// --------------------------------------------------------------------- + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + namespace GridTools + \{ + template std::tuple>, + std::vector>, + SubCellData> + get_coarse_mesh_description( + const Triangulation &tria); + + template void + delete_unused_vertices(std::vector> &, + std::vector> &, + SubCellData &); + + template void + delete_duplicated_vertices(std::vector> &, + std::vector> &, + SubCellData &, + std::vector &, + double); + +# if deal_II_dimension == deal_II_space_dimension + template void + delete_duplicated_vertices(std::vector> &, + const double); +# endif + + template void + invert_all_negative_measure_cells( + const std::vector> &, + std::vector> &); + + template std::size_t + invert_cells_with_negative_measure( + const std::vector> &, + std::vector> &); + +# if deal_II_dimension == deal_II_space_dimension + template void + consistently_order_cells(std::vector> &); +# endif + + template std::map> + get_all_vertices_at_boundary( + const Triangulation &tria); + + template std::vector< + std::set::active_cell_iterator>> + vertex_to_cell_map( + const Triangulation + &triangulation); + + template std::map> + extract_used_vertices( + const Triangulation &mesh, + const Mapping &mapping); + + template void + get_face_connectivity_of_cells( + const Triangulation + &triangulation, + DynamicSparsityPattern &cell_connectivity); + + template void + get_vertex_connectivity_of_cells( + const Triangulation + &triangulation, + DynamicSparsityPattern &cell_connectivity); + + template void + get_vertex_connectivity_of_cells_on_level( + const Triangulation + &triangulation, + const unsigned int level, + DynamicSparsityPattern &cell_connectivity); + + \} +#endif + } + +for (deal_II_dimension : DIMENSIONS) + { + template void + GridTools::remove_hanging_nodes( + Triangulation & tria, bool, unsigned int); + + template void + GridTools::remove_anisotropy( + Triangulation &, double, unsigned int); + +#if deal_II_dimension < 3 + template void + GridTools::remove_hanging_nodes( + Triangulation & tria, + bool, + unsigned int); + + template void + GridTools::remove_anisotropy( + Triangulation &, + double, + unsigned int); +#endif + } -- 2.39.5