From ff24cca5dd5a236c63b32296762cb5b07dc1e284 Mon Sep 17 00:00:00 2001
From: Luca Heltai <luca.heltai@sissa.it>
Date: Tue, 19 Aug 2014 14:37:08 +0200
Subject: [PATCH] GridGenerator::flatten_triangulation.

Added GridTools::flatten_triangulation + one test.

Fixed WB comments.

Moved flatten triangulation to GridGenerator.

Fixed TH and WB comments.

Fixed indentation issue.
---
 doc/news/changes.h                    |  18 ++
 include/deal.II/grid/grid_generator.h |  33 ++++
 include/deal.II/grid/grid_tools.h     |   1 +
 source/grid/grid_generator.cc         |  83 +++++++++
 source/grid/grid_generator.inst.in    |  13 ++
 source/grid/grid_out.inst.in          |   3 +
 source/grid/grid_tools.cc             |   2 -
 source/grid/grid_tools.inst.in        |   1 +
 tests/grid/flatten_grid_01.cc         |  63 +++++++
 tests/grid/flatten_grid_01.output     | 231 ++++++++++++++++++++++++++
 10 files changed, 446 insertions(+), 2 deletions(-)
 create mode 100644 tests/grid/flatten_grid_01.cc
 create mode 100644 tests/grid/flatten_grid_01.output

diff --git a/doc/news/changes.h b/doc/news/changes.h
index 2e6d6f72b8..859da6936a 100644
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -97,6 +97,24 @@ inconvenience this causes.
 
 
 <ol>
+  <li> New: There is now a GridGenerator::flatten_triangulation()
+  taking a Triangulation<dim, spacedim_1> as input and returning
+  a Triangulation<dim, spacedim_2> as output. The output
+  triangulation will contain a single level with all active
+  cells of the input triangulation, and will be topologically
+  equivalent to the input triangulation. If the two spacedimensions
+  are equal, then this function will copy the triangulation
+  removing all levels, e.g., flattening it. If the two spacedimensions
+  are different, then this function will copy the vertices only
+  up to the smallest spacedimension parameter. <br>
+  Using this function, you can create a Triangulation<2,3> from
+  a Triangulation<2,2> or project to the plane z=0 your
+  Triangulation<2,3>. No checks are performed on the validity of
+  the resulting Triangulation.
+  <br>
+  (Luca Heltai, 2014/08/19)
+  </li>
+
   <li> Fixed: Newer versions of GCC (e.g. 4.9.x) are no longer compatible
   with BOOST 1.46. Consequently, the CMake scripts now require at least
   BOOST 1.48 in order to use a BOOST version found on the system. If no
diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h
index 44a45c0edf..14ea1a5924 100644
--- a/include/deal.II/grid/grid_generator.h
+++ b/include/deal.II/grid/grid_generator.h
@@ -690,6 +690,39 @@ namespace GridGenerator
                          const double               height,
                          Triangulation<3,3>        &result);
 
+  /**
+   * Given an input triangulation @p in_tria, this function makes a new
+   * flat triangulation @p out_tria which contains a single level with
+   * all active cells of the input triangulation. If spacedim1 and
+   * spacedim2 are different, only the smallest spacedim components of
+   * the vertices are copied over. This is useful to create a
+   * Triangulation<2,3> out of a Triangulation<2,2>, or to project a
+   * Triangulation<2,3> into a Triangulation<2,2>, by neglecting the z
+   * components of the vertices.
+   *
+   * No internal checks are performed on the vertices, which are
+   * assumed to make sense topologically in the target #spacedim2
+   * dimensional space. If this is not the case, you will encounter
+   * problems when using the triangulation later on.
+   *
+   * All informations about cell manifold_ids and material ids are
+   * copied from one triangulation to the other, and only the boundary
+   * manifold_ids and boundary_ids are copied over from the faces of
+   * @p in_tria to the faces of @p out_tria. If you need to specify
+   * manifold ids on interior faces, they have to be specified
+   * manually after the triangulation is created.
+   *
+   * This function will fail if the input Triangulation is of type
+   * parallel::distributed::Triangulation, as well as when the input
+   * Triangulation contains hanging nodes.
+   *
+   * @author Luca Heltai, 2014
+   */
+  template <int dim, int spacedim1, int spacedim2>
+  void flatten_triangulation(const Triangulation<dim,spacedim1> &in_tria,
+			     Triangulation<dim,spacedim2> &out_tria);
+
+
   /**
    * This function transforms the @p Triangulation @p tria smoothly to a
    * domain that is described by the boundary points in the map @p
diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h
index 76dbe802e3..4355d8984a 100644
--- a/include/deal.II/grid/grid_tools.h
+++ b/include/deal.II/grid/grid_tools.h
@@ -1001,6 +1001,7 @@ namespace GridTools
   std::vector<typename Container::active_cell_iterator>
   get_patch_around_cell(const typename Container::active_cell_iterator &cell);
 
+
   /*@}*/
   /**
    *  @name Lower-dimensional meshes for parts of higher-dimensional meshes 
diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc
index 4f4cf102bf..3eca78a3f7 100644
--- a/source/grid/grid_generator.cc
+++ b/source/grid/grid_generator.cc
@@ -38,6 +38,8 @@
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/numerics/matrix_tools.h>
 
+#include <deal.II/distributed/tria.h>
+
 #include <iostream>
 #include <cmath>
 #include <limits>
@@ -3866,6 +3868,87 @@ namespace GridGenerator
             }
       }
   }
+  
+  template <int dim, int spacedim1, int spacedim2>
+  void flatten_triangulation(const Triangulation<dim, spacedim1> &in_tria,
+			     Triangulation<dim,spacedim2> &out_tria) 
+  {
+    const parallel::distributed::Triangulation<dim, spacedim1> * pt =
+      dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim1> *>(&in_tria);
+    
+    Assert (pt == NULL, 
+	    ExcMessage("Cannot use this function on parallel::distributed::Triangulation."));
+
+    std::vector<Point<spacedim2> > v;
+    std::vector<CellData<dim> > cells;
+    SubCellData subcelldata;
+  
+    const unsigned int spacedim = std::min(spacedim1,spacedim2);
+    const std::vector<Point<spacedim1> > &in_vertices = in_tria.get_vertices();
+  
+    v.resize(in_vertices.size());
+    for(unsigned int i=0; i<in_vertices.size(); ++i)
+      for(unsigned int d=0; d<spacedim; ++d)
+	v[i][d] = in_vertices[i][d];
+
+    cells.resize(in_tria.n_active_cells());
+    typename Triangulation<dim,spacedim1>::active_cell_iterator 
+      cell = in_tria.begin_active(),
+      endc = in_tria.end();
+  
+    for(unsigned int id=0; cell != endc; ++cell, ++id) 
+      {
+	for(unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i) 
+	  cells[id].vertices[i] = cell->vertex_index(i);
+	cells[id].material_id = cell->material_id();
+	cells[id].manifold_id = cell->manifold_id();
+      }
+
+    if(dim>1) {
+      typename Triangulation<dim,spacedim1>::active_face_iterator
+	face = in_tria.begin_active_face(),
+	endf = in_tria.end_face();
+      
+      // Face counter for both dim == 2 and dim == 3
+      unsigned int f=0;
+      switch(dim) {
+      case 2:
+	{
+	  subcelldata.boundary_lines.resize(in_tria.n_active_faces());
+	  for(; face != endf; ++face)
+	    if(face->at_boundary()) 
+	      {
+		for(unsigned int i=0; i<GeometryInfo<dim>::vertices_per_face; ++i) 
+		  subcelldata.boundary_lines[f].vertices[i] = face->vertex_index(i);
+		subcelldata.boundary_lines[f].boundary_id = face->boundary_indicator();
+		subcelldata.boundary_lines[f].manifold_id = face->manifold_id();
+		++f;
+	      }
+	  subcelldata.boundary_lines.resize(f);
+	}
+	break;
+      case 3:
+	{
+	  subcelldata.boundary_quads.resize(in_tria.n_active_faces());
+	  for(; face != endf; ++face) 
+	    if(face->at_boundary()) 
+	      {
+		for(unsigned int i=0; i<GeometryInfo<dim>::vertices_per_face; ++i) 
+		  subcelldata.boundary_quads[f].vertices[i] = face->vertex_index(i);
+		subcelldata.boundary_quads[f].boundary_id = face->boundary_indicator();
+		subcelldata.boundary_quads[f].manifold_id = face->manifold_id();
+		++f;
+	      }
+	  subcelldata.boundary_quads.resize(f);
+	}
+	break;
+      default:
+	Assert(false, ExcInternalError());
+      }
+    }
+    out_tria.create_triangulation(v, cells, subcelldata);
+  }
+
 }
 
 // explicit instantiations
diff --git a/source/grid/grid_generator.inst.in b/source/grid/grid_generator.inst.in
index 4101f47e36..74600bcdda 100644
--- a/source/grid/grid_generator.inst.in
+++ b/source/grid/grid_generator.inst.in
@@ -99,3 +99,16 @@ namespace GridGenerator \{
 \}  
  }
 
+
+for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS; deal_II_space_dimension_2 : SPACE_DIMENSIONS)
+{
+namespace GridGenerator \{
+#if (deal_II_dimension <= deal_II_space_dimension) && (deal_II_dimension <= deal_II_space_dimension_2)
+     template 
+     void 
+     flatten_triangulation<>(const Triangulation<deal_II_dimension,deal_II_space_dimension> &,
+                             Triangulation<deal_II_dimension,deal_II_space_dimension_2>&);
+#endif
+\}    
+}
+
diff --git a/source/grid/grid_out.inst.in b/source/grid/grid_out.inst.in
index c9cbe7f4bb..bb3390e29e 100644
--- a/source/grid/grid_out.inst.in
+++ b/source/grid/grid_out.inst.in
@@ -96,6 +96,9 @@ for (deal_II_dimension : DIMENSIONS)
    template void GridOut::write_ucd
      (const Triangulation<1, 3> &,
       std::ostream &) const;
+   template void GridOut::write_msh
+     (const Triangulation<1, 3> &,
+      std::ostream &) const;
 #endif
   }
 
diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index c2ba9e2d4a..01d3e4e9fc 100644
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -2595,8 +2595,6 @@ next_cell:
 #endif
   }
 
-
-
 } /* namespace GridTools */
 
 
diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in
index 6191916d47..6198f016cb 100644
--- a/source/grid/grid_tools.inst.in
+++ b/source/grid/grid_tools.inst.in
@@ -373,3 +373,4 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS
 #endif
 }
 
+
diff --git a/tests/grid/flatten_grid_01.cc b/tests/grid/flatten_grid_01.cc
new file mode 100644
index 0000000000..c795eecd81
--- /dev/null
+++ b/tests/grid/flatten_grid_01.cc
@@ -0,0 +1,63 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2004 - 2013 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.
+//
+// ---------------------------------------------------------------------
+
+// Generate a grid, refine it once, flatten it and output the result.
+
+#include "../tests.h"
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_generator.h>
+
+template <int dim, int spacedim1, int spacedim2> 
+void test() {
+  deallog << "Testing <" << dim << "," << spacedim1
+	  << "> VS <" << dim << "," << spacedim2 
+	  << ">" << std::endl;
+  
+  Triangulation<dim, spacedim1> tria1;
+  GridGenerator::hyper_cube(tria1);
+  tria1.refine_global(1);
+  
+  Triangulation<dim, spacedim2> tria2;
+  GridGenerator::flatten_triangulation(tria1, tria2);
+  GridOut go;
+  go.write_msh(tria2, deallog.get_file_stream());
+}
+
+int main()
+{
+  initlog();
+  test<1,1,1>();
+  test<1,1,2>();
+  test<1,1,3>();
+  // 
+  test<1,2,1>();
+  test<1,2,2>();
+  test<1,2,3>();
+  // 
+  test<1,3,1>();
+  test<1,3,2>();
+  test<1,3,3>();
+  //
+  test<2,2,2>();
+  test<2,2,3>();
+  //
+  test<2,3,2>();
+  test<2,3,3>();
+  //
+  test<3,3,3>();
+}
diff --git a/tests/grid/flatten_grid_01.output b/tests/grid/flatten_grid_01.output
new file mode 100644
index 0000000000..7304f8a15a
--- /dev/null
+++ b/tests/grid/flatten_grid_01.output
@@ -0,0 +1,231 @@
+
+DEAL::Testing <1,1> VS <1,1>
+$NOD
+3
+1  0.00000 0 0
+2  1.00000 0 0
+3  0.500000 0 0
+$ENDNOD
+$ELM
+2
+1 1 0 0 2 1 3 
+2 1 0 0 2 3 2 
+$ENDELM
+DEAL::Testing <1,1> VS <1,2>
+$NOD
+3
+1  0.00000 0.00000 0
+2  1.00000 0.00000 0
+3  0.500000 0.00000 0
+$ENDNOD
+$ELM
+2
+1 1 0 0 2 1 3 
+2 1 0 0 2 3 2 
+$ENDELM
+DEAL::Testing <1,1> VS <1,3>
+$NOD
+3
+1  0.00000 0.00000 0.00000
+2  1.00000 0.00000 0.00000
+3  0.500000 0.00000 0.00000
+$ENDNOD
+$ELM
+2
+1 1 0 0 2 1 3 
+2 1 0 0 2 3 2 
+$ENDELM
+DEAL::Testing <1,2> VS <1,1>
+$NOD
+3
+1  0.00000 0 0
+2  1.00000 0 0
+3  0.500000 0 0
+$ENDNOD
+$ELM
+2
+1 1 0 0 2 1 3 
+2 1 0 0 2 3 2 
+$ENDELM
+DEAL::Testing <1,2> VS <1,2>
+$NOD
+3
+1  0.00000 0.00000 0
+2  1.00000 0.00000 0
+3  0.500000 0.00000 0
+$ENDNOD
+$ELM
+2
+1 1 0 0 2 1 3 
+2 1 0 0 2 3 2 
+$ENDELM
+DEAL::Testing <1,2> VS <1,3>
+$NOD
+3
+1  0.00000 0.00000 0.00000
+2  1.00000 0.00000 0.00000
+3  0.500000 0.00000 0.00000
+$ENDNOD
+$ELM
+2
+1 1 0 0 2 1 3 
+2 1 0 0 2 3 2 
+$ENDELM
+DEAL::Testing <1,3> VS <1,1>
+$NOD
+3
+1  0.00000 0 0
+2  1.00000 0 0
+3  0.500000 0 0
+$ENDNOD
+$ELM
+2
+1 1 0 0 2 1 3 
+2 1 0 0 2 3 2 
+$ENDELM
+DEAL::Testing <1,3> VS <1,2>
+$NOD
+3
+1  0.00000 0.00000 0
+2  1.00000 0.00000 0
+3  0.500000 0.00000 0
+$ENDNOD
+$ELM
+2
+1 1 0 0 2 1 3 
+2 1 0 0 2 3 2 
+$ENDELM
+DEAL::Testing <1,3> VS <1,3>
+$NOD
+3
+1  0.00000 0.00000 0.00000
+2  1.00000 0.00000 0.00000
+3  0.500000 0.00000 0.00000
+$ENDNOD
+$ELM
+2
+1 1 0 0 2 1 3 
+2 1 0 0 2 3 2 
+$ENDELM
+DEAL::Testing <2,2> VS <2,2>
+$NOD
+9
+1  0.00000 0.00000 0
+2  1.00000 0.00000 0
+3  0.00000 1.00000 0
+4  1.00000 1.00000 0
+5  0.500000 0.00000 0
+6  0.00000 0.500000 0
+7  1.00000 0.500000 0
+8  0.500000 1.00000 0
+9  0.500000 0.500000 0
+$ENDNOD
+$ELM
+4
+1 3 0 0 4 1 5 9 6 
+2 3 0 0 4 5 2 7 9 
+3 3 0 0 4 6 9 8 3 
+4 3 0 0 4 9 7 4 8 
+$ENDELM
+DEAL::Testing <2,2> VS <2,3>
+$NOD
+9
+1  0.00000 0.00000 0.00000
+2  1.00000 0.00000 0.00000
+3  0.00000 1.00000 0.00000
+4  1.00000 1.00000 0.00000
+5  0.500000 0.00000 0.00000
+6  0.00000 0.500000 0.00000
+7  1.00000 0.500000 0.00000
+8  0.500000 1.00000 0.00000
+9  0.500000 0.500000 0.00000
+$ENDNOD
+$ELM
+4
+1 3 0 0 4 1 5 9 6 
+2 3 0 0 4 5 2 7 9 
+3 3 0 0 4 6 9 8 3 
+4 3 0 0 4 9 7 4 8 
+$ENDELM
+DEAL::Testing <2,3> VS <2,2>
+$NOD
+9
+1  0.00000 0.00000 0
+2  1.00000 0.00000 0
+3  0.00000 1.00000 0
+4  1.00000 1.00000 0
+5  0.500000 0.00000 0
+6  0.00000 0.500000 0
+7  1.00000 0.500000 0
+8  0.500000 1.00000 0
+9  0.500000 0.500000 0
+$ENDNOD
+$ELM
+4
+1 3 0 0 4 1 5 9 6 
+2 3 0 0 4 5 2 7 9 
+3 3 0 0 4 6 9 8 3 
+4 3 0 0 4 9 7 4 8 
+$ENDELM
+DEAL::Testing <2,3> VS <2,3>
+$NOD
+9
+1  0.00000 0.00000 0.00000
+2  1.00000 0.00000 0.00000
+3  0.00000 1.00000 0.00000
+4  1.00000 1.00000 0.00000
+5  0.500000 0.00000 0.00000
+6  0.00000 0.500000 0.00000
+7  1.00000 0.500000 0.00000
+8  0.500000 1.00000 0.00000
+9  0.500000 0.500000 0.00000
+$ENDNOD
+$ELM
+4
+1 3 0 0 4 1 5 9 6 
+2 3 0 0 4 5 2 7 9 
+3 3 0 0 4 6 9 8 3 
+4 3 0 0 4 9 7 4 8 
+$ENDELM
+DEAL::Testing <3,3> VS <3,3>
+$NOD
+27
+1  0.00000 0.00000 0.00000
+2  1.00000 0.00000 0.00000
+3  0.00000 1.00000 0.00000
+4  1.00000 1.00000 0.00000
+5  0.00000 0.00000 1.00000
+6  1.00000 0.00000 1.00000
+7  0.00000 1.00000 1.00000
+8  1.00000 1.00000 1.00000
+9  0.500000 0.00000 0.00000
+10  0.00000 0.500000 0.00000
+11  0.00000 0.00000 0.500000
+12  1.00000 0.500000 0.00000
+13  1.00000 0.00000 0.500000
+14  0.500000 1.00000 0.00000
+15  0.00000 1.00000 0.500000
+16  1.00000 1.00000 0.500000
+17  0.500000 0.00000 1.00000
+18  0.00000 0.500000 1.00000
+19  1.00000 0.500000 1.00000
+20  0.500000 1.00000 1.00000
+21  0.500000 0.00000 0.500000
+22  0.500000 0.500000 0.00000
+23  0.00000 0.500000 0.500000
+24  1.00000 0.500000 0.500000
+25  0.500000 1.00000 0.500000
+26  0.500000 0.500000 1.00000
+27  0.500000 0.500000 0.500000
+$ENDNOD
+$ELM
+8
+1 5 0 0 8 1 9 21 11 10 22 27 23 
+2 5 0 0 8 9 2 13 21 22 12 24 27 
+3 5 0 0 8 10 22 27 23 3 14 25 15 
+4 5 0 0 8 22 12 24 27 14 4 16 25 
+5 5 0 0 8 11 21 17 5 23 27 26 18 
+6 5 0 0 8 21 13 6 17 27 24 19 26 
+7 5 0 0 8 23 27 26 18 15 25 20 7 
+8 5 0 0 8 27 24 19 26 25 16 8 20 
+$ENDELM
-- 
2.39.5