From 67e6946f84a56c67daa2b14ddcfbcd9bb62ad320 Mon Sep 17 00:00:00 2001 From: heister Date: Sun, 17 Feb 2013 05:34:31 +0000 Subject: [PATCH] introduce GridGenerator::extrude_triangulation git-svn-id: https://svn.dealii.org/trunk@28431 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 5 + deal.II/include/deal.II/grid/grid_generator.h | 15 +++ deal.II/source/grid/grid_generator.cc | 95 +++++++++++++++++++ 3 files changed, 115 insertions(+) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index dc0f680338..527ad31d64 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -193,6 +193,11 @@ DoFHandler, in particular removal of specializations.

Specific improvements

    +
  1. New: GridGenerator::extrude_triangulation() allows +you to extrude a 2d mesh to turn it into a 3d mesh. +
    +(Timo Heister, 2013/02/16) +
  2. PETScWrappers::MPI::Vector with ghost entries are read-only now.
    diff --git a/deal.II/include/deal.II/grid/grid_generator.h b/deal.II/include/deal.II/grid/grid_generator.h index 68fa62da86..b92327f0bf 100644 --- 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 &triangulation_2, Triangulation &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 6e43558f2b..0425941ce5 100644 --- a/deal.II/source/grid/grid_generator.cc +++ b/deal.II/source/grid/grid_generator.cc @@ -3376,6 +3376,101 @@ merge_triangulations (const Triangulation &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 > points(n_slices*input.n_vertices()); + std::vector > cells; + cells.reserve((n_slices-1)*input.n_active_cells()); + + for (unsigned int slice=0;slice & 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 this_cell; + for (unsigned int v=0; v::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;sliceface(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 -- 2.39.5