From: David Wells <drwells@email.unc.edu>
Date: Sun, 11 Nov 2018 22:53:27 +0000 (-0500)
Subject: Add GridTools::get_coarse_mesh_description.
X-Git-Tag: v9.1.0-rc1~532^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9e6ba48d9597c0b6922555cc99b4598f18adbf40;p=dealii.git

Add GridTools::get_coarse_mesh_description.
---

diff --git a/doc/news/changes/minor/20181111DavidWells b/doc/news/changes/minor/20181111DavidWells
new file mode 100644
index 0000000000..9b95b3b3fa
--- /dev/null
+++ b/doc/news/changes/minor/20181111DavidWells
@@ -0,0 +1,5 @@
+New: A new function GridTools::get_coarse_mesh_description has been added. This
+function extracts the vertices, CellData, and SubCellData objects from a
+Triangulation (i.e., it is the inverse of Triangulation::create_triangulation).
+<br>
+(David Wells, 2018/11/11)
diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h
index 2762529c6e..9cdb1543d1 100644
--- a/include/deal.II/grid/grid_tools.h
+++ b/include/deal.II/grid/grid_tools.h
@@ -263,6 +263,19 @@ namespace GridTools
     const Iterator &                                      object,
     const Point<Iterator::AccessorType::space_dimension> &trial_point);
 
