From 549e39982f28ea49ca6f9a8a2d06f346dbbdad11 Mon Sep 17 00:00:00 2001 From: Kyle Schwiebert Date: Wed, 9 Apr 2025 20:15:11 -0400 Subject: [PATCH] Add basic periodicity for tetrahedra with a mesh generator and a test. --- doc/news/changes/minor/20250404KyleSchwiebert | 9 + include/deal.II/fe/fe_simplex_p.h | 17 + include/deal.II/grid/grid_generator.h | 28 +- source/fe/fe_simplex_p.cc | 137 +++++ source/grid/grid_generator.cc | 476 ++++++++++++------ source/grid/grid_generator.inst.in | 6 +- source/grid/grid_generator_from_name.cc | 3 +- source/grid/grid_tools_dof_handlers.cc | 17 +- ...rid_generator_from_name_and_argument_02.cc | 10 +- ...generator_from_name_and_argument_02.output | 8 +- tests/simplex/periodic_tet.cc | 300 +++++++++++ tests/simplex/periodic_tet.output | 385 ++++++++++++++ 12 files changed, 1218 insertions(+), 178 deletions(-) create mode 100644 doc/news/changes/minor/20250404KyleSchwiebert create mode 100644 tests/simplex/periodic_tet.cc create mode 100644 tests/simplex/periodic_tet.output diff --git a/doc/news/changes/minor/20250404KyleSchwiebert b/doc/news/changes/minor/20250404KyleSchwiebert new file mode 100644 index 0000000000..f25ad59d72 --- /dev/null +++ b/doc/news/changes/minor/20250404KyleSchwiebert @@ -0,0 +1,9 @@ +New: FE_Simplex_Poly now implements FiniteElement::face_to_cell_index(), +enabling periodicity for derived classes, with restrictions similar to the +implementation in FE_Q_Base. Additionally, +GridGenerator::hyper_rectangle_with_simplicies() now accepts an additional +boolean argument which, in 3D, breaks each cubical cell into six tetrahedra +instead of five. This mesh meets the restrictions of +FE_Simplex_Poly::face_to_cell_index(). +
+(Kyle Schwiebert, 2025/04/04) diff --git a/include/deal.II/fe/fe_simplex_p.h b/include/deal.II/fe/fe_simplex_p.h index f493c16649..7ce1999318 100644 --- a/include/deal.II/fe/fe_simplex_p.h +++ b/include/deal.II/fe/fe_simplex_p.h @@ -19,6 +19,7 @@ #include #include +#include #include @@ -65,6 +66,22 @@ public: const RefinementCase &refinement_case = RefinementCase::isotropic_refinement) const override; + /** + * @see FiniteElement::face_to_cell_index() + * + * @note This function has some limitations in 3D. If @p dim is 3, this + * function works if either: + * 1. @p combined_orientation is the default orientation. + * 2. @p combined_orientation is `(false, false, false)` and the triangular + * face has no more than one degree of freedom in its interior. + * Otherwise this function throws an error. + */ + virtual unsigned int + face_to_cell_index(const unsigned int face_dof_index, + const unsigned int face, + const types::geometric_orientation combined_orientation = + numbers::default_geometric_orientation) const override; + /** * @copydoc dealii::FiniteElement::get_restriction_matrix() * diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h index 73a4415d6e..c081235abc 100644 --- a/include/deal.II/grid/grid_generator.h +++ b/include/deal.II/grid/grid_generator.h @@ -2651,9 +2651,14 @@ namespace GridGenerator * The number of vertices in coordinate * direction @p i is given by repetitions[i]+1. * - * @note This function takes the mesh produced by subdivided_hyper_rectangle() - * and further subdivides each cell into 2 triangles (for @p dim 2) or - * 5 tetrahedra (for @p dim 3), respectively. + * This function takes the mesh produced by subdivided_hyper_rectangle() and + * further subdivides each cell. For @p dim 2, it subdivides each cell into 2 + * triangles. For @p dim 3, it subdivides each cell into 5 or 6 tetrahedra + * based on the value of @p periodic. If @p periodic is true, then we split + * each cell into 6 cells so that each face of the rectangular prism has + * the same stencil, enabling periodicity. If @p periodic is false, we + * instead subdivide each hexahedral cell into 5 tetrahedra. If @p dim is not + * 3, then @p periodic has no effect. * * @note Currently, this function only works for `dim==spacedim`. * @@ -2667,7 +2672,8 @@ namespace GridGenerator const std::vector &repetitions, const Point &p1, const Point &p2, - const bool colorize = false); + const bool colorize = false, + const bool periodic = false); /** * Initialize the given triangulation with a hypercube (square in 2d and @@ -2679,9 +2685,14 @@ namespace GridGenerator * the limits are given as arguments. They default to zero and unity, then * producing the unit hypercube. * - * @note This function takes the mesh produced by subdivided_hyper_cube() - * and further subdivides each cell into 2 triangles (for @p dim 2) or - * 5 tetrahedra (for @p dim 3), respectively. + * This function takes the mesh produced by subdivided_hyper_cube() and + * further subdivides each cell. For @p dim 2, it subdivides each cell into 2 + * triangles. For @p dim 3, it subdivides each cell into 5 or 6 tetrahedra + * based on the value of @p periodic. If @p periodic is true, then we split + * each cell into 6 cells so that each face of the rectangular prism has + * the same stencil, enabling periodicity. If @p periodic is false, we + * instead subdivide each hexahedral cell into 5 tetrahedra. If @p dim is not + * 3, then @p periodic has no effect. * * Also see * @ref simplex "Simplex support". @@ -2692,7 +2703,8 @@ namespace GridGenerator const unsigned int repetitions, const double p1 = 0.0, const double p2 = 1.0, - const bool colorize = false); + const bool colorize = false, + const bool periodic = false); /** @} */ diff --git a/source/fe/fe_simplex_p.cc b/source/fe/fe_simplex_p.cc index e958760b5c..5c2a65e1fc 100644 --- a/source/fe/fe_simplex_p.cc +++ b/source/fe/fe_simplex_p.cc @@ -14,8 +14,10 @@ #include +#include #include #include +#include #include #include @@ -24,6 +26,7 @@ #include #include +#include "deal.II/grid/tria_orientation.h" #include DEAL_II_NAMESPACE_OPEN @@ -436,6 +439,140 @@ FE_SimplexPoly::get_prolongation_matrix( +template +unsigned int +FE_SimplexPoly::face_to_cell_index( + const unsigned int face_dof_index, + const unsigned int face, + const types::geometric_orientation combined_orientation) const +{ + // Function largely lifted from FE_Q_Base::face_to_cell_index() + AssertIndexRange(face_dof_index, this->n_dofs_per_face(face)); + AssertIndexRange(face, this->reference_cell().n_faces()); + + // TODO: once the default orientation is switched to 0 then we can remove this + // special case for 1D. + if (dim == 1) + Assert(combined_orientation == numbers::default_geometric_orientation, + ExcMessage("In 1D, all faces must have the default orientation.")); + else + AssertIndexRange(combined_orientation, + this->reference_cell().n_face_orientations(face)); + + // we need to distinguish between DoFs on vertices, lines and in 3d quads. + // do so in a sequence of if-else statements + if (face_dof_index < this->get_first_face_line_index(face)) + // DoF is on a vertex + { + // get the number of the vertex on the face that corresponds to this DoF, + // along with the number of the DoF on this vertex + const unsigned int face_vertex = + face_dof_index / this->n_dofs_per_vertex(); + const unsigned int dof_index_on_vertex = + face_dof_index % this->n_dofs_per_vertex(); + + // then get the number of this vertex on the cell and translate + // this to a DoF number on the cell + return this->reference_cell().face_to_cell_vertices( + face, face_vertex, combined_orientation) * + this->n_dofs_per_vertex() + + dof_index_on_vertex; + } + else if (face_dof_index < this->get_first_face_quad_index(face)) + // DoF is on a line + { + // do the same kind of translation as before. we need to only consider + // DoFs on the lines, i.e., ignoring those on the vertices + const unsigned int index = + face_dof_index - this->get_first_face_line_index(face); + + const unsigned int face_line = index / this->n_dofs_per_line(); + const unsigned int dof_index_on_line = index % this->n_dofs_per_line(); + + // we now also need to adjust the line index for the case of + // face orientation, face flips, etc + unsigned int adjusted_dof_index_on_line = 0; + switch (dim) + { + case 1: + DEAL_II_ASSERT_UNREACHABLE(); + break; + + case 2: + if (combined_orientation == numbers::default_geometric_orientation) + adjusted_dof_index_on_line = dof_index_on_line; + else + adjusted_dof_index_on_line = + this->n_dofs_per_line() - dof_index_on_line - 1; + break; + + case 3: + // in 3d, things are difficult. someone will have to think + // about how this code here should look like, by drawing a bunch + // of pictures of how all the faces can look like with the various + // flips and rotations. + // + // that said, we can implement a couple other situations easily: + // if the face orientation is + // numbers::default_combined_face_orientation then things are + // simple. Likewise if the face orientation is + // internal::combined_face_orientation(false,false,false) then we + // just know we need to reverse the DoF order on each line. + // + // this at least allows for tetrahedral meshes where periodic face + // pairs are each in one of the two above configurations. Doing so + // may turn out to be always possible or it may not, but at least + // this is less restrictive. + Assert((this->n_dofs_per_line() <= 1) || + combined_orientation == + numbers::default_geometric_orientation || + combined_orientation == + internal::combined_face_orientation(false, false, false), + ExcNotImplemented()); + + if (combined_orientation == numbers::default_geometric_orientation) + { + adjusted_dof_index_on_line = dof_index_on_line; + } + else + { + adjusted_dof_index_on_line = + this->n_dofs_per_line() - dof_index_on_line - 1; + } + break; + + default: + DEAL_II_ASSERT_UNREACHABLE(); + } + + return (this->get_first_line_index() + + this->reference_cell().face_to_cell_lines(face, + face_line, + combined_orientation) * + this->n_dofs_per_line() + + adjusted_dof_index_on_line); + } + else + // DoF is on a quad + { + Assert(dim >= 3, ExcInternalError()); + + // ignore vertex and line dofs + const unsigned int index = + face_dof_index - this->get_first_face_quad_index(face); + + // the same is true here as above for the 3d case -- someone will + // just have to draw a bunch of pictures. in the meantime, + // we can implement the degree <= 3 case in which it is simple + Assert((this->n_dofs_per_quad(face) <= 1) || + combined_orientation == numbers::default_geometric_orientation, + ExcNotImplemented()); + return (this->get_first_quad_index(face) + index); + } +} + + + template const FullMatrix & FE_SimplexPoly::get_restriction_matrix( diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index 5acf678de8..40d87b3614 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -22,6 +22,7 @@ #include #include #include +#include #include #include #include @@ -8940,179 +8941,348 @@ namespace GridGenerator - template - void - subdivided_hyper_rectangle_with_simplices( - Triangulation &tria, - const std::vector &repetitions, - const Point &p1, - const Point &p2, - const bool colorize) + // Hide the implementation for two cases of + // subdivided_hyper_rectangle_with_simplices in an anonymous namespace. + namespace { - AssertDimension(dim, spacedim); + template + void + subdivided_hyper_rectangle_with_simplices_no_periodic( + Triangulation &tria, + const std::vector &repetitions, + const Point &p1, + const Point &p2, + const bool colorize) + { + AssertDimension(dim, spacedim); - std::vector> vertices; - std::vector> cells; + std::vector> vertices; + std::vector> cells; - if (dim == 2) - { - // determine cell sizes - const Point dx((p2[0] - p1[0]) / repetitions[0], - (p2[1] - p1[1]) / repetitions[1]); + if (dim == 2) + { + // determine cell sizes + const Point dx((p2[0] - p1[0]) / repetitions[0], + (p2[1] - p1[1]) / repetitions[1]); + + // create vertices + for (unsigned int j = 0; j <= repetitions[1]; ++j) + for (unsigned int i = 0; i <= repetitions[0]; ++i) + vertices.push_back( + Point(p1[0] + dx[0] * i, p1[1] + dx[1] * j)); + + // create cells + for (unsigned int j = 0; j < repetitions[1]; ++j) + for (unsigned int i = 0; i < repetitions[0]; ++i) + { + // create reference QUAD cell + std::array quad{{ + (j + 0) * (repetitions[0] + 1) + i + 0, // + (j + 0) * (repetitions[0] + 1) + i + 1, // + (j + 1) * (repetitions[0] + 1) + i + 0, // + (j + 1) * (repetitions[0] + 1) + i + 1 // + }}; // + + // TRI cell 0 + { + CellData tri; + tri.vertices = {quad[0], quad[1], quad[2]}; + cells.push_back(tri); + } + + // TRI cell 1 + { + CellData tri; + tri.vertices = {quad[3], quad[2], quad[1]}; + cells.push_back(tri); + } + } + } + else if (dim == 3) + { + // determine cell sizes + const Point dx((p2[0] - p1[0]) / repetitions[0], + (p2[1] - p1[1]) / repetitions[1], + (p2[2] - p1[2]) / repetitions[2]); + + // create vertices + for (unsigned int k = 0; k <= repetitions[2]; ++k) + for (unsigned int j = 0; j <= repetitions[1]; ++j) + for (unsigned int i = 0; i <= repetitions[0]; ++i) + vertices.push_back(Point(p1[0] + dx[0] * i, + p1[1] + dx[1] * j, + p1[2] + dx[2] * k)); + + // create cells + for (unsigned int k = 0; k < repetitions[2]; ++k) + for (unsigned int j = 0; j < repetitions[1]; ++j) + for (unsigned int i = 0; i < repetitions[0]; ++i) + { + // create reference HEX cell + std::array quad{ + {(k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 0, + (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 1, + (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 0, + (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 1, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 0, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 1, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 0, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 1}}; + + // TET cell 0 + { + CellData cell; + if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) + cell.vertices = {{quad[0], quad[1], quad[2], quad[4]}}; + else + cell.vertices = {{quad[0], quad[1], quad[3], quad[5]}}; + + cells.push_back(cell); + } - // create vertices + // TET cell 1 + { + CellData cell; + if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) + cell.vertices = {{quad[2], quad[1], quad[3], quad[7]}}; + else + cell.vertices = {{quad[0], quad[3], quad[2], quad[6]}}; + cells.push_back(cell); + } + + // TET cell 2 + { + CellData cell; + if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) + cell.vertices = {{quad[1], quad[4], quad[5], quad[7]}}; + else + cell.vertices = {{quad[0], quad[4], quad[5], quad[6]}}; + cells.push_back(cell); + } + + // TET cell 3 + { + CellData cell; + if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) + cell.vertices = {{quad[2], quad[4], quad[7], quad[6]}}; + else + cell.vertices = {{quad[3], quad[5], quad[7], quad[6]}}; + cells.push_back(cell); + } + + // TET cell 4 + { + CellData cell; + if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) + cell.vertices = {{quad[1], quad[2], quad[4], quad[7]}}; + else + cell.vertices = {{quad[0], quad[3], quad[6], quad[5]}}; + cells.push_back(cell); + } + } + } + else + { + AssertThrow(false, ExcNotImplemented()); + } + + // actually create triangulation + tria.create_triangulation(vertices, cells, SubCellData()); + + if (colorize) + { + // to colorize, run through all + // faces of all cells and set + // boundary indicator to the + // correct value if it was 0. + + // use a large epsilon to + // compare numbers to avoid + // roundoff problems. + double epsilon = std::numeric_limits::max(); + for (unsigned int i = 0; i < dim; ++i) + epsilon = + std::min(epsilon, + 0.01 * (std::abs(p2[i] - p1[i]) / repetitions[i])); + Assert(epsilon > 0, + ExcMessage( + "The distance between corner points must be positive.")); + + // actual code is external since + // 1-D is different from 2/3d. + colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon); + } + } + + + + // This function is only needed in 3D. + template + void + subdivided_hyper_rectangle_with_simplices_periodic( + Triangulation &tria, + const std::vector &repetitions, + const Point &p1, + const Point &p2, + const bool colorize) + { + // This function is only needed in 3D (and hypothetically in higher + // dimension), so library internals should ensure it is never called + // unless dim == 3. + Assert(dim == 3, ExcInternalError()); + AssertDimension(dim, spacedim); + + std::vector> vertices; + std::vector> cells; + + // determine cell sizes + const Point dx((p2[0] - p1[0]) / repetitions[0], + (p2[1] - p1[1]) / repetitions[1], + (p2[2] - p1[2]) / repetitions[2]); + + // create vertices + for (unsigned int k = 0; k <= repetitions[2]; ++k) for (unsigned int j = 0; j <= repetitions[1]; ++j) for (unsigned int i = 0; i <= repetitions[0]; ++i) - vertices.push_back( - Point(p1[0] + dx[0] * i, p1[1] + dx[1] * j)); + vertices.push_back(Point(p1[0] + dx[0] * i, + p1[1] + dx[1] * j, + p1[2] + dx[2] * k)); - // create cells + // create cells + for (unsigned int k = 0; k < repetitions[2]; ++k) for (unsigned int j = 0; j < repetitions[1]; ++j) for (unsigned int i = 0; i < repetitions[0]; ++i) { - // create reference QUAD cell - std::array quad{{ - (j + 0) * (repetitions[0] + 1) + i + 0, // - (j + 0) * (repetitions[0] + 1) + i + 1, // - (j + 1) * (repetitions[0] + 1) + i + 0, // - (j + 1) * (repetitions[0] + 1) + i + 1 // - }}; // - - // TRI cell 0 + // create reference HEX cell + std::array quad{ + {(k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 0, + (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 1, + (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 0, + (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 1, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 0, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 1, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 0, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 1}}; + + // TET cell 0 { - CellData tri; - tri.vertices = {quad[0], quad[1], quad[2]}; - cells.push_back(tri); + CellData cell; + cell.vertices = {{quad[0], quad[1], quad[3], quad[7]}}; + cells.push_back(cell); } - // TRI cell 1 + // TET cell 1 { - CellData tri; - tri.vertices = {quad[3], quad[2], quad[1]}; - cells.push_back(tri); + CellData cell; + cell.vertices = {{quad[0], quad[1], quad[7], quad[5]}}; + cells.push_back(cell); } - } - } - else if (dim == 3) - { - // determine cell sizes - const Point dx((p2[0] - p1[0]) / repetitions[0], - (p2[1] - p1[1]) / repetitions[1], - (p2[2] - p1[2]) / repetitions[2]); - // create vertices - for (unsigned int k = 0; k <= repetitions[2]; ++k) - for (unsigned int j = 0; j <= repetitions[1]; ++j) - for (unsigned int i = 0; i <= repetitions[0]; ++i) - vertices.push_back(Point(p1[0] + dx[0] * i, - p1[1] + dx[1] * j, - p1[2] + dx[2] * k)); + // TET cell 2 + { + CellData cell; + cell.vertices = {{quad[0], quad[7], quad[3], quad[2]}}; + cells.push_back(cell); + } - // create cells - for (unsigned int k = 0; k < repetitions[2]; ++k) - for (unsigned int j = 0; j < repetitions[1]; ++j) - for (unsigned int i = 0; i < repetitions[0]; ++i) + // TET cell 3 { - // create reference HEX cell - std::array quad{ - {(k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + - (j + 0) * (repetitions[0] + 1) + i + 0, - (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + - (j + 0) * (repetitions[0] + 1) + i + 1, - (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + - (j + 1) * (repetitions[0] + 1) + i + 0, - (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + - (j + 1) * (repetitions[0] + 1) + i + 1, - (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + - (j + 0) * (repetitions[0] + 1) + i + 0, - (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + - (j + 0) * (repetitions[0] + 1) + i + 1, - (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + - (j + 1) * (repetitions[0] + 1) + i + 0, - (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + - (j + 1) * (repetitions[0] + 1) + i + 1}}; - - // TET cell 0 - { - CellData cell; - if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) - cell.vertices = {{quad[0], quad[1], quad[2], quad[4]}}; - else - cell.vertices = {{quad[0], quad[1], quad[3], quad[5]}}; + CellData cell; + cell.vertices = {{quad[2], quad[6], quad[0], quad[7]}}; + cells.push_back(cell); + } - cells.push_back(cell); - } + // TET cell 4 + { + CellData cell; + cell.vertices = {{quad[4], quad[7], quad[5], quad[0]}}; + cells.push_back(cell); + } - // TET cell 1 - { - CellData cell; - if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) - cell.vertices = {{quad[2], quad[1], quad[3], quad[7]}}; - else - cell.vertices = {{quad[0], quad[3], quad[2], quad[6]}}; - cells.push_back(cell); - } + // TET cell 5 + { + CellData cell; + cell.vertices = {{quad[4], quad[6], quad[7], quad[0]}}; + cells.push_back(cell); + } + } - // TET cell 2 - { - CellData cell; - if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) - cell.vertices = {{quad[1], quad[4], quad[5], quad[7]}}; - else - cell.vertices = {{quad[0], quad[4], quad[5], quad[6]}}; - cells.push_back(cell); - } + // actually create triangulation + tria.create_triangulation(vertices, cells, SubCellData()); - // TET cell 3 - { - CellData cell; - if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) - cell.vertices = {{quad[2], quad[4], quad[7], quad[6]}}; - else - cell.vertices = {{quad[3], quad[5], quad[7], quad[6]}}; - cells.push_back(cell); - } + if (colorize) + { + // to colorize, run through all + // faces of all cells and set + // boundary indicator to the + // correct value if it was 0. + + // use a large epsilon to + // compare numbers to avoid + // roundoff problems. + double epsilon = std::numeric_limits::max(); + for (unsigned int i = 0; i < dim; ++i) + epsilon = + std::min(epsilon, + 0.01 * (std::abs(p2[i] - p1[i]) / repetitions[i])); + Assert(epsilon > 0, + ExcMessage( + "The distance between corner points must be positive.")); - // TET cell 4 - { - CellData cell; - if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) - cell.vertices = {{quad[1], quad[2], quad[4], quad[7]}}; - else - cell.vertices = {{quad[0], quad[3], quad[6], quad[5]}}; - cells.push_back(cell); - } - } + // actual code is external since + // 1-D is different from 2/3d. + colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon); + } + } + } // namespace + + + + template + void + subdivided_hyper_rectangle_with_simplices( + Triangulation &tria, + const std::vector &repetitions, + const Point &p1, + const Point &p2, + const bool colorize, + const bool periodic) + { + // We only need to call the "periodic" variant if it was requested and we + // are in 3D. + if (dim != 3) + { + subdivided_hyper_rectangle_with_simplices_no_periodic( + tria, repetitions, p1, p2, colorize); + return; } - else + else if (!periodic) { - AssertThrow(false, ExcNotImplemented()); + subdivided_hyper_rectangle_with_simplices_no_periodic( + tria, repetitions, p1, p2, colorize); + return; } - - // actually create triangulation - tria.create_triangulation(vertices, cells, SubCellData()); - - if (colorize) + else { - // to colorize, run through all - // faces of all cells and set - // boundary indicator to the - // correct value if it was 0. - - // use a large epsilon to - // compare numbers to avoid - // roundoff problems. - double epsilon = std::numeric_limits::max(); - for (unsigned int i = 0; i < dim; ++i) - epsilon = std::min(epsilon, - 0.01 * (std::abs(p2[i] - p1[i]) / repetitions[i])); - Assert(epsilon > 0, - ExcMessage( - "The distance between corner points must be positive.")); - - // actual code is external since - // 1-D is different from 2/3d. - colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon); + subdivided_hyper_rectangle_with_simplices_periodic( + tria, repetitions, p1, p2, colorize); + return; } } @@ -9124,12 +9294,17 @@ namespace GridGenerator const unsigned int repetitions, const double p1, const double p2, - const bool colorize) + const bool colorize, + const bool periodic) { if (dim == 2) { - subdivided_hyper_rectangle_with_simplices( - tria, {{repetitions, repetitions}}, {p1, p1}, {p2, p2}, colorize); + subdivided_hyper_rectangle_with_simplices(tria, + {{repetitions, repetitions}}, + {p1, p1}, + {p2, p2}, + colorize, + periodic); } else if (dim == 3) { @@ -9138,7 +9313,8 @@ namespace GridGenerator {{repetitions, repetitions, repetitions}}, {p1, p1, p1}, {p2, p2, p2}, - colorize); + colorize, + periodic); } else { diff --git a/source/grid/grid_generator.inst.in b/source/grid/grid_generator.inst.in index 189ce57a84..90d15de511 100644 --- a/source/grid/grid_generator.inst.in +++ b/source/grid/grid_generator.inst.in @@ -302,12 +302,14 @@ for (deal_II_dimension : DIMENSIONS) const std::vector &repetitions, const Point &p1, const Point &p2, - const bool colorize); + const bool colorize, + const bool periodic); template void GridGenerator::subdivided_hyper_cube_with_simplices( Triangulation & tria, const unsigned int repetitions, const double p1, const double p2, - const bool colorize); + const bool colorize, + const bool periodic); } diff --git a/source/grid/grid_generator_from_name.cc b/source/grid/grid_generator_from_name.cc index 361915d5f2..bc0b89a38f 100644 --- a/source/grid/grid_generator_from_name.cc +++ b/source/grid/grid_generator_from_name.cc @@ -232,7 +232,7 @@ namespace GridGenerator bool>(concentric_hyper_shells, arguments, tria); else if (name == "subdivided_hyper_cube_with_simplices") - parse_and_create( + parse_and_create( subdivided_hyper_cube_with_simplices, arguments, tria); else if (name == "subdivided_hyper_rectangle_with_simplices") @@ -241,6 +241,7 @@ namespace GridGenerator const std::vector &, const Point &, const Point &, + bool, bool>(subdivided_hyper_rectangle_with_simplices, arguments, tria); diff --git a/source/grid/grid_tools_dof_handlers.cc b/source/grid/grid_tools_dof_handlers.cc index 009e8d3274..3f7da7e914 100644 --- a/source/grid/grid_tools_dof_handlers.cc +++ b/source/grid/grid_tools_dof_handlers.cc @@ -15,6 +15,7 @@ #include #include #include +#include #include #include @@ -2432,19 +2433,17 @@ namespace GridTools { Assert(matrix.m() == matrix.n(), ExcMessage("The supplied matrix must be a square matrix")); - - static const int dim = FaceIterator::AccessorType::dimension; + Assert(face1->reference_cell() == face2->reference_cell(), + ExcMessage( + "The faces to be matched must have the same reference cell.")); // Do a full matching of the face vertices: - - std::array::vertices_per_face> - face1_vertices, face2_vertices; - - face1_vertices.fill(numbers::invalid_unsigned_int); - face2_vertices.fill(numbers::invalid_unsigned_int); - AssertDimension(face1->n_vertices(), face2->n_vertices()); + std::vector face1_vertices(face1->n_vertices(), + numbers::invalid_unsigned_int), + face2_vertices(face2->n_vertices(), numbers::invalid_unsigned_int); + std::set face2_vertices_set; for (unsigned int i = 0; i < face1->n_vertices(); ++i) face2_vertices_set.insert(i); diff --git a/tests/grid/grid_generator_from_name_and_argument_02.cc b/tests/grid/grid_generator_from_name_and_argument_02.cc index 2883b09d6e..40d9b19297 100644 --- a/tests/grid/grid_generator_from_name_and_argument_02.cc +++ b/tests/grid/grid_generator_from_name_and_argument_02.cc @@ -54,13 +54,15 @@ main() test<2, 2>("hyper_ball_balanced", "0,0 : 1"); test<3, 3>("hyper_ball_balanced", "0,0,0 : 1"); - test<2, 2>("subdivided_hyper_cube_with_simplices", "2 : 0.0 : 1.0 : false"); - test<3, 3>("subdivided_hyper_cube_with_simplices", "2 : 0.0 : 1.0 : false"); + test<2, 2>("subdivided_hyper_cube_with_simplices", + "2 : 0.0 : 1.0 : false : false"); + test<3, 3>("subdivided_hyper_cube_with_simplices", + "2 : 0.0 : 1.0 : false : false"); test<2, 2>("subdivided_hyper_rectangle_with_simplices", - "2, 2 : 0.0, 0.0 : 1.0, 2.0 : false"); + "2, 2 : 0.0, 0.0 : 1.0, 2.0 : false : false"); test<3, 3>("subdivided_hyper_rectangle_with_simplices", - "2, 2, 3 : 0.0, 0.0, 1.0 : 1.0, 2.0, 3.0 : false"); + "2, 2, 3 : 0.0, 0.0, 1.0 : 1.0, 2.0, 3.0 : false : false"); test<2, 2>("subdivided_hyper_L", "5, 5 : 0, 0 : 1, 1 : 2, 3"); test<3, 3>("subdivided_hyper_L", "5, 5, 5 : 0, 0, 0 : 1, 1, 1 : 2, 2, 3"); diff --git a/tests/grid/grid_generator_from_name_and_argument_02.output b/tests/grid/grid_generator_from_name_and_argument_02.output index e5a39f8d22..31b233c3f5 100644 --- a/tests/grid/grid_generator_from_name_and_argument_02.output +++ b/tests/grid/grid_generator_from_name_and_argument_02.output @@ -191,7 +191,7 @@ $ELM 31 5 0 0 8 18 16 38 37 27 25 52 53 32 5 0 0 8 43 36 48 50 37 38 52 53 $ENDELM -DEAL::Generating Triangulation<2, 2> : subdivided_hyper_cube_with_simplices(2 : 0.0 : 1.0 : false) +DEAL::Generating Triangulation<2, 2> : subdivided_hyper_cube_with_simplices(2 : 0.0 : 1.0 : false : false) $NOD 9 1 0.00000 0.00000 0 @@ -215,7 +215,7 @@ $ELM 7 2 0 0 3 5 6 8 8 2 0 0 3 9 8 6 $ENDELM -DEAL::Generating Triangulation<3, 3> : subdivided_hyper_cube_with_simplices(2 : 0.0 : 1.0 : false) +DEAL::Generating Triangulation<3, 3> : subdivided_hyper_cube_with_simplices(2 : 0.0 : 1.0 : false : false) $NOD 27 1 0.00000 0.00000 0.00000 @@ -289,7 +289,7 @@ $ELM 39 4 0 0 4 18 24 27 26 40 4 0 0 4 14 18 26 24 $ENDELM -DEAL::Generating Triangulation<2, 2> : subdivided_hyper_rectangle_with_simplices(2, 2 : 0.0, 0.0 : 1.0, 2.0 : false) +DEAL::Generating Triangulation<2, 2> : subdivided_hyper_rectangle_with_simplices(2, 2 : 0.0, 0.0 : 1.0, 2.0 : false : false) $NOD 9 1 0.00000 0.00000 0 @@ -313,7 +313,7 @@ $ELM 7 2 0 0 3 5 6 8 8 2 0 0 3 9 8 6 $ENDELM -DEAL::Generating Triangulation<3, 3> : subdivided_hyper_rectangle_with_simplices(2, 2, 3 : 0.0, 0.0, 1.0 : 1.0, 2.0, 3.0 : false) +DEAL::Generating Triangulation<3, 3> : subdivided_hyper_rectangle_with_simplices(2, 2, 3 : 0.0, 0.0, 1.0 : 1.0, 2.0, 3.0 : false : false) $NOD 36 1 0.00000 0.00000 1.00000 diff --git a/tests/simplex/periodic_tet.cc b/tests/simplex/periodic_tet.cc new file mode 100644 index 0000000000..d3745afddb --- /dev/null +++ b/tests/simplex/periodic_tet.cc @@ -0,0 +1,300 @@ +/* ------------------------------------------------------------------------ + * + * SPDX-License-Identifier: LGPL-2.1-or-later + * Copyright (C) 2020 - 2024 by the deal.II authors + * + * This file is part of the deal.II library. + * + * Part of the source code is dual licensed under Apache-2.0 WITH + * LLVM-exception OR LGPL-2.1-or-later. Detailed license information + * governing the source code and code contributions can be found in + * LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. + * + * ------------------------------------------------------------------------ + */ + +// Test that the function +// GridGenerator::subdivided_hyper_rectable_with_simplices() can give a mesh +// which is suitible for periodicity with tetrahedra. + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "../tests.h" + +void +check() +{ + Triangulation<3> tria; + Point<3> corner_p1 = {-0.5, 1.0, 0.0}; + Point<3> corner_p2 = {0.5, 2.0, 1.2}; + std::vector sub_divisions = {2, 2, 2}; + GridGenerator::subdivided_hyper_rectangle_with_simplices( + tria, sub_divisions, corner_p1, corner_p2, true, true); + + // Check that we have made a grid with the expected number of cells and only + // one level of refinement. + if (tria.n_cells() != 48) + { + deallog << "FAIL! Wrong number of cells." << std::endl; + } + + if (tria.n_levels() != 1) + { + deallog << "FAIL! Mesh has more than one refinement level." << std::endl; + } + + // Check that each cell has a positive measure. + { + for (const auto &cell : tria.active_cell_iterators()) + { + double measure = cell->measure(); + if (measure <= 0.0) + { + deallog << "FAIL! Found a cell with nonpositive measure: " + << cell->active_cell_index() << std::endl; + } + } + } + + // Collect periodic face pairs and ensure that they make sense and we have + // the right numbers of them. + std::array::cell_iterator>>, + 3> + pairs; + GridTools::collect_periodic_faces(tria, 0, 1, 0, pairs[0]); + GridTools::collect_periodic_faces(tria, 2, 3, 1, pairs[1]); + GridTools::collect_periodic_faces(tria, 4, 5, 2, pairs[2]); + tria.add_periodicity(pairs[0]); + tria.add_periodicity(pairs[1]); + tria.add_periodicity(pairs[2]); + + // To enable periodicity the easy way, we need each boundary pair to be in + // default orientation. Check that this holds. + { + IteratorFilters::AtBoundary filter; + const types::geometric_orientation default_orientation = + numbers::default_geometric_orientation; + const types::geometric_orientation default_mirror_orientation = + internal::combined_face_orientation(false, false, false); + for (const auto &pair : pairs) + { + for (const auto &match : pair) + { + if (match.orientation != default_orientation && + match.orientation != default_mirror_orientation) + { + deallog << "FAIL! Found a face pair which is not in default " + << "orientation." << std::endl; + deallog << "Orientation is " + << static_cast(match.orientation) + << " vs default " + << static_cast(default_orientation) + << std::endl; + deallog << "Cells " << match.cell[0]->active_cell_index() + << " and " << match.cell[1]->active_cell_index() + << std::endl; + deallog << "Faces " << match.face_idx[0] << " and " + << match.face_idx[1] << std::endl; + return; + } + } + } + } + + { + // Check that we have gotten all the periodicity. + bool mismatched_faces = false; + unsigned int changing_direction = 0; + for (const auto &pair : pairs) + { + if (pair.size() != 8) + { + deallog << "FAIL! There should be 8 matches on each face." + << std::endl; + return; + } + for (int face = 0; face < static_cast(pair.size()); ++face) + { + Point<3> p1( + pair[face].cell[0]->face(pair[face].face_idx[0])->center()); + Point<3> p2( + pair[face].cell[1]->face(pair[face].face_idx[1])->center()); + const std::array results = { + {std::abs(p1[0] - p2[0]) > 1e-10, + std::abs(p1[1] - p2[1]) > 1e-10, + std::abs(p1[2] - p2[2]) > 1e-10}}; + bool fail = false; + + for (unsigned int i = 0; i < 3; ++i) + { + if (i == changing_direction) + { + continue; + } + else + { + fail = fail || results[i]; + } + } + + if (fail) + { + deallog << "FAIL! Found non matching faces." << std::endl; + return; + } + } + ++changing_direction; + } + } + + // write 2 outputs (total mesh and only surface mesh) + const auto grid_out = [](const auto &tria, const bool surface_mesh_only) { + GridOutFlags::Vtk flags; + + if (surface_mesh_only) + { + flags.output_cells = false; + flags.output_faces = true; + flags.output_edges = false; + flags.output_only_relevant = false; + } + + GridOut grid_out; + grid_out.set_flags(flags); + + grid_out.write_vtk(tria, deallog.get_file_stream()); + }; + + grid_out(tria, false); // total mesh + grid_out(tria, true); // only surface mesh + + // If we got this far, we have a mesh which supports periodicity. + // Let's see if we can proscribe periodicity with actual DoFs now. + // Check for C0 tetrahedra up to degree 3. (Highest implemented at the time of + // writing.) + for (const unsigned int degree : {1, 2, 3}) + { + DoFHandler<3> dof_handler(tria); + FE_SimplexP<3> fe(degree); + dof_handler.distribute_dofs(fe); + std::array, 3> constraints; + + // Periodic in the x direction: + DoFTools::make_periodicity_constraints( + dof_handler, 0, 1, 0, constraints[0]); + constraints[0].close(); + + // Periodic in the y direction: + DoFTools::make_periodicity_constraints( + dof_handler, 2, 3, 1, constraints[1]); + constraints[1].close(); + + // Periodic in the z direction: + DoFTools::make_periodicity_constraints( + dof_handler, 4, 5, 2, constraints[2]); + constraints[2].close(); + + // Some information for checking constraints. + MappingFE<3> mapping(FE_SimplexP<3>(1)); + auto map = DoFTools::map_dofs_to_support_points(mapping, dof_handler); + unsigned int changing_direction = 0; + + for (const auto &cons : constraints) + { + // Check the total number of constraints. + if (cons.n_constraints() != (2 * degree + 1) * (2 * degree + 1)) + { + deallog << "FAIL! Wrong number of constraints!" << std::endl; + } + + // Loop through the constraints and verify that the correct components + // of the DoF's support points are equal. + for (const auto &constraint : cons.get_lines()) + { + if (constraint.entries.size() == 1) + { + if (std::abs(constraint.entries[0].second - 1.0) > 1e-10) + { + deallog << "FAIL! All constraint entries should be one." + << std::endl; + continue; + } + + const types::global_dof_index dof_1 = constraint.index; + const types::global_dof_index dof_2 = + constraint.entries[0].first; + const std::array results = { + {std::abs(map[dof_1][0] - map[dof_2][0]) > 1e-10, + std::abs(map[dof_1][1] - map[dof_2][1]) > 1e-10, + std::abs(map[dof_1][2] - map[dof_2][2]) > 1e-10}}; + bool fail = false; + + for (unsigned int i = 0; i < 3; ++i) + { + if (i == changing_direction) + { + continue; + } + else + { + fail = fail || results[i]; + } + } + + if (fail) + { + deallog << "FAIL! Found non matching constraint." + << std::endl; + deallog << "(" << map[dof_1][0] << "," << map[dof_1][1] + << "," << map[dof_1][2] << ")" << std::endl; + deallog << "(" << map[dof_2][0] << "," << map[dof_2][1] + << "," << map[dof_2][2] << ")" << std::endl; + } + } + else + { + deallog + << "FAIL! " + << "Found constraint with the wrong number of entries." + << std::endl; + return; + } + } + ++changing_direction; + + // If degree is 3, put down a little more data for the output. + if (degree == 3) + { + cons.print(deallog.get_file_stream()); + } + } + } + + deallog << "OK!" << std::endl; +} + +int +main() +{ + initlog(); + + check(); +} diff --git a/tests/simplex/periodic_tet.output b/tests/simplex/periodic_tet.output new file mode 100644 index 0000000000..ca3a39b96d --- /dev/null +++ b/tests/simplex/periodic_tet.output @@ -0,0 +1,385 @@ + +# vtk DataFile Version 3.0 +Triangulation generated with deal.II +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 27 double +-0.500000 1.00000 0.00000 +0.00000 1.00000 0.00000 +0.500000 1.00000 0.00000 +-0.500000 1.50000 0.00000 +0.00000 1.50000 0.00000 +0.500000 1.50000 0.00000 +-0.500000 2.00000 0.00000 +0.00000 2.00000 0.00000 +0.500000 2.00000 0.00000 +-0.500000 1.00000 0.600000 +0.00000 1.00000 0.600000 +0.500000 1.00000 0.600000 +-0.500000 1.50000 0.600000 +0.00000 1.50000 0.600000 +0.500000 1.50000 0.600000 +-0.500000 2.00000 0.600000 +0.00000 2.00000 0.600000 +0.500000 2.00000 0.600000 +-0.500000 1.00000 1.20000 +0.00000 1.00000 1.20000 +0.500000 1.00000 1.20000 +-0.500000 1.50000 1.20000 +0.00000 1.50000 1.20000 +0.500000 1.50000 1.20000 +-0.500000 2.00000 1.20000 +0.00000 2.00000 1.20000 +0.500000 2.00000 1.20000 + +CELLS 88 400 +4 0 1 4 13 +4 0 1 13 10 +4 0 13 4 3 +4 3 12 0 13 +4 9 13 10 0 +4 9 12 13 0 +4 1 2 5 14 +4 1 2 14 11 +4 1 14 5 4 +4 4 13 1 14 +4 10 14 11 1 +4 10 13 14 1 +4 3 4 7 16 +4 3 4 16 13 +4 3 16 7 6 +4 6 15 3 16 +4 12 16 13 3 +4 12 15 16 3 +4 4 5 8 17 +4 4 5 17 14 +4 4 17 8 7 +4 7 16 4 17 +4 13 17 14 4 +4 13 16 17 4 +4 9 10 13 22 +4 9 10 22 19 +4 9 22 13 12 +4 12 21 9 22 +4 18 22 19 9 +4 18 21 22 9 +4 10 11 14 23 +4 10 11 23 20 +4 10 23 14 13 +4 13 22 10 23 +4 19 23 20 10 +4 19 22 23 10 +4 12 13 16 25 +4 12 13 25 22 +4 12 25 16 15 +4 15 24 12 25 +4 21 25 22 12 +4 21 24 25 12 +4 13 14 17 26 +4 13 14 26 23 +4 13 26 17 16 +4 16 25 13 26 +4 22 26 23 13 +4 22 25 26 13 +3 1 0 10 +3 0 1 4 +3 2 1 11 +3 1 2 5 +3 5 2 14 +3 14 2 11 +3 0 4 3 +3 3 4 7 +3 1 5 4 +3 4 5 8 +3 8 5 17 +3 17 5 14 +3 15 6 16 +3 3 7 6 +3 7 16 6 +3 16 7 17 +3 4 8 7 +3 8 17 7 +3 9 10 0 +3 10 9 19 +3 10 11 1 +3 11 10 20 +3 23 11 20 +3 14 11 23 +3 26 14 23 +3 24 15 25 +3 16 25 15 +3 25 16 26 +3 17 14 26 +3 17 26 16 +3 18 21 22 +3 18 22 19 +3 18 19 9 +3 19 23 20 +3 19 20 10 +3 21 24 25 +3 21 25 22 +3 19 22 23 +3 22 26 23 +3 22 25 26 + +CELL_TYPES 88 +10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 +5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + + +CELL_DATA 88 +SCALARS MaterialID int 1 +LOOKUP_TABLE default +6 6 6 6 6 6 7 7 7 7 7 7 6 6 6 6 6 6 7 7 7 7 7 7 6 6 6 6 6 6 7 7 7 7 7 7 6 6 6 6 6 6 7 7 7 7 7 7 +2 4 2 4 1 1 4 4 4 4 1 1 3 4 3 3 4 3 2 2 2 2 1 1 1 3 3 3 1 3 5 5 2 5 2 5 5 5 5 5 + + +SCALARS ManifoldID int 1 +LOOKUP_TABLE default +-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 +-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + +# vtk DataFile Version 3.0 +Triangulation generated with deal.II +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 27 double +-0.500000 1.00000 0.00000 +0.00000 1.00000 0.00000 +0.500000 1.00000 0.00000 +-0.500000 1.50000 0.00000 +0.00000 1.50000 0.00000 +0.500000 1.50000 0.00000 +-0.500000 2.00000 0.00000 +0.00000 2.00000 0.00000 +0.500000 2.00000 0.00000 +-0.500000 1.00000 0.600000 +0.00000 1.00000 0.600000 +0.500000 1.00000 0.600000 +-0.500000 1.50000 0.600000 +0.00000 1.50000 0.600000 +0.500000 1.50000 0.600000 +-0.500000 2.00000 0.600000 +0.00000 2.00000 0.600000 +0.500000 2.00000 0.600000 +-0.500000 1.00000 1.20000 +0.00000 1.00000 1.20000 +0.500000 1.00000 1.20000 +-0.500000 1.50000 1.20000 +0.00000 1.50000 1.20000 +0.500000 1.50000 1.20000 +-0.500000 2.00000 1.20000 +0.00000 2.00000 1.20000 +0.500000 2.00000 1.20000 + +CELLS 48 192 +3 1 0 10 +3 0 1 4 +3 2 1 11 +3 1 2 5 +3 5 2 14 +3 14 2 11 +3 3 12 0 +3 0 4 3 +3 3 4 7 +3 1 5 4 +3 4 5 8 +3 8 5 17 +3 17 5 14 +3 6 15 3 +3 15 6 16 +3 3 7 6 +3 7 16 6 +3 16 7 17 +3 4 8 7 +3 8 17 7 +3 12 9 0 +3 9 10 0 +3 10 9 19 +3 10 11 1 +3 11 10 20 +3 23 11 20 +3 15 12 3 +3 12 21 9 +3 14 11 23 +3 26 14 23 +3 15 24 12 +3 24 15 25 +3 16 25 15 +3 25 16 26 +3 17 14 26 +3 17 26 16 +3 21 18 9 +3 18 21 22 +3 18 22 19 +3 18 19 9 +3 19 23 20 +3 19 20 10 +3 24 21 12 +3 21 24 25 +3 21 25 22 +3 19 22 23 +3 22 26 23 +3 22 25 26 + +CELL_TYPES 48 +5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + + +CELL_DATA 48 +SCALARS MaterialID int 1 +LOOKUP_TABLE default +2 4 2 4 1 1 0 4 4 4 4 1 1 0 3 4 3 3 4 3 0 2 2 2 2 1 0 0 1 1 0 3 3 3 1 3 0 5 5 2 5 2 0 5 5 5 5 5 + + +SCALARS ManifoldID int 1 +LOOKUP_TABLE default +-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + 64 0: 1.00000 + 65 30: 1.00000 + 66 40: 1.00000 + 69 31: 1.00000 + 70 32: 1.00000 + 75 44: 1.00000 + 76 43: 1.00000 + 77 41: 1.00000 + 78 42: 1.00000 + 82 47: 1.00000 + 83 50: 1.00000 + 86 56: 1.00000 + 87 55: 1.00000 + 88 61: 1.00000 + 89 60: 1.00000 + 92 63: 1.00000 + 160 132: 1.00000 + 161 142: 1.00000 + 162 133: 1.00000 + 163 134: 1.00000 + 168 146: 1.00000 + 169 145: 1.00000 + 170 143: 1.00000 + 171 144: 1.00000 + 175 149: 1.00000 + 176 157: 1.00000 + 177 156: 1.00000 + 179 159: 1.00000 + 244 220: 1.00000 + 247 224: 1.00000 + 248 223: 1.00000 + 249 221: 1.00000 + 250 222: 1.00000 + 253 227: 1.00000 + 254 230: 1.00000 + 257 236: 1.00000 + 258 235: 1.00000 + 259 241: 1.00000 + 260 240: 1.00000 + 263 243: 1.00000 + 316 298: 1.00000 + 319 302: 1.00000 + 320 301: 1.00000 + 321 299: 1.00000 + 322 300: 1.00000 + 325 305: 1.00000 + 326 313: 1.00000 + 327 312: 1.00000 + 329 315: 1.00000 + 112 1: 1.00000 + 113 20: 1.00000 + 122 23: 1.00000 + 123 24: 1.00000 + 132 0: 1.00000 + 135 22: 1.00000 + 136 21: 1.00000 + 137 5: 1.00000 + 138 4: 1.00000 + 141 27: 1.00000 + 142 50: 1.00000 + 143 56: 1.00000 + 144 55: 1.00000 + 147 54: 1.00000 + 148 53: 1.00000 + 150 59: 1.00000 + 160 64: 1.00000 + 161 83: 1.00000 + 170 86: 1.00000 + 171 87: 1.00000 + 180 85: 1.00000 + 181 84: 1.00000 + 182 68: 1.00000 + 183 67: 1.00000 + 186 90: 1.00000 + 187 107: 1.00000 + 188 106: 1.00000 + 189 110: 1.00000 + 280 206: 1.00000 + 285 209: 1.00000 + 286 210: 1.00000 + 294 208: 1.00000 + 295 207: 1.00000 + 297 213: 1.00000 + 298 230: 1.00000 + 299 236: 1.00000 + 300 235: 1.00000 + 303 234: 1.00000 + 304 233: 1.00000 + 306 239: 1.00000 + 316 254: 1.00000 + 321 257: 1.00000 + 322 258: 1.00000 + 330 256: 1.00000 + 331 255: 1.00000 + 333 261: 1.00000 + 334 275: 1.00000 + 335 274: 1.00000 + 336 278: 1.00000 + 196 2: 1.00000 + 206 1: 1.00000 + 211 7: 1.00000 + 212 6: 1.00000 + 220 30: 1.00000 + 225 36: 1.00000 + 226 35: 1.00000 + 230 0: 1.00000 + 231 9: 1.00000 + 232 8: 1.00000 + 233 5: 1.00000 + 234 4: 1.00000 + 237 16: 1.00000 + 240 31: 1.00000 + 241 32: 1.00000 + 242 38: 1.00000 + 244 65: 1.00000 + 254 64: 1.00000 + 259 70: 1.00000 + 260 69: 1.00000 + 268 96: 1.00000 + 269 95: 1.00000 + 272 72: 1.00000 + 273 71: 1.00000 + 274 68: 1.00000 + 275 67: 1.00000 + 276 79: 1.00000 + 279 98: 1.00000 + 280 112: 1.00000 + 290 115: 1.00000 + 291 114: 1.00000 + 298 132: 1.00000 + 303 138: 1.00000 + 304 137: 1.00000 + 308 117: 1.00000 + 309 116: 1.00000 + 310 124: 1.00000 + 312 133: 1.00000 + 313 134: 1.00000 + 314 140: 1.00000 + 316 160: 1.00000 + 326 163: 1.00000 + 327 162: 1.00000 + 334 183: 1.00000 + 335 182: 1.00000 + 338 165: 1.00000 + 339 164: 1.00000 + 340 172: 1.00000 + 342 185: 1.00000 +DEAL::OK! -- 2.39.5