#include <deal.II/fe/mapping.h>
#include <deal.II/grid/grid_tools_geometry.h>
+#include <deal.II/grid/grid_tools_topology.h>
#include <deal.II/grid/manifold.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <boost/serialization/array.hpp>
#include <boost/serialization/vector.hpp>
-#include <optional>
-
#ifdef DEAL_II_WITH_ZLIB
# include <boost/iostreams/device/back_inserter.hpp>
# include <boost/iostreams/filter/gzip.hpp>
#endif
#include <bitset>
-#include <list>
+#include <optional>
#include <set>
#ifdef DEAL_II_HAVE_CXX20
*/
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 <tt>GridIn::read_*</tt> 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 <tt>GridIn::read_*</tt> 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 <int dim, int spacedim>
- void
- delete_unused_vertices(std::vector<Point<spacedim>> &vertices,
- std::vector<CellData<dim>> &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 <tt>GridIn::read_*</tt> 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 <int dim, int spacedim>
- void
- delete_duplicated_vertices(std::vector<Point<spacedim>> &all_vertices,
- std::vector<CellData<dim>> &cells,
- SubCellData &subcelldata,
- std::vector<unsigned int> &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 <int dim>
- void
- delete_duplicated_vertices(std::vector<Point<dim>> &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 <int dim, int spacedim>
- void
- invert_all_negative_measure_cells(
- const std::vector<Point<spacedim>> &all_vertices,
- std::vector<CellData<dim>> &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 <int dim, int spacedim>
- std::size_t
- invert_cells_with_negative_measure(
- const std::vector<Point<spacedim>> &all_vertices,
- std::vector<CellData<dim>> &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 <int dim>
- void
- consistently_order_cells(std::vector<CellData<dim>> &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 <int dim, int spacedim>
- std::
- tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
- get_coarse_mesh_description(const Triangulation<dim, spacedim> &tria);
/**
* Transform the vertices of the given triangulation by applying the
const Function<dim, double> *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 <int dim, int spacedim>
- std::map<unsigned int, Point<spacedim>>
- get_all_vertices_at_boundary(const Triangulation<dim, spacedim> &tria);
-
/**
* Scale the entire triangulation by the given factor. To preserve the
* orientation of the triangulation, the factor must be positive.
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 <int dim, int spacedim>
- void
- remove_hanging_nodes(Triangulation<dim, spacedim> &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 <int dim, int spacedim>
- void
- remove_anisotropy(Triangulation<dim, spacedim> &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
} // namespace internal
- /**
- * Return a map `vertex index -> Point<spacedim>` 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 <int dim, int spacedim>
- std::map<unsigned int, Point<spacedim>>
- extract_used_vertices(
- const Triangulation<dim, spacedim> &container,
- const Mapping<dim, spacedim> &mapping =
- (ReferenceCells::get_hypercube<dim>()
-#ifndef _MSC_VER
- .template get_default_linear_mapping<dim, spacedim>()
-#else
- .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
-#endif
- ));
-
/**
* Find and return the index of the closest vertex to a given point in the
* map of vertices passed as the first argument.
const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
const std::vector<Point<spacedim>> &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 <code>cell-@>vertex_index()</code>. Notice that only the indices
- * marked in the array returned by
- * Triangulation<dim,spacedim>::get_used_vertices() are used.
- */
- template <int dim, int spacedim>
- std::vector<
- std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
- vertex_to_cell_map(const Triangulation<dim, spacedim> &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
*/
/** @{ */
- /**
- * 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 <int dim, int spacedim>
- void
- get_face_connectivity_of_cells(
- const Triangulation<dim, spacedim> &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 <int dim, int spacedim>
- void
- get_vertex_connectivity_of_cells(
- const Triangulation<dim, spacedim> &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 <int dim, int spacedim>
- void
- get_vertex_connectivity_of_cells_on_level(
- const Triangulation<dim, spacedim> &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
std::vector<bool>
get_locally_owned_vertices(const Triangulation<dim, spacedim> &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<MeshType>}
- */
- template <typename MeshType>
- DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
- std::list<std::pair<
- typename MeshType::cell_iterator,
- typename MeshType::cell_iterator>> 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 <int dim, int spacedim>
- bool
- have_same_coarse_mesh(const Triangulation<dim, spacedim> &mesh_1,
- const Triangulation<dim, spacedim> &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<MeshType>}
- */
- template <typename MeshType>
- DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
- bool have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2);
-
/** @} */
/**
* @name Dealing with distorted cells
<< "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*/
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/config.h>
+
+#include <deal.II/base/point.h>
+#include <deal.II/base/template_constraints.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_description.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+
+#include <list>
+#include <tuple>
+#include <vector>
+
+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 <int dim, int spacedim>
+ std::
+ tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
+ get_coarse_mesh_description(const Triangulation<dim, spacedim> &tria);
+
+ /**
+ * Remove vertices that are not referenced by any of the cells. This
+ * function is called by all <tt>GridIn::read_*</tt> 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 <tt>GridIn::read_*</tt> 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 <int dim, int spacedim>
+ void
+ delete_unused_vertices(std::vector<Point<spacedim>> &vertices,
+ std::vector<CellData<dim>> &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 <tt>GridIn::read_*</tt> 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 <int dim, int spacedim>
+ void
+ delete_duplicated_vertices(std::vector<Point<spacedim>> &all_vertices,
+ std::vector<CellData<dim>> &cells,
+ SubCellData &subcelldata,
+ std::vector<unsigned int> &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 <int dim>
+ void
+ delete_duplicated_vertices(std::vector<Point<dim>> &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 <int dim, int spacedim>
+ void
+ invert_all_negative_measure_cells(
+ const std::vector<Point<spacedim>> &all_vertices,
+ std::vector<CellData<dim>> &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 <int dim, int spacedim>
+ std::size_t
+ invert_cells_with_negative_measure(
+ const std::vector<Point<spacedim>> &all_vertices,
+ std::vector<CellData<dim>> &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 <int dim>
+ void
+ consistently_order_cells(std::vector<CellData<dim>> &cells);
+
+ /**
+ * Return a std::map with all vertices of faces located in the boundary
+ *
+ * @param[in] tria The Triangulation object.
+ */
+ template <int dim, int spacedim>
+ std::map<unsigned int, Point<spacedim>>
+ get_all_vertices_at_boundary(const Triangulation<dim, spacedim> &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 <int dim, int spacedim>
+ void
+ remove_hanging_nodes(Triangulation<dim, spacedim> &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 <int dim, int spacedim>
+ void
+ remove_anisotropy(Triangulation<dim, spacedim> &tria,
+ const double max_ratio = 1.6180339887,
+ const unsigned int max_iterations = 5);
+
+ /**
+ * Return a map `vertex index -> Point<spacedim>` 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 <int dim, int spacedim>
+ std::map<unsigned int, Point<spacedim>>
+ extract_used_vertices(
+ const Triangulation<dim, spacedim> &container,
+ const Mapping<dim, spacedim> &mapping =
+ (ReferenceCells::get_hypercube<dim>()
+#ifndef _MSC_VER
+ .template get_default_linear_mapping<dim, spacedim>()
+#else
+ .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
+#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 <code>cell-@>vertex_index()</code>. Notice that only the indices
+ * marked in the array returned by
+ * Triangulation<dim,spacedim>::get_used_vertices() are used.
+ */
+ template <int dim, int spacedim>
+ std::vector<
+ std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
+ vertex_to_cell_map(const Triangulation<dim, spacedim> &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 <int dim, int spacedim>
+ void
+ get_face_connectivity_of_cells(
+ const Triangulation<dim, spacedim> &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 <int dim, int spacedim>
+ void
+ get_vertex_connectivity_of_cells(
+ const Triangulation<dim, spacedim> &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 <int dim, int spacedim>
+ void
+ get_vertex_connectivity_of_cells_on_level(
+ const Triangulation<dim, spacedim> &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<MeshType>}
+ */
+ template <typename MeshType>
+ DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
+ std::list<std::pair<
+ typename MeshType::cell_iterator,
+ typename MeshType::cell_iterator>> 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 <int dim, int spacedim>
+ bool
+ have_same_coarse_mesh(const Triangulation<dim, spacedim> &mesh_1,
+ const Triangulation<dim, spacedim> &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<MeshType>}
+ */
+ template <typename MeshType>
+ DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
+ 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
grid_tools_dof_handlers.cc
grid_tools_geometry.cc
grid_tools_nontemplates.cc
+ grid_tools_topology.cc
tria.cc
)
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
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 <int structdim>
- struct CellDataComparator
- {
- bool
- operator()(const CellData<structdim> &a,
- const CellData<structdim> &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<dim-1>, 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 <int dim>
- class FaceDataHelper
- {
- public:
- /**
- * Store the data about the given face @p face.
- */
- template <typename FaceIteratorType>
- void
- insert_face_data(const FaceIteratorType &face)
- {
- CellData<dim - 1> 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<dim - 1> &face_cell_data : face_data)
- internal::append_face_data(face_cell_data, subcell_data);
- return subcell_data;
- }
-
-
- private:
- std::set<CellData<dim - 1>, internal::CellDataComparator<dim - 1>>
- face_data;
- };
-
-
- // Do nothing for dim=1:
- template <>
- class FaceDataHelper<1>
- {
- public:
- template <typename FaceIteratorType>
- void
- insert_face_data(const FaceIteratorType &)
- {}
-
- SubCellData
- get()
- {
- return SubCellData();
- }
- };
- } // namespace internal
-
-
-
- template <int dim, int spacedim>
- std::
- tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
- get_coarse_mesh_description(const Triangulation<dim, spacedim> &tria)
- {
- Assert(tria.n_levels() >= 1,
- ExcMessage("The input triangulation must be non-empty."));
-
- std::vector<Point<spacedim>> vertices = tria.get_vertices();
- std::vector<CellData<dim>> cells;
-
- internal::FaceDataHelper<dim> face_data;
- std::set<CellData<1>, internal::CellDataComparator<1>>
- line_data; // only used in 3d
-
- for (const auto &cell : tria.cell_iterators_on_level(0))
- {
- // Save cell data
- CellData<dim> 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<Point<spacedim>>,
- std::vector<CellData<dim>>,
- SubCellData>(std::move(vertices),
- std::move(cells),
- std::move(subcell_data));
- }
-
-
-
- template <int dim, int spacedim>
- void
- delete_unused_vertices(std::vector<Point<spacedim>> &vertices,
- std::vector<CellData<dim>> &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<bool> 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<unsigned int> 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<Point<spacedim>> 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 <int dim, int spacedim>
- void
- delete_duplicated_vertices(std::vector<Point<spacedim>> &vertices,
- std::vector<CellData<dim>> &cells,
- SubCellData &subcelldata,
- std::vector<unsigned int> &considered_vertices,
- const double tol)
- {
- if (tol == 0.0)
- return; // nothing to do per definition
-
- AssertIndexRange(2, vertices.size());
- std::vector<unsigned int> 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<spacedim> 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<std::pair<unsigned int, Point<spacedim>>> 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<unsigned int, Point<spacedim>> &a,
- const std::pair<unsigned int, Point<spacedim>> &b) {
- return a.second[longest_coordinate_direction] <
- b.second[longest_coordinate_direction];
- });
-
- auto within_tolerance = [=](const Point<spacedim> &a,
- const Point<spacedim> &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<unsigned int, Point<spacedim>> &a,
- const std::pair<unsigned int, Point<spacedim>> &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 <int dim>
- void
- delete_duplicated_vertices(std::vector<Point<dim>> &vertices,
- const double tol)
- {
- if (vertices.empty())
- return;
-
- // 1) map point to local vertex index
- std::map<Point<dim>, unsigned int, FloatingPointComparator<double>>
- map_point_to_local_vertex_index{FloatingPointComparator<double>(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 <int dim, int spacedim>
- std::size_t
- invert_cells_with_negative_measure(
- const std::vector<Point<spacedim>> &all_vertices,
- std::vector<CellData<dim>> &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<const unsigned int> 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 <int dim, int spacedim>
- void
- invert_all_negative_measure_cells(
- const std::vector<Point<spacedim>> &all_vertices,
- std::vector<CellData<dim>> &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 <int dim>
- bool
- is_consistent(const std::vector<CellData<dim>> &cells)
- {
- std::set<CheapEdge> edges;
-
- for (typename std::vector<CellData<dim>>::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<dim>::lines_per_cell; ++l)
- {
- const CheapEdge reverse_edge(
- c->vertices[GeometryInfo<dim>::line_to_cell_vertices(l, 1)],
- c->vertices[GeometryInfo<dim>::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<dim>::line_to_cell_vertices(l, 0)],
- c->vertices[GeometryInfo<dim>::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 <int dim>
- 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<dim>::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 <int dim>
- 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<AdjacentCell>, 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<AdjacentCell>
- {};
-
-
- /**
- * 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 <int dim>
- 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<dim> &cell, const unsigned int edge_number)
- : orientation_status(not_oriented)
- {
- Assert(edge_number < GeometryInfo<dim>::lines_per_cell,
- ExcInternalError());
-
- // copy vertices for this particular line
- vertex_indices[0] =
- cell
- .vertices[GeometryInfo<dim>::line_to_cell_vertices(edge_number, 0)];
- vertex_indices[1] =
- cell
- .vertices[GeometryInfo<dim>::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<dim> &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<dim> &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<dim> adjacent_cells;
- };
-
-
-
- /**
- * A data structure that represents a cell with all of its vertices
- * and edges.
- */
- template <int dim>
- 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<dim> &c, const std::vector<Edge<dim>> &edge_list)
- {
- for (const unsigned int i : GeometryInfo<dim>::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<dim>::lines_per_cell; ++l)
- {
- const Edge<dim> 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<dim>::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<dim>::lines_per_cell];
- };
-
-
-
- template <int dim>
- 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<unsigned int>
- {};
-
-
-
- /**
- * From a list of cells, build a sorted vector that contains all of the
- * edges that exist in the mesh.
- */
- template <int dim>
- std::vector<Edge<dim>>
- build_edges(const std::vector<CellData<dim>> &cells)
- {
- // build the edge list for all cells. because each cell has
- // GeometryInfo<dim>::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<dim>> edge_list;
- edge_list.reserve(cells.size() * GeometryInfo<dim>::lines_per_cell);
- for (unsigned int i = 0; i < cells.size(); ++i)
- for (unsigned int l = 0; l < GeometryInfo<dim>::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 <int dim>
- std::vector<Cell<dim>>
- build_cells_and_connect_edges(const std::vector<CellData<dim>> &cells,
- std::vector<Edge<dim>> &edges)
- {
- std::vector<Cell<dim>> 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<dim>::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 <int dim>
- unsigned int
- get_next_unoriented_cell(const std::vector<Cell<dim>> &cells,
- const std::vector<Edge<dim>> &edges,
- const unsigned int current_cell)
- {
- for (unsigned int c = current_cell; c < cells.size(); ++c)
- for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
- if (edges[cells[c].edge_indices[l]].orientation_status ==
- Edge<dim>::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 <int dim>
- void
- orient_one_set_of_parallel_edges(const std::vector<Cell<dim>> &cells,
- std::vector<Edge<dim>> &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<dim>::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<dim>::backward : Edge<dim>::forward);
- else
- {
- Assert(
- edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] ==
- cells[cell].vertex_indices
- [GeometryInfo<dim>::line_to_cell_vertices(local_edge, 1)],
- ExcInternalError());
- Assert(
- edges[cells[cell].edge_indices[local_edge]].vertex_indices[1] ==
- cells[cell].vertex_indices
- [GeometryInfo<dim>::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<dim>::forward : Edge<dim>::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<dim> Delta_k;
- EdgeDeltaSet<dim> 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<dim>::const_iterator delta =
- Delta_k_minus_1.begin();
- delta != Delta_k_minus_1.end();
- ++delta)
- {
- Assert(edges[*delta].orientation_status !=
- Edge<dim>::not_oriented,
- ExcInternalError());
-
- // now go through the cells adjacent to this edge
- for (typename AdjacentCells<dim>::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<dim>::forward ?
- edges[*delta].vertex_indices[0] :
- edges[*delta].vertex_indices[1]);
- const unsigned int first_edge_vertex_in_K =
- cells[K]
- .vertex_indices[GeometryInfo<dim>::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<dim>::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<dim>::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<dim>::OrientationStatus
- opposite_edge_orientation =
- (edges[opposite_edge].vertex_indices[0] ==
- first_opposite_edge_vertex ?
- Edge<dim>::forward :
- Edge<dim>::backward);
-
- // see if the opposite edge (there is only one in 2d) has
- // already been oriented.
- if (edges[opposite_edge].orientation_status ==
- Edge<dim>::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 <code>raw_cells[cell_index]</code>.
- */
- template <int dim>
- void
- rotate_cell(const std::vector<Cell<dim>> &cell_list,
- const std::vector<Edge<dim>> &edge_list,
- const unsigned int cell_index,
- std::vector<CellData<dim>> &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<dim>::lines_per_cell];
- for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++e)
- {
- Assert(edge_list[cell_list[cell_index].edge_indices[e]]
- .orientation_status != Edge<dim>::not_oriented,
- ExcInternalError());
- if (edge_list[cell_list[cell_index].edge_indices[e]]
- .orientation_status == Edge<dim>::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<dim>::vertices_per_cell;
- ++origin_vertex_of_cell)
- if (std::count(starting_vertex_of_edge,
- starting_vertex_of_edge +
- GeometryInfo<dim>::lines_per_cell,
- cell_list[cell_index]
- .vertex_indices[origin_vertex_of_cell]) == dim)
- break;
- Assert(origin_vertex_of_cell <
- GeometryInfo<dim>::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<dim>::vertices_per_cell];
- for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
- temp_vertex_indices[v] =
- raw_cells[cell_index]
- .vertices[cube_permutations[origin_vertex_of_cell][v]];
- for (const unsigned int v : GeometryInfo<dim>::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 <int dim>
- void
- reorient(std::vector<CellData<dim>> &cells)
- {
- // first build the arrays that connect cells to edges and the other
- // way around
- std::vector<Edge<dim>> edge_list = build_edges(cells);
- std::vector<Cell<dim>> 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<dim>::starter_edges[l]]]
- .orientation_status == Edge<dim>::not_oriented)
- orient_one_set_of_parallel_edges(
- cell_list,
- edge_list,
- next_cell_with_unoriented_edge,
- ParallelEdges<dim>::starter_edges[l]);
-
- // ensure that we have really oriented all edges now, not just
- // the starter edges
- for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
- Assert(edge_list[cell_list[next_cell_with_unoriented_edge]
- .edge_indices[l]]
- .orientation_status != Edge<dim>::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<CellData<1>> &)
- {}
- } // namespace
-
- template <int dim>
- void
- consistently_order_cells(std::vector<CellData<dim>> &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
{
}
}
- template <int dim, int spacedim>
- std::map<unsigned int, Point<spacedim>>
- get_all_vertices_at_boundary(const Triangulation<dim, spacedim> &tria)
- {
- std::map<unsigned int, Point<spacedim>> vertex_map;
- typename Triangulation<dim, spacedim>::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<dim, spacedim>::face_iterator &face =
- cell->face(i);
- if (face->at_boundary())
- {
- for (unsigned j = 0; j < face->n_vertices(); ++j)
- {
- const Point<spacedim> &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.
}
- template <int dim, int spacedim>
- std::vector<
- std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
- vertex_to_cell_map(const Triangulation<dim, spacedim> &triangulation)
- {
- std::vector<
- std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
- vertex_to_cell_map(triangulation.n_vertices());
- typename Triangulation<dim, spacedim>::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<dim, spacedim>::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<dim, spacedim>::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 <int dim, int spacedim>
std::map<unsigned int, types::global_vertex_index>
- template <int dim, int spacedim>
- void
- get_face_connectivity_of_cells(
- const Triangulation<dim, spacedim> &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 <int dim, int spacedim>
- void
- get_vertex_connectivity_of_cells(
- const Triangulation<dim, spacedim> &triangulation,
- DynamicSparsityPattern &cell_connectivity)
- {
- std::vector<std::vector<unsigned int>> 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 <int dim, int spacedim>
- void
- get_vertex_connectivity_of_cells_on_level(
- const Triangulation<dim, spacedim> &triangulation,
- const unsigned int level,
- DynamicSparsityPattern &cell_connectivity)
- {
- std::vector<std::vector<unsigned int>> vertex_to_cell(
- triangulation.n_vertices());
- for (typename Triangulation<dim, spacedim>::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<dim, spacedim>::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 <int dim, int spacedim>
void
partition_triangulation(const unsigned int n_partitions,
- template <int dim, int spacedim>
- void
- remove_hanging_nodes(Triangulation<dim, spacedim> &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 <int dim, int spacedim>
- void
- remove_anisotropy(Triangulation<dim, spacedim> &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<unsigned int, double> info =
- GridTools::get_longest_direction<dim, spacedim>(cell);
- if (info.second > max_ratio)
- {
- cell->set_refine_flag(
- RefinementCase<dim>::cut_axis(info.first));
- continue_refinement = true;
- }
- }
- tria.execute_coarsening_and_refinement();
- }
- }
-
-
template <int dim, int spacedim>
void
regularize_corner_cells(Triangulation<dim, spacedim> &tria,
- template <int dim, int spacedim>
- std::map<unsigned int, Point<spacedim>>
- extract_used_vertices(const Triangulation<dim, spacedim> &container,
- const Mapping<dim, spacedim> &mapping)
- {
- std::map<unsigned int, Point<spacedim>> 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 <int spacedim>
unsigned int
find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
deal_II_space_dimension>
&triangulation);
- template std::map<unsigned int, Point<deal_II_space_dimension>>
- extract_used_vertices(
- const Triangulation<deal_II_dimension, deal_II_space_dimension> &mesh,
- const Mapping<deal_II_dimension, deal_II_space_dimension> &mapping);
-
template std::pair<
typename Triangulation<deal_II_dimension,
deal_II_space_dimension>::active_cell_iterator,
#if deal_II_dimension <= deal_II_space_dimension
namespace GridTools
\{
- template std::tuple<std::vector<Point<deal_II_space_dimension>>,
- std::vector<CellData<deal_II_dimension>>,
- SubCellData>
- get_coarse_mesh_description(
- const Triangulation<deal_II_dimension, deal_II_space_dimension> &tria);
-
- template void
- delete_unused_vertices(std::vector<Point<deal_II_space_dimension>> &,
- std::vector<CellData<deal_II_dimension>> &,
- SubCellData &);
-
- template void
- delete_duplicated_vertices(std::vector<Point<deal_II_space_dimension>> &,
- std::vector<CellData<deal_II_dimension>> &,
- SubCellData &,
- std::vector<unsigned int> &,
- double);
-
- template void
- invert_all_negative_measure_cells(
- const std::vector<Point<deal_II_space_dimension>> &,
- std::vector<CellData<deal_II_dimension>> &);
-
- template std::size_t
- invert_cells_with_negative_measure(
- const std::vector<Point<deal_II_space_dimension>> &,
- std::vector<CellData<deal_II_dimension>> &);
-
-# if deal_II_dimension == deal_II_space_dimension
- template void
- consistently_order_cells(std::vector<CellData<deal_II_dimension>> &);
-# endif
-
template void
shift<deal_II_dimension>(
const Tensor<1, deal_II_space_dimension> &,
const bool,
const unsigned int);
- template void
- get_face_connectivity_of_cells(
- const Triangulation<deal_II_dimension, deal_II_space_dimension>
- &triangulation,
- DynamicSparsityPattern &cell_connectivity);
-
- template void
- get_vertex_connectivity_of_cells(
- const Triangulation<deal_II_dimension, deal_II_space_dimension>
- &triangulation,
- DynamicSparsityPattern &cell_connectivity);
-
- template void
- get_vertex_connectivity_of_cells_on_level(
- const Triangulation<deal_II_dimension, deal_II_space_dimension>
- &triangulation,
- const unsigned int level,
- DynamicSparsityPattern &cell_connectivity);
-
template void
partition_triangulation(
const unsigned int,
get_locally_owned_vertices(
const Triangulation<deal_II_dimension, deal_II_space_dimension> &);
- template std::map<unsigned int, Point<deal_II_space_dimension>>
- get_all_vertices_at_boundary(
- const Triangulation<deal_II_dimension, deal_II_space_dimension> &tria);
-
- template std::vector<
- std::set<Triangulation<deal_II_dimension,
- deal_II_space_dimension>::active_cell_iterator>>
- vertex_to_cell_map(
- const Triangulation<deal_II_dimension, deal_II_space_dimension>
- &triangulation);
-
template std::vector<std::vector<Tensor<1, deal_II_space_dimension>>>
vertex_to_cell_centers_directions(
const Triangulation<deal_II_dimension, deal_II_space_dimension> &mesh,
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/bounding_box.h>
+#include <deal.II/base/floating_point_comparator.h>
+#include <deal.II/base/geometry_info.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/grid/grid_tools_geometry.h>
+#include <deal.II/grid/grid_tools_topology.h>
+
+#include <algorithm>
+#include <map>
+#include <set>
+#include <vector>
+
+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 <int structdim>
+ struct CellDataComparator
+ {
+ bool
+ operator()(const CellData<structdim> &a,
+ const CellData<structdim> &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<dim-1>, 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 <int dim>
+ class FaceDataHelper
+ {
+ public:
+ /**
+ * Store the data about the given face @p face.
+ */
+ template <typename FaceIteratorType>
+ void
+ insert_face_data(const FaceIteratorType &face)
+ {
+ CellData<dim - 1> 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<dim - 1> &face_cell_data : face_data)
+ internal::append_face_data(face_cell_data, subcell_data);
+ return subcell_data;
+ }
+
+
+ private:
+ std::set<CellData<dim - 1>, internal::CellDataComparator<dim - 1>>
+ face_data;
+ };
+
+
+ // Do nothing for dim=1:
+ template <>
+ class FaceDataHelper<1>
+ {
+ public:
+ template <typename FaceIteratorType>
+ void
+ insert_face_data(const FaceIteratorType &)
+ {}
+
+ SubCellData
+ get()
+ {
+ return SubCellData();
+ }
+ };
+ } // namespace internal
+
+
+
+ template <int dim, int spacedim>
+ std::
+ tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
+ get_coarse_mesh_description(const Triangulation<dim, spacedim> &tria)
+ {
+ Assert(tria.n_levels() >= 1,
+ ExcMessage("The input triangulation must be non-empty."));
+
+ std::vector<Point<spacedim>> vertices = tria.get_vertices();
+ std::vector<CellData<dim>> cells;
+
+ internal::FaceDataHelper<dim> face_data;
+ std::set<CellData<1>, internal::CellDataComparator<1>>
+ line_data; // only used in 3d
+
+ for (const auto &cell : tria.cell_iterators_on_level(0))
+ {
+ // Save cell data
+ CellData<dim> 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<Point<spacedim>>,
+ std::vector<CellData<dim>>,
+ SubCellData>(std::move(vertices),
+ std::move(cells),
+ std::move(subcell_data));
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ delete_unused_vertices(std::vector<Point<spacedim>> &vertices,
+ std::vector<CellData<dim>> &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<bool> 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<unsigned int> 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<Point<spacedim>> 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 <int dim, int spacedim>
+ void
+ delete_duplicated_vertices(std::vector<Point<spacedim>> &vertices,
+ std::vector<CellData<dim>> &cells,
+ SubCellData &subcelldata,
+ std::vector<unsigned int> &considered_vertices,
+ const double tol)
+ {
+ if (tol == 0.0)
+ return; // nothing to do per definition
+
+ AssertIndexRange(2, vertices.size());
+ std::vector<unsigned int> 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<spacedim> 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<std::pair<unsigned int, Point<spacedim>>> 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<unsigned int, Point<spacedim>> &a,
+ const std::pair<unsigned int, Point<spacedim>> &b) {
+ return a.second[longest_coordinate_direction] <
+ b.second[longest_coordinate_direction];
+ });
+
+ auto within_tolerance = [=](const Point<spacedim> &a,
+ const Point<spacedim> &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<unsigned int, Point<spacedim>> &a,
+ const std::pair<unsigned int, Point<spacedim>> &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 <int dim>
+ void
+ delete_duplicated_vertices(std::vector<Point<dim>> &vertices,
+ const double tol)
+ {
+ if (vertices.empty())
+ return;
+
+ // 1) map point to local vertex index
+ std::map<Point<dim>, unsigned int, FloatingPointComparator<double>>
+ map_point_to_local_vertex_index{FloatingPointComparator<double>(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 <int dim, int spacedim>
+ std::size_t
+ invert_cells_with_negative_measure(
+ const std::vector<Point<spacedim>> &all_vertices,
+ std::vector<CellData<dim>> &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<const unsigned int> 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 <int dim, int spacedim>
+ void
+ invert_all_negative_measure_cells(
+ const std::vector<Point<spacedim>> &all_vertices,
+ std::vector<CellData<dim>> &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 <int dim>
+ bool
+ is_consistent(const std::vector<CellData<dim>> &cells)
+ {
+ std::set<CheapEdge> edges;
+
+ for (typename std::vector<CellData<dim>>::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<dim>::lines_per_cell; ++l)
+ {
+ const CheapEdge reverse_edge(
+ c->vertices[GeometryInfo<dim>::line_to_cell_vertices(l, 1)],
+ c->vertices[GeometryInfo<dim>::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<dim>::line_to_cell_vertices(l, 0)],
+ c->vertices[GeometryInfo<dim>::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 <int dim>
+ 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<dim>::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 <int dim>
+ 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<AdjacentCell>, 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<AdjacentCell>
+ {};
+
+
+ /**
+ * 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 <int dim>
+ 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<dim> &cell, const unsigned int edge_number)
+ : orientation_status(not_oriented)
+ {
+ Assert(edge_number < GeometryInfo<dim>::lines_per_cell,
+ ExcInternalError());
+
+ // copy vertices for this particular line
+ vertex_indices[0] =
+ cell
+ .vertices[GeometryInfo<dim>::line_to_cell_vertices(edge_number, 0)];
+ vertex_indices[1] =
+ cell
+ .vertices[GeometryInfo<dim>::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<dim> &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<dim> &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<dim> adjacent_cells;
+ };
+
+
+
+ /**
+ * A data structure that represents a cell with all of its vertices
+ * and edges.
+ */
+ template <int dim>
+ 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<dim> &c, const std::vector<Edge<dim>> &edge_list)
+ {
+ for (const unsigned int i : GeometryInfo<dim>::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<dim>::lines_per_cell; ++l)
+ {
+ const Edge<dim> 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<dim>::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<dim>::lines_per_cell];
+ };
+
+
+
+ template <int dim>
+ 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<unsigned int>
+ {};
+
+
+
+ /**
+ * From a list of cells, build a sorted vector that contains all of the
+ * edges that exist in the mesh.
+ */
+ template <int dim>
+ std::vector<Edge<dim>>
+ build_edges(const std::vector<CellData<dim>> &cells)
+ {
+ // build the edge list for all cells. because each cell has
+ // GeometryInfo<dim>::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<dim>> edge_list;
+ edge_list.reserve(cells.size() * GeometryInfo<dim>::lines_per_cell);
+ for (unsigned int i = 0; i < cells.size(); ++i)
+ for (unsigned int l = 0; l < GeometryInfo<dim>::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 <int dim>
+ std::vector<Cell<dim>>
+ build_cells_and_connect_edges(const std::vector<CellData<dim>> &cells,
+ std::vector<Edge<dim>> &edges)
+ {
+ std::vector<Cell<dim>> 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<dim>::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 <int dim>
+ unsigned int
+ get_next_unoriented_cell(const std::vector<Cell<dim>> &cells,
+ const std::vector<Edge<dim>> &edges,
+ const unsigned int current_cell)
+ {
+ for (unsigned int c = current_cell; c < cells.size(); ++c)
+ for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
+ if (edges[cells[c].edge_indices[l]].orientation_status ==
+ Edge<dim>::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 <int dim>
+ void
+ orient_one_set_of_parallel_edges(const std::vector<Cell<dim>> &cells,
+ std::vector<Edge<dim>> &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<dim>::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<dim>::backward : Edge<dim>::forward);
+ else
+ {
+ Assert(
+ edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] ==
+ cells[cell].vertex_indices
+ [GeometryInfo<dim>::line_to_cell_vertices(local_edge, 1)],
+ ExcInternalError());
+ Assert(
+ edges[cells[cell].edge_indices[local_edge]].vertex_indices[1] ==
+ cells[cell].vertex_indices
+ [GeometryInfo<dim>::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<dim>::forward : Edge<dim>::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<dim> Delta_k;
+ EdgeDeltaSet<dim> 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<dim>::const_iterator delta =
+ Delta_k_minus_1.begin();
+ delta != Delta_k_minus_1.end();
+ ++delta)
+ {
+ Assert(edges[*delta].orientation_status !=
+ Edge<dim>::not_oriented,
+ ExcInternalError());
+
+ // now go through the cells adjacent to this edge
+ for (typename AdjacentCells<dim>::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<dim>::forward ?
+ edges[*delta].vertex_indices[0] :
+ edges[*delta].vertex_indices[1]);
+ const unsigned int first_edge_vertex_in_K =
+ cells[K]
+ .vertex_indices[GeometryInfo<dim>::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<dim>::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<dim>::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<dim>::OrientationStatus
+ opposite_edge_orientation =
+ (edges[opposite_edge].vertex_indices[0] ==
+ first_opposite_edge_vertex ?
+ Edge<dim>::forward :
+ Edge<dim>::backward);
+
+ // see if the opposite edge (there is only one in 2d) has
+ // already been oriented.
+ if (edges[opposite_edge].orientation_status ==
+ Edge<dim>::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 <code>raw_cells[cell_index]</code>.
+ */
+ template <int dim>
+ void
+ rotate_cell(const std::vector<Cell<dim>> &cell_list,
+ const std::vector<Edge<dim>> &edge_list,
+ const unsigned int cell_index,
+ std::vector<CellData<dim>> &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<dim>::lines_per_cell];
+ for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++e)
+ {
+ Assert(edge_list[cell_list[cell_index].edge_indices[e]]
+ .orientation_status != Edge<dim>::not_oriented,
+ ExcInternalError());
+ if (edge_list[cell_list[cell_index].edge_indices[e]]
+ .orientation_status == Edge<dim>::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<dim>::vertices_per_cell;
+ ++origin_vertex_of_cell)
+ if (std::count(starting_vertex_of_edge,
+ starting_vertex_of_edge +
+ GeometryInfo<dim>::lines_per_cell,
+ cell_list[cell_index]
+ .vertex_indices[origin_vertex_of_cell]) == dim)
+ break;
+ Assert(origin_vertex_of_cell <
+ GeometryInfo<dim>::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<dim>::vertices_per_cell];
+ for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
+ temp_vertex_indices[v] =
+ raw_cells[cell_index]
+ .vertices[cube_permutations[origin_vertex_of_cell][v]];
+ for (const unsigned int v : GeometryInfo<dim>::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 <int dim>
+ void
+ reorient(std::vector<CellData<dim>> &cells)
+ {
+ // first build the arrays that connect cells to edges and the other
+ // way around
+ std::vector<Edge<dim>> edge_list = build_edges(cells);
+ std::vector<Cell<dim>> 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<dim>::starter_edges[l]]]
+ .orientation_status == Edge<dim>::not_oriented)
+ orient_one_set_of_parallel_edges(
+ cell_list,
+ edge_list,
+ next_cell_with_unoriented_edge,
+ ParallelEdges<dim>::starter_edges[l]);
+
+ // ensure that we have really oriented all edges now, not just
+ // the starter edges
+ for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
+ Assert(edge_list[cell_list[next_cell_with_unoriented_edge]
+ .edge_indices[l]]
+ .orientation_status != Edge<dim>::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<CellData<1>> &)
+ {}
+ } // namespace
+
+
+
+ template <int dim>
+ void
+ consistently_order_cells(std::vector<CellData<dim>> &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 <int dim, int spacedim>
+ std::map<unsigned int, Point<spacedim>>
+ get_all_vertices_at_boundary(const Triangulation<dim, spacedim> &tria)
+ {
+ std::map<unsigned int, Point<spacedim>> vertex_map;
+ typename Triangulation<dim, spacedim>::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<dim, spacedim>::face_iterator &face =
+ cell->face(i);
+ if (face->at_boundary())
+ {
+ for (unsigned j = 0; j < face->n_vertices(); ++j)
+ {
+ const Point<spacedim> &vertex = face->vertex(j);
+ const unsigned int vertex_index = face->vertex_index(j);
+ vertex_map[vertex_index] = vertex;
+ }
+ }
+ }
+ }
+ return vertex_map;
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ remove_hanging_nodes(Triangulation<dim, spacedim> &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 <int dim, int spacedim>
+ void
+ remove_anisotropy(Triangulation<dim, spacedim> &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<unsigned int, double> info =
+ GridTools::get_longest_direction<dim, spacedim>(cell);
+ if (info.second > max_ratio)
+ {
+ cell->set_refine_flag(
+ RefinementCase<dim>::cut_axis(info.first));
+ continue_refinement = true;
+ }
+ }
+ tria.execute_coarsening_and_refinement();
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::map<unsigned int, Point<spacedim>>
+ extract_used_vertices(const Triangulation<dim, spacedim> &container,
+ const Mapping<dim, spacedim> &mapping)
+ {
+ std::map<unsigned int, Point<spacedim>> 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 <int dim, int spacedim>
+ std::vector<
+ std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
+ vertex_to_cell_map(const Triangulation<dim, spacedim> &triangulation)
+ {
+ std::vector<
+ std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
+ vertex_to_cell_map(triangulation.n_vertices());
+ typename Triangulation<dim, spacedim>::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<dim, spacedim>::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<dim, spacedim>::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 <int dim, int spacedim>
+ void
+ get_face_connectivity_of_cells(
+ const Triangulation<dim, spacedim> &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 <int dim, int spacedim>
+ void
+ get_vertex_connectivity_of_cells(
+ const Triangulation<dim, spacedim> &triangulation,
+ DynamicSparsityPattern &cell_connectivity)
+ {
+ std::vector<std::vector<unsigned int>> 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 <int dim, int spacedim>
+ void
+ get_vertex_connectivity_of_cells_on_level(
+ const Triangulation<dim, spacedim> &triangulation,
+ const unsigned int level,
+ DynamicSparsityPattern &cell_connectivity)
+ {
+ std::vector<std::vector<unsigned int>> vertex_to_cell(
+ triangulation.n_vertices());
+ for (typename Triangulation<dim, spacedim>::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<dim, spacedim>::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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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<Point<deal_II_space_dimension>>,
+ std::vector<CellData<deal_II_dimension>>,
+ SubCellData>
+ get_coarse_mesh_description(
+ const Triangulation<deal_II_dimension, deal_II_space_dimension> &tria);
+
+ template void
+ delete_unused_vertices(std::vector<Point<deal_II_space_dimension>> &,
+ std::vector<CellData<deal_II_dimension>> &,
+ SubCellData &);
+
+ template void
+ delete_duplicated_vertices(std::vector<Point<deal_II_space_dimension>> &,
+ std::vector<CellData<deal_II_dimension>> &,
+ SubCellData &,
+ std::vector<unsigned int> &,
+ double);
+
+# if deal_II_dimension == deal_II_space_dimension
+ template void
+ delete_duplicated_vertices(std::vector<Point<deal_II_space_dimension>> &,
+ const double);
+# endif
+
+ template void
+ invert_all_negative_measure_cells(
+ const std::vector<Point<deal_II_space_dimension>> &,
+ std::vector<CellData<deal_II_dimension>> &);
+
+ template std::size_t
+ invert_cells_with_negative_measure(
+ const std::vector<Point<deal_II_space_dimension>> &,
+ std::vector<CellData<deal_II_dimension>> &);
+
+# if deal_II_dimension == deal_II_space_dimension
+ template void
+ consistently_order_cells(std::vector<CellData<deal_II_dimension>> &);
+# endif
+
+ template std::map<unsigned int, Point<deal_II_space_dimension>>
+ get_all_vertices_at_boundary(
+ const Triangulation<deal_II_dimension, deal_II_space_dimension> &tria);
+
+ template std::vector<
+ std::set<Triangulation<deal_II_dimension,
+ deal_II_space_dimension>::active_cell_iterator>>
+ vertex_to_cell_map(
+ const Triangulation<deal_II_dimension, deal_II_space_dimension>
+ &triangulation);
+
+ template std::map<unsigned int, Point<deal_II_space_dimension>>
+ extract_used_vertices(
+ const Triangulation<deal_II_dimension, deal_II_space_dimension> &mesh,
+ const Mapping<deal_II_dimension, deal_II_space_dimension> &mapping);
+
+ template void
+ get_face_connectivity_of_cells(
+ const Triangulation<deal_II_dimension, deal_II_space_dimension>
+ &triangulation,
+ DynamicSparsityPattern &cell_connectivity);
+
+ template void
+ get_vertex_connectivity_of_cells(
+ const Triangulation<deal_II_dimension, deal_II_space_dimension>
+ &triangulation,
+ DynamicSparsityPattern &cell_connectivity);
+
+ template void
+ get_vertex_connectivity_of_cells_on_level(
+ const Triangulation<deal_II_dimension, deal_II_space_dimension>
+ &triangulation,
+ const unsigned int level,
+ DynamicSparsityPattern &cell_connectivity);
+
+ \}
+#endif
+ }
+
+for (deal_II_dimension : DIMENSIONS)
+ {
+ template void
+ GridTools::remove_hanging_nodes<deal_II_dimension, deal_II_dimension>(
+ Triangulation<deal_II_dimension> & tria, bool, unsigned int);
+
+ template void
+ GridTools::remove_anisotropy<deal_II_dimension, deal_II_dimension>(
+ Triangulation<deal_II_dimension> &, double, unsigned int);
+
+#if deal_II_dimension < 3
+ template void
+ GridTools::remove_hanging_nodes<deal_II_dimension, deal_II_dimension + 1>(
+ Triangulation<deal_II_dimension, deal_II_dimension + 1> & tria,
+ bool,
+ unsigned int);
+
+ template void
+ GridTools::remove_anisotropy<deal_II_dimension, deal_II_dimension + 1>(
+ Triangulation<deal_II_dimension, deal_II_dimension + 1> &,
+ double,
+ unsigned int);
+#endif
+ }
projects / dealii.git / commitdiff