+  /**
+   * Return the arrays that define the coarse mesh of a Triangulation. This
+   * function is the inverse of Triangulation::create_triangulation().
+   *
+   * This function is useful in cases where one needs to deconstruct a
+   * Triangulation or manipulate the numbering of the vertices in some way: an
+   * example is GridGenerator::merge_triangulations().
+   */
+  template <int dim, int spacedim>
+  std::
+    tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
+    get_coarse_mesh_description(const Triangulation<dim, spacedim> &tria);
+
   /*@}*/
   /**
    * @name Functions supporting the creation of meshes
diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index 48f77b7a03..24def8abd2 100644
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -429,6 +429,184 @@ namespace GridTools
 
 
 
+  // Generic functions for appending face data in 2D or 3D. TODO: we can
+  // remove these once we have 'if constexpr'.
+  namespace
+  {
+    void
+    append_face_data(const CellData<0> & /*face_data*/,
+                     SubCellData & /*subcell_data*/)
+    {
+      Assert(false, ExcInternalError());
+    }
+
+
+
+    void
+    append_face_data(const CellData<1> &face_data, SubCellData &subcell_data)
+    {
+      subcell_data.boundary_lines.push_back(face_data);
+    }
+
+
+
+    void
+    append_face_data(const CellData<2> &face_data, SubCellData &subcell_data)
+    {
+      subcell_data.boundary_quads.push_back(face_data);
+    }
+
+
+
+    // Lexical comparison for sorting CellData objects.
+    template <int structdim>
+    struct CellDataComparator
+    {
+      bool
+      operator()(const CellData<structdim> &a,
+                 const CellData<structdim> &b) const
+      {
+        // Check vertices:
+        if (std::lexicographical_compare(std::begin(a.vertices),
+                                         std::end(a.vertices),
+                                         std::begin(b.vertices),
+                                         std::end(b.vertices)))
+          return true;
+          // it should never be necessary to check the material or manifold
+          // ids as a 'tiebreaker' (since they must be equal if the vertex
+          // indices are equal). Assert it anyway:
+#ifdef DEBUG
+        if (std::equal(std::begin(a.vertices),
+                       std::end(a.vertices),
+                       std::begin(b.vertices)))
+          {
+            Assert(a.material_id == b.material_id &&
+                     a.manifold_id == b.manifold_id,
+                   ExcMessage(
+                     "Two CellData objects with equal vertices must "
+                     "have the same material/boundary ids and manifold "
+                     "ids."));
+          }
+#endif
+        return false;
+      }
+    };
+  } // namespace
+
+
+
+  template <int dim, int spacedim>
+  std::
+    tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
+    get_coarse_mesh_description(const Triangulation<dim, spacedim> &tria)
+  {
+    Assert(1 <= tria.n_levels(),
+           ExcMessage("The input triangulation must be non-empty."));
+
+    std::vector<Point<spacedim>> vertices;
+    std::vector<CellData<dim>>   cells;
+    SubCellData                  subcell_data;
+
+    unsigned int max_level_0_vertex_n = 0;
+    for (const auto &cell : tria.cell_iterators_on_level(0))
+      for (unsigned int cell_vertex_n = 0;
+           cell_vertex_n < GeometryInfo<dim>::vertices_per_cell;
+           ++cell_vertex_n)
+        max_level_0_vertex_n =
+          std::max(cell->vertex_index(cell_vertex_n), max_level_0_vertex_n);
+    vertices.resize(max_level_0_vertex_n + 1);
+    std::set<CellData<dim - 1>, CellDataComparator<dim - 1>> face_data;
+    std::set<CellData<1>, CellDataComparator<1>> line_data; // only used in 3D
+    for (const auto &cell : tria.cell_iterators_on_level(0))
+      {
+        // Save cell data
+        CellData<dim> cell_data;
+        for (unsigned int cell_vertex_n = 0;
+             cell_vertex_n < GeometryInfo<dim>::vertices_per_cell;
+             ++cell_vertex_n)
+          {
+            Assert(cell->vertex_index(cell_vertex_n) < vertices.size(),
+                   ExcInternalError());
+            vertices[cell->vertex_index(cell_vertex_n)] =
+              cell->vertex(cell_vertex_n);
+            cell_data.vertices[cell_vertex_n] =
+              cell->vertex_index(cell_vertex_n);
+          }
+        cell_data.material_id = cell->material_id();
+        cell_data.manifold_id = cell->manifold_id();
+        cells.push_back(cell_data);
+
+        // Save face data
+        if (dim != 1)
+          {
+            for (unsigned int face_n = 0;
+                 face_n < GeometryInfo<dim>::faces_per_cell;
+                 ++face_n)
+              {
+                const auto        face = cell->face(face_n);
+                CellData<dim - 1> face_cell_data;
+                for (unsigned int vertex_n = 0;
+                     vertex_n < GeometryInfo<dim>::vertices_per_face;
+                     ++vertex_n)
+                  face_cell_data.vertices[vertex_n] =
+                    face->vertex_index(vertex_n);
+                face_cell_data.boundary_id = face->boundary_id();
+                face_cell_data.manifold_id = face->manifold_id();
+
+                face_data.insert(face_cell_data);
+              }
+          }
+        // Save line data
+        if (dim == 3)
+          {
+            for (unsigned int line_n = 0;
+                 line_n < GeometryInfo<dim>::lines_per_cell;
+                 ++line_n)
+              {
+                const auto  line = cell->line(line_n);
+                CellData<1> line_cell_data;
+                for (unsigned int vertex_n = 0;
+                     vertex_n < GeometryInfo<2>::vertices_per_face;
+                     ++vertex_n)
+                  line_cell_data.vertices[vertex_n] =
+                    line->vertex_index(vertex_n);
+                line_cell_data.boundary_id = line->boundary_id();
+                line_cell_data.manifold_id = line->manifold_id();
+
+                line_data.insert(line_cell_data);
+              }
+          }
+      }
+      // Double-check that there are no unused vertices:
+#ifdef DEBUG
+    {
+      std::vector<bool> used_vertices(vertices.size());
+      for (const CellData<dim> &cell_data : cells)
+        for (unsigned int cell_vertex_n = 0;
+             cell_vertex_n < GeometryInfo<dim>::vertices_per_cell;
+             ++cell_vertex_n)
+          used_vertices[cell_data.vertices[cell_vertex_n]] = true;
+      Assert(std::find(used_vertices.begin(), used_vertices.end(), false) ==
+               used_vertices.end(),
+             ExcMessage("The level zero vertices should form a contiguous "
+                        "range."));
+    }
+#endif
+
+    for (const CellData<dim - 1> &face_cell_data : face_data)
+      append_face_data(face_cell_data, subcell_data);
+    if (dim == 3)
+      for (const CellData<1> &face_line_data : line_data)
+        subcell_data.boundary_lines.push_back(face_line_data);
+    return std::tuple<std::vector<Point<spacedim>>,
+                      std::vector<CellData<dim>>,
+                      SubCellData>(std::move(vertices),
+                                   std::move(cells),
+                                   std::move(subcell_data));
+  }
+
+
+
   template <int dim, int spacedim>
   void
   delete_unused_vertices(std::vector<Point<spacedim>> &vertices,
diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in
index 673d937c91..03bb9c6d38 100644
--- a/source/grid/grid_tools.inst.in
+++ b/source/grid/grid_tools.inst.in
@@ -169,6 +169,12 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
       compute_bounding_box(
         const Triangulation<deal_II_dimension, deal_II_space_dimension> &);
 
+      template std::tuple<std::vector<Point<deal_II_space_dimension>>,
+                          std::vector<CellData<deal_II_dimension>>,
+                          SubCellData>
+      get_coarse_mesh_description(
+        const Triangulation<deal_II_dimension, deal_II_space_dimension> &tria);
+
       template void
       delete_unused_vertices(std::vector<Point<deal_II_space_dimension>> &,
                              std::vector<CellData<deal_II_dimension>> &,
diff --git a/tests/grid/get_coarse_mesh_description_01.cc b/tests/grid/get_coarse_mesh_description_01.cc
new file mode 100644
index 0000000000..23980a44a1
--- /dev/null
+++ b/tests/grid/get_coarse_mesh_description_01.cc
@@ -0,0 +1,154 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Test GridTools::get_coarse_mesh_description()
+
+#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 "../tests.h"
+
+
+
+void setup_tria(Triangulation<1> &tria)
+{
+  GridGenerator::subdivided_hyper_rectangle(
+    tria, std::vector<unsigned int>(5u, 1), Point<1>(), Point<1>(1.5), true);
+  for (auto cell : tria.active_cell_iterators())
+    cell->set_material_id(
+      static_cast<types::material_id>(10.0 * (3.0 + cell->center()[0])));
+  for (auto face : tria.active_face_iterators())
+    if (0.5 < face->center()[0])
+      face->set_all_manifold_ids(42);
+  for (auto cell : tria.active_cell_iterators())
+    cell->set_manifold_id(
+      static_cast<types::material_id>(10.0 * (2.0 + cell->center()[0])));
+}
+
+
+
+template <int dim>
+void
+setup_tria(Triangulation<dim> &tria)
+{
+  static_assert(dim != 1, "not implemented for dim == 1");
+  GridGenerator::hyper_ball(tria);
+  // Make the boundary/material id check more robust
+  for (auto face : tria.active_face_iterators())
+    if (face->at_boundary() && 0.0 < face->center()[0])
+      face->set_boundary_id(42);
+  for (auto cell : tria.active_cell_iterators())
+    if (cell->center() != Point<dim>())
+      {
+        const double angle = std::atan2(cell->center()[0], cell->center()[1]);
+        cell->set_material_id(
+          static_cast<types::material_id>(360 + angle * 180.0 / numbers::PI));
+      }
+  // Make the manifold id check more robust. This relies on the implementation
+  // detail that unassigned manifolds are flat.
+  for (auto face : tria.active_face_iterators())
+    if (-0.1 < face->center()[0])
+      face->set_all_manifold_ids(42);
+  for (auto cell : tria.active_cell_iterators())
+    if (cell->center() != Point<dim>())
+      {
+        const double angle = std::atan2(cell->center()[0], cell->center()[1]);
+        cell->set_manifold_id(
+          static_cast<types::material_id>(361 - angle * 180.0 / numbers::PI));
+      }
+}
+
+
+template <int dim>
+void
+test()
+{
+  {
+    // Check that we can copy a coarse Triangulation
+    Triangulation<dim> tria;
+    setup_tria(tria);
+
+    std::vector<Point<dim>>    vertices;
+    std::vector<CellData<dim>> cells;
+    SubCellData                subcell_data;
+    std::tie(vertices, cells, subcell_data) =
+      GridTools::get_coarse_mesh_description(tria);
+
+    Triangulation<dim> tria_2;
+    tria_2.create_triangulation(vertices, cells, subcell_data);
+    Assert(GridTools::have_same_coarse_mesh(tria, tria_2), ExcInternalError());
+
+    GridOut grid_out;
+    deallog << "Original Triangulation:" << std::endl;
+    grid_out.write_vtk(tria, deallog.get_file_stream());
+    deallog << "Triangulation constructed from an unrefined Triangulation:"
+            << std::endl;
+    grid_out.write_vtk(tria_2, deallog.get_file_stream());
+  }
+
+  {
+    // Check that we can copy a refined Triangulation
+    Triangulation<dim> tria;
+    setup_tria(tria);
+
+    tria.refine_global(2);
+    unsigned int cell_n = 0;
+    for (const auto &cell : tria.active_cell_iterators())
+      {
+        if (cell_n % 3 == 0)
+          cell->set_refine_flag();
+        else if (cell_n % 5 == 0)
+          cell->set_coarsen_flag();
+        ++cell_n;
+      }
+    tria.execute_coarsening_and_refinement();
+
+    std::vector<Point<dim>>    vertices;
+    std::vector<CellData<dim>> cells;
+    SubCellData                subcell_data;
+    std::tie(vertices, cells, subcell_data) =
+      GridTools::get_coarse_mesh_description(tria);
+
+    Triangulation<dim> tria_2;
+    tria_2.create_triangulation(vertices, cells, subcell_data);
+    Assert(GridTools::have_same_coarse_mesh(tria, tria_2), ExcInternalError());
+
+    GridOut grid_out;
+    deallog << "Triangulation constructed from a refined Triangulation:"
+            << std::endl;
+    grid_out.write_vtk(tria_2, deallog.get_file_stream());
+  }
+}
+
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("1d");
+  test<2>();
+  deallog.pop();
+
+  deallog.push("2d");
+  test<2>();
+  deallog.pop();
+
+  deallog.push("3d");
+  test<3>();
+  deallog.pop();
+}
diff --git a/tests/grid/get_coarse_mesh_description_01.output b/tests/grid/get_coarse_mesh_description_01.output
new file mode 100644
index 0000000000..f9dd4c1309
--- /dev/null
+++ b/tests/grid/get_coarse_mesh_description_01.output
@@ -0,0 +1,589 @@
+
+DEAL:1d::Original Triangulation:
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 8 double
+-0.707107 -0.707107 0
+0.707107 -0.707107 0
+-0.292893 -0.292893 0
+0.292893 -0.292893 0
+-0.292893 0.292893 0
+0.292893 0.292893 0
+-0.707107 0.707107 0
+0.707107 0.707107 0
+
+CELLS 14 52
+4 0 1 3 2
+4 0 2 4 6
+4 2 3 5 4
+4 1 7 5 3
+4 6 4 5 7
+2 0 1
+2 0 6
+2 1 3
+2 1 7
+2 2 3
+2 3 5
+2 4 5
+2 6 7
+2 7 5
+
+CELL_TYPES 14
+9 9 9 9 9
+3 3 3 3 3 3 3 3 3
+
+
+CELL_DATA 14
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+540 270 0 450 360
+0 0 -1 42 -1 -1 -1 0 -1
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+181 451 1 271 361
+42 0 42 42 42 42 42 42 42
+
+DEAL:1d::Triangulation constructed from an unrefined Triangulation:
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 8 double
+-0.707107 -0.707107 0
+0.707107 -0.707107 0
+-0.292893 -0.292893 0
+0.292893 -0.292893 0
+-0.292893 0.292893 0
+0.292893 0.292893 0
+-0.707107 0.707107 0
+0.707107 0.707107 0
+
+CELLS 14 52
+4 0 1 3 2
+4 0 2 4 6
+4 2 3 5 4
+4 1 7 5 3
+4 6 4 5 7
+2 0 1
+2 0 6
+2 1 3
+2 1 7
+2 2 3
+2 3 5
+2 4 5
+2 6 7
+2 7 5
+
+CELL_TYPES 14
+9 9 9 9 9
+3 3 3 3 3 3 3 3 3
+
+
+CELL_DATA 14
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+540 270 0 450 360
+0 0 -1 42 -1 -1 -1 0 -1
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+181 451 1 271 361
+42 0 42 42 42 42 42 42 42
+
+DEAL:1d::Triangulation constructed from a refined Triangulation:
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 8 double
+-0.707107 -0.707107 0
+0.707107 -0.707107 0
+-0.292893 -0.292893 0
+0.292893 -0.292893 0
+-0.292893 0.292893 0
+0.292893 0.292893 0
+-0.707107 0.707107 0
+0.707107 0.707107 0
+
+CELLS 14 52
+4 0 1 3 2
+4 0 2 4 6
+4 2 3 5 4
+4 1 7 5 3
+4 6 4 5 7
+2 0 1
+2 0 6
+2 1 3
+2 1 7
+2 2 3
+2 3 5
+2 4 5
+2 6 7
+2 7 5
+
+CELL_TYPES 14
+9 9 9 9 9
+3 3 3 3 3 3 3 3 3
+
+
+CELL_DATA 14
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+540 270 0 450 360
+0 0 -1 42 -1 -1 -1 0 -1
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+181 451 1 271 361
+42 0 42 42 42 42 42 42 42
+
+DEAL:2d::Original Triangulation:
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 8 double
+-0.707107 -0.707107 0
+0.707107 -0.707107 0
+-0.292893 -0.292893 0
+0.292893 -0.292893 0
+-0.292893 0.292893 0
+0.292893 0.292893 0
+-0.707107 0.707107 0
+0.707107 0.707107 0
+
+CELLS 14 52
+4 0 1 3 2
+4 0 2 4 6
+4 2 3 5 4
+4 1 7 5 3
+4 6 4 5 7
+2 0 1
+2 0 6
+2 1 3
+2 1 7
+2 2 3
+2 3 5
+2 4 5
+2 6 7
+2 7 5
+
+CELL_TYPES 14
+9 9 9 9 9
+3 3 3 3 3 3 3 3 3
+
+
+CELL_DATA 14
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+540 270 0 450 360
+0 0 -1 42 -1 -1 -1 0 -1
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+181 451 1 271 361
+42 0 42 42 42 42 42 42 42
+
+DEAL:2d::Triangulation constructed from an unrefined Triangulation:
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 8 double
+-0.707107 -0.707107 0
+0.707107 -0.707107 0
+-0.292893 -0.292893 0
+0.292893 -0.292893 0
+-0.292893 0.292893 0
+0.292893 0.292893 0
+-0.707107 0.707107 0
+0.707107 0.707107 0
+
+CELLS 14 52
+4 0 1 3 2
+4 0 2 4 6
+4 2 3 5 4
+4 1 7 5 3
+4 6 4 5 7
+2 0 1
+2 0 6
+2 1 3
+2 1 7
+2 2 3
+2 3 5
+2 4 5
+2 6 7
+2 7 5
+
+CELL_TYPES 14
+9 9 9 9 9
+3 3 3 3 3 3 3 3 3
+
+
+CELL_DATA 14
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+540 270 0 450 360
+0 0 -1 42 -1 -1 -1 0 -1
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+181 451 1 271 361
+42 0 42 42 42 42 42 42 42
+
+DEAL:2d::Triangulation constructed from a refined Triangulation:
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 8 double
+-0.707107 -0.707107 0
+0.707107 -0.707107 0
+-0.292893 -0.292893 0
+0.292893 -0.292893 0
+-0.292893 0.292893 0
+0.292893 0.292893 0
+-0.707107 0.707107 0
+0.707107 0.707107 0
+
+CELLS 14 52
+4 0 1 3 2
+4 0 2 4 6
+4 2 3 5 4
+4 1 7 5 3
+4 6 4 5 7
+2 0 1
+2 0 6
+2 1 3
+2 1 7
+2 2 3
+2 3 5
+2 4 5
+2 6 7
+2 7 5
+
+CELL_TYPES 14
+9 9 9 9 9
+3 3 3 3 3 3 3 3 3
+
+
+CELL_DATA 14
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+540 270 0 450 360
+0 0 -1 42 -1 -1 -1 0 -1
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+181 451 1 271 361
+42 0 42 42 42 42 42 42 42
+
+DEAL:3d::Original Triangulation:
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 16 double
+-0.211325 -0.211325 -0.211325
+0.211325 -0.211325 -0.211325
+0.211325 -0.211325 0.211325
+-0.211325 -0.211325 0.211325
+-0.211325 0.211325 -0.211325
+0.211325 0.211325 -0.211325
+0.211325 0.211325 0.211325
+-0.211325 0.211325 0.211325
+-0.577350 -0.577350 -0.577350
+0.577350 -0.577350 -0.577350
+0.577350 -0.577350 0.577350
+-0.577350 -0.577350 0.577350
+-0.577350 0.577350 -0.577350
+0.577350 0.577350 -0.577350
+0.577350 0.577350 0.577350
+-0.577350 0.577350 0.577350
+
+CELLS 58 254
+8 0 1 2 3 4 5 6 7
+8 8 9 1 0 12 13 5 4
+8 9 13 14 10 1 5 6 2
+8 11 10 14 15 3 2 6 7
+8 8 0 3 11 12 4 7 15
+8 8 9 10 11 0 1 2 3
+8 12 4 7 15 13 5 6 14
+4 0 3 2 1
+4 0 1 5 4
+4 1 5 6 2
+4 3 2 6 7
+4 4 7 6 5
+4 8 0 1 9
+4 8 11 10 9
+4 8 12 15 11
+4 8 9 13 12
+4 9 13 5 1
+4 9 1 2 10
+4 9 10 14 13
+4 10 14 6 2
+4 11 10 2 3
+4 11 15 14 10
+4 12 4 5 13
+4 12 13 14 15
+4 13 5 6 14
+4 15 14 6 7
+2 0 1
+2 3 2
+2 0 3
+2 1 2
+2 4 7
+2 0 4
+2 3 7
+2 1 5
+2 4 5
+2 5 6
+2 2 6
+2 7 6
+2 8 0
+2 11 3
+2 8 11
+2 12 4
+2 8 12
+2 8 9
+2 9 1
+2 11 10
+2 9 10
+2 12 15
+2 11 15
+2 9 13
+2 12 13
+2 13 5
+2 10 2
+2 10 14
+2 13 14
+2 14 6
+2 15 14
+2 15 7
+
+CELL_TYPES 58
+12 12 12 12 12 12 12
+9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
+3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
+
+CELL_DATA 58
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+360 360 450 360 270 540 360
+-1 -1 -1 -1 -1 -1 0 0 0 -1 -1 42 -1 -1 0 -1 0 -1 -1
+-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 0 -1 0 0 -1 0 0 0 0 0 0 -1 -1 0 0 -1 0 -1
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+361 361 271 361 451 181 361
+42 42 42 42 42 42 42 0 42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42
+DEAL:3d::Triangulation constructed from an unrefined Triangulation:
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 16 double
+-0.211325 -0.211325 -0.211325
+0.211325 -0.211325 -0.211325
+0.211325 -0.211325 0.211325
+-0.211325 -0.211325 0.211325
+-0.211325 0.211325 -0.211325
+0.211325 0.211325 -0.211325
+0.211325 0.211325 0.211325
+-0.211325 0.211325 0.211325
+-0.577350 -0.577350 -0.577350
+0.577350 -0.577350 -0.577350
+0.577350 -0.577350 0.577350
+-0.577350 -0.577350 0.577350
+-0.577350 0.577350 -0.577350
+0.577350 0.577350 -0.577350
+0.577350 0.577350 0.577350
+-0.577350 0.577350 0.577350
+
+CELLS 58 254
+8 0 1 2 3 4 5 6 7
+8 8 9 1 0 12 13 5 4
+8 9 13 14 10 1 5 6 2
+8 11 10 14 15 3 2 6 7
+8 8 0 3 11 12 4 7 15
+8 8 9 10 11 0 1 2 3
+8 12 4 7 15 13 5 6 14
+4 0 3 2 1
+4 0 1 5 4
+4 1 5 6 2
+4 3 2 6 7
+4 4 7 6 5
+4 8 0 1 9
+4 8 11 10 9
+4 8 12 15 11
+4 8 9 13 12
+4 9 13 5 1
+4 9 1 2 10
+4 9 10 14 13
+4 10 14 6 2
+4 11 10 2 3
+4 11 15 14 10
+4 12 4 5 13
+4 12 13 14 15
+4 13 5 6 14
+4 15 14 6 7
+2 0 1
+2 3 2
+2 0 3
+2 1 2
+2 4 7
+2 0 4
+2 3 7
+2 1 5
+2 4 5
+2 5 6
+2 2 6
+2 7 6
+2 8 0
+2 11 3
+2 8 11
+2 12 4
+2 8 12
+2 8 9
+2 9 1
+2 11 10
+2 9 10
+2 12 15
+2 11 15
+2 9 13
+2 12 13
+2 13 5
+2 10 2
+2 10 14
+2 13 14
+2 14 6
+2 15 14
+2 15 7
+
+CELL_TYPES 58
+12 12 12 12 12 12 12
+9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
+3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
+
+CELL_DATA 58
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+360 360 450 360 270 540 360
+-1 -1 -1 -1 -1 -1 0 0 0 -1 -1 42 -1 -1 0 -1 0 -1 -1
+-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 0 -1 0 0 -1 0 0 0 0 0 0 -1 -1 0 0 -1 0 -1
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+361 361 271 361 451 181 361
+42 42 42 42 42 42 42 0 42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42
+DEAL:3d::Triangulation constructed from a refined Triangulation:
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 16 double
+-0.211325 -0.211325 -0.211325
+0.211325 -0.211325 -0.211325
+0.211325 -0.211325 0.211325
+-0.211325 -0.211325 0.211325
+-0.211325 0.211325 -0.211325
+0.211325 0.211325 -0.211325
+0.211325 0.211325 0.211325
+-0.211325 0.211325 0.211325
+-0.577350 -0.577350 -0.577350
+0.577350 -0.577350 -0.577350
+0.577350 -0.577350 0.577350
+-0.577350 -0.577350 0.577350
+-0.577350 0.577350 -0.577350
+0.577350 0.577350 -0.577350
+0.577350 0.577350 0.577350
+-0.577350 0.577350 0.577350
+
+CELLS 58 254
+8 0 1 2 3 4 5 6 7
+8 8 9 1 0 12 13 5 4
+8 9 13 14 10 1 5 6 2
+8 11 10 14 15 3 2 6 7
+8 8 0 3 11 12 4 7 15
+8 8 9 10 11 0 1 2 3
+8 12 4 7 15 13 5 6 14
+4 0 3 2 1
+4 0 1 5 4
+4 1 5 6 2
+4 3 2 6 7
+4 4 7 6 5
+4 8 0 1 9
+4 8 11 10 9
+4 8 12 15 11
+4 8 9 13 12
+4 9 13 5 1
+4 9 1 2 10
+4 9 10 14 13
+4 10 14 6 2
+4 11 10 2 3
+4 11 15 14 10
+4 12 4 5 13
+4 12 13 14 15
+4 13 5 6 14
+4 15 14 6 7
+2 0 1
+2 3 2
+2 0 3
+2 1 2
+2 4 7
+2 0 4
+2 3 7
+2 1 5
+2 4 5
+2 5 6
+2 2 6
+2 7 6
+2 8 0
+2 11 3
+2 8 11
+2 12 4
+2 8 12
+2 8 9
+2 9 1
+2 11 10
+2 9 10
+2 12 15
+2 11 15
+2 9 13
+2 12 13
+2 13 5
+2 10 2
+2 10 14
+2 13 14
+2 14 6
+2 15 14
+2 15 7
+
+CELL_TYPES 58
+12 12 12 12 12 12 12
+9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
+3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
+
+CELL_DATA 58
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+360 360 450 360 270 540 360
+-1 -1 -1 -1 -1 -1 0 0 0 -1 -1 42 -1 -1 0 -1 0 -1 -1
+-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 0 -1 0 0 -1 0 0 0 0 0 0 -1 -1 0 0 -1 0 -1
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+361 361 271 361 451 181 361
+42 42 42 42 42 42 42 0 42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42 42