--- /dev/null
+New: Added a new function GridGenerator::replicate_triangulation() for creating
+a new triangulation by copying a given triangulation repeatedly along the
+coordinate axes.
+<br>
+(David Wells, Victor Zheng, 2020/01/13)
const double duplicated_vertex_tolerance = 1.0e-12,
const bool copy_manifold_ids = false);
+ /**
+ * \brief Replicate a given triangulation in multiple coordinate axes
+ *
+ * @param input The triangulation which will be replicated along the
+ * coordinate axes.
+ *
+ * @param extents A vector with <tt>dim</tt> entries specifying how many
+ * copies of a triangulation should be present along each coordinate axis.
+ *
+ * @param result The triangulation to be created. It needs to be empty upon
+ * calling this function.
+ *
+ * This function creates a new Triangulation equal to a
+ * <tt>dim</tt>-dimensional array of copies of @p input. Copies of @p input
+ * are created by translating @p input along the coordinate axes. Boundary
+ * ids of faces (but not lines in 3D) and all manifold ids are copied but
+ * Manifold objects are not since most Manifold objects do not work
+ * correctly when a Triangulation has been translated.
+ *
+ * To see how this works, consider the following code:
+ * @code
+ * Triangulation<2> input;
+ * GridGenerator::hyper_cube_with_cylindrical_hole(input);
+ * Triangulation<2> output;
+ * GridGenerator::replicate_triangulation(input, {3, 2}, output);
+ * @endcode
+ * results in
+ *
+ * @image html replicated_tria_2d.png
+ *
+ * And, similarly, in 3D:
+ * @code
+ * Triangulation<3> input;
+ * GridGenerator::hyper_cross(1, 1, 1, 2, 1, 2);
+ * Triangulation<3> output;
+ * GridGenerator::replicate_triangulation(input, {3, 2, 1}, output);
+ * @endcode
+ * results in
+ *
+ * @image html replicated_tria_3d.png
+ *
+ * @note This function determines the spacing of the copies of @p input
+ * based on the BoundingBox of @p input. If the boundary faces of @p input
+ * are not aligned with the coordinate axes then the copies might not share
+ * common faces; i.e., this function is intended for simple geometries with
+ * boundary faces aligned along the coordinate axes.
+ */
+ template <int dim, int spacedim = dim>
+ void
+ replicate_triangulation(const Triangulation<dim, spacedim> &input,
+ const std::vector<unsigned int> & extents,
+ Triangulation<dim, spacedim> & result);
+
/**
* Given the two triangulations specified as the first two arguments, create
* the triangulation that contains the finest cells of both triangulation
+ namespace
+ {
+ /**
+ * Merging or replicating triangulations usually results in duplicated
+ * boundary objects - in particular, faces that used to be boundary faces
+ * will now be internal faces.
+ *
+ * This function modifies @p subcell_data by detecting duplicated objects,
+ * marking said duplicates as internal faces, and then deleting all but
+ * one of the duplicates.
+ *
+ * This function relies on some implementation details of
+ * create_triangulation to uniquify objects in SubCellData - in
+ * particular, quadrilaterals are only identified by their lines and not
+ * by their orientation or volume, so rotating and flipping a
+ * quadrilateral doesn't effect the way said quadrilateral is read by that
+ * function.
+ *
+ * @warning even though this function is implemented for structdim 1 and
+ * structdim 2, it will produce <em>wrong</em> results when called for
+ * boundary lines in 3D in most cases since a boundary line can be shared
+ * by an arbitrary number of cells in 3D.
+ */
+ template <int structdim>
+ void
+ delete_duplicated_objects(std::vector<CellData<structdim>> &subcell_data)
+ {
+ static_assert(structdim == 1 || structdim == 2,
+ "This function is only implemented for lines and "
+ "quadrilaterals.");
+ // start by making sure that all objects representing the same vertices
+ // are numbered in the same way by canonicalizing the numberings. This
+ // makes it possible to detect duplicates.
+ for (CellData<structdim> &cell_data : subcell_data)
+ {
+ if (structdim == 1)
+ std::sort(std::begin(cell_data.vertices),
+ std::end(cell_data.vertices));
+ else if (structdim == 2)
+ {
+ // rotate the vertex numbers so that the lowest one is first
+ std::array<unsigned int, 4> renumbering;
+ std::copy(std::begin(cell_data.vertices),
+ std::end(cell_data.vertices),
+ renumbering.begin());
+
+ // convert to old style vertex numbering. This makes the
+ // permutations easy since the valid configurations are
+ //
+ // 3 2 2 1 1 0 0 3
+ // 0 1 3 0 2 3 1 2
+ // (0123) (3012) (2310) (1230)
+ //
+ // rather than the lexical ordering which is harder to permute
+ // by rotation.
+ std::swap(renumbering[2], renumbering[3]);
+ std::rotate(renumbering.begin(),
+ std::min_element(renumbering.begin(),
+ renumbering.end()),
+ renumbering.end());
+ // convert to new style
+ std::swap(renumbering[2], renumbering[3]);
+ // deal with cases where we might have
+ //
+ // 3 2 1 2
+ // 0 1 0 3
+ //
+ // by forcing the second vertex (in lexical ordering) to be
+ // smaller than the third
+ if (renumbering[1] > renumbering[2])
+ std::swap(renumbering[1], renumbering[2]);
+ std::copy(renumbering.begin(),
+ renumbering.end(),
+ std::begin(cell_data.vertices));
+ }
+ }
+
+ // Now that all cell objects have been canonicalized they can be sorted:
+ auto compare = [](const CellData<structdim> &a,
+ const CellData<structdim> &b) {
+ return std::lexicographical_compare(std::begin(a.vertices),
+ std::end(a.vertices),
+ std::begin(b.vertices),
+ std::end(b.vertices));
+ };
+ std::sort(subcell_data.begin(), subcell_data.end(), compare);
+
+ // Finally, determine which objects are duplicates. Duplicates are
+ // assumed to be interior objects, so delete all but one and change the
+ // boundary id:
+ auto left = subcell_data.begin();
+ while (left != subcell_data.end())
+ {
+ const auto right =
+ std::upper_bound(left, subcell_data.end(), *left, compare);
+ // if the range has more than one item, then there are duplicates -
+ // set all boundary ids in the range to the internal boundary id
+ if (left + 1 != right)
+ for (auto it = left; it != right; ++it)
+ {
+ it->boundary_id = numbers::internal_face_boundary_id;
+ Assert(it->manifold_id == left->manifold_id,
+ ExcMessage(
+ "In the process of grid generation a single "
+ "line or quadrilateral has been assigned two "
+ "different manifold ids. This can happen when "
+ "a Triangulation is copied, e.g., via "
+ "GridGenerator::replicate_triangulation() and "
+ "not all external boundary faces have the same "
+ "manifold id. Double check that all faces "
+ "which you expect to be merged together have "
+ "the same manifold id."));
+ }
+ left = right;
+ }
+
+ subcell_data.erase(std::unique(subcell_data.begin(), subcell_data.end()),
+ subcell_data.end());
+ }
+ } // namespace
+
+
+
+ template <int dim, int spacedim>
+ void
+ replicate_triangulation(const Triangulation<dim, spacedim> &input,
+ const std::vector<unsigned int> & extents,
+ Triangulation<dim, spacedim> & result)
+ {
+ AssertDimension(dim, extents.size());
+#ifdef DEBUG
+ for (const auto &extent : extents)
+ Assert(0 < extent,
+ ExcMessage("The Triangulation must be copied at least one time in "
+ "each coordinate dimension."));
+#endif
+ const BoundingBox<spacedim> bbox(input.get_vertices());
+ const auto & min = bbox.get_boundary_points().first;
+ const auto & max = bbox.get_boundary_points().second;
+
+ std::array<Tensor<1, spacedim>, dim> offsets;
+ for (unsigned int d = 0; d < dim; ++d)
+ offsets[d][d] = max[d] - min[d];
+
+ Triangulation<dim, spacedim> tria_to_replicate;
+ tria_to_replicate.copy_triangulation(input);
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ std::vector<Point<spacedim>> input_vertices;
+ std::vector<CellData<dim>> input_cell_data;
+ SubCellData input_subcell_data;
+ std::tie(input_vertices, input_cell_data, input_subcell_data) =
+ GridTools::get_coarse_mesh_description(tria_to_replicate);
+ std::vector<Point<spacedim>> output_vertices = input_vertices;
+ std::vector<CellData<dim>> output_cell_data = input_cell_data;
+ SubCellData output_subcell_data = input_subcell_data;
+
+ for (unsigned int k = 1; k < extents[d]; ++k)
+ {
+ const std::size_t vertex_offset = k * input_vertices.size();
+ // vertices
+ for (const Point<spacedim> &point : input_vertices)
+ output_vertices.push_back(point + double(k) * offsets[d]);
+ // cell data
+ for (const CellData<dim> &cell_data : input_cell_data)
+ {
+ output_cell_data.push_back(cell_data);
+ for (unsigned int &vertex : output_cell_data.back().vertices)
+ vertex += vertex_offset;
+ }
+ // subcell data
+ for (const CellData<1> &boundary_line :
+ input_subcell_data.boundary_lines)
+ {
+ output_subcell_data.boundary_lines.push_back(boundary_line);
+ for (unsigned int &vertex :
+ output_subcell_data.boundary_lines.back().vertices)
+ vertex += vertex_offset;
+ }
+ for (const CellData<2> &boundary_quad :
+ input_subcell_data.boundary_quads)
+ {
+ output_subcell_data.boundary_quads.push_back(boundary_quad);
+ for (unsigned int &vertex :
+ output_subcell_data.boundary_quads.back().vertices)
+ vertex += vertex_offset;
+ }
+ }
+ // check all vertices: since the grid is coarse, most will be on the
+ // boundary anyway
+ std::vector<unsigned int> boundary_vertices;
+ GridTools::delete_duplicated_vertices(
+ output_vertices,
+ output_cell_data,
+ output_subcell_data,
+ boundary_vertices,
+ 1e-6 * input.begin_active()->diameter());
+ // delete_duplicated_vertices also deletes any unused vertices
+ // deal with any reordering issues created by delete_duplicated_vertices
+ GridReordering<dim>::reorder_cells(output_cell_data, true);
+ // clean up the boundary ids of the boundary objects: note that we
+ // have to do this after delete_duplicated_vertices so that boundary
+ // objects are actually duplicated at this point
+ if (dim == 2)
+ delete_duplicated_objects(output_subcell_data.boundary_lines);
+ else if (dim == 3)
+ {
+ delete_duplicated_objects(output_subcell_data.boundary_quads);
+ for (CellData<1> &boundary_line :
+ output_subcell_data.boundary_lines)
+ // set boundary lines to the default value - let
+ // create_triangulation figure out the rest.
+ boundary_line.boundary_id = numbers::internal_face_boundary_id;
+ }
+
+ tria_to_replicate.clear();
+ tria_to_replicate.create_triangulation(output_vertices,
+ output_cell_data,
+ output_subcell_data);
+ }
+
+ result.copy_triangulation(tria_to_replicate);
+ }
+
+
+
template <int dim, int spacedim>
void
create_union_triangulation(
const double duplicated_vertex_tolerance,
const bool copy_manifold_ids);
+ template void
+ replicate_triangulation(
+ const Triangulation<deal_II_dimension, deal_II_space_dimension> &input,
+ const std::vector<unsigned int> & extents,
+ Triangulation<deal_II_dimension, deal_II_space_dimension> &result);
+
template void
create_union_triangulation(
const Triangulation<deal_II_dimension, deal_II_space_dimension>
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 - 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// test GridGenerator::replicate_triangulation
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include "../tests.h"
+
+void
+test()
+{
+ {
+ Triangulation<1, 1> tr1;
+ GridGenerator::hyper_cube(tr1);
+ std::vector<unsigned int> reps1 = {2};
+ Triangulation<1> res1;
+ GridGenerator::replicate_triangulation(tr1, reps1, res1);
+ deallog << "1d triangulation has the following cell centers:" << std::endl;
+ for (const auto &cell : res1.active_cell_iterators())
+ deallog << '(' << cell->center() << ')' << std::endl;
+
+ deallog << "Number of global vertices: " << res1.n_vertices() << std::endl;
+ }
+
+ {
+ Triangulation<2, 2> tr2;
+ GridGenerator::hyper_cube_with_cylindrical_hole(tr2, 0.3, 0.5);
+ const std::vector<unsigned int> reps2 = {{2, 2}};
+ Triangulation<2, 2> res2;
+ GridGenerator::replicate_triangulation(tr2, reps2, res2);
+ deallog << "2d triangulation has the following cell centers:" << std::endl;
+ for (const auto &cell : res2.active_cell_iterators())
+ deallog << '(' << cell->center() << ')' << std::endl;
+
+ deallog << "Number of global vertices: " << res2.n_vertices() << std::endl;
+ }
+
+ {
+ Triangulation<3, 3> tr3;
+ std::vector<unsigned int> reps3 = {2, 2, 2};
+ GridGenerator::subdivided_hyper_rectangle(tr3,
+ reps3,
+ Point<3>(-0.7, -0.5, -1.5),
+ Point<3>(0.5, 0.5, 2.));
+ Triangulation<3, 3> res3;
+ GridGenerator::replicate_triangulation(tr3, reps3, res3);
+ deallog << "3d triangulation has the following cell centers:" << std::endl;
+ for (const auto &cell : res3.active_cell_iterators())
+ deallog << '(' << cell->center() << ')' << std::endl;
+
+ deallog << "Number of global vertices: " << res3.n_vertices() << std::endl;
+ }
+}
+
+int
+main()
+{
+ initlog();
+ test();
+}
--- /dev/null
+
+DEAL::1d triangulation has the following cell centers:
+DEAL::(0.500000)
+DEAL::(1.50000)
+DEAL::Number of global vertices: 3
+DEAL::2d triangulation has the following cell centers:
+DEAL::(0.378033 0.178033)
+DEAL::(0.178033 0.378033)
+DEAL::(-0.178033 0.378033)
+DEAL::(-0.378033 0.178033)
+DEAL::(-0.378033 -0.178033)
+DEAL::(-0.178033 -0.378033)
+DEAL::(0.178033 -0.378033)
+DEAL::(0.378033 -0.178033)
+DEAL::(1.37803 0.178033)
+DEAL::(1.17803 0.378033)
+DEAL::(0.821967 0.378033)
+DEAL::(0.621967 0.178033)
+DEAL::(0.621967 -0.178033)
+DEAL::(0.821967 -0.378033)
+DEAL::(1.17803 -0.378033)
+DEAL::(1.37803 -0.178033)
+DEAL::(0.378033 1.17803)
+DEAL::(0.178033 1.37803)
+DEAL::(-0.178033 1.37803)
+DEAL::(-0.378033 1.17803)
+DEAL::(-0.378033 0.821967)
+DEAL::(-0.178033 0.621967)
+DEAL::(0.178033 0.621967)
+DEAL::(0.378033 0.821967)
+DEAL::(1.37803 1.17803)
+DEAL::(1.17803 1.37803)
+DEAL::(0.821967 1.37803)
+DEAL::(0.621967 1.17803)
+DEAL::(0.621967 0.821967)
+DEAL::(0.821967 0.621967)
+DEAL::(1.17803 0.621967)
+DEAL::(1.37803 0.821967)
+DEAL::Number of global vertices: 53
+DEAL::3d triangulation has the following cell centers:
+DEAL::(-0.400000 -0.250000 -0.625000)
+DEAL::(0.200000 -0.250000 -0.625000)
+DEAL::(-0.400000 0.250000 -0.625000)
+DEAL::(0.200000 0.250000 -0.625000)
+DEAL::(-0.400000 -0.250000 1.12500)
+DEAL::(0.200000 -0.250000 1.12500)
+DEAL::(-0.400000 0.250000 1.12500)
+DEAL::(0.200000 0.250000 1.12500)
+DEAL::(0.800000 -0.250000 -0.625000)
+DEAL::(1.40000 -0.250000 -0.625000)
+DEAL::(0.800000 0.250000 -0.625000)
+DEAL::(1.40000 0.250000 -0.625000)
+DEAL::(0.800000 -0.250000 1.12500)
+DEAL::(1.40000 -0.250000 1.12500)
+DEAL::(0.800000 0.250000 1.12500)
+DEAL::(1.40000 0.250000 1.12500)
+DEAL::(-0.400000 0.750000 -0.625000)
+DEAL::(0.200000 0.750000 -0.625000)
+DEAL::(-0.400000 1.25000 -0.625000)
+DEAL::(0.200000 1.25000 -0.625000)
+DEAL::(-0.400000 0.750000 1.12500)
+DEAL::(0.200000 0.750000 1.12500)
+DEAL::(-0.400000 1.25000 1.12500)
+DEAL::(0.200000 1.25000 1.12500)
+DEAL::(0.800000 0.750000 -0.625000)
+DEAL::(1.40000 0.750000 -0.625000)
+DEAL::(0.800000 1.25000 -0.625000)
+DEAL::(1.40000 1.25000 -0.625000)
+DEAL::(0.800000 0.750000 1.12500)
+DEAL::(1.40000 0.750000 1.12500)
+DEAL::(0.800000 1.25000 1.12500)
+DEAL::(1.40000 1.25000 1.12500)
+DEAL::(-0.400000 -0.250000 2.87500)
+DEAL::(0.200000 -0.250000 2.87500)
+DEAL::(-0.400000 0.250000 2.87500)
+DEAL::(0.200000 0.250000 2.87500)
+DEAL::(-0.400000 -0.250000 4.62500)
+DEAL::(0.200000 -0.250000 4.62500)
+DEAL::(-0.400000 0.250000 4.62500)
+DEAL::(0.200000 0.250000 4.62500)
+DEAL::(0.800000 -0.250000 2.87500)
+DEAL::(1.40000 -0.250000 2.87500)
+DEAL::(0.800000 0.250000 2.87500)
+DEAL::(1.40000 0.250000 2.87500)
+DEAL::(0.800000 -0.250000 4.62500)
+DEAL::(1.40000 -0.250000 4.62500)
+DEAL::(0.800000 0.250000 4.62500)
+DEAL::(1.40000 0.250000 4.62500)
+DEAL::(-0.400000 0.750000 2.87500)
+DEAL::(0.200000 0.750000 2.87500)
+DEAL::(-0.400000 1.25000 2.87500)
+DEAL::(0.200000 1.25000 2.87500)
+DEAL::(-0.400000 0.750000 4.62500)
+DEAL::(0.200000 0.750000 4.62500)
+DEAL::(-0.400000 1.25000 4.62500)
+DEAL::(0.200000 1.25000 4.62500)
+DEAL::(0.800000 0.750000 2.87500)
+DEAL::(1.40000 0.750000 2.87500)
+DEAL::(0.800000 1.25000 2.87500)
+DEAL::(1.40000 1.25000 2.87500)
+DEAL::(0.800000 0.750000 4.62500)
+DEAL::(1.40000 0.750000 4.62500)
+DEAL::(0.800000 1.25000 4.62500)
+DEAL::(1.40000 1.25000 4.62500)
+DEAL::Number of global vertices: 125
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// test GridGenerator::replicate_triangulation's ability to copy boundary ids
+// correctly
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include "../tests.h"
+
+void
+test()
+{
+ // no subcell data in 1D so nothing to test
+
+ {
+ Triangulation<2, 2> tr2;
+ GridGenerator::hyper_cube(tr2, 0.0, 1.0, true);
+ const std::vector<unsigned int> reps2 = {{2, 2}};
+ Triangulation<2, 2> res2;
+ GridGenerator::replicate_triangulation(tr2, reps2, res2);
+ for (const auto &face : res2.active_face_iterators())
+ if (face->at_boundary())
+ deallog << "face: " << face->center()
+ << " boundary id: " << face->boundary_id() << std::endl;
+ }
+
+ {
+ Triangulation<3, 3> tr3;
+ GridGenerator::hyper_cube(tr3, 0.0, 1.0, true);
+ // actually set some boundary line boundary ids to something besides zero
+ // for testing purposes
+ for (auto &face : tr3.active_face_iterators())
+ {
+ if (face->boundary_id() == 4)
+ face->set_all_boundary_ids(4);
+ else if (face->boundary_id() == 5)
+ face->set_all_boundary_ids(5);
+ }
+ // set all (i.e., overwrite two set above) lines on the x = 1 face to
+ // boundary id 1
+ for (auto &face : tr3.active_face_iterators())
+ if (face->boundary_id() == 1)
+ face->set_all_boundary_ids(1);
+ std::vector<unsigned int> reps3 = {2, 1, 1};
+ Triangulation<3, 3> res3;
+ GridGenerator::replicate_triangulation(tr3, reps3, res3);
+ for (const auto &face : res3.active_face_iterators())
+ if (face->at_boundary())
+ {
+ deallog << "face: " << face->center()
+ << " boundary id: " << face->boundary_id() << std::endl;
+ for (unsigned int l = 0; l < GeometryInfo<2>::lines_per_cell; ++l)
+ deallog << "line: " << face->line(l)->center()
+ << " boundary id: " << face->line(l)->boundary_id()
+ << std::endl;
+ }
+ }
+}
+
+int
+main()
+{
+ initlog();
+ test();
+}
--- /dev/null
+
+DEAL::face: 0.500000 0.00000 boundary id: 2
+DEAL::face: 0.00000 0.500000 boundary id: 0
+DEAL::face: 1.50000 0.00000 boundary id: 2
+DEAL::face: 0.00000 1.50000 boundary id: 0
+DEAL::face: 2.00000 0.500000 boundary id: 1
+DEAL::face: 2.00000 1.50000 boundary id: 1
+DEAL::face: 0.500000 2.00000 boundary id: 3
+DEAL::face: 1.50000 2.00000 boundary id: 3
+DEAL::face: 0.500000 0.00000 0.500000 boundary id: 2
+DEAL::line: 0.500000 0.00000 0.00000 boundary id: 0
+DEAL::line: 0.500000 0.00000 1.00000 boundary id: 0
+DEAL::line: 0.00000 0.00000 0.500000 boundary id: 0
+DEAL::line: 1.00000 0.00000 0.500000 boundary id: 0
+DEAL::face: 0.500000 0.500000 0.00000 boundary id: 4
+DEAL::line: 0.00000 0.500000 0.00000 boundary id: 0
+DEAL::line: 1.00000 0.500000 0.00000 boundary id: 0
+DEAL::line: 0.500000 0.00000 0.00000 boundary id: 0
+DEAL::line: 0.500000 1.00000 0.00000 boundary id: 0
+DEAL::face: 0.00000 0.500000 0.500000 boundary id: 0
+DEAL::line: 0.00000 0.00000 0.500000 boundary id: 0
+DEAL::line: 0.00000 1.00000 0.500000 boundary id: 0
+DEAL::line: 0.00000 0.500000 0.00000 boundary id: 0
+DEAL::line: 0.00000 0.500000 1.00000 boundary id: 0
+DEAL::face: 1.50000 0.500000 0.00000 boundary id: 4
+DEAL::line: 1.00000 0.500000 0.00000 boundary id: 0
+DEAL::line: 2.00000 0.500000 0.00000 boundary id: 0
+DEAL::line: 1.50000 0.00000 0.00000 boundary id: 0
+DEAL::line: 1.50000 1.00000 0.00000 boundary id: 0
+DEAL::face: 1.50000 0.00000 0.500000 boundary id: 2
+DEAL::line: 1.50000 0.00000 0.00000 boundary id: 0
+DEAL::line: 1.50000 0.00000 1.00000 boundary id: 0
+DEAL::line: 1.00000 0.00000 0.500000 boundary id: 0
+DEAL::line: 2.00000 0.00000 0.500000 boundary id: 0
+DEAL::face: 0.500000 1.00000 0.500000 boundary id: 3
+DEAL::line: 0.500000 1.00000 0.00000 boundary id: 0
+DEAL::line: 0.500000 1.00000 1.00000 boundary id: 0
+DEAL::line: 0.00000 1.00000 0.500000 boundary id: 0
+DEAL::line: 1.00000 1.00000 0.500000 boundary id: 0
+DEAL::face: 1.50000 1.00000 0.500000 boundary id: 3
+DEAL::line: 1.50000 1.00000 0.00000 boundary id: 0
+DEAL::line: 1.50000 1.00000 1.00000 boundary id: 0
+DEAL::line: 1.00000 1.00000 0.500000 boundary id: 0
+DEAL::line: 2.00000 1.00000 0.500000 boundary id: 0
+DEAL::face: 0.500000 0.500000 1.00000 boundary id: 5
+DEAL::line: 0.00000 0.500000 1.00000 boundary id: 0
+DEAL::line: 1.00000 0.500000 1.00000 boundary id: 0
+DEAL::line: 0.500000 0.00000 1.00000 boundary id: 0
+DEAL::line: 0.500000 1.00000 1.00000 boundary id: 0
+DEAL::face: 1.50000 0.500000 1.00000 boundary id: 5
+DEAL::line: 1.00000 0.500000 1.00000 boundary id: 0
+DEAL::line: 2.00000 0.500000 1.00000 boundary id: 0
+DEAL::line: 1.50000 0.00000 1.00000 boundary id: 0
+DEAL::line: 1.50000 1.00000 1.00000 boundary id: 0
+DEAL::face: 2.00000 0.500000 0.500000 boundary id: 1
+DEAL::line: 2.00000 0.00000 0.500000 boundary id: 0
+DEAL::line: 2.00000 1.00000 0.500000 boundary id: 0
+DEAL::line: 2.00000 0.500000 0.00000 boundary id: 0
+DEAL::line: 2.00000 0.500000 1.00000 boundary id: 0