--- /dev/null
+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().
+<br>
+(Kyle Schwiebert, 2025/04/04)
#include <deal.II/base/mutex.h>
#include <deal.II/base/polynomials_barycentric.h>
+#include <deal.II/base/types.h>
#include <deal.II/fe/fe_poly.h>
const RefinementCase<dim> &refinement_case =
RefinementCase<dim>::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()
*
* The number of vertices in coordinate
* direction @p i is given by <tt>repetitions[i]+1</tt>.
*
- * @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`.
*
const std::vector<unsigned int> &repetitions,
const Point<dim> &p1,
const Point<dim> &p2,
- const bool colorize = false);
+ const bool colorize = false,
+ const bool periodic = false);
/**
* Initialize the given triangulation with a hypercube (square in 2d and
* 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".
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);
/** @} */
#include <deal.II/base/config.h>
+#include <deal.II/base/exceptions.h>
#include <deal.II/base/polynomials_barycentric.h>
#include <deal.II/base/qprojector.h>
+#include <deal.II/base/types.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/mapping.h>
+#include "deal.II/grid/tria_orientation.h"
#include <deal.II/grid/grid_generator.h>
DEAL_II_NAMESPACE_OPEN
+template <int dim, int spacedim>
+unsigned int
+FE_SimplexPoly<dim, spacedim>::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 <int dim, int spacedim>
const FullMatrix<double> &
FE_SimplexPoly<dim, spacedim>::get_restriction_matrix(
#include <deal.II/grid/filtered_iterator.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_tools_topology.h>
#include <deal.II/grid/intergrid_map.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/tria.h>
- template <int dim, int spacedim>
- void
- subdivided_hyper_rectangle_with_simplices(
- Triangulation<dim, spacedim> &tria,
- const std::vector<unsigned int> &repetitions,
- const Point<dim> &p1,
- const Point<dim> &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 <int dim, int spacedim>
+ void
+ subdivided_hyper_rectangle_with_simplices_no_periodic(
+ Triangulation<dim, spacedim> &tria,
+ const std::vector<unsigned int> &repetitions,
+ const Point<dim> &p1,
+ const Point<dim> &p2,
+ const bool colorize)
+ {
+ AssertDimension(dim, spacedim);
- std::vector<Point<spacedim>> vertices;
- std::vector<CellData<dim>> cells;
+ std::vector<Point<spacedim>> vertices;
+ std::vector<CellData<dim>> cells;
- if (dim == 2)
- {
- // determine cell sizes
- const Point<dim> dx((p2[0] - p1[0]) / repetitions[0],
- (p2[1] - p1[1]) / repetitions[1]);
+ if (dim == 2)
+ {
+ // determine cell sizes
+ const Point<dim> 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<spacedim>(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<unsigned int, 4> 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<dim> tri;
+ tri.vertices = {quad[0], quad[1], quad[2]};
+ cells.push_back(tri);
+ }
+
+ // TRI cell 1
+ {
+ CellData<dim> tri;
+ tri.vertices = {quad[3], quad[2], quad[1]};
+ cells.push_back(tri);
+ }
+ }
+ }
+ else if (dim == 3)
+ {
+ // determine cell sizes
+ const Point<dim> 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<spacedim>(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<unsigned int, 8> 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<dim> 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<dim> 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<dim> 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<dim> 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<dim> 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<double>::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 <int dim, int spacedim>
+ void
+ subdivided_hyper_rectangle_with_simplices_periodic(
+ Triangulation<dim, spacedim> &tria,
+ const std::vector<unsigned int> &repetitions,
+ const Point<dim> &p1,
+ const Point<dim> &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<Point<spacedim>> vertices;
+ std::vector<CellData<dim>> cells;
+
+ // determine cell sizes
+ const Point<dim> 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<spacedim>(p1[0] + dx[0] * i, p1[1] + dx[1] * j));
+ vertices.push_back(Point<spacedim>(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<unsigned int, 4> 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<unsigned int, 8> 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<dim> tri;
- tri.vertices = {quad[0], quad[1], quad[2]};
- cells.push_back(tri);
+ CellData<dim> cell;
+ cell.vertices = {{quad[0], quad[1], quad[3], quad[7]}};
+ cells.push_back(cell);
}
- // TRI cell 1
+ // TET cell 1
{
- CellData<dim> tri;
- tri.vertices = {quad[3], quad[2], quad[1]};
- cells.push_back(tri);
+ CellData<dim> cell;
+ cell.vertices = {{quad[0], quad[1], quad[7], quad[5]}};
+ cells.push_back(cell);
}
- }
- }
- else if (dim == 3)
- {
- // determine cell sizes
- const Point<dim> 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<spacedim>(p1[0] + dx[0] * i,
- p1[1] + dx[1] * j,
- p1[2] + dx[2] * k));
+ // TET cell 2
+ {
+ CellData<dim> 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<unsigned int, 8> 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<dim> 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<dim> cell;
+ cell.vertices = {{quad[2], quad[6], quad[0], quad[7]}};
+ cells.push_back(cell);
+ }
- cells.push_back(cell);
- }
+ // TET cell 4
+ {
+ CellData<dim> cell;
+ cell.vertices = {{quad[4], quad[7], quad[5], quad[0]}};
+ cells.push_back(cell);
+ }
- // TET cell 1
- {
- CellData<dim> 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<dim> cell;
+ cell.vertices = {{quad[4], quad[6], quad[7], quad[0]}};
+ cells.push_back(cell);
+ }
+ }
- // TET cell 2
- {
- CellData<dim> 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<dim> 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<double>::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<dim> 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 <int dim, int spacedim>
+ void
+ subdivided_hyper_rectangle_with_simplices(
+ Triangulation<dim, spacedim> &tria,
+ const std::vector<unsigned int> &repetitions,
+ const Point<dim> &p1,
+ const Point<dim> &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<double>::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;
}
}
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)
{
{{repetitions, repetitions, repetitions}},
{p1, p1, p1},
{p2, p2, p2},
- colorize);
+ colorize,
+ periodic);
}
else
{
const std::vector<unsigned int> &repetitions,
const Point<deal_II_dimension> &p1,
const Point<deal_II_dimension> &p2,
- const bool colorize);
+ const bool colorize,
+ const bool periodic);
template void GridGenerator::subdivided_hyper_cube_with_simplices(
Triangulation<deal_II_dimension> & tria,
const unsigned int repetitions,
const double p1,
const double p2,
- const bool colorize);
+ const bool colorize,
+ const bool periodic);
}
bool>(concentric_hyper_shells, arguments, tria);
else if (name == "subdivided_hyper_cube_with_simplices")
- parse_and_create<dim, dim, unsigned int, double, double, bool>(
+ parse_and_create<dim, dim, unsigned int, double, double, bool, bool>(
subdivided_hyper_cube_with_simplices, arguments, tria);
else if (name == "subdivided_hyper_rectangle_with_simplices")
const std::vector<unsigned int> &,
const Point<dim> &,
const Point<dim> &,
+ bool,
bool>(subdivided_hyper_rectangle_with_simplices,
arguments,
tria);
#include <deal.II/base/geometry_info.h>
#include <deal.II/base/point.h>
#include <deal.II/base/tensor.h>
+#include <deal.II/base/types.h>
#include <deal.II/distributed/fully_distributed_tria.h>
#include <deal.II/distributed/shared_tria.h>
{
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<unsigned int, GeometryInfo<dim>::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<unsigned int> face1_vertices(face1->n_vertices(),
+ numbers::invalid_unsigned_int),
+ face2_vertices(face2->n_vertices(), numbers::invalid_unsigned_int);
+
std::set<unsigned int> face2_vertices_set;
for (unsigned int i = 0; i < face1->n_vertices(); ++i)
face2_vertices_set.insert(i);
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");
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
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
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
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
--- /dev/null
+/* ------------------------------------------------------------------------
+ *
+ * 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 <deal.II/base/types.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_orientation.h>
+
+#include <deal.II/lac/affine_constraints.h>
+
+#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<unsigned int> 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<std::vector<GridTools::PeriodicFacePair<
+ typename Triangulation<3>::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<unsigned int>(match.orientation)
+ << " vs default "
+ << static_cast<unsigned int>(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<int>(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<bool, 3> 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<AffineConstraints<double>, 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<bool, 3> 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();
+}
--- /dev/null
+
+# 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!