--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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_tria_connectivity_h
+#define dealii_tria_connectivity_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/array_view.h>
+
+#include <deal.II/grid/reference_cell.h>
+#include <deal.II/grid/tria_description.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace internal
+{
+ namespace TriangulationImplementation
+ {
+ /**
+ * Interface of geometric cell entities with the focus on creating a
+ * reduced connectivity table.
+ */
+ struct CellTypeBase
+ {
+ /**
+ * Default destructor.
+ */
+ virtual ~CellTypeBase() = default;
+
+ /**
+ * Number of sub-entities of dimension @p d.
+ */
+ virtual unsigned int
+ n_entities(const unsigned int d) const
+ {
+ Assert(false, ExcNotImplemented());
+ (void)d;
+
+ return 0;
+ }
+
+ /**
+ * Number of vertices of the @p e-th sub-entity of dimension @p d.
+ */
+ virtual dealii::ArrayView<const unsigned int>
+ vertices_of_entity(const unsigned int d, const unsigned int e) const
+ {
+ Assert(false, ExcNotImplemented());
+ (void)d;
+ (void)e;
+
+ return {};
+ }
+
+ /**
+ * Geometric entity type of the @p e-th sub-entity of dimension @p d.
+ */
+ virtual ReferenceCell::Type
+ type_of_entity(const unsigned int d, const unsigned int e) const
+ {
+ Assert(false, ExcNotImplemented());
+ (void)d;
+ (void)e;
+
+ return ReferenceCell::Type::Vertex;
+ }
+
+ /**
+ * Number of lines of @p face-th surface.
+ */
+ virtual unsigned int
+ n_lines_of_surface(const unsigned int face) const
+ {
+ Assert(false, ExcNotImplemented());
+ (void)face;
+
+ return 0;
+ }
+
+ /**
+ * Index of the @p line-th lines of @p face-th surface.
+ */
+ virtual unsigned int
+ nth_line_of_surface(const unsigned int line,
+ const unsigned int face) const
+ {
+ Assert(false, ExcNotImplemented());
+ (void)line;
+ (void)face;
+
+ return 0;
+ }
+
+ /**
+ * Vertex indices of the @p line-th lines of @p face-th surface.
+ */
+ virtual const std::array<unsigned int, 2> &
+ vertices_of_nth_line_of_surface(const unsigned int line,
+ const unsigned int face) const
+ {
+ Assert(false, ExcNotImplemented());
+ (void)line;
+ (void)face;
+
+ const static std::array<unsigned int, 2> table = {};
+
+ return table;
+ }
+ };
+
+
+
+ /**
+ * Implementation for lines.
+ */
+ struct CellTypeLine : public CellTypeBase
+ {
+ dealii::ArrayView<const unsigned int>
+ vertices_of_entity(const unsigned int d,
+ const unsigned int e) const override
+ {
+ (void)e;
+
+ if (d == 1)
+ {
+ static const std::array<unsigned int, 2> table = {{0, 1}};
+
+ AssertDimension(e, 0);
+
+ return {table};
+ }
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
+
+ ReferenceCell::Type
+ type_of_entity(const unsigned int d, const unsigned int e) const override
+ {
+ (void)e;
+
+ if (d == 1)
+ return ReferenceCell::Type::Line;
+
+ Assert(false, ExcNotImplemented());
+
+ return ReferenceCell::Type::Vertex;
+ }
+
+ unsigned int
+ n_entities(const unsigned int d) const override
+ {
+ static std::array<unsigned int, 3> table = {2, 1};
+ return table[d];
+ }
+ };
+
+
+
+ /**
+ * Implementation for triangles.
+ */
+ struct CellTypeTri : public CellTypeBase
+ {
+ dealii::ArrayView<const unsigned int>
+ vertices_of_entity(const unsigned int d,
+ const unsigned int e) const override
+ {
+ if (d == 2)
+ {
+ static const std::array<unsigned int, 3> table = {0, 1, 2};
+
+ AssertDimension(e, 0);
+
+ return {table};
+ }
+
+ if (d == 1)
+ {
+ static const std::array<std::array<unsigned int, 2>, 3> table = {
+ {{0, 1}, {1, 2}, {2, 0}}};
+
+ return {table[e]};
+ }
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
+
+ virtual ReferenceCell::Type
+ type_of_entity(const unsigned int d, const unsigned int e) const override
+ {
+ (void)e;
+
+ if (d == 2)
+ return ReferenceCell::Type::Tri;
+
+ if (d == 1)
+ return ReferenceCell::Type::Line;
+
+ Assert(false, ExcNotImplemented());
+
+ return ReferenceCell::Type::Vertex;
+ }
+
+ unsigned int
+ n_entities(const unsigned int d) const override
+ {
+ static std::array<unsigned int, 3> table = {3, 3, 1};
+ return table[d];
+ }
+ };
+
+
+
+ /**
+ * Implementation for quadrilaterals.
+ */
+ struct CellTypeQuad : public CellTypeBase
+ {
+ dealii::ArrayView<const unsigned int>
+ vertices_of_entity(const unsigned int d,
+ const unsigned int e) const override
+ {
+ if (d == 2)
+ {
+ static const std::array<unsigned int, 4> table = {0, 1, 2, 3};
+
+ AssertDimension(e, 0);
+
+ return {table};
+ }
+
+ if (d == 1)
+ {
+ static const std::array<std::array<unsigned int, 2>, 4> table = {
+ {{0, 2}, {1, 3}, {0, 1}, {2, 3}}};
+
+ return {table[e]};
+ }
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
+
+ virtual ReferenceCell::Type
+ type_of_entity(const unsigned int d, const unsigned int e) const override
+ {
+ (void)e;
+
+ if (d == 2)
+ return ReferenceCell::Type::Quad;
+
+ if (d == 1)
+ return ReferenceCell::Type::Line;
+
+ Assert(false, ExcNotImplemented());
+
+ return ReferenceCell::Type::Vertex;
+ }
+
+ unsigned int
+ n_entities(const unsigned int d) const override
+ {
+ static std::array<unsigned int, 3> table = {4, 4, 1};
+ return table[d];
+ }
+ };
+
+
+
+ /**
+ * Implementation for tetrahedrons.
+ */
+ struct CellTypeTet : public CellTypeBase
+ {
+ dealii::ArrayView<const unsigned int>
+ vertices_of_entity(const unsigned int d,
+ const unsigned int e) const override
+ {
+ if (d == 3)
+ {
+ static const std::array<unsigned int, 4> table = {0, 1, 2, 3};
+
+ AssertDimension(e, 0);
+
+ return {table};
+ }
+
+ if (d == 2)
+ {
+ static const std::array<std::array<unsigned int, 3>, 4> table = {
+ {{0, 1, 2}, {1, 0, 3}, {0, 2, 3}, {2, 1, 3}}};
+
+ return {table[e]};
+ }
+
+ if (d == 1)
+ {
+ static const std::array<std::array<unsigned int, 2>, 6> table = {
+ {{0, 1}, {1, 2}, {2, 0}, {0, 3}, {1, 3}, {2, 3}}};
+
+ return {table[e]};
+ }
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
+
+ virtual ReferenceCell::Type
+ type_of_entity(const unsigned int d, const unsigned int e) const override
+ {
+ (void)e;
+
+ if (d == 3)
+ return ReferenceCell::Type::Tet;
+
+ if (d == 2)
+ return ReferenceCell::Type::Tri;
+
+ if (d == 1)
+ return ReferenceCell::Type::Line;
+
+ Assert(false, ExcNotImplemented());
+
+ return ReferenceCell::Type::Vertex;
+ }
+
+ unsigned int
+ n_entities(const unsigned int d) const override
+ {
+ static std::array<unsigned int, 4> table = {4, 6, 4, 1};
+ return table[d];
+ }
+
+ unsigned int
+ n_lines_of_surface(const unsigned int line) const override
+ {
+ (void)line;
+ return 3;
+ }
+
+ unsigned int
+ nth_line_of_surface(const unsigned int line,
+ const unsigned int face) const override
+ {
+ const static std::array<std::array<unsigned int, 3>, 4> table = {
+ {{0, 1, 2}, {0, 3, 4}, {2, 5, 3}, {1, 4, 5}}};
+
+ return table[face][line];
+ }
+
+ const std::array<unsigned int, 2> &
+ vertices_of_nth_line_of_surface(const unsigned int line,
+ const unsigned int face) const override
+ {
+ const static std::array<std::array<std::array<unsigned int, 2>, 3>, 4>
+ table = {{{{{0, 1}, {1, 2}, {2, 0}}},
+ {{{1, 0}, {0, 3}, {3, 1}}},
+ {{{0, 2}, {2, 3}, {3, 0}}},
+ {{{2, 1}, {1, 3}, {3, 2}}}}};
+
+ return table[face][line];
+ }
+ };
+
+
+ /**
+ * Implementation for pyramids.
+ */
+
+ struct CellTypePyramid : public CellTypeBase
+ {
+ dealii::ArrayView<const unsigned int>
+ vertices_of_entity(const unsigned int d,
+ const unsigned int e) const override
+ {
+ if (d == 3)
+ {
+ static const std::array<unsigned int, 5> table = {0, 1, 2, 3, 4};
+
+ AssertDimension(e, 0);
+
+ return {table};
+ }
+
+ if (d == 2)
+ {
+ if (e == 0)
+ {
+ static const std::array<unsigned int, 4> table = {0, 1, 2, 3};
+ return {table};
+ }
+
+ static const std::array<std::array<unsigned int, 3>, 4> table = {
+ {{0, 2, 4}, {3, 1, 4}, {1, 0, 4}, {2, 3, 4}}};
+
+ return {table[e - 1]};
+ }
+
+ if (d == 1)
+ {
+ static const std::array<std::array<unsigned int, 2>, 8> table = {
+ {{0, 2}, {1, 3}, {0, 1}, {2, 3}, {0, 4}, {1, 4}, {2, 4}, {3, 4}}};
+
+ return {table[e]};
+ }
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
+
+ virtual ReferenceCell::Type
+ type_of_entity(const unsigned int d, const unsigned int e) const override
+ {
+ (void)e;
+
+ if (d == 3)
+ return ReferenceCell::Type::Pyramid;
+
+ if (d == 2 && e == 0)
+ return ReferenceCell::Type::Quad;
+ else if (d == 2)
+ return ReferenceCell::Type::Tri;
+
+ if (d == 1)
+ return ReferenceCell::Type::Line;
+
+ Assert(false, ExcNotImplemented());
+
+ return ReferenceCell::Type::Vertex;
+ }
+
+ unsigned int
+ n_entities(const unsigned int d) const override
+ {
+ static std::array<unsigned int, 4> table = {5, 8, 5, 1};
+ return table[d];
+ }
+
+ unsigned int
+ n_lines_of_surface(const unsigned int surface) const override
+ {
+ if (surface == 0)
+ return 4;
+
+ return 3;
+ }
+
+ unsigned int
+ nth_line_of_surface(const unsigned int line,
+ const unsigned int face) const override
+ {
+ const static std::array<std::array<unsigned int, 4>, 5> table = {
+ {{0, 1, 2, 3},
+ {0, 6, 4, numbers::invalid_unsigned_int},
+ {1, 5, 7, numbers::invalid_unsigned_int},
+ {2, 4, 5, numbers::invalid_unsigned_int},
+ {3, 7, 6, numbers::invalid_unsigned_int}}};
+
+ return table[face][line];
+ }
+
+ const std::array<unsigned int, 2> &
+ vertices_of_nth_line_of_surface(const unsigned int line,
+ const unsigned int face) const override
+ {
+ static const unsigned int X = static_cast<unsigned int>(-1);
+
+ const static std::array<std::array<std::array<unsigned int, 2>, 4>, 5>
+ table = {{{{{0, 2}, {1, 3}, {0, 1}, {2, 3}}},
+ {{{0, 2}, {2, 4}, {4, 0}, {X, X}}},
+ {{{3, 1}, {1, 4}, {4, 3}, {X, X}}},
+ {{{1, 0}, {0, 4}, {4, 1}, {X, X}}},
+ {{{2, 3}, {3, 4}, {4, 2}, {X, X}}}}};
+
+ return table[face][line];
+ }
+ };
+
+
+
+ /**
+ * Implementation for wedges.
+ */
+ struct CellTypeWedge : public CellTypeBase
+ {
+ dealii::ArrayView<const unsigned int>
+ vertices_of_entity(const unsigned int d,
+ const unsigned int e) const override
+ {
+ if (d == 3)
+ {
+ static const std::array<unsigned int, 6> table = {0, 1, 2, 3, 4, 5};
+
+ AssertDimension(e, 0);
+
+ return {table};
+ }
+
+ if (d == 2)
+ {
+ if (e == 0 || e == 1)
+ {
+ static const std::array<std::array<unsigned int, 3>, 2> table =
+ {{{1, 0, 2}, {3, 4, 5}}};
+
+ return {table[e]};
+ }
+
+ static const std::array<std::array<unsigned int, 4>, 3> table = {
+ {{0, 1, 3, 4}, {1, 2, 4, 5}, {2, 0, 5, 3}}};
+
+ return {table[e - 2]};
+ }
+
+ if (d == 1)
+ {
+ static const std::array<std::array<unsigned int, 2>, 9> table = {
+ {{0, 1},
+ {1, 2},
+ {2, 0},
+ {3, 4},
+ {4, 5},
+ {5, 3},
+ {0, 3},
+ {1, 4},
+ {2, 5}}};
+
+ return {table[e]};
+ }
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
+
+ virtual ReferenceCell::Type
+ type_of_entity(const unsigned int d, const unsigned int e) const override
+ {
+ (void)e;
+
+ if (d == 3)
+ return ReferenceCell::Type::Wedge;
+
+ if (d == 2 && e > 1)
+ return ReferenceCell::Type::Quad;
+ else if (d == 2)
+ return ReferenceCell::Type::Tri;
+
+ if (d == 1)
+ return ReferenceCell::Type::Line;
+
+ Assert(false, ExcNotImplemented());
+
+ return ReferenceCell::Type::Vertex;
+ }
+
+ unsigned int
+ n_entities(const unsigned int d) const override
+ {
+ static std::array<unsigned int, 4> table = {6, 9, 5, 1};
+ return table[d];
+ }
+
+ unsigned int
+ n_lines_of_surface(const unsigned int surface) const override
+ {
+ if (surface > 1)
+ return 4;
+
+ return 3;
+ }
+
+ unsigned int
+ nth_line_of_surface(const unsigned int line,
+ const unsigned int face) const override
+ {
+ static const unsigned int X = static_cast<unsigned int>(-1);
+
+ const static std::array<std::array<unsigned int, 4>, 5> table = {
+ {{0, 2, 1, X},
+ {3, 4, 5, X},
+ {6, 7, 0, 3},
+ {7, 8, 1, 4},
+ {8, 6, 5, 2}}};
+
+ return table[face][line];
+ }
+
+ const std::array<unsigned int, 2> &
+ vertices_of_nth_line_of_surface(const unsigned int line,
+ const unsigned int face) const override
+ {
+ static const unsigned int X = static_cast<unsigned int>(-1);
+
+ const static std::array<std::array<std::array<unsigned int, 2>, 4>, 5>
+ table = {{{{{1, 0}, {0, 2}, {2, 1}, {X, X}}},
+ {{{3, 4}, {4, 5}, {5, 3}, {X, X}}},
+ {{{0, 3}, {1, 4}, {0, 1}, {3, 4}}},
+ {{{1, 4}, {2, 5}, {1, 2}, {4, 5}}},
+ {{{2, 5}, {0, 3}, {2, 0}, {5, 3}}}}};
+
+ return table[face][line];
+ }
+ };
+
+
+
+ /**
+ * Implementation for hexahedrons.
+ */
+ struct CellTypeHex : public CellTypeBase
+ {
+ dealii::ArrayView<const unsigned int>
+ vertices_of_entity(const unsigned int d,
+ const unsigned int e) const override
+ {
+ if (d == 3)
+ {
+ static const std::array<unsigned int, 8> table = {
+ 0, 1, 2, 3, 4, 5, 6, 7};
+
+ AssertDimension(e, 0);
+
+ return {table};
+ }
+
+ if (d == 2)
+ {
+ static const std::array<std::array<unsigned int, 4>, 6> table = {
+ {{0, 2, 4, 6},
+ {1, 3, 5, 7},
+ {0, 4, 1, 5},
+ {2, 6, 3, 7},
+ {0, 1, 2, 3},
+ {4, 5, 6, 7}}};
+
+ return {table[e]};
+ }
+
+ if (d == 1)
+ {
+ static const std::array<std::array<unsigned int, 2>, 12> table = {
+ {{0, 2},
+ {1, 3},
+ {0, 1},
+ {2, 3},
+ {4, 6},
+ {5, 7},
+ {4, 5},
+ {6, 7},
+ {0, 4},
+ {1, 5},
+ {2, 6},
+ {3, 7}}};
+
+ return {table[e]};
+ }
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
+
+ virtual ReferenceCell::Type
+ type_of_entity(const unsigned int d, const unsigned int e) const override
+ {
+ (void)e;
+
+ if (d == 3)
+ return ReferenceCell::Type::Hex;
+
+ if (d == 2)
+ return ReferenceCell::Type::Quad;
+
+ if (d == 1)
+ return ReferenceCell::Type::Line;
+
+ Assert(false, ExcNotImplemented());
+
+ return ReferenceCell::Type::Vertex;
+ }
+
+ unsigned int
+ n_entities(const unsigned int d) const override
+ {
+ static std::array<unsigned int, 4> table = {8, 12, 6, 1};
+ return table[d];
+ }
+
+ unsigned int
+ n_lines_of_surface(const unsigned int surface) const override
+ {
+ (void)surface;
+ return 4;
+ }
+
+ unsigned int
+ nth_line_of_surface(const unsigned int line,
+ const unsigned int face) const override
+ {
+ const static std::array<std::array<unsigned int, 4>, 6> table = {
+ {{8, 10, 0, 4},
+ {9, 11, 1, 5},
+ {2, 6, 8, 9},
+ {3, 7, 10, 11},
+ {0, 1, 2, 3},
+ {4, 5, 6, 7}}};
+
+ return table[face][line];
+ }
+
+ const std::array<unsigned int, 2> &
+ vertices_of_nth_line_of_surface(const unsigned int line,
+ const unsigned int face) const override
+ {
+ const static std::array<std::array<std::array<unsigned int, 2>, 4>, 6>
+ table = {{{{{0, 4}, {2, 6}, {0, 2}, {4, 6}}},
+ {{{1, 5}, {3, 7}, {1, 3}, {5, 7}}},
+ {{{0, 1}, {4, 5}, {0, 4}, {1, 5}}},
+ {{{2, 3}, {6, 7}, {2, 6}, {3, 7}}},
+ {{{0, 2}, {1, 3}, {0, 1}, {2, 3}}},
+ {{{4, 6}, {5, 7}, {4, 5}, {6, 7}}}}};
+
+ return table[face][line];
+ }
+ };
+
+
+
+ /**
+ * Compressed row storage sparse matrix. This class is similar to
+ * SparsityPattern but reduced to the bare minimum as needed here - in the
+ * context of setting up the connectivity - and allowing direct simplified
+ * access to the entries.
+ */
+ template <typename T = unsigned int>
+ struct CRS
+ {
+ /**
+ * Default constructor.
+ */
+ CRS()
+ : ptr{0} {};
+
+ /**
+ * Constructor which allows to set the internal fields directly.
+ */
+ CRS(const std::vector<std::size_t> &ptr, const std::vector<T> &col)
+ : ptr(ptr)
+ , col(col)
+ {}
+
+ // row index
+ std::vector<std::size_t> ptr;
+
+ // column index
+ std::vector<T> col;
+ };
+
+
+
+ /**
+ * Class for storing the reduced connectivity table.
+ *
+ * A full connectivity table contains all possible connectivities of
+ * entities of dimension d and entities of dimension d' with 0<=d,d'<=dim.
+ * However, in the library we only need the following types of
+ * connectivities:
+ * - dim-dimensional neighbors of dim-dimensional entities (connected via
+ * faces)
+ * - d-dimensional entity to it's (d-1)-dimension bounding entities
+ * - quad (2 - 3D), line (1 - 2D/3D) to vertices (0) to be able to process
+ * the user provided SubCellData during
+ * Triangulation::create_triangulation().
+ * We call a table, which computes the corresponding entries of a full
+ * connectivity table a reduced table.
+ *
+ * The entries of the reduced table are as follows for 1D-3D:
+ *
+ * 1D : | 0 1 2D: | 0 1 2 3D: | 0 1 2 3
+ * ---+----- ---+------- ---+--------
+ * 0 | 0 | 0 |
+ * 1 | x n 1 | x 1 | x
+ * 2 | s x n 2 | s x
+ * 3 | x n
+ *
+ * with markers highlighting the reason for the entry x:=bounding entities;
+ * n:= neighboring entities; s:=sub-cell data
+ */
+ template <typename T = unsigned int>
+ struct Connectivity
+ {
+ Connectivity(const unsigned int dim,
+ const std::vector<ReferenceCell::Type> &cell_types)
+ : dim(dim)
+ , cell_types(cell_types)
+ {}
+
+ inline std::vector<unsigned char> &
+ entity_orientations(const unsigned int structdim)
+ {
+ if (structdim == 1)
+ return line_orientation;
+
+ AssertDimension(structdim, 2);
+
+ return quad_orientation;
+ }
+
+ inline const std::vector<unsigned char> &
+ entity_orientations(const unsigned int structdim) const
+ {
+ if (structdim == 1)
+ return line_orientation;
+
+ AssertDimension(structdim, 2);
+
+ return quad_orientation;
+ }
+
+ inline std::vector<ReferenceCell::Type> &
+ entity_types(const unsigned int structdim)
+ {
+ if (structdim == dim)
+ return cell_types;
+
+ // for vertices/lines the entity types are clear (0/1)
+ AssertDimension(structdim, 2);
+ AssertDimension(dim, 3);
+
+ return quad_types;
+ }
+
+ inline const std::vector<ReferenceCell::Type> &
+ entity_types(const unsigned int structdim) const
+ {
+ if (structdim == dim)
+ return cell_types;
+
+ // for vertices/lines the entity types are clear (0/1)
+ AssertDimension(structdim, 2);
+ AssertDimension(dim, 3);
+
+ return quad_types;
+ }
+
+ inline CRS<T> &
+ entity_to_entities(const unsigned int from, const unsigned int to)
+ {
+ if (from == dim && to == dim)
+ return neighbors;
+ else if (from == dim && to == dim - 1)
+ return cell_entities;
+ else if (dim == 3 && from == 2 && to == 0)
+ return quad_vertices;
+ else if (dim == 3 && from == 2 && to == 1)
+ return quad_lines;
+ else if (from == 1 && to == 0)
+ return line_vertices;
+
+ Assert(false, ExcNotImplemented());
+
+ return cell_entities;
+ }
+
+ inline const CRS<T> &
+ entity_to_entities(const unsigned int from, const unsigned int to) const
+ {
+ if (from == dim && to == dim)
+ return neighbors;
+ else if (from == dim && to == dim - 1)
+ return cell_entities;
+ else if (dim == 3 && from == 2 && to == 0)
+ return quad_vertices;
+ else if (dim == 3 && from == 2 && to == 1)
+ return quad_lines;
+ else if (from == 1 && to == 0)
+ return line_vertices;
+
+ Assert(false, ExcNotImplemented());
+
+ return cell_entities;
+ }
+
+ private:
+ const unsigned int dim;
+ std::vector<ReferenceCell::Type> cell_types;
+
+ CRS<T> line_vertices;
+
+ std::vector<unsigned char> line_orientation;
+
+ CRS<T> quad_vertices;
+ CRS<T> quad_lines;
+
+ std::vector<unsigned char> quad_orientation;
+
+ CRS<T> cell_entities;
+ CRS<T> neighbors;
+
+ std::vector<ReferenceCell::Type> quad_types;
+ };
+
+
+
+ /**
+ * Determine the neighbors of all cells.
+ *
+ * @p con_cf connectivity cell-face
+ * @p con_cc connectivity cell-cell (for each cell-face it contains the
+ * the index of the neighboring cell or -1 for boundary face)
+ */
+ template <typename T>
+ void
+ determine_neighbors(const CRS<T> &con_cf, CRS<T> &con_cc)
+ {
+ const auto &col_cf = con_cf.col;
+ const auto &ptr_cf = con_cf.ptr;
+
+ auto &col_cc = con_cc.col;
+ auto &ptr_cc = con_cc.ptr;
+
+ const unsigned int n_faces =
+ *std::max_element(col_cf.begin(), col_cf.end()) + 1;
+
+ // clear and initialize with -1 (assume that all faces are at the
+ // boundary)
+ col_cc = std::vector<T>(col_cf.size(), -1);
+ ptr_cc = ptr_cf;
+
+ std::vector<std::pair<T, unsigned int>> neighbors(n_faces, {-1, -1});
+
+ // loop over all cells
+ for (unsigned int i_0 = 0; i_0 < ptr_cf.size() - 1; i_0++)
+ {
+ // ... and all its faces
+ for (std::size_t j_0 = ptr_cf[i_0]; j_0 < ptr_cf[i_0 + 1]; j_0++)
+ {
+ if (neighbors[col_cf[j_0]].first == static_cast<unsigned int>(-1))
+ {
+ // face is visited the first time -> save the visiting cell
+ // and the face pointer
+ neighbors[col_cf[j_0]] = std::pair<T, unsigned int>(i_0, j_0);
+ }
+ else
+ {
+ // face is visited the second time -> now we know the cells
+ // on both sides of the face and we can determine for both
+ // cells the neigbor
+ col_cc[j_0] = neighbors[col_cf[j_0]].first;
+ col_cc[neighbors[col_cf[j_0]].second] = i_0;
+ }
+ }
+ }
+ }
+
+
+
+ /**
+ * Determine the orientation of an entity of @p type described by its
+ * vertices @p var_1 relative to an entity described by @p var_0.
+ */
+ template <typename T, std::size_t N>
+ inline unsigned char
+ compute_orientation(const ReferenceCell::Type entity_type,
+ const std::array<T, N> & vertices_0,
+ const std::array<T, N> & vertices_1)
+ {
+ if (entity_type == ReferenceCell::Type::Line)
+ {
+ const std::array<T, 2> i{vertices_0[0], vertices_0[1]};
+ const std::array<T, 3> j{vertices_1[0], vertices_1[1]};
+
+ // line_orientation=true
+ if (i == std::array<T, 2>{{j[0], j[1]}})
+ return 1;
+
+ // line_orientation=false
+ if (i == std::array<T, 2>{{j[1], j[0]}})
+ return 0;
+ }
+ else if (entity_type == ReferenceCell::Type::Tri)
+ {
+ const std::array<T, 3> i{vertices_0[0], vertices_0[1], vertices_0[2]};
+ const std::array<T, 3> j{vertices_1[0], vertices_1[1], vertices_1[2]};
+
+ // face_orientation=true, face_rotation=false, face_flip=false
+ if (i == std::array<T, 3>{{j[0], j[1], j[2]}})
+ return 1;
+
+ // face_orientation=true, face_rotation=true, face_flip=false
+ if (i == std::array<T, 3>{{j[1], j[0], j[2]}})
+ return 3;
+
+ // face_orientation=true, face_rotation=false, face_flip=true
+ if (i == std::array<T, 3>{{j[2], j[0], j[1]}})
+ return 5;
+
+ // face_orientation=false, face_rotation=false, face_flip=false
+ if (i == std::array<T, 3>{{j[0], j[2], j[1]}})
+ return 0;
+
+ // face_orientation=false, face_rotation=true, face_flip=false
+ if (i == std::array<T, 3>{{j[1], j[2], j[0]}})
+ return 2;
+
+ // face_orientation=false, face_rotation=false, face_flip=true
+ if (i == std::array<T, 3>{{j[2], j[1], j[0]}})
+ return 4;
+ }
+ else if (entity_type == ReferenceCell::Type::Quad)
+ {
+ const std::array<T, 4> i{vertices_0[0],
+ vertices_0[1],
+ vertices_0[2],
+ vertices_0[3]};
+ const std::array<T, 4> j{vertices_1[0],
+ vertices_1[1],
+ vertices_1[2],
+ vertices_1[3]};
+
+ // face_orientation=true, face_rotation=false, face_flip=false
+ if (i == std::array<T, 4>{{j[0], j[1], j[2], j[3]}})
+ return 1;
+
+ // face_orientation=true, face_rotation=true, face_flip=false
+ if (i == std::array<T, 4>{{j[1], j[3], j[0], j[2]}})
+ return 3;
+
+ // face_orientation=true, face_rotation=false, face_flip=true
+ if (i == std::array<T, 4>{{j[3], j[2], j[1], j[0]}})
+ return 5;
+
+ // face_orientation=true, face_rotation=true, face_flip=true
+ if (i == std::array<T, 4>{{j[2], j[0], j[3], j[1]}})
+ return 7;
+
+ // face_orientation=false, face_rotation=false, face_flip=false
+ if (i == std::array<T, 4>{{j[0], j[2], j[1], j[3]}})
+ return 0;
+
+ // face_orientation=false, face_rotation=true, face_flip=false
+ if (i == std::array<T, 4>{{j[2], j[3], j[0], j[1]}})
+ return 2;
+
+ // face_orientation=false, face_rotation=false, face_flip=true
+ if (i == std::array<T, 4>{{j[3], j[1], j[2], j[0]}})
+ return 4;
+
+ // face_orientation=false, face_rotation=true, face_flip=true
+ if (i == std::array<T, 4>{{j[1], j[0], j[3], j[2]}})
+ return 6;
+ }
+
+ AssertThrow(false, dealii::StandardExceptions::ExcNotImplemented());
+
+ return -1;
+ }
+
+
+
+ /**
+ * Build entities of dimension d (with 0<d<dim). Entities are described by
+ * a set of vertices.
+ *
+ * Furthermore, the function determines for each cell of which d-dimensional
+ * entity it consists of and its orientation relative to the cell.
+ */
+ template <int key_length, typename FU>
+ void
+ build_entity_templated(
+ const unsigned int d,
+ const std::vector<std::shared_ptr<CellTypeBase>> &cell_types,
+ const std::vector<ReferenceCell::Type> & cell_types_index,
+ const CRS<unsigned int> & crs,
+ CRS<unsigned int> & crs_d, // result
+ CRS<unsigned int> & crs_0, // result
+ std::vector<unsigned char> & orientations, // result
+ const FU & second_key_function)
+ {
+ const bool compatibility_mode = true;
+
+ const std::vector<std::size_t> & cell_ptr = crs.ptr;
+ const std::vector<unsigned int> &cell_vertices = crs.col;
+ std::vector<std::size_t> & ptr_d = crs_d.ptr;
+ std::vector<unsigned int> & col_d = crs_d.col;
+
+ // note: we do not pre-allocate memory for these arrays because it turned
+ // out that counting unique entities is more expensive than push_back().
+ std::vector<std::size_t> & ptr_0 = crs_0.ptr;
+ std::vector<unsigned int> &col_0 = crs_0.col;
+
+ // clear
+ ptr_0 = {};
+ col_0 = {};
+
+ unsigned int n_entities = 0;
+
+ for (const auto &c : cell_types_index)
+ n_entities +=
+ cell_types[static_cast<types::geometric_entity_type>(c)]->n_entities(
+ d);
+
+ // step 1: store each d-dimensional entity of a cell (described by their
+ // vertices) into a vector and create a key for them
+ //
+ // note: it turned out to be more efficient to have a vector of tuples
+ // than to have two vectors (sorting becomes inefficient)
+ std::vector<
+ std::tuple<std::array<unsigned int, key_length>, unsigned int>>
+ keys; // key (sorted vertices), cell-entity index
+
+ std::vector<std::array<unsigned int, key_length>> ad_entity_vertices;
+ std::vector<ReferenceCell::Type> ad_entity_types;
+ std::vector<std::array<unsigned int, key_length>> ad_compatibility;
+
+ keys.reserve(n_entities);
+ ad_entity_vertices.reserve(n_entities);
+ ad_entity_types.reserve(n_entities);
+ ad_compatibility.reserve(n_entities);
+
+ ptr_d.resize(cell_types_index.size() + 1);
+ ptr_d[0] = 0;
+
+ static const unsigned int offset = 1;
+
+ // loop over all cells
+ for (unsigned int c = 0, counter = 0; c < cell_types_index.size(); c++)
+ {
+ const auto &cell_type =
+ cell_types[static_cast<types::geometric_entity_type>(
+ cell_types_index[c])];
+ ptr_d[c + 1] = ptr_d[c] + cell_type->n_entities(d);
+
+ // ... collect vertices of cell
+ const dealii::ArrayView<const unsigned int> cell_vertice(
+ cell_vertices.data() + cell_ptr[c], cell_ptr[c + 1] - cell_ptr[c]);
+
+ // ... loop over all its entities
+ for (unsigned int e = 0; e < cell_type->n_entities(d); e++)
+ {
+ // ... determine global entity vertices
+ const auto &local_entity_vertices =
+ cell_type->vertices_of_entity(d, e);
+
+ std::array<unsigned int, key_length> entity_vertices;
+ std::fill(entity_vertices.begin(), entity_vertices.end(), 0);
+
+ for (unsigned int i = 0; i < local_entity_vertices.size(); i++)
+ entity_vertices[i] =
+ cell_vertice[local_entity_vertices[i]] + offset;
+
+ // ... create key
+ std::array<unsigned int, key_length> key = entity_vertices;
+ std::sort(key.begin(), key.end());
+ keys.emplace_back(key, counter++);
+
+ ad_entity_vertices.emplace_back(entity_vertices);
+
+ ad_entity_types.emplace_back(cell_type->type_of_entity(d, e));
+
+ if (compatibility_mode)
+ ad_compatibility.emplace_back(
+ second_key_function(entity_vertices, cell_type, c, e));
+ }
+ }
+
+ col_d.resize(keys.size());
+ orientations.resize(keys.size());
+
+ // step 2: sort according to key so that entities with same key can be
+ // merged
+ std::sort(keys.begin(), keys.end());
+
+
+ if (compatibility_mode)
+ {
+ unsigned int n_unique_entities = 0;
+ unsigned int n_unique_entity_vertices = 0;
+
+ std::array<unsigned int, key_length> ref_key, new_key;
+ std::fill(ref_key.begin(), ref_key.end(), 0);
+ for (unsigned int i = 0; i < keys.size(); ++i)
+ {
+ const auto offset_i = std::get<1>(keys[i]);
+
+ if (ref_key != std::get<0>(keys[i]))
+ {
+ ref_key = std::get<0>(keys[i]);
+
+ n_unique_entities++;
+ n_unique_entity_vertices +=
+ cell_types[static_cast<types::geometric_entity_type>(
+ ad_entity_types[offset_i])]
+ ->n_entities(0);
+
+ new_key = ad_compatibility[offset_i];
+ }
+
+ std::get<0>(keys[i]) = new_key;
+ }
+
+ std::sort(keys.begin(), keys.end());
+
+ ptr_0.reserve(n_unique_entities);
+ col_0.reserve(n_unique_entity_vertices);
+ }
+
+
+ std::array<unsigned int, key_length> ref_key;
+ std::array<unsigned int, key_length> ref_indices;
+ std::fill(ref_key.begin(), ref_key.end(), 0);
+
+ for (unsigned int i = 0, counter = dealii::numbers::invalid_unsigned_int;
+ i < keys.size();
+ i++)
+ {
+ const auto offset_i = std::get<1>(keys[i]);
+
+ if (ref_key != std::get<0>(keys[i]))
+ {
+ // new key
+ counter++;
+ ref_key = std::get<0>(keys[i]);
+ ref_indices = ad_entity_vertices[offset_i];
+
+ ptr_0.push_back(col_0.size());
+ for (const auto j : ad_entity_vertices[offset_i])
+ if (j != 0)
+ col_0.push_back(j - offset);
+
+ // take its orientation as default
+ col_d[offset_i] = counter;
+ orientations[offset_i] = 1;
+ }
+ else
+ {
+ col_d[offset_i] = counter;
+ orientations[offset_i] =
+ compute_orientation(ad_entity_types[offset_i],
+ ref_indices,
+ ad_entity_vertices[offset_i]);
+ }
+ }
+ ptr_0.push_back(col_0.size());
+ }
+
+
+
+ /**
+ * Call the right templated function to be able to use std::array instead
+ * of std::vector.
+ */
+ template <typename FU>
+ void
+ build_entity(const unsigned int d,
+ const std::vector<std::shared_ptr<CellTypeBase>> &cell_types,
+ const std::vector<ReferenceCell::Type> &cell_types_index,
+ const CRS<unsigned int> & crs,
+ CRS<unsigned int> & crs_d,
+ CRS<unsigned int> & crs_0,
+ std::vector<unsigned char> & orientations,
+ const FU & second_key_function)
+ {
+ std::size_t key_length = 0;
+
+ for (const auto &c : cell_types_index)
+ {
+ const auto &cell_type =
+ cell_types[static_cast<types::geometric_entity_type>(c)];
+ for (unsigned int e = 0; e < cell_type->n_entities(d); e++)
+ key_length =
+ std::max(key_length, cell_type->vertices_of_entity(d, e).size());
+ }
+
+ if (key_length == 2)
+ build_entity_templated<2>(d,
+ cell_types,
+ cell_types_index,
+ crs,
+ crs_d,
+ crs_0,
+ orientations,
+ second_key_function);
+ else if (key_length == 3)
+ build_entity_templated<3>(d,
+ cell_types,
+ cell_types_index,
+ crs,
+ crs_d,
+ crs_0,
+ orientations,
+ second_key_function);
+ else if (key_length == 4)
+ build_entity_templated<4>(d,
+ cell_types,
+ cell_types_index,
+ crs,
+ crs_d,
+ crs_0,
+ orientations,
+ second_key_function);
+ else
+ AssertThrow(false, dealii::StandardExceptions::ExcNotImplemented());
+ }
+
+
+
+ /**
+ * Build surface lines described by:
+ * - connectivity quad -> line
+ * - orientation of line relative to the quad
+ *
+ * Furthermore, the type of the quad is determined.
+ */
+ void
+ build_intersection(
+ const std::vector<std::shared_ptr<CellTypeBase>> &cell_types,
+ const std::vector<ReferenceCell::Type> & cell_types_index,
+ const CRS<unsigned int> & con_cv,
+ const CRS<unsigned int> & con_cl,
+ const CRS<unsigned int> & con_lv,
+ const CRS<unsigned int> & con_cq,
+ const CRS<unsigned int> & con_qv,
+ const std::vector<unsigned char> & ori_cq,
+ CRS<unsigned int> & con_ql, // result
+ std::vector<unsigned char> & ori_ql, // result
+ std::vector<ReferenceCell::Type> & quad_t_id // result
+ )
+ {
+ // reset output
+ ori_ql = {};
+ con_ql.ptr = {};
+ con_ql.col = {};
+
+ con_ql.ptr.resize(con_qv.ptr.size());
+ con_ql.ptr[0] = 0;
+
+ quad_t_id.resize(con_qv.ptr.size() - 1);
+
+ // count the number of lines of each face
+ for (unsigned int c = 0; c < con_cq.ptr.size() - 1; ++c)
+ {
+ const auto &cell_type =
+ cell_types[static_cast<types::geometric_entity_type>(
+ cell_types_index[c])];
+
+ // loop over faces
+ for (unsigned int f_ = con_cq.ptr[c], f_index = 0;
+ f_ < con_cq.ptr[c + 1];
+ ++f_, ++f_index)
+ {
+ const unsigned int f = con_cq.col[f_];
+
+ con_ql.ptr[f + 1] = cell_type->n_lines_of_surface(f_index);
+ }
+ }
+
+ // use the counts to determine the offsets -> prefix sum
+ for (unsigned int i = 0; i < con_ql.ptr.size() - 1; ++i)
+ con_ql.ptr[i + 1] += con_ql.ptr[i];
+
+ // allocate memory
+ con_ql.col.resize(con_ql.ptr.back());
+ ori_ql.resize(con_ql.ptr.back());
+
+ // loop over cells
+ for (unsigned int c = 0; c < con_cq.ptr.size() - 1; ++c)
+ {
+ const auto &cell_type =
+ cell_types[static_cast<types::geometric_entity_type>(
+ cell_types_index[c])];
+
+ // loop over faces
+ for (unsigned int f_ = con_cq.ptr[c], f_index = 0;
+ f_ < con_cq.ptr[c + 1];
+ ++f_, ++f_index)
+ {
+ const unsigned int f = con_cq.col[f_];
+
+ // only faces with default orientation have to do something
+ if (ori_cq[f_] != 1)
+ continue;
+
+ // determine entity type of face
+ quad_t_id[f] = cell_type->type_of_entity(2, f_index);
+
+ // loop over lines
+ for (unsigned int l = 0;
+ l < cell_type->n_lines_of_surface(f_index);
+ ++l)
+ {
+ // determine global index of line
+ const unsigned int local_line_index =
+ cell_type->nth_line_of_surface(l, f_index);
+ const unsigned int global_line_index =
+ con_cl.col[con_cl.ptr[c] + local_line_index];
+ con_ql.col[con_ql.ptr[f] + l] = global_line_index;
+
+ // determine orientation of line
+ const auto line_vertices_1_ref =
+ cell_type->vertices_of_nth_line_of_surface(l, f_index);
+
+ bool same = true;
+ for (unsigned int v = 0; v < line_vertices_1_ref.size(); ++v)
+ if (con_cv.col[con_cv.ptr[c] + line_vertices_1_ref[v]] !=
+ con_lv.col[con_lv.ptr[global_line_index] + v])
+ {
+ same = false;
+ break;
+ }
+
+ // ... comparison gives orientation
+ ori_ql[con_ql.ptr[f] + l] = (same ? 1 : 0);
+ }
+ }
+ }
+ }
+
+
+
+ /**
+ * Build the reduced connectivity table for the given dimension @p dim.
+ *
+ * This function is inspired by the publication Anders Logg "Efficient
+ * Representation of Computational Meshes" and the FEniCS's DOLFIN mesh
+ * implementation. It has been strongly adjusted to efficiently solely meet
+ * our connectivity needs while sacrificing some of the flexibility there.
+ */
+ template <typename T>
+ Connectivity<T>
+ build_connectivity(const unsigned int dim,
+ const std::vector<std::shared_ptr<CellTypeBase>> &cell_t,
+ const std::vector<ReferenceCell::Type> &cell_t_id,
+ const CRS<T> & con_cv)
+ {
+ Connectivity<T> connectivity(dim, cell_t_id);
+
+ CRS<T> temp1; // needed for 3D
+
+ if (dim == 1)
+ connectivity.entity_to_entities(1, 0) = con_cv;
+
+ if (dim == 2 || dim == 3) // build lines
+ {
+ std::vector<unsigned char> dummy;
+
+ build_entity(1,
+ cell_t,
+ connectivity.entity_types(dim),
+ con_cv,
+ dim == 2 ? connectivity.entity_to_entities(2, 1) : temp1,
+ connectivity.entity_to_entities(1, 0),
+ dim == 2 ? connectivity.entity_orientations(1) : dummy,
+ [](auto key, const auto &, const auto &, const auto &) {
+ // to ensure same enumeration as in deal.II
+ return key;
+ });
+ }
+
+ if (dim == 3) // build quads
+ {
+ build_entity(
+ 2,
+ cell_t,
+ connectivity.entity_types(3),
+ con_cv,
+ connectivity.entity_to_entities(3, 2),
+ connectivity.entity_to_entities(2, 0),
+ connectivity.entity_orientations(2),
+ [&](auto key, const auto &cell_type, const auto &c, const auto &f) {
+ // to ensure same enumeration as in deal.II
+ AssertIndexRange(cell_type->n_lines_of_surface(f),
+ key.size() + 1);
+
+ unsigned int l = 0;
+
+ for (; l < cell_type->n_lines_of_surface(f); ++l)
+ key[l] =
+ temp1
+ .col[temp1.ptr[c] + cell_type->nth_line_of_surface(l, f)] +
+ 1 /*offset!*/;
+
+ for (; l < key.size(); ++l)
+ key[l] = 0;
+
+ return key;
+ });
+
+ // create connectivity: quad -> line
+ build_intersection(cell_t,
+ connectivity.entity_types(3),
+ con_cv,
+ temp1,
+ connectivity.entity_to_entities(1, 0),
+ connectivity.entity_to_entities(3, 2),
+ connectivity.entity_to_entities(2, 0),
+ connectivity.entity_orientations(2),
+ connectivity.entity_to_entities(2, 1),
+ connectivity.entity_orientations(1),
+ connectivity.entity_types(2));
+ }
+
+ // determine neighbors
+ determine_neighbors(connectivity.entity_to_entities(dim, dim - 1),
+ connectivity.entity_to_entities(dim, dim));
+
+ return connectivity;
+ }
+
+
+
+ /**
+ * Preprocessing step to remove the template argument dim.
+ */
+ template <typename T, int dim>
+ Connectivity<T>
+ build_connectivity(const std::vector<CellData<dim>> &cells)
+ {
+ // vector of possible cell entity types
+ std::vector<std::shared_ptr<CellTypeBase>> cell_types_impl(8);
+
+ cell_types_impl[static_cast<types::geometric_entity_type>(
+ ReferenceCell::Type::Line)]
+ .reset(new CellTypeLine());
+ cell_types_impl[static_cast<types::geometric_entity_type>(
+ ReferenceCell::Type::Tri)]
+ .reset(new CellTypeTri());
+ cell_types_impl[static_cast<types::geometric_entity_type>(
+ ReferenceCell::Type::Quad)]
+ .reset(new CellTypeQuad());
+ cell_types_impl[static_cast<types::geometric_entity_type>(
+ ReferenceCell::Type::Tet)]
+ .reset(new CellTypeTet());
+ cell_types_impl[static_cast<types::geometric_entity_type>(
+ ReferenceCell::Type::Pyramid)]
+ .reset(new CellTypePyramid());
+ cell_types_impl[static_cast<types::geometric_entity_type>(
+ ReferenceCell::Type::Wedge)]
+ .reset(new CellTypeWedge());
+ cell_types_impl[static_cast<types::geometric_entity_type>(
+ ReferenceCell::Type::Hex)]
+ .reset(new CellTypeHex());
+
+ // jump table to pick the right entity type
+ static const ReferenceCell::Type X = ReferenceCell::Type::Invalid;
+ static const std::array<const std::array<ReferenceCell::Type, 9>, 4>
+ table = {{{X, ReferenceCell::Type::Vertex, X, X, X, X, X, X, X},
+ {X, X, ReferenceCell::Type::Line, X, X, X, X, X, X},
+ {X,
+ X,
+ X,
+ ReferenceCell::Type::Tri,
+ ReferenceCell::Type::Quad,
+ X,
+ X,
+ X,
+ X},
+ {X,
+ X,
+ X,
+ X,
+ ReferenceCell::Type::Tet,
+ ReferenceCell::Type::Pyramid,
+ ReferenceCell::Type::Wedge,
+ X,
+ ReferenceCell::Type::Hex}}};
+
+ // determine cell types and process vertices
+ std::vector<T> cell_vertices;
+ cell_vertices.reserve(
+ std::accumulate(cells.begin(),
+ cells.end(),
+ 0,
+ [](const auto &result, const auto &cell) {
+ return result + cell.vertices.size();
+ }));
+
+ std::vector<std::size_t> cell_vertices_ptr;
+ cell_vertices_ptr.reserve(cells.size() + 1);
+ cell_vertices_ptr.push_back(0);
+
+ std::vector<ReferenceCell::Type> cell_types_indices;
+ cell_types_indices.reserve(cells.size());
+
+ // loop over cells and create CRS
+ for (const auto &cell : cells)
+ {
+ // determine cell type
+ const ReferenceCell::Type cell_type =
+ table[dim][cell.vertices.size()];
+
+ Assert(cell_type != ReferenceCell::Type::Invalid,
+ ExcNotImplemented());
+ AssertIndexRange(static_cast<types::geometric_entity_type>(cell_type),
+ cell_types_impl.size());
+ Assert(cell_types_impl[static_cast<types::geometric_entity_type>(
+ cell_type)]
+ .get() != nullptr,
+ ExcNotImplemented());
+
+ cell_types_indices.push_back(cell_type);
+
+ // create CRS of vertices (to remove template argument dim)
+ for (const auto &vertex : cell.vertices)
+ cell_vertices.push_back(vertex);
+
+ cell_vertices_ptr.push_back(cell_vertices.size());
+ }
+
+ // do the actual work
+ return build_connectivity<T>(dim,
+ cell_types_impl,
+ cell_types_indices,
+ {cell_vertices_ptr, cell_vertices});
+ }
+ } // namespace TriangulationImplementation
+} // namespace internal
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/grid/connectivity.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/magic_numbers.h>
#include <deal.II/grid/manifold.h>
/**
- * Create a triangulation from
- * given data. This function does
- * this work for 1-dimensional
- * triangulations independently
- * of the actual space dimension.
+ * Create a triangulation from given data.
*/
- template <int spacedim>
+ template <int dim, int spacedim>
static void
- create_triangulation(const std::vector<Point<spacedim>> &v,
- const std::vector<CellData<1>> & cells,
- const SubCellData & /*subcelldata*/,
- Triangulation<1, spacedim> &triangulation)
+ create_triangulation(const std::vector<Point<spacedim>> &vertices,
+ const std::vector<CellData<dim>> & cells,
+ const SubCellData & subcelldata,
+ Triangulation<dim, spacedim> & tria)
{
- AssertThrow(v.size() > 0, ExcMessage("No vertices given"));
- AssertThrow(cells.size() > 0, ExcMessage("No cells given"));
-
- // note: since no boundary
- // information can be given in one
- // dimension, the @p{subcelldata}
- // field is ignored. (only used for
- // error checking, which is a good
- // idea in any case)
- const unsigned int dim = 1;
-
- // copy vertices
- triangulation.vertices = v;
- triangulation.vertices_used = std::vector<bool>(v.size(), true);
-
- // Check that all cells have positive volume. This check is not run in
- // the codimension one or two cases since cell_measure is not
- // implemented for those.
-#ifndef _MSC_VER
- // TODO: The following code does not compile with MSVC. Find a way
- // around it
- if (dim == spacedim)
- {
- for (unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
- {
- // If we should check for distorted cells, then we permit them
- // to exist. If a cell has negative measure, then it must be
- // distorted (the converse is not necessarily true); hence
- // throw an exception if no such cells should exist.
- if (!triangulation.check_for_distorted_cells)
- {
- unsigned int vertices[GeometryInfo<1>::vertices_per_cell];
-
- for (unsigned int i = 0;
- i < GeometryInfo<1>::vertices_per_cell;
- ++i)
- vertices[i] = cells[cell_no].vertices[i];
-
- const double cell_measure =
- GridTools::cell_measure<1>(triangulation.vertices,
- vertices);
- AssertThrow(cell_measure > 0,
- ExcGridHasInvalidCell(cell_no));
- }
- }
- }
-#endif
+ // clear old content
+ tria.levels.clear();
+ tria.levels.push_back(
+ std::make_unique<
+ dealii::internal::TriangulationImplementation::TriaLevel>(dim));
+ if (dim > 1)
+ tria.faces = std::make_unique<
+ dealii::internal::TriangulationImplementation::TriaFaces>(dim);
- // store the indices of the lines
- // which are adjacent to a given
- // vertex
- std::vector<std::vector<int>> lines_at_vertex(v.size());
+ // copy vertices
+ tria.vertices = vertices;
+ tria.vertices_used.assign(vertices.size(), true);
- // reserve enough space
- triangulation.levels.push_back(
- std::make_unique<internal::TriangulationImplementation::TriaLevel>(
- dim));
- reserve_space(*triangulation.levels[0], cells.size(), dim, spacedim);
- reserve_space(triangulation.levels[0]->cells, 0, cells.size());
+ // compute connectivity
+ const auto connectivity = build_connectivity<unsigned int>(cells);
+ const unsigned int n_cell = cells.size();
- // make up cells
- typename Triangulation<dim, spacedim>::raw_line_iterator
- next_free_line = triangulation.begin_raw_line();
- for (unsigned int cell = 0; cell < cells.size(); ++cell)
+ // TriaObjects: lines
+ if (dim >= 2)
{
- while (next_free_line->used())
- ++next_free_line;
-
- next_free_line->set_bounding_object_indices(
- {cells[cell].vertices[0], cells[cell].vertices[1]});
- next_free_line->set_used_flag();
- next_free_line->set_material_id(cells[cell].material_id);
- next_free_line->set_manifold_id(cells[cell].manifold_id);
- next_free_line->clear_user_data();
- next_free_line->set_subdomain_id(0);
-
- // note that this cell is
- // adjacent to these vertices
- lines_at_vertex[cells[cell].vertices[0]].push_back(cell);
- lines_at_vertex[cells[cell].vertices[1]].push_back(cell);
- }
-
+ auto &lines_0 = tria.faces->lines; // data structure to be filled
- // some security tests
- {
- unsigned int boundary_nodes = 0;
- for (const auto &line : lines_at_vertex)
- switch (line.size())
- {
- case 1:
- // this vertex has only
- // one adjacent line
- ++boundary_nodes;
- break;
- case 2:
- break;
- default:
- AssertThrow(
- false,
- ExcMessage(
- "You have a vertex in your triangulation "
- "at which more than two cells come together. "
- "(For one dimensional triangulation, cells are "
- "line segments.)"
- "\n\n"
- "This is not currently supported because the "
- "Triangulation class makes the assumption that "
- "every cell has zero or one neighbors behind "
- "each face (here, behind each vertex), but in your "
- "situation there would be more than one."
- "\n\n"
- "Support for this is not currently implemented. "
- "If you need to work with triangulations where "
- "more than two cells come together at a vertex, "
- "duplicate the vertices once per cell (i.e., put "
- "multiple vertices at the same physical location, "
- "but using different vertex indices for each) "
- "and then ensure continuity of the solution by "
- "explicitly creating constraints that the degrees "
- "of freedom at these vertices have the same "
- "value, using the AffineConstraints class."));
- }
- }
-
-
-
- // update neighborship info
- typename Triangulation<dim, spacedim>::active_line_iterator line =
- triangulation.begin_active_line();
- // for all lines
- for (; line != triangulation.end(); ++line)
- // for each of the two vertices
- for (const unsigned int vertex : GeometryInfo<dim>::vertex_indices())
- // if first cell adjacent to
- // this vertex is the present
- // one, then the neighbor is
- // the second adjacent cell and
- // vice versa
- if (lines_at_vertex[line->vertex_index(vertex)][0] == line->index())
- if (lines_at_vertex[line->vertex_index(vertex)].size() == 2)
- {
- const typename Triangulation<dim, spacedim>::cell_iterator
- neighbor(&triangulation,
- 0, // level
- lines_at_vertex[line->vertex_index(vertex)][1]);
- line->set_neighbor(vertex, neighbor);
- }
- else
- // no second adjacent cell
- // entered -> cell at
- // boundary
- line->set_neighbor(vertex, triangulation.end());
- else
- // present line is not first
- // adjacent one -> first
- // adjacent one is neighbor
- {
- const typename Triangulation<dim, spacedim>::cell_iterator
- neighbor(&triangulation,
- 0, // level
- lines_at_vertex[line->vertex_index(vertex)][0]);
- line->set_neighbor(vertex, neighbor);
- }
+ // get connectivity between quads and lines
+ const auto & crs = connectivity.entity_to_entities(1, 0);
+ const unsigned int n_lines = crs.ptr.size() - 1;
- // finally set the
- // vertex_to_boundary_id_map_1d
- // and vertex_to_manifold_id_map_1d
- // maps
- triangulation.vertex_to_boundary_id_map_1d->clear();
- triangulation.vertex_to_manifold_id_map_1d->clear();
- for (const auto &cell : triangulation.active_cell_iterators())
- for (auto f : GeometryInfo<dim>::face_indices())
- {
- (*triangulation.vertex_to_manifold_id_map_1d)
- [cell->face(f)->vertex_index()] = numbers::flat_manifold_id;
+ // allocate memory
+ reserve_space_(lines_0, n_lines);
- if (cell->at_boundary(f))
- (*triangulation.vertex_to_boundary_id_map_1d)
- [cell->face(f)->vertex_index()] = f;
- }
- }
+ // loop over lines
+ for (unsigned int line = 0; line < n_lines; ++line)
+ for (unsigned int i = crs.ptr[line], j = 0; i < crs.ptr[line + 1];
+ ++i, ++j)
+ lines_0.cells[line * GeometryInfo<1>::faces_per_cell + j] =
+ crs.col[i]; // set vertex indices
+ }
+ // TriaObjects: quads
+ if (dim == 3)
+ {
+ auto &quads_0 = tria.faces->quads; // data structures to be filled
+ auto &faces = *tria.faces;
- /**
- * Create a triangulation from
- * given data. This function does
- * this work for 2-dimensional
- * triangulations independently
- * of the actual space dimension.
- */
- template <int spacedim>
- static void
- create_triangulation(const std::vector<Point<spacedim>> &v,
- const std::vector<CellData<2>> & cells,
- const SubCellData & subcelldata,
- Triangulation<2, spacedim> & triangulation)
- {
- AssertThrow(v.size() > 0, ExcMessage("No vertices given"));
- AssertThrow(cells.size() > 0, ExcMessage("No cells given"));
+ // get connectivity between quads and lines
+ const auto & crs = connectivity.entity_to_entities(2, 1);
+ const unsigned int n_quads = crs.ptr.size() - 1;
- const unsigned int dim = 2;
+ // allocate memory
+ reserve_space_(quads_0, n_quads);
+ reserve_space_(faces, 2 /*structdim*/, n_quads);
- // copy vertices
- triangulation.vertices = v;
- triangulation.vertices_used = std::vector<bool>(v.size(), true);
-
- // Check that all cells have positive volume. This check is not run in
- // the codimension one or two cases since cell_measure is not
- // implemented for those.
-#ifndef _MSC_VER
- // TODO: The following code does not compile with MSVC. Find a way
- // around it
- if (dim == spacedim)
- {
- for (unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
+ // loop over all quads -> entity type, line indices/orientations
+ for (unsigned int q = 0, k = 0; q < n_quads; ++q)
{
- // See the note in the 1D function on this if statement.
- if (!triangulation.check_for_distorted_cells)
+ // set entity type of quads
+ faces.quad_reference_cell_type[q] =
+ connectivity.entity_types(2)[q];
+
+ // loop over all its lines
+ for (unsigned int i = crs.ptr[q], j = 0; i < crs.ptr[q + 1];
+ ++i, ++j, ++k)
{
- unsigned int vertices[GeometryInfo<2>::vertices_per_cell];
-
- for (unsigned int i = 0;
- i < GeometryInfo<2>::vertices_per_cell;
- ++i)
- vertices[i] = cells[cell_no].vertices[i];
-
- const double cell_measure =
- GridTools::cell_measure<2>(triangulation.vertices,
- vertices);
- AssertThrow(cell_measure > 0,
- ExcGridHasInvalidCell(cell_no));
+ // set line index
+ quads_0.cells[q * GeometryInfo<2>::faces_per_cell + j] =
+ crs.col[i];
+
+ // set line orientations
+ faces.quads_line_orientations
+ [q * GeometryInfo<2>::faces_per_cell + j] =
+ connectivity.entity_orientations(1)[k];
}
}
}
-#endif
-
- // make up a list of the needed
- // lines each line is a pair of
- // vertices. The list is kept
- // sorted and it is guaranteed that
- // each line is inserted only once.
- // While the key of such an entry
- // is the pair of vertices, the
- // thing it points to is an
- // iterator pointing to the line
- // object itself. In the first run,
- // these iterators are all invalid
- // ones, but they are filled
- // afterwards
- std::map<std::pair<int, int>,
- typename Triangulation<dim, spacedim>::line_iterator>
- needed_lines;
- for (unsigned int cell = 0; cell < cells.size(); ++cell)
- {
- for (const auto vertex : cells[cell].vertices)
- AssertThrow(vertex < triangulation.vertices.size(),
- ExcInvalidVertexIndex(cell,
- vertex,
- triangulation.vertices.size()));
-
- for (const unsigned int line : GeometryInfo<dim>::face_indices())
- {
- // given a line vertex number (0,1) on a specific line
- // we get the cell vertex number (0-4) through the
- // line_to_cell_vertices function
- std::pair<int, int> line_vertices(
- cells[cell].vertices[GeometryInfo<dim>::line_to_cell_vertices(
- line, 0)],
- cells[cell].vertices[GeometryInfo<dim>::line_to_cell_vertices(
- line, 1)]);
-
- // assert that the line was not already inserted in
- // reverse order. This happens in spite of the vertex
- // rotation above, if the sense of the cell was
- // incorrect.
- //
- // Here is what usually happened when this exception
- // is thrown: consider these two cells and the
- // vertices
- // 3---4---5
- // | | |
- // 0---1---2
- // If in the input vector the two cells are given with
- // vertices <0 1 3 4> and <4 1 5 2>, in the first cell
- // the middle line would have direction 1->4, while in
- // the second it would be 4->1. This will cause the
- // exception.
- AssertThrow(needed_lines.find(std::make_pair(
- line_vertices.second, line_vertices.first)) ==
- needed_lines.end(),
- ExcGridHasInvalidCell(cell));
-
- // insert line, with
- // invalid iterator if line
- // already exists, then
- // nothing bad happens here
- needed_lines[line_vertices] = triangulation.end_line();
- }
- }
-
- // check that every vertex has at
- // least two adjacent lines
+ // TriaObjects/TriaLevel: cell
{
- std::vector<unsigned short int> vertex_touch_count(v.size(), 0);
- typename std::map<
- std::pair<int, int>,
- typename Triangulation<dim, spacedim>::line_iterator>::iterator i;
- for (i = needed_lines.begin(); i != needed_lines.end(); ++i)
+ auto &cells_0 = tria.levels[0]->cells; // data structure to be filled
+ auto &level = *tria.levels[0];
+
+ // get connectivity between cells/faces and cells/cells
+ const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
+ const auto &nei = connectivity.entity_to_entities(dim, dim);
+
+ // in 2D optional: since in in pure QUAD meshes same line
+ // orientations can be guaranteed
+ const bool orientation_needed =
+ dim == 3 ||
+ (dim == 2 &&
+ std::any_of(connectivity.entity_orientations(1).begin(),
+ connectivity.entity_orientations(1).end(),
+ [](const auto &i) { return i == 0; }));
+
+ // allocate memory
+ reserve_space_(cells_0, n_cell);
+ reserve_space_(level, spacedim, n_cell, orientation_needed);
+
+ // loop over all cells
+ for (unsigned int cell = 0; cell < n_cell; ++cell)
{
- // touch the vertices of
- // this line
- ++vertex_touch_count[i->first.first];
- ++vertex_touch_count[i->first.second];
- }
+ // set material ids
+ cells_0.boundary_or_material_id[cell].material_id =
+ cells[cell].material_id;
- // assert minimum touch count
- // is at least two. if not so,
- // then clean triangulation and
- // exit with an exception
- AssertThrow(*(std::min_element(vertex_touch_count.begin(),
- vertex_touch_count.end())) >= 2,
- ExcMessage(
- "During creation of a triangulation, a part of the "
- "algorithm encountered a vertex that is part of only "
- "a single adjacent line. However, in 2d, every vertex "
- "needs to be at least part of two lines."));
- }
+ // set manifold ids
+ cells_0.manifold_id[cell] = cells[cell].manifold_id;
- // reserve enough space
- triangulation.levels.push_back(
- std::make_unique<internal::TriangulationImplementation::TriaLevel>(
- dim));
- triangulation.faces =
- std::make_unique<internal::TriangulationImplementation::TriaFaces>(
- dim);
- reserve_space(*triangulation.levels[0], cells.size(), dim, spacedim);
- reserve_space(triangulation.faces->lines, 0, needed_lines.size());
- reserve_space(triangulation.levels[0]->cells, 0, cells.size());
-
- // make up lines
- {
- typename Triangulation<dim, spacedim>::raw_line_iterator line =
- triangulation.begin_raw_line();
- typename std::map<
- std::pair<int, int>,
- typename Triangulation<dim, spacedim>::line_iterator>::iterator i;
- for (i = needed_lines.begin(); line != triangulation.end_line();
- ++line, ++i)
- {
- line->set_bounding_object_indices(
- {i->first.first, i->first.second});
- line->set_used_flag();
- line->clear_user_flag();
- line->clear_user_data();
- i->second = line;
+ // set entity types
+ level.reference_cell_type[cell] =
+ connectivity.entity_types(dim)[cell];
+
+ // loop over faces
+ for (unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
+ ++i, ++j)
+ {
+ // set neighbor if not at boundary
+ if (nei.col[i] != static_cast<unsigned int>(-1))
+ level.neighbors[cell * GeometryInfo<dim>::faces_per_cell +
+ j] = {0, nei.col[i]};
+
+ // set face indices
+ cells_0.cells[cell * GeometryInfo<dim>::faces_per_cell + j] =
+ crs.col[i];
+
+ // set face orientation if needed
+ if (orientation_needed)
+ level.face_orientations
+ [cell * GeometryInfo<dim>::faces_per_cell + j] =
+ connectivity.entity_orientations(dim - 1)[i];
+ }
}
}
+ // TriaFaces: boundary id of boundary faces
+ if (dim > 1)
+ {
+ auto &bids_face = dim == 3 ?
+ tria.faces->quads.boundary_or_material_id :
+ tria.faces->lines.boundary_or_material_id;
- // store for each line index
- // the adjacent cells
- std::map<
- int,
- std::vector<typename Triangulation<dim, spacedim>::cell_iterator>>
- adjacent_cells;
+ // count number of cells a face is belonging to
+ std::vector<unsigned int> count(bids_face.size(), 0);
- // finally make up cells
- {
- typename Triangulation<dim, spacedim>::raw_cell_iterator cell =
- triangulation.begin_raw_quad();
- for (unsigned int c = 0; c < cells.size(); ++c, ++cell)
- {
- typename Triangulation<dim, spacedim>::line_iterator
- lines[GeometryInfo<dim>::lines_per_cell];
- for (unsigned int line = 0;
- line < GeometryInfo<dim>::lines_per_cell;
- ++line)
- lines[line] = needed_lines[std::make_pair(
- cells[c].vertices[GeometryInfo<dim>::line_to_cell_vertices(
- line, 0)],
- cells[c].vertices[GeometryInfo<dim>::line_to_cell_vertices(
- line, 1)])];
-
- cell->set_bounding_object_indices({lines[0]->index(),
- lines[1]->index(),
- lines[2]->index(),
- lines[3]->index()});
-
- cell->set_used_flag();
- cell->set_material_id(cells[c].material_id);
- cell->set_manifold_id(cells[c].manifold_id);
- cell->clear_user_data();
- cell->set_subdomain_id(0);
-
- // note that this cell is
- // adjacent to the four
- // lines
- for (const auto &line : lines)
- adjacent_cells[line->index()].push_back(cell);
- }
- }
+ // get connectivity between cells/faces
+ const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
+ // count how many cells are adjacent to the same face
+ for (unsigned int cell = 0; cell < cells.size(); ++cell)
+ for (unsigned int i = crs.ptr[cell]; i < crs.ptr[cell + 1]; ++i)
+ count[crs.col[i]]++;
- for (typename Triangulation<dim, spacedim>::line_iterator line =
- triangulation.begin_line();
- line != triangulation.end_line();
- ++line)
- {
- const unsigned int n_adj_cells =
- adjacent_cells[line->index()].size();
-
- // assert that every line has one or two adjacent cells.
- // this has to be the case for 2d triangulations in 2d.
- // in higher dimensions, this may happen but is not
- // implemented
- if (spacedim == 2)
+ // loop over all faces
+ for (unsigned int face = 0; face < count.size(); ++face)
{
- AssertThrow((n_adj_cells >= 1) && (n_adj_cells <= 2),
- ExcInternalError());
- }
- else
- {
- AssertThrow(
- (n_adj_cells >= 1) && (n_adj_cells <= 2),
- ExcMessage("You have a line in your triangulation at which "
- "more than two cells come together."
- "\n\n"
- "This is not currently supported because the "
- "Triangulation class makes the assumption that "
- "every cell has zero or one neighbors behind each "
- "face (here, behind each line), but in your "
- "situation there would be more than one."
- "\n\n"
- "Support for this is not currently implemented. "
- "If you need to work with triangulations where "
- "more than two cells come together at a line, "
- "duplicate the vertices once per cell (i.e., put "
- "multiple vertices at the same physical location, "
- "but using different vertex indices for each) "
- "and then ensure continuity of the solution by "
- "explicitly creating constraints that the degrees "
- "of freedom at these lines have the same "
- "value, using the AffineConstraints class."));
- }
+ if (count[face] != 1) // inner face
+ continue;
- // if only one cell: line is at boundary -> give it the boundary
- // indicator zero by default
- line->set_boundary_id_internal(
- (n_adj_cells == 1) ? 0 : numbers::internal_face_boundary_id);
- line->set_manifold_id(numbers::flat_manifold_id);
- }
+ // boundary faces ...
+ bids_face[face].boundary_id = 0;
- // set boundary indicators where given
- for (const auto &subcell_line : subcelldata.boundary_lines)
- {
- typename Triangulation<dim, spacedim>::line_iterator line;
- std::pair<int, int> line_vertices(
- std::make_pair(subcell_line.vertices[0],
- subcell_line.vertices[1]));
- if (needed_lines.find(line_vertices) != needed_lines.end())
- // line found in this direction
- line = needed_lines[line_vertices];
- else
- {
- // look whether it exists in reverse direction
- std::swap(line_vertices.first, line_vertices.second);
- if (needed_lines.find(line_vertices) != needed_lines.end())
- line = needed_lines[line_vertices];
- else
- // line does not exist
- AssertThrow(false,
- ExcLineInexistant(line_vertices.first,
- line_vertices.second));
- }
+ if (dim != 3)
+ continue;
- // assert that we only set boundary info once
- AssertThrow(!(line->boundary_id() != 0 &&
- line->boundary_id() !=
- numbers::internal_face_boundary_id),
- ExcMultiplySetLineInfoOfLine(line_vertices.first,
- line_vertices.second));
-
- // assert that the manifold id is not yet set or consistent
- // with the previous id
- AssertThrow(line->manifold_id() == numbers::flat_manifold_id ||
- line->manifold_id() == subcell_line.manifold_id,
- ExcInconsistentLineInfoOfLine(line_vertices.first,
- line_vertices.second,
- "manifold ids"));
- line->set_manifold_id(subcell_line.manifold_id);
-
- // assert that only exterior lines are given a boundary
- // indicator
- if (subcell_line.boundary_id != numbers::internal_face_boundary_id)
- {
- AssertThrow(
- line->boundary_id() != numbers::internal_face_boundary_id,
- ExcInteriorLineCantBeBoundary(line->vertex_index(0),
- line->vertex_index(1),
- subcell_line.boundary_id));
- line->set_boundary_id_internal(subcell_line.boundary_id);
+ // ... and the lines of quads in 3D
+ const auto &crs = connectivity.entity_to_entities(2, 1);
+ for (unsigned int i = crs.ptr[face]; i < crs.ptr[face + 1]; ++i)
+ tria.faces->lines.boundary_or_material_id[crs.col[i]]
+ .boundary_id = 0;
}
}
+ else // 1D
+ {
+ static const unsigned int t_tba = static_cast<unsigned int>(-1);
+ static const unsigned int t_inner = static_cast<unsigned int>(-2);
+ std::vector<unsigned int> type(vertices.size(), t_tba);
- // finally update neighborship info
- for (const auto &cell : triangulation.cell_iterators())
- for (unsigned int side = 0; side < 4; ++side)
- if (adjacent_cells[cell->line(side)->index()][0] == cell)
- // first adjacent cell is
- // this one
+ const auto &crs = connectivity.entity_to_entities(1, 0);
+
+ for (unsigned int cell = 0; cell < cells.size(); ++cell)
+ for (unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
+ ++i, ++j)
+ if (type[crs.col[i]] != t_inner)
+ type[crs.col[i]] = type[crs.col[i]] == t_tba ? j : t_inner;
+
+ for (unsigned int face = 0; face < type.size(); ++face)
{
- if (adjacent_cells[cell->line(side)->index()].size() == 2)
- // there is another
- // adjacent cell
- cell->set_neighbor(
- side, adjacent_cells[cell->line(side)->index()][1]);
+ // note: we also treat manifolds here!?
+ (*tria.vertex_to_manifold_id_map_1d)[face] =
+ numbers::flat_manifold_id;
+ if (type[face] != t_inner && type[face] != t_tba)
+ (*tria.vertex_to_boundary_id_map_1d)[face] = type[face];
}
- // first adjacent cell is not this
- // one, -> it must be the neighbor
- // we are looking for
- else
- cell->set_neighbor(side,
- adjacent_cells[cell->line(side)->index()][0]);
- }
+ }
+ // SubCellData: line
+ if (dim >= 2)
+ process_subcelldata(connectivity.entity_to_entities(1, 0),
+ tria.faces->lines,
+ subcelldata.boundary_lines);
- /**
- * Invent an object which compares two std::vector objects against each
- * other. This comparison is needed in order to establish a map of vertex
- * index tuples to iterators in the
- * Triangulation<3,3>::create_triangulation function.
- */
- struct QuadComparator
- {
- inline bool
- operator()(const std::vector<int> &q1, const std::vector<int> &q2) const
- {
- Assert(q1.size() == 4, ExcInternalError());
- Assert(q2.size() == 4, ExcInternalError());
-
- // here is room to
- // optimize the repeated
- // equality test of the
- // previous lines; the
- // compiler will probably
- // take care of most of
- // it anyway
- if ((q1[0] < q2[0]) || ((q1[0] == q2[0]) && (q1[1] < q2[1])) ||
- ((q1[0] == q2[0]) && (q1[1] == q2[1]) && (q1[2] < q2[2])) ||
- ((q1[0] == q2[0]) && (q1[1] == q2[1]) && (q1[2] == q2[2]) &&
- (q1[3] < q2[3])))
- return true;
- else
- return false;
- }
- };
+ // SubCellData: quad
+ if (dim == 3)
+ process_subcelldata(connectivity.entity_to_entities(2, 0),
+ tria.faces->quads,
+ subcelldata.boundary_quads);
+ }
- /**
- * Create a triangulation from
- * given data. This function does
- * this work for 3-dimensional
- * triangulations independently
- * of the actual space dimension.
- */
- template <int spacedim>
+ template <int structdim, typename T>
static void
- create_triangulation(const std::vector<Point<spacedim>> &v,
- const std::vector<CellData<3>> & cells,
- const SubCellData & subcelldata,
- Triangulation<3, spacedim> & triangulation)
+ process_subcelldata(
+ const CRS<T> & crs,
+ TriaObjects & obj,
+ const std::vector<CellData<structdim>> &boundary_objects_in)
{
- AssertThrow(v.size() > 0, ExcMessage("No vertices given"));
- AssertThrow(cells.size() > 0, ExcMessage("No cells given"));
+ AssertDimension(obj.structdim, structdim);
- const unsigned int dim = 3;
+ if (boundary_objects_in.size() == 0)
+ return; // empty subcelldata -> nothing to do
- // copy vertices
- triangulation.vertices = v;
- triangulation.vertices_used = std::vector<bool>(v.size(), true);
-
- // Check that all cells have positive volume.
-#ifndef _MSC_VER
- // TODO: The following code does not compile with MSVC. Find a way
- // around it
- for (unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
- {
- // See the note in the 1D function on this if statement.
- if (!triangulation.check_for_distorted_cells)
- {
- unsigned int vertices[GeometryInfo<3>::vertices_per_cell];
+ // pre-sort subcelldata
+ auto boundary_objects = boundary_objects_in;
- for (unsigned int i = 0; i < GeometryInfo<3>::vertices_per_cell;
- ++i)
- vertices[i] = cells[cell_no].vertices[i];
+ // ... sort vertices
+ for (auto &boundary_object : boundary_objects)
+ std::sort(boundary_object.vertices.begin(),
+ boundary_object.vertices.end());
- const double cell_measure =
- GridTools::cell_measure<3>(triangulation.vertices, vertices);
- AssertThrow(cell_measure > 0, ExcGridHasInvalidCell(cell_no));
- }
- }
-#endif
+ // ... sort cells
+ std::sort(boundary_objects.begin(),
+ boundary_objects.end(),
+ [](const auto &a, const auto &b) {
+ return a.vertices < b.vertices;
+ });
- ///////////////////////////////////////
- // first set up some collections of data
- //
- // make up a list of the needed
- // lines
- //
- // each line is a pair of
- // vertices. The list is kept
- // sorted and it is guaranteed that
- // each line is inserted only once.
- // While the key of such an entry
- // is the pair of vertices, the
- // thing it points to is an
- // iterator pointing to the line
- // object itself. In the first run,
- // these iterators are all invalid
- // ones, but they are filled
- // afterwards same applies for the
- // quads
- typename std::map<std::pair<int, int>,
- typename Triangulation<dim, spacedim>::line_iterator>
- needed_lines;
- for (unsigned int cell = 0; cell < cells.size(); ++cell)
+ unsigned int counter = 0;
+
+ std::vector<unsigned int> key;
+ key.reserve(GeometryInfo<structdim>::vertices_per_cell);
+
+ for (unsigned int o = 0; o < obj.n_objects(); ++o)
{
- // check whether vertex indices
- // are valid ones
- for (const auto vertex : cells[cell].vertices)
- AssertThrow(vertex < triangulation.vertices.size(),
- ExcInvalidVertexIndex(cell,
- vertex,
- triangulation.vertices.size()));
+ auto &boundary_id = obj.boundary_or_material_id[o].boundary_id;
+ auto &manifold_id = obj.manifold_id[o];
+
+ // assert that object has not been visited yet and its value
+ // has not been modified yet
+ AssertThrow(boundary_id == 0 ||
+ boundary_id == numbers::internal_face_boundary_id,
+ ExcNotImplemented());
+ AssertThrow(manifold_id == numbers::flat_manifold_id,
+ ExcNotImplemented());
+
+ // create key
+ key.assign(crs.col.data() + crs.ptr[o],
+ crs.col.data() + crs.ptr[o + 1]);
+ std::sort(key.begin(), key.end());
+
+ // is subcelldata provided? -> binary search
+ const auto subcell_object =
+ std::lower_bound(boundary_objects.begin(),
+ boundary_objects.end(),
+ key,
+ [&](const auto &cell, const auto &key) {
+ return cell.vertices < key;
+ });
+
+ // no subcelldata provided for this object
+ if (subcell_object == boundary_objects.end() ||
+ subcell_object->vertices != key)
+ continue;
- for (unsigned int line = 0;
- line < GeometryInfo<dim>::lines_per_cell;
- ++line)
+ counter++;
+
+ // set manifold id
+ manifold_id = subcell_object->manifold_id;
+
+ // set boundary id
+ if (subcell_object->boundary_id !=
+ numbers::internal_face_boundary_id)
{
- // given a line vertex number
- // (0,1) on a specific line we
- // get the cell vertex number
- // (0-7) through the
- // line_to_cell_vertices
- // function
- std::pair<int, int> line_vertices(
- cells[cell].vertices[GeometryInfo<dim>::line_to_cell_vertices(
- line, 0)],
- cells[cell].vertices[GeometryInfo<dim>::line_to_cell_vertices(
- line, 1)]);
-
- // if that line was already inserted
- // in reverse order do nothing, else
- // insert the line
- if ((needed_lines.find(std::make_pair(line_vertices.second,
- line_vertices.first)) ==
- needed_lines.end()))
- {
- // insert line, with
- // invalid iterator. if line
- // already exists, then
- // nothing bad happens here
- needed_lines[line_vertices] = triangulation.end_line();
- }
+ AssertThrow(boundary_id != numbers::internal_face_boundary_id,
+ ExcNotImplemented());
+ boundary_id = subcell_object->boundary_id;
}
}
+ // make sure that all subcelldata entries have been processed
+ // TODO: this is not guaranteed, why?
+ // AssertDimension(counter, boundary_objects_in.size());
+ }
- /////////////////////////////////
- // now for some sanity-checks:
- //
- // check that every vertex has at
- // least tree adjacent lines
- {
- std::vector<unsigned short int> vertex_touch_count(v.size(), 0);
- typename std::map<
- std::pair<int, int>,
- typename Triangulation<dim, spacedim>::line_iterator>::iterator i;
- for (i = needed_lines.begin(); i != needed_lines.end(); ++i)
- {
- // touch the vertices of
- // this line
- ++vertex_touch_count[i->first.first];
- ++vertex_touch_count[i->first.second];
- }
-
- // assert minimum touch count
- // is at least three. if not so,
- // then clean triangulation and
- // exit with an exception
- AssertThrow(
- *(std::min_element(vertex_touch_count.begin(),
- vertex_touch_count.end())) >= 3,
- ExcMessage(
- "During creation of a triangulation, a part of the "
- "algorithm encountered a vertex that is part of only "
- "one or two adjacent lines. However, in 3d, every vertex "
- "needs to be at least part of three lines."));
- }
- ///////////////////////////////////
- // actually set up data structures
- // for the lines
- // reserve enough space
- triangulation.levels.push_back(
- std::make_unique<internal::TriangulationImplementation::TriaLevel>(
- dim));
- triangulation.faces =
- std::make_unique<internal::TriangulationImplementation::TriaFaces>(
- dim);
- reserve_space(*triangulation.levels[0], cells.size(), dim, spacedim);
- reserve_space(triangulation.faces->lines, 0, needed_lines.size());
-
- // make up lines
- {
- typename Triangulation<dim, spacedim>::raw_line_iterator line =
- triangulation.begin_raw_line();
- typename std::map<
- std::pair<int, int>,
- typename Triangulation<dim, spacedim>::line_iterator>::iterator i;
- for (i = needed_lines.begin(); line != triangulation.end_line();
- ++line, ++i)
- {
- line->set_bounding_object_indices(
- {i->first.first, i->first.second});
- line->set_used_flag();
- line->clear_user_flag();
- line->clear_user_data();
-
- // now set the iterator for
- // this line
- i->second = line;
- }
- }
+ static void
+ reserve_space_(TriaFaces & faces,
+ const unsigned structdim,
+ const unsigned int size)
+ {
+ const unsigned int dim = faces.dim;
+ const unsigned int faces_per_cell =
+ structdim == 1 ? GeometryInfo<1>::faces_per_cell :
+ (structdim == 2 ? GeometryInfo<2>::faces_per_cell :
+ GeometryInfo<3>::faces_per_cell);
- ///////////////////////////////////////////
- // make up the quads of this triangulation
- //
- // same thing: the iterators are
- // set to the invalid value at
- // first, we only collect the data
- // now
-
- // the bool array stores, whether the lines
- // are in the standard orientation or not
-
- // note that QuadComparator is a
- // class declared and defined in
- // this file
- std::map<std::vector<int>,
- std::pair<typename Triangulation<dim, spacedim>::quad_iterator,
- std::array<bool, GeometryInfo<dim>::lines_per_face>>,
- QuadComparator>
- needed_quads;
- for (const auto &cell : cells)
+ if (dim == 3 && structdim == 2)
{
- // the faces are quads which
- // consist of four numbers
- // denoting the index of the
- // four lines bounding the
- // quad. we can get this index
- // by asking @p{needed_lines}
- // for an iterator to this
- // line, dereferencing it and
- // thus return an iterator into
- // the @p{lines} array of the
- // triangulation, which is
- // already set up. we can then
- // ask this iterator for its
- // index within the present
- // level (the level is zero, of
- // course)
- //
- // to make things easier, we
- // don't create the lines
- // (pairs of their vertex
- // indices) in place, but
- // before they are really
- // needed.
- std::pair<int, int> line_list[GeometryInfo<dim>::lines_per_cell],
- inverse_line_list[GeometryInfo<dim>::lines_per_cell];
- unsigned int face_line_list[GeometryInfo<dim>::lines_per_face];
- std::array<bool, GeometryInfo<dim>::lines_per_face> orientation;
-
- for (unsigned int line = 0;
- line < GeometryInfo<dim>::lines_per_cell;
- ++line)
- {
- line_list[line] = std::pair<int, int>(
- cell.vertices[GeometryInfo<dim>::line_to_cell_vertices(line,
- 0)],
- cell.vertices[GeometryInfo<dim>::line_to_cell_vertices(line,
- 1)]);
- inverse_line_list[line] = std::pair<int, int>(
- cell.vertices[GeometryInfo<dim>::line_to_cell_vertices(line,
- 1)],
- cell.vertices[GeometryInfo<dim>::line_to_cell_vertices(line,
- 0)]);
- }
-
- for (const unsigned int face : GeometryInfo<dim>::face_indices())
- {
- // set up a list of the lines to be
- // used for this face. check the
- // direction for each line
- //
- // given a face line number (0-3) on
- // a specific face we get the cell
- // line number (0-11) through the
- // face_to_cell_lines function
- for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_face;
- ++l)
- if (needed_lines.find(
- inverse_line_list[GeometryInfo<dim>::face_to_cell_lines(
- face, l)]) == needed_lines.end())
- {
- face_line_list[l] =
- needed_lines[line_list[GeometryInfo<
- dim>::face_to_cell_lines(face, l)]]
- ->index();
- orientation[l] = true;
- }
- else
- {
- face_line_list[l] =
- needed_lines[inverse_line_list[GeometryInfo<
- dim>::face_to_cell_lines(face, l)]]
- ->index();
- orientation[l] = false;
- }
+ // quad entity types
+ faces.quad_reference_cell_type.assign(size,
+ ReferenceCell::Type::Invalid);
-
- const std::vector<int> quad(
- {static_cast<int>(face_line_list[0]),
- static_cast<int>(face_line_list[1]),
- static_cast<int>(face_line_list[2]),
- static_cast<int>(face_line_list[3])});
-
- // insert quad, with invalid iterator
- //
- // if quad already exists, then nothing bad happens here, as
- // this will then simply become an interior face of the
- // triangulation. however, we will run into major trouble if the
- // face was already inserted in the opposite direction. there
- // are really only two orientations for a face to be in, since
- // the edge directions are already set. thus, vertex 0 is the
- // one from which two edges originate, and vertex 3 is the one
- // to which they converge. we are then left with orientations
- // 0-1-2-3 and 2-3-0-1 for the order of lines. the corresponding
- // quad can be easily constructed by exchanging lines. we do so
- // here, just to check that that flipped quad isn't already in
- // the triangulation. if it is, then don't insert the new one
- // and instead later set the face_orientation flag
-
- // face_orientation=false, face_flip=false, face_rotation=false
- const std::vector<int> test_quad_1(
- {quad[2], quad[3], quad[0], quad[1]}),
- // face_orientation=false, face_flip=false, face_rotation=true
- test_quad_2({quad[0], quad[1], quad[3], quad[2]}),
- // face_orientation=false, face_flip=true, face_rotation=false
- test_quad_3({quad[3], quad[2], quad[1], quad[0]}),
- // face_orientation=false, face_flip=true, face_rotation=true
- test_quad_4({quad[1], quad[0], quad[2], quad[3]}),
- // face_orientation=true, face_flip=false, face_rotation=true
- test_quad_5({quad[2], quad[3], quad[1], quad[0]}),
- // face_orientation=true, face_flip=true, face_rotation=false
- test_quad_6({quad[1], quad[0], quad[3], quad[2]}),
- // face_orientation=true, face_flip=true, face_rotation=true
- test_quad_7({quad[3], quad[2], quad[0], quad[1]});
-
- if (needed_quads.find(test_quad_1) == needed_quads.end() &&
- needed_quads.find(test_quad_2) == needed_quads.end() &&
- needed_quads.find(test_quad_3) == needed_quads.end() &&
- needed_quads.find(test_quad_4) == needed_quads.end() &&
- needed_quads.find(test_quad_5) == needed_quads.end() &&
- needed_quads.find(test_quad_6) == needed_quads.end() &&
- needed_quads.find(test_quad_7) == needed_quads.end())
- needed_quads[quad] =
- std::make_pair(triangulation.end_quad(), orientation);
- }
+ // quad line orientations
+ faces.quads_line_orientations.assign(size * faces_per_cell, -1);
}
+ }
- /////////////////////////////////
- // enter the resulting quads into
- // the arrays of the Triangulation
- //
- // first reserve enough space
- reserve_space(*triangulation.faces, 0, needed_quads.size());
- reserve_space(triangulation.faces->quads, 0, needed_quads.size());
-
- {
- typename Triangulation<dim, spacedim>::raw_quad_iterator quad =
- triangulation.begin_raw_quad();
- typename std::map<
- std::vector<int>,
- std::pair<typename Triangulation<dim, spacedim>::quad_iterator,
- std::array<bool, GeometryInfo<dim>::lines_per_face>>,
- QuadComparator>::iterator q;
- for (q = needed_quads.begin(); quad != triangulation.end_quad();
- ++quad, ++q)
- {
- quad->set_bounding_object_indices(
- {q->first[0], q->first[1], q->first[2], q->first[3]});
- quad->set_used_flag();
- quad->clear_user_flag();
- quad->clear_user_data();
- // set the line orientation
- quad->set_line_orientation(0, q->second.second[0]);
- quad->set_line_orientation(1, q->second.second[1]);
- quad->set_line_orientation(2, q->second.second[2]);
- quad->set_line_orientation(3, q->second.second[3]);
-
-
- // now set the iterator for
- // this quad
- q->second.first = quad;
- }
- }
-
- /////////////////////////////////
- // finally create the cells
- reserve_space(triangulation.levels[0]->cells, cells.size());
- // store for each quad index the
- // adjacent cells
- std::map<
- int,
- std::vector<typename Triangulation<dim, spacedim>::cell_iterator>>
- adjacent_cells;
+ static void
+ reserve_space_(TriaLevel & level,
+ const unsigned int spacedim,
+ const unsigned int size,
+ const bool orientation_needed)
+ {
+ const unsigned int dim = level.dim;
- // finally make up cells
- {
- typename Triangulation<dim, spacedim>::raw_cell_iterator cell =
- triangulation.begin_raw_hex();
- for (unsigned int c = 0; c < cells.size(); ++c, ++cell)
- {
- // first find for each of
- // the cells the quad
- // iterator of the
- // respective faces.
- //
- // to this end, set up the
- // lines of this cell and
- // find the quads that are
- // bounded by these lines;
- // these are then the faces
- // of the present cell
- std::pair<int, int> line_list[GeometryInfo<dim>::lines_per_cell],
- inverse_line_list[GeometryInfo<dim>::lines_per_cell];
- unsigned int face_line_list[4];
- for (unsigned int line = 0;
- line < GeometryInfo<dim>::lines_per_cell;
- ++line)
- {
- line_list[line] = std::make_pair(
- cells[c].vertices[GeometryInfo<dim>::line_to_cell_vertices(
- line, 0)],
- cells[c].vertices[GeometryInfo<dim>::line_to_cell_vertices(
- line, 1)]);
- inverse_line_list[line] = std::pair<int, int>(
- cells[c].vertices[GeometryInfo<dim>::line_to_cell_vertices(
- line, 1)],
- cells[c].vertices[GeometryInfo<dim>::line_to_cell_vertices(
- line, 0)]);
- }
+ const unsigned int faces_per_cell =
+ dim == 1 ? GeometryInfo<1>::faces_per_cell :
+ (dim == 2 ? GeometryInfo<2>::faces_per_cell :
+ GeometryInfo<3>::faces_per_cell);
- // get the iterators
- // corresponding to the
- // faces. also store
- // whether they are
- // reversed or not
- typename Triangulation<dim, spacedim>::quad_iterator
- face_iterator[GeometryInfo<dim>::faces_per_cell];
- bool face_orientation[GeometryInfo<dim>::faces_per_cell];
- bool face_flip[GeometryInfo<dim>::faces_per_cell];
- bool face_rotation[GeometryInfo<dim>::faces_per_cell];
- for (const unsigned int face : GeometryInfo<dim>::face_indices())
- {
- for (unsigned int l = 0;
- l < GeometryInfo<dim>::lines_per_face;
- ++l)
- if (needed_lines.find(inverse_line_list[GeometryInfo<
- dim>::face_to_cell_lines(face, l)]) ==
- needed_lines.end())
- face_line_list[l] =
- needed_lines[line_list[GeometryInfo<
- dim>::face_to_cell_lines(face, l)]]
- ->index();
- else
- face_line_list[l] =
- needed_lines[inverse_line_list[GeometryInfo<
- dim>::face_to_cell_lines(face, l)]]
- ->index();
-
- const std::vector<int> quad(
- {static_cast<int>(face_line_list[0]),
- static_cast<int>(face_line_list[1]),
- static_cast<int>(face_line_list[2]),
- static_cast<int>(face_line_list[3])});
-
- if (needed_quads.find(quad) != needed_quads.end())
- {
- // face is in standard
- // orientation (and not
- // flipped or rotated). this
- // must be true for at least
- // one of the two cells
- // containing this face
- // (i.e. for the cell which
- // originally inserted the
- // face)
- face_iterator[face] = needed_quads[quad].first;
- face_orientation[face] = true;
- face_flip[face] = false;
- face_rotation[face] = false;
- }
- else
- {
- // face must be available in reverse order then. construct
- // all possibilities and check them one after the other
-
- // face_orientation=false, face_flip=false,
- // face_rotation=false
- const std::vector<int> test_quad_1(
- {quad[2], quad[3], quad[0], quad[1]}),
- // face_orientation=false, face_flip=false,
- // face_rotation=true
- test_quad_2({quad[0], quad[1], quad[3], quad[2]}),
- // face_orientation=false, face_flip=true,
- // face_rotation=false
- test_quad_3({quad[3], quad[2], quad[1], quad[0]}),
- // face_orientation=false, face_flip=true,
- // face_rotation=true
- test_quad_4({quad[1], quad[0], quad[2], quad[3]}),
- // face_orientation=true, face_flip=false,
- // face_rotation=true
- test_quad_5({quad[2], quad[3], quad[1], quad[0]}),
- // face_orientation=true, face_flip=true,
- // face_rotation=false
- test_quad_6({quad[1], quad[0], quad[3], quad[2]}),
- // face_orientation=true, face_flip=true,
- // face_rotation=true
- test_quad_7({quad[3], quad[2], quad[0], quad[1]});
-
- if (needed_quads.find(test_quad_1) != needed_quads.end())
- {
- face_iterator[face] = needed_quads[test_quad_1].first;
- face_orientation[face] = false;
- face_flip[face] = false;
- face_rotation[face] = false;
- }
- else if (needed_quads.find(test_quad_2) !=
- needed_quads.end())
- {
- face_iterator[face] = needed_quads[test_quad_2].first;
- face_orientation[face] = false;
- face_flip[face] = false;
- face_rotation[face] = true;
- }
- else if (needed_quads.find(test_quad_3) !=
- needed_quads.end())
- {
- face_iterator[face] = needed_quads[test_quad_3].first;
- face_orientation[face] = false;
- face_flip[face] = true;
- face_rotation[face] = false;
- }
- else if (needed_quads.find(test_quad_4) !=
- needed_quads.end())
- {
- face_iterator[face] = needed_quads[test_quad_4].first;
- face_orientation[face] = false;
- face_flip[face] = true;
- face_rotation[face] = true;
- }
- else if (needed_quads.find(test_quad_5) !=
- needed_quads.end())
- {
- face_iterator[face] = needed_quads[test_quad_5].first;
- face_orientation[face] = true;
- face_flip[face] = false;
- face_rotation[face] = true;
- }
- else if (needed_quads.find(test_quad_6) !=
- needed_quads.end())
- {
- face_iterator[face] = needed_quads[test_quad_6].first;
- face_orientation[face] = true;
- face_flip[face] = true;
- face_rotation[face] = false;
- }
- else if (needed_quads.find(test_quad_7) !=
- needed_quads.end())
- {
- face_iterator[face] = needed_quads[test_quad_7].first;
- face_orientation[face] = true;
- face_flip[face] = true;
- face_rotation[face] = true;
- }
+ level.active_cell_indices.assign(size, -1);
+ level.subdomain_ids.assign(size, 0);
+ level.level_subdomain_ids.assign(size, 0);
- else
- // we didn't find the
- // face in any direction,
- // so something went
- // wrong above
- Assert(false, ExcInternalError());
- }
- } // for all faces
-
- // make the cell out of
- // these iterators
- cell->set_bounding_object_indices({face_iterator[0]->index(),
- face_iterator[1]->index(),
- face_iterator[2]->index(),
- face_iterator[3]->index(),
- face_iterator[4]->index(),
- face_iterator[5]->index()});
-
- cell->set_used_flag();
- cell->set_material_id(cells[c].material_id);
- cell->set_manifold_id(cells[c].manifold_id);
- cell->clear_user_flag();
- cell->clear_user_data();
- cell->set_subdomain_id(0);
-
- // set orientation flag for
- // each of the faces
- for (const unsigned int quad : GeometryInfo<dim>::face_indices())
- {
- cell->set_face_orientation(quad, face_orientation[quad]);
- cell->set_face_flip(quad, face_flip[quad]);
- cell->set_face_rotation(quad, face_rotation[quad]);
- }
+ level.refine_flags.assign(size, false);
+ level.coarsen_flags.assign(size, false);
+ level.parents.assign((size + 1) / 2, -1);
- // note that this cell is
- // adjacent to the six
- // quads
- for (const auto &quad : face_iterator)
- adjacent_cells[quad->index()].push_back(cell);
+ if (dim < spacedim)
+ level.direction_flags.assign(size, true);
-#ifdef DEBUG
- // make some checks on the
- // lines and their
- // ordering
-
- // first map all cell lines
- // to the two face lines
- // which should
- // coincide. all face lines
- // are included with a cell
- // line number (0-11)
- // key. At the end all keys
- // will be included twice
- // (for each of the two
- // coinciding lines once)
- std::multimap<unsigned int, std::pair<unsigned int, unsigned int>>
- cell_to_face_lines;
- for (const unsigned int face : GeometryInfo<dim>::face_indices())
- for (unsigned int line = 0;
- line < GeometryInfo<dim>::lines_per_face;
- ++line)
- cell_to_face_lines.insert(
- std::pair<unsigned int,
- std::pair<unsigned int, unsigned int>>(
- GeometryInfo<dim>::face_to_cell_lines(face, line),
- std::pair<unsigned int, unsigned int>(face, line)));
- std::multimap<unsigned int,
- std::pair<unsigned int, unsigned int>>::
- const_iterator map_iter = cell_to_face_lines.begin();
-
- for (; map_iter != cell_to_face_lines.end(); ++map_iter)
- {
- const unsigned int cell_line = map_iter->first;
- const unsigned int face1 = map_iter->second.first;
- const unsigned int line1 = map_iter->second.second;
- ++map_iter;
- Assert(map_iter != cell_to_face_lines.end(),
- ExcInternalErrorOnCell(c));
- Assert(map_iter->first == cell_line,
- ExcInternalErrorOnCell(c));
- const unsigned int face2 = map_iter->second.first;
- const unsigned int line2 = map_iter->second.second;
-
- // check that the pair
- // of lines really
- // coincide. Take care
- // about the face
- // orientation;
- Assert(face_iterator[face1]->line(
- GeometryInfo<dim>::standard_to_real_face_line(
- line1,
- face_orientation[face1],
- face_flip[face1],
- face_rotation[face1])) ==
- face_iterator[face2]->line(
- GeometryInfo<dim>::standard_to_real_face_line(
- line2,
- face_orientation[face2],
- face_flip[face2],
- face_rotation[face2])),
- ExcInternalErrorOnCell(c));
- }
-#endif
- }
- }
+ level.neighbors.assign(size * faces_per_cell, {-1, -1});
+ level.reference_cell_type.assign(size, ReferenceCell::Type::Invalid);
- /////////////////////////////////////////
- // find those quads which are at the
- // boundary and mark them appropriately
- for (typename Triangulation<dim, spacedim>::quad_iterator quad =
- triangulation.begin_quad();
- quad != triangulation.end_quad();
- ++quad)
- {
- const unsigned int n_adj_cells =
- adjacent_cells[quad->index()].size();
- // assert that every quad has
- // one or two adjacent cells
- AssertThrow((n_adj_cells >= 1) && (n_adj_cells <= 2),
- ExcInternalError());
-
- // if only one cell: quad is at boundary -> give it the boundary
- // indicator zero by default
- quad->set_boundary_id_internal(
- (n_adj_cells == 1) ? 0 : numbers::internal_face_boundary_id);
-
- // Manifold ids are set independently of where they are
- quad->set_manifold_id(numbers::flat_manifold_id);
- }
+ if (orientation_needed)
+ level.face_orientations.assign(size * faces_per_cell, -1);
+ }
- /////////////////////////////////////////
- // next find those lines which are at
- // the boundary and mark all others as
- // interior ones
- //
- // for this: first mark all lines as interior. use this loop
- // to also set all manifold ids of all lines
- for (typename Triangulation<dim, spacedim>::line_iterator line =
- triangulation.begin_line();
- line != triangulation.end_line();
- ++line)
- {
- line->set_boundary_id_internal(numbers::internal_face_boundary_id);
- line->set_manifold_id(numbers::flat_manifold_id);
- }
- // next reset all lines bounding
- // boundary quads as on the
- // boundary also. note that since
- // we are in 3d, there are cases
- // where one or more lines of a
- // quad that is not on the
- // boundary, are actually boundary
- // lines. they will not be marked
- // when visiting this
- // face. however, since we do not
- // support dim-2 dimensional
- // boundaries (i.e. internal lines
- // constituting boundaries), every
- // such line is also part of a face
- // that is actually on the
- // boundary, so sooner or later we
- // get to mark that line for being
- // on the boundary
- for (typename Triangulation<dim, spacedim>::quad_iterator quad =
- triangulation.begin_quad();
- quad != triangulation.end_quad();
- ++quad)
- if (quad->at_boundary())
- {
- for (unsigned int l = 0; l < 4; ++l)
- {
- typename Triangulation<dim, spacedim>::line_iterator line =
- quad->line(l);
- line->set_boundary_id_internal(0);
- }
- }
- ///////////////////////////////////////
- // now set boundary indicators
- // where given
- //
- // first do so for lines
- for (const auto &subcell_line : subcelldata.boundary_lines)
+ static void
+ reserve_space_(TriaObjects &obj, const unsigned int size)
+ {
+ const unsigned int structdim = obj.structdim;
+
+ const unsigned int max_children_per_cell =
+ structdim == 1 ?
+ GeometryInfo<1>::max_children_per_cell :
+ (structdim == 2 ? GeometryInfo<2>::max_children_per_cell :
+ GeometryInfo<3>::max_children_per_cell);
+ const unsigned int faces_per_cell =
+ structdim == 1 ? GeometryInfo<1>::faces_per_cell :
+ (structdim == 2 ? GeometryInfo<2>::faces_per_cell :
+ GeometryInfo<3>::faces_per_cell);
+
+ obj.used.assign(size, true);
+ obj.boundary_or_material_id.assign(
+ size,
+ internal::TriangulationImplementation::TriaObjects::
+ BoundaryOrMaterialId());
+ obj.manifold_id.assign(size, -1);
+ obj.user_flags.assign(size, false);
+ obj.user_data.resize(size);
+
+ if (structdim > 1) // TODO: why?
+ obj.refinement_cases.assign(size, 0);
+
+ obj.children.assign(max_children_per_cell / 2 * size, -1);
+
+ obj.cells.assign(faces_per_cell * size, -1);
+
+ if (structdim <= 2)
{
- typename Triangulation<dim, spacedim>::line_iterator line;
- std::pair<int, int> line_vertices(
- std::make_pair(subcell_line.vertices[0],
- subcell_line.vertices[1]));
- if (needed_lines.find(line_vertices) != needed_lines.end())
- // line found in this
- // direction
- line = needed_lines[line_vertices];
-
- else
- {
- // look whether it exists in
- // reverse direction
- std::swap(line_vertices.first, line_vertices.second);
- if (needed_lines.find(line_vertices) != needed_lines.end())
- line = needed_lines[line_vertices];
- else
- // line does not exist
- AssertThrow(false,
- ExcLineInexistant(line_vertices.first,
- line_vertices.second));
- }
- // Only exterior lines can be given a boundary indicator
- if (line->at_boundary())
- {
- // make sure that we don't attempt to reset the boundary
- // indicator to a different than the previously set value
- AssertThrow(line->boundary_id() == 0 ||
- line->boundary_id() == subcell_line.boundary_id,
- ExcInconsistentLineInfoOfLine(line_vertices.first,
- line_vertices.second,
- "boundary ids"));
- // If the boundary id provided in subcell_line
- // is anything other than the default
- // (internal_face_boundary_id), then set it in the new
- // triangulation.
- if (subcell_line.boundary_id !=
- numbers::internal_face_boundary_id)
- line->set_boundary_id(subcell_line.boundary_id);
- }
- // Set manifold id if given
- AssertThrow(line->manifold_id() == numbers::flat_manifold_id ||
- line->manifold_id() == subcell_line.manifold_id,
- ExcInconsistentLineInfoOfLine(line_vertices.first,
- line_vertices.second,
- "manifold ids"));
- line->set_manifold_id(subcell_line.manifold_id);
+ obj.next_free_single = size - 1;
+ obj.next_free_pair = 0;
+ obj.reverse_order_next_free_single = true;
}
-
-
- // now go on with the faces
- for (const auto &subcell_quad : subcelldata.boundary_quads)
+ else
{
- typename Triangulation<dim, spacedim>::quad_iterator quad;
- typename Triangulation<dim, spacedim>::line_iterator line[4];
-
- // first find the lines that
- // are made up of the given
- // vertices, then build up a
- // quad from these lines
- // finally use the find
- // function of the map template
- // to find the quad
- for (unsigned int i = 0; i < 4; ++i)
- {
- std::pair<int, int> line_vertices(
- subcell_quad
- .vertices[GeometryInfo<dim - 1>::line_to_cell_vertices(i,
- 0)],
- subcell_quad
- .vertices[GeometryInfo<dim - 1>::line_to_cell_vertices(i,
- 1)]);
-
- // check whether line
- // already exists
- if (needed_lines.find(line_vertices) != needed_lines.end())
- line[i] = needed_lines[line_vertices];
- else
- // look whether it exists
- // in reverse direction
- {
- std::swap(line_vertices.first, line_vertices.second);
- if (needed_lines.find(line_vertices) != needed_lines.end())
- line[i] = needed_lines[line_vertices];
- else
- // line does
- // not exist
- AssertThrow(false,
- ExcLineInexistant(line_vertices.first,
- line_vertices.second));
- }
- }
-
-
- // Set up 2 quads that are
- // built up from the lines for
- // reasons of comparison to
- // needed_quads. The second
- // quad is the reversed version
- // of the first quad in order
- // find the quad regardless of
- // its orientation. This is
- // introduced for convenience
- // and because boundary quad
- // orientation does not carry
- // any information.
- std::vector<int> quad_compare_1({line[0]->index(),
- line[1]->index(),
- line[2]->index(),
- line[3]->index()});
- std::vector<int> quad_compare_2({line[2]->index(),
- line[3]->index(),
- line[0]->index(),
- line[1]->index()});
-
- // try to find the quad with
- // lines situated as
- // constructed above. if it
- // could not be found, rotate
- // the boundary lines 3 times
- // until it is found or it does
- // not exist.
-
- // mapping from counterclock to
- // lexicographic ordering of
- // quad lines
- static const unsigned int lex2cclock[4] = {3, 1, 0, 2};
- // copy lines from
- // lexicographic to
- // counterclock ordering, as
- // rotation is much simpler in
- // counterclock ordering
- typename Triangulation<dim, spacedim>::line_iterator
- line_counterclock[4];
- for (unsigned int i = 0; i < 4; ++i)
- line_counterclock[lex2cclock[i]] = line[i];
- unsigned int n_rotations = 0;
- bool not_found_quad_1;
- while ((not_found_quad_1 = (needed_quads.find(quad_compare_1) ==
- needed_quads.end())) &&
- (needed_quads.find(quad_compare_2) == needed_quads.end()) &&
- (n_rotations < 4))
- {
- // use the rotate defined
- // in <algorithms>
- std::rotate(line_counterclock,
- line_counterclock + 1,
- line_counterclock + 4);
- // update the quads with
- // rotated lines (i runs in
- // lexicographic ordering)
- for (unsigned int i = 0; i < 4; ++i)
- {
- quad_compare_1[i] =
- line_counterclock[lex2cclock[i]]->index();
- quad_compare_2[(i + 2) % 4] =
- line_counterclock[lex2cclock[i]]->index();
- }
-
- ++n_rotations;
- }
-
- AssertThrow(n_rotations != 4,
- ExcQuadInexistant(line[0]->index(),
- line[1]->index(),
- line[2]->index(),
- line[3]->index()));
-
- if (not_found_quad_1)
- quad = needed_quads[quad_compare_2].first;
- else
- quad = needed_quads[quad_compare_1].first;
-
- // check whether this face is
- // really an exterior one
- if (quad->at_boundary())
- {
- // and make sure that we don't attempt to reset the boundary
- // indicator to a different than the previously set value
- AssertThrow(quad->boundary_id() == 0 ||
- quad->boundary_id() == subcell_quad.boundary_id,
- ExcInconsistentQuadInfoOfQuad(line[0]->index(),
- line[1]->index(),
- line[2]->index(),
- line[3]->index(),
- "boundary ids"));
- // If the boundary id provided in subcell_line
- // is anything other than the default
- // (internal_face_boundary_id), then set it in the new
- // triangulation.
- if (subcell_quad.boundary_id !=
- numbers::internal_face_boundary_id)
- quad->set_boundary_id(subcell_quad.boundary_id);
- }
- // Set manifold id if given
- if (quad->manifold_id() != numbers::flat_manifold_id)
- AssertThrow(quad->manifold_id() == subcell_quad.manifold_id,
- ExcInconsistentQuadInfoOfQuad(line[0]->index(),
- line[1]->index(),
- line[2]->index(),
- line[3]->index(),
- "manifold ids"));
-
- quad->set_manifold_id(subcell_quad.manifold_id);
+ obj.next_free_single = obj.next_free_pair = 0;
}
-
-
- /////////////////////////////////////////
- // finally update neighborship info
- for (const auto &cell : triangulation.cell_iterators())
- for (unsigned int face = 0; face < 6; ++face)
- if (adjacent_cells[cell->quad(face)->index()][0] == cell)
- // first adjacent cell is
- // this one
- {
- if (adjacent_cells[cell->quad(face)->index()].size() == 2)
- // there is another
- // adjacent cell
- cell->set_neighbor(
- face, adjacent_cells[cell->quad(face)->index()][1]);
- }
- // first adjacent cell is not this
- // one, -> it must be the neighbor
- // we are looking for
- else
- cell->set_neighbor(face,
- adjacent_cells[cell->quad(face)->index()][0]);
}
// because sometimes other objects are already attached to it:
try
{
- const bool arbitray_mesh_provided =
- std::any_of(cells.begin(), cells.end(), [](const auto &cell) {
- return cell.vertices.size() != GeometryInfo<dim>::vertices_per_cell;
- });
+#ifndef DEAL_II_WITH_SIMPLEX_SUPPORT
+ AssertThrow(
+ std::any_of(cells.begin(),
+ cells.end(),
+ [](const auto &cell) {
+ return cell.vertices.size() !=
+ GeometryInfo<dim>::vertices_per_cell;
+ }) == false,
+ ExcMessage(
+ "A cell with invalid number of vertices has been provided."));
+#endif
- if (arbitray_mesh_provided == false)
- {
- internal::TriangulationImplementation::Implementation::
- create_triangulation(v, cells, subcelldata, *this);
- }
- else
- {
- AssertThrow(
- false,
- ExcMessage(
- "A cell with invalid number of vertices has been provided."));
- }
+ internal::TriangulationImplementation::Implementation::
+ create_triangulation(v, cells, subcelldata, *this);
}
catch (...)
{
projects / dealii.git / commitdiff