introduce GridGenerator::extrude_triangulation

git-svn-id: https://svn.dealii.org/trunk@28431 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index dc0f680338028f777837025052601134e1372bd5..527ad31d6464bb7b4d5b63aa858cea0b7124af6f 100644 (file)
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -193,6 +193,11 @@ DoFHandler, in particular removal of specializations.
 <h3>Specific improvements</h3>
 
 <ol>
+<li> New: GridGenerator::extrude_triangulation() allows
+you to extrude a 2d mesh to turn it into a 3d mesh.
+<br>
+(Timo Heister, 2013/02/16)
+
 <li> PETScWrappers::MPI::Vector with ghost entries are read-only
 now.
 <br>

diff --git a/deal.II/include/deal.II/grid/grid_generator.h b/deal.II/include/deal.II/grid/grid_generator.h
index 68fa62da86efb5df044f01e3ab6ca0a1cc238767..b92327f0bfea56b7522cb265880fa3042bfb3603 100644 (file)
--- a/deal.II/include/deal.II/grid/grid_generator.h
+++ b/deal.II/include/deal.II/grid/grid_generator.h
@@ -673,6 +673,21 @@ public:
                         const Triangulation<dim, spacedim> &triangulation_2,
                         Triangulation<dim, spacedim>       &result);
 
+
+  /**
+   * Take a 2d Triangulation that is being extruded in z direction
+   * by the total height of @p height using @p n_slices slices (minimum is 2).
+   * The boundary indicators of the faces of @p input are going to be assigned
+   * to the corresponding side walls in z direction. The bottom and top
+   * get the next two free boundary indicators.
+   */
+  static
+  void
+  extrude_triangulation(const Triangulation<2, 2> & input,
+                        const unsigned int n_slices,
+                        const double height,
+                        Triangulation<3,3> &result);
+
   /**
    * This function transformes the @p Triangulation @p tria smoothly to a
    * domain that is described by the boundary points in the map @p

diff --git a/deal.II/source/grid/grid_generator.cc b/deal.II/source/grid/grid_generator.cc
index 6e43558f2bb4791420c17cecd333daa51dbb6472..0425941ce551095edaa8bc637b83ec42f52e3c1b 100644 (file)
--- a/deal.II/source/grid/grid_generator.cc
+++ b/deal.II/source/grid/grid_generator.cc
@@ -3376,6 +3376,101 @@ merge_triangulations (const Triangulation<dim, spacedim> &triangulation_1,
   result.create_triangulation (vertices, cells, subcell_data);
 }
 
+void
+GridGenerator::
+extrude_triangulation(const Triangulation<2, 2> & input,
+                      const unsigned int n_slices,
+                      const double height,
+                      Triangulation<3,3> &result)
+{
+  Assert (input.n_levels() == 1,
+          ExcMessage ("The input triangulations must be coarse meshes."));
+  Assert(result.n_cells()==0, ExcMessage("resultin Triangulation need to be empty upon calling extrude_triangulation."));
+  Assert(height>0, ExcMessage("The height in extrude_triangulation needs to be positive."));
+  Assert(n_slices>=2, ExcMessage("The number of slices in extrude_triangulation needs to be at least 2."));
+
+  std::vector<Point<3> > points(n_slices*input.n_vertices());
+  std::vector<CellData<3> > cells;
+  cells.reserve((n_slices-1)*input.n_active_cells());
+
+  for (unsigned int slice=0;slice<n_slices;++slice)
+    {
+      for (unsigned int i=0;i<input.n_vertices();++i)
+
+      {
+        const Point<2> & v = input.get_vertices()[i];
+        points[i+slice*input.n_vertices()](0) = v(0);
+        points[i+slice*input.n_vertices()](1) = v(1);
+        points[i+slice*input.n_vertices()](2) = height * slice / n_slices;
+      }
+    }
+
+  for (typename Triangulation<2,2>::cell_iterator
+      cell = input.begin(); cell != input.end(); ++cell)
+    {
+      for (unsigned int slice=0;slice<n_slices-1;++slice)
+        {
+          CellData<3> this_cell;
+          for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
+            {
+              this_cell.vertices[v]
+                                 = cell->vertex_index(v)+slice*input.n_vertices();
+              this_cell.vertices[v+GeometryInfo<2>::vertices_per_cell]
+                                 = cell->vertex_index(v)+(slice+1)*input.n_vertices();
+            }
+
+          this_cell.material_id = cell->material_id();
+          cells.push_back(this_cell);
+        }
+    }
+
+  SubCellData s;
+  types::boundary_id bid=0;
+  s.boundary_quads.reserve(input.n_active_lines()*(n_slices-1) + input.n_active_cells()*2);
+  for (typename Triangulation<2,2>::cell_iterator
+      cell = input.begin(); cell != input.end(); ++cell)
+    {
+      CellData<2> quad;
+      for (unsigned int f=0; f<4;++f)
+        if (cell->at_boundary(f))
+          {
+            quad.boundary_id = cell->face(f)->boundary_indicator();
+            bid = std::max(bid, quad.boundary_id);
+          for (unsigned int slice=0;slice<n_slices-1;++slice)
+            {
+              quad.vertices[0] = cell->face(f)->vertex_index(0)+slice*input.n_vertices();
+              quad.vertices[1] = cell->face(f)->vertex_index(1)+slice*input.n_vertices();
+              quad.vertices[2] = cell->face(f)->vertex_index(0)+(slice+1)*input.n_vertices();
+              quad.vertices[3] = cell->face(f)->vertex_index(1)+(slice+1)*input.n_vertices();
+              s.boundary_quads.push_back(quad);
+            }
+          }
+    }
+
+  for (typename Triangulation<2,2>::cell_iterator
+      cell = input.begin(); cell != input.end(); ++cell)
+    {
+      CellData<2> quad;
+      quad.boundary_id = bid + 1;
+      quad.vertices[0] = cell->vertex_index(0);
+      quad.vertices[1] = cell->vertex_index(1);
+      quad.vertices[2] = cell->vertex_index(2);
+      quad.vertices[3] = cell->vertex_index(3);
+      s.boundary_quads.push_back(quad);
+
+      quad.boundary_id = bid + 2;
+      for (int i=0;i<4;++i)
+        quad.vertices[i] += (n_slices-1)*input.n_vertices();
+      s.boundary_quads.push_back(quad);
+    }
+
+
+
+  result.create_triangulation (
+    points,
+    cells,
+    s);
+}
 
 
 // make the following function inline. this is so that it becomes marked