From 74545664bfc7ff483728d87f082e7e1c53c5253a Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Mon, 29 Jun 2020 22:31:52 +0200 Subject: [PATCH] Refactor Triangulation::create_triangulation() --- include/deal.II/base/types.h | 5 + include/deal.II/grid/connectivity.h | 1645 ++++++++++++++++++++++++ source/grid/tria.cc | 1798 +++++---------------------- 3 files changed, 1984 insertions(+), 1464 deletions(-) create mode 100644 include/deal.II/grid/connectivity.h diff --git a/include/deal.II/base/types.h b/include/deal.II/base/types.h index 53310129e9..ac262ad545 100644 --- a/include/deal.II/base/types.h +++ b/include/deal.II/base/types.h @@ -150,6 +150,11 @@ namespace types * @ref GlossMaterialId "Glossary entry on material indicators" */ using material_id = unsigned int; + + /** + * The type used to denote geometric entity types. + */ + using geometric_entity_type = std::uint8_t; } // namespace types /** diff --git a/include/deal.II/grid/connectivity.h b/include/deal.II/grid/connectivity.h new file mode 100644 index 0000000000..8ef923c927 --- /dev/null +++ b/include/deal.II/grid/connectivity.h @@ -0,0 +1,1645 @@ +// --------------------------------------------------------------------- +// +// 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 + +#include + +#include +#include + + +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 + 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 & + 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 table = {}; + + return table; + } + }; + + + + /** + * Implementation for lines. + */ + struct CellTypeLine : public CellTypeBase + { + dealii::ArrayView + vertices_of_entity(const unsigned int d, + const unsigned int e) const override + { + (void)e; + + if (d == 1) + { + static const std::array 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 table = {2, 1}; + return table[d]; + } + }; + + + + /** + * Implementation for triangles. + */ + struct CellTypeTri : public CellTypeBase + { + dealii::ArrayView + vertices_of_entity(const unsigned int d, + const unsigned int e) const override + { + if (d == 2) + { + static const std::array table = {0, 1, 2}; + + AssertDimension(e, 0); + + return {table}; + } + + if (d == 1) + { + static const std::array, 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 table = {3, 3, 1}; + return table[d]; + } + }; + + + + /** + * Implementation for quadrilaterals. + */ + struct CellTypeQuad : public CellTypeBase + { + dealii::ArrayView + vertices_of_entity(const unsigned int d, + const unsigned int e) const override + { + if (d == 2) + { + static const std::array table = {0, 1, 2, 3}; + + AssertDimension(e, 0); + + return {table}; + } + + if (d == 1) + { + static const std::array, 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 table = {4, 4, 1}; + return table[d]; + } + }; + + + + /** + * Implementation for tetrahedrons. + */ + struct CellTypeTet : public CellTypeBase + { + dealii::ArrayView + vertices_of_entity(const unsigned int d, + const unsigned int e) const override + { + if (d == 3) + { + static const std::array table = {0, 1, 2, 3}; + + AssertDimension(e, 0); + + return {table}; + } + + if (d == 2) + { + static const std::array, 4> table = { + {{0, 1, 2}, {1, 0, 3}, {0, 2, 3}, {2, 1, 3}}}; + + return {table[e]}; + } + + if (d == 1) + { + static const std::array, 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 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, 4> table = { + {{0, 1, 2}, {0, 3, 4}, {2, 5, 3}, {1, 4, 5}}}; + + return table[face][line]; + } + + const std::array & + vertices_of_nth_line_of_surface(const unsigned int line, + const unsigned int face) const override + { + const static std::array, 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 + vertices_of_entity(const unsigned int d, + const unsigned int e) const override + { + if (d == 3) + { + static const std::array table = {0, 1, 2, 3, 4}; + + AssertDimension(e, 0); + + return {table}; + } + + if (d == 2) + { + if (e == 0) + { + static const std::array table = {0, 1, 2, 3}; + return {table}; + } + + static const std::array, 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, 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 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, 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 & + vertices_of_nth_line_of_surface(const unsigned int line, + const unsigned int face) const override + { + static const unsigned int X = static_cast(-1); + + const static std::array, 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 + vertices_of_entity(const unsigned int d, + const unsigned int e) const override + { + if (d == 3) + { + static const std::array table = {0, 1, 2, 3, 4, 5}; + + AssertDimension(e, 0); + + return {table}; + } + + if (d == 2) + { + if (e == 0 || e == 1) + { + static const std::array, 2> table = + {{{1, 0, 2}, {3, 4, 5}}}; + + return {table[e]}; + } + + static const std::array, 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, 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 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(-1); + + const static std::array, 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 & + vertices_of_nth_line_of_surface(const unsigned int line, + const unsigned int face) const override + { + static const unsigned int X = static_cast(-1); + + const static std::array, 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 + vertices_of_entity(const unsigned int d, + const unsigned int e) const override + { + if (d == 3) + { + static const std::array table = { + 0, 1, 2, 3, 4, 5, 6, 7}; + + AssertDimension(e, 0); + + return {table}; + } + + if (d == 2) + { + static const std::array, 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, 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 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, 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 & + vertices_of_nth_line_of_surface(const unsigned int line, + const unsigned int face) const override + { + const static std::array, 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 + struct CRS + { + /** + * Default constructor. + */ + CRS() + : ptr{0} {}; + + /** + * Constructor which allows to set the internal fields directly. + */ + CRS(const std::vector &ptr, const std::vector &col) + : ptr(ptr) + , col(col) + {} + + // row index + std::vector ptr; + + // column index + std::vector 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 + struct Connectivity + { + Connectivity(const unsigned int dim, + const std::vector &cell_types) + : dim(dim) + , cell_types(cell_types) + {} + + inline std::vector & + entity_orientations(const unsigned int structdim) + { + if (structdim == 1) + return line_orientation; + + AssertDimension(structdim, 2); + + return quad_orientation; + } + + inline const std::vector & + entity_orientations(const unsigned int structdim) const + { + if (structdim == 1) + return line_orientation; + + AssertDimension(structdim, 2); + + return quad_orientation; + } + + inline std::vector & + 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 & + 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 & + 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 & + 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 cell_types; + + CRS line_vertices; + + std::vector line_orientation; + + CRS quad_vertices; + CRS quad_lines; + + std::vector quad_orientation; + + CRS cell_entities; + CRS neighbors; + + std::vector 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 + void + determine_neighbors(const CRS &con_cf, CRS &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(col_cf.size(), -1); + ptr_cc = ptr_cf; + + std::vector> 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(-1)) + { + // face is visited the first time -> save the visiting cell + // and the face pointer + neighbors[col_cf[j_0]] = std::pair(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 + inline unsigned char + compute_orientation(const ReferenceCell::Type entity_type, + const std::array & vertices_0, + const std::array & vertices_1) + { + if (entity_type == ReferenceCell::Type::Line) + { + const std::array i{vertices_0[0], vertices_0[1]}; + const std::array j{vertices_1[0], vertices_1[1]}; + + // line_orientation=true + if (i == std::array{{j[0], j[1]}}) + return 1; + + // line_orientation=false + if (i == std::array{{j[1], j[0]}}) + return 0; + } + else if (entity_type == ReferenceCell::Type::Tri) + { + const std::array i{vertices_0[0], vertices_0[1], vertices_0[2]}; + const std::array j{vertices_1[0], vertices_1[1], vertices_1[2]}; + + // face_orientation=true, face_rotation=false, face_flip=false + if (i == std::array{{j[0], j[1], j[2]}}) + return 1; + + // face_orientation=true, face_rotation=true, face_flip=false + if (i == std::array{{j[1], j[0], j[2]}}) + return 3; + + // face_orientation=true, face_rotation=false, face_flip=true + if (i == std::array{{j[2], j[0], j[1]}}) + return 5; + + // face_orientation=false, face_rotation=false, face_flip=false + if (i == std::array{{j[0], j[2], j[1]}}) + return 0; + + // face_orientation=false, face_rotation=true, face_flip=false + if (i == std::array{{j[1], j[2], j[0]}}) + return 2; + + // face_orientation=false, face_rotation=false, face_flip=true + if (i == std::array{{j[2], j[1], j[0]}}) + return 4; + } + else if (entity_type == ReferenceCell::Type::Quad) + { + const std::array i{vertices_0[0], + vertices_0[1], + vertices_0[2], + vertices_0[3]}; + const std::array 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{{j[0], j[1], j[2], j[3]}}) + return 1; + + // face_orientation=true, face_rotation=true, face_flip=false + if (i == std::array{{j[1], j[3], j[0], j[2]}}) + return 3; + + // face_orientation=true, face_rotation=false, face_flip=true + if (i == std::array{{j[3], j[2], j[1], j[0]}}) + return 5; + + // face_orientation=true, face_rotation=true, face_flip=true + if (i == std::array{{j[2], j[0], j[3], j[1]}}) + return 7; + + // face_orientation=false, face_rotation=false, face_flip=false + if (i == std::array{{j[0], j[2], j[1], j[3]}}) + return 0; + + // face_orientation=false, face_rotation=true, face_flip=false + if (i == std::array{{j[2], j[3], j[0], j[1]}}) + return 2; + + // face_orientation=false, face_rotation=false, face_flip=true + if (i == std::array{{j[3], j[1], j[2], j[0]}}) + return 4; + + // face_orientation=false, face_rotation=true, face_flip=true + if (i == std::array{{j[1], j[0], j[3], j[2]}}) + return 6; + } + + AssertThrow(false, dealii::StandardExceptions::ExcNotImplemented()); + + return -1; + } + + + + /** + * Build entities of dimension d (with 0 + void + build_entity_templated( + const unsigned int d, + const std::vector> &cell_types, + const std::vector & cell_types_index, + const CRS & crs, + CRS & crs_d, // result + CRS & crs_0, // result + std::vector & orientations, // result + const FU & second_key_function) + { + const bool compatibility_mode = true; + + const std::vector & cell_ptr = crs.ptr; + const std::vector &cell_vertices = crs.col; + std::vector & ptr_d = crs_d.ptr; + std::vector & 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 & ptr_0 = crs_0.ptr; + std::vector &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(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, unsigned int>> + keys; // key (sorted vertices), cell-entity index + + std::vector> ad_entity_vertices; + std::vector ad_entity_types; + std::vector> 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( + cell_types_index[c])]; + ptr_d[c + 1] = ptr_d[c] + cell_type->n_entities(d); + + // ... collect vertices of cell + const dealii::ArrayView 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 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 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 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( + 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 ref_key; + std::array 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 + void + build_entity(const unsigned int d, + const std::vector> &cell_types, + const std::vector &cell_types_index, + const CRS & crs, + CRS & crs_d, + CRS & crs_0, + std::vector & 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(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> &cell_types, + const std::vector & cell_types_index, + const CRS & con_cv, + const CRS & con_cl, + const CRS & con_lv, + const CRS & con_cq, + const CRS & con_qv, + const std::vector & ori_cq, + CRS & con_ql, // result + std::vector & ori_ql, // result + std::vector & 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( + 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( + 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 + Connectivity + build_connectivity(const unsigned int dim, + const std::vector> &cell_t, + const std::vector &cell_t_id, + const CRS & con_cv) + { + Connectivity connectivity(dim, cell_t_id); + + CRS temp1; // needed for 3D + + if (dim == 1) + connectivity.entity_to_entities(1, 0) = con_cv; + + if (dim == 2 || dim == 3) // build lines + { + std::vector 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 + Connectivity + build_connectivity(const std::vector> &cells) + { + // vector of possible cell entity types + std::vector> cell_types_impl(8); + + cell_types_impl[static_cast( + ReferenceCell::Type::Line)] + .reset(new CellTypeLine()); + cell_types_impl[static_cast( + ReferenceCell::Type::Tri)] + .reset(new CellTypeTri()); + cell_types_impl[static_cast( + ReferenceCell::Type::Quad)] + .reset(new CellTypeQuad()); + cell_types_impl[static_cast( + ReferenceCell::Type::Tet)] + .reset(new CellTypeTet()); + cell_types_impl[static_cast( + ReferenceCell::Type::Pyramid)] + .reset(new CellTypePyramid()); + cell_types_impl[static_cast( + ReferenceCell::Type::Wedge)] + .reset(new CellTypeWedge()); + cell_types_impl[static_cast( + 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, 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 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 cell_vertices_ptr; + cell_vertices_ptr.reserve(cells.size() + 1); + cell_vertices_ptr.push_back(0); + + std::vector 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(cell_type), + cell_types_impl.size()); + Assert(cell_types_impl[static_cast( + 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(dim, + cell_types_impl, + cell_types_indices, + {cell_vertices_ptr, cell_vertices}); + } + } // namespace TriangulationImplementation +} // namespace internal + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/grid/tria.cc b/source/grid/tria.cc index 21f9946bde..eb0b44e86c 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -19,6 +19,7 @@ #include +#include #include #include #include @@ -2173,1539 +2174,411 @@ namespace internal /** - * 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 + template static void - create_triangulation(const std::vector> &v, - const std::vector> & cells, - const SubCellData & /*subcelldata*/, - Triangulation<1, spacedim> &triangulation) + create_triangulation(const std::vector> &vertices, + const std::vector> & cells, + const SubCellData & subcelldata, + Triangulation & 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(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> 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( - 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(cells); + const unsigned int n_cell = cells.size(); - // make up cells - typename Triangulation::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::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::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::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::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::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 - static void - create_triangulation(const std::vector> &v, - const std::vector> & 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(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, - typename Triangulation::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::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 line_vertices( - cells[cell].vertices[GeometryInfo::line_to_cell_vertices( - line, 0)], - cells[cell].vertices[GeometryInfo::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 vertex_touch_count(v.size(), 0); - typename std::map< - std::pair, - typename Triangulation::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( - dim)); - triangulation.faces = - std::make_unique( - 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::raw_line_iterator line = - triangulation.begin_raw_line(); - typename std::map< - std::pair, - typename Triangulation::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(-1)) + level.neighbors[cell * GeometryInfo::faces_per_cell + + j] = {0, nei.col[i]}; + + // set face indices + cells_0.cells[cell * GeometryInfo::faces_per_cell + j] = + crs.col[i]; + + // set face orientation if needed + if (orientation_needed) + level.face_orientations + [cell * GeometryInfo::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::cell_iterator>> - adjacent_cells; + // count number of cells a face is belonging to + std::vector count(bids_face.size(), 0); - // finally make up cells - { - typename Triangulation::raw_cell_iterator cell = - triangulation.begin_raw_quad(); - for (unsigned int c = 0; c < cells.size(); ++c, ++cell) - { - typename Triangulation::line_iterator - lines[GeometryInfo::lines_per_cell]; - for (unsigned int line = 0; - line < GeometryInfo::lines_per_cell; - ++line) - lines[line] = needed_lines[std::make_pair( - cells[c].vertices[GeometryInfo::line_to_cell_vertices( - line, 0)], - cells[c].vertices[GeometryInfo::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::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::line_iterator line; - std::pair 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(-1); + static const unsigned int t_inner = static_cast(-2); + std::vector 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 &q1, const std::vector &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 + template static void - create_triangulation(const std::vector> &v, - const std::vector> & cells, - const SubCellData & subcelldata, - Triangulation<3, spacedim> & triangulation) + process_subcelldata( + const CRS & crs, + TriaObjects & obj, + const std::vector> &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(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, - typename Triangulation::line_iterator> - needed_lines; - for (unsigned int cell = 0; cell < cells.size(); ++cell) + unsigned int counter = 0; + + std::vector key; + key.reserve(GeometryInfo::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::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 line_vertices( - cells[cell].vertices[GeometryInfo::line_to_cell_vertices( - line, 0)], - cells[cell].vertices[GeometryInfo::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 vertex_touch_count(v.size(), 0); - typename std::map< - std::pair, - typename Triangulation::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( - dim)); - triangulation.faces = - std::make_unique( - dim); - reserve_space(*triangulation.levels[0], cells.size(), dim, spacedim); - reserve_space(triangulation.faces->lines, 0, needed_lines.size()); - - // make up lines - { - typename Triangulation::raw_line_iterator line = - triangulation.begin_raw_line(); - typename std::map< - std::pair, - typename Triangulation::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::pair::quad_iterator, - std::array::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 line_list[GeometryInfo::lines_per_cell], - inverse_line_list[GeometryInfo::lines_per_cell]; - unsigned int face_line_list[GeometryInfo::lines_per_face]; - std::array::lines_per_face> orientation; - - for (unsigned int line = 0; - line < GeometryInfo::lines_per_cell; - ++line) - { - line_list[line] = std::pair( - cell.vertices[GeometryInfo::line_to_cell_vertices(line, - 0)], - cell.vertices[GeometryInfo::line_to_cell_vertices(line, - 1)]); - inverse_line_list[line] = std::pair( - cell.vertices[GeometryInfo::line_to_cell_vertices(line, - 1)], - cell.vertices[GeometryInfo::line_to_cell_vertices(line, - 0)]); - } - - for (const unsigned int face : GeometryInfo::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::lines_per_face; - ++l) - if (needed_lines.find( - inverse_line_list[GeometryInfo::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 quad( - {static_cast(face_line_list[0]), - static_cast(face_line_list[1]), - static_cast(face_line_list[2]), - static_cast(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 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::raw_quad_iterator quad = - triangulation.begin_raw_quad(); - typename std::map< - std::vector, - std::pair::quad_iterator, - std::array::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::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::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 line_list[GeometryInfo::lines_per_cell], - inverse_line_list[GeometryInfo::lines_per_cell]; - unsigned int face_line_list[4]; - for (unsigned int line = 0; - line < GeometryInfo::lines_per_cell; - ++line) - { - line_list[line] = std::make_pair( - cells[c].vertices[GeometryInfo::line_to_cell_vertices( - line, 0)], - cells[c].vertices[GeometryInfo::line_to_cell_vertices( - line, 1)]); - inverse_line_list[line] = std::pair( - cells[c].vertices[GeometryInfo::line_to_cell_vertices( - line, 1)], - cells[c].vertices[GeometryInfo::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::quad_iterator - face_iterator[GeometryInfo::faces_per_cell]; - bool face_orientation[GeometryInfo::faces_per_cell]; - bool face_flip[GeometryInfo::faces_per_cell]; - bool face_rotation[GeometryInfo::faces_per_cell]; - for (const unsigned int face : GeometryInfo::face_indices()) - { - for (unsigned int l = 0; - l < GeometryInfo::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 quad( - {static_cast(face_line_list[0]), - static_cast(face_line_list[1]), - static_cast(face_line_list[2]), - static_cast(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 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::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> - cell_to_face_lines; - for (const unsigned int face : GeometryInfo::face_indices()) - for (unsigned int line = 0; - line < GeometryInfo::lines_per_face; - ++line) - cell_to_face_lines.insert( - std::pair>( - GeometryInfo::face_to_cell_lines(face, line), - std::pair(face, line))); - std::multimap>:: - 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::standard_to_real_face_line( - line1, - face_orientation[face1], - face_flip[face1], - face_rotation[face1])) == - face_iterator[face2]->line( - GeometryInfo::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::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::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::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::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::line_iterator line; - std::pair 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::quad_iterator quad; - typename Triangulation::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 line_vertices( - subcell_quad - .vertices[GeometryInfo::line_to_cell_vertices(i, - 0)], - subcell_quad - .vertices[GeometryInfo::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 quad_compare_1({line[0]->index(), - line[1]->index(), - line[2]->index(), - line[3]->index()}); - std::vector 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::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 - 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]); } @@ -10757,23 +9630,20 @@ Triangulation::create_triangulation( // 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::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::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 (...) { -- 2.39.5