From: pauletti Date: Mon, 15 Nov 2010 17:09:16 +0000 (+0000) Subject: Added the function GridTools::extract_boundary_mesh (...). X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=686e5fe34d3d4d01e92ade39986f2ec707994941;p=dealii-svn.git Added the function GridTools::extract_boundary_mesh (...). git-svn-id: https://svn.dealii.org/trunk@22747 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/grid/grid_tools.h b/deal.II/include/deal.II/grid/grid_tools.h index 606709fc68..15237f9475 100644 --- a/deal.II/include/deal.II/grid/grid_tools.h +++ b/deal.II/include/deal.II/grid/grid_tools.h @@ -795,6 +795,28 @@ class GridTools fix_up_distorted_child_cells (const typename Triangulation::DistortedCellList &distorted_cells, Triangulation &triangulation); + /** + This function implements a boundary subgrid extraction. + Given a -Triangulation (the "volume mesh") + the function extracts a subset of its boundary (the "surface mesh"). + The boundary to be extracted is specified by a list of boundary_ids. + If none is specified the whole boundary will be extracted. + + It also builds a mapping linking the cells on the surface mesh + to the corresponding faces on the volume one. + + + */ + template + static void + extract_boundary_mesh (const Triangulation &volume_mesh, + Triangulation &surface_mesh, + std::map::cell_iterator, + typename Triangulation::face_iterator> + &surface_to_volume_mapping, + const std::set &boundary_ids + = std::set()); + /** * Exception */ @@ -842,6 +864,8 @@ class GridTools unsigned int, << "The given vertex " << arg1 << " is not used in the given triangulation"); + + }; @@ -986,6 +1010,11 @@ GridTools::cell_measure<2>(const std::vector > &all_vertices, // const unsigned int (&vertex_indices) [GeometryInfo<2>::vertices_per_cell]); + + + + + DEAL_II_NAMESPACE_CLOSE /*---------------------------- grid_tools.h ---------------------------*/ diff --git a/deal.II/source/grid/grid_tools.cc b/deal.II/source/grid/grid_tools.cc index 1cd240ce10..862fe7a252 100644 --- a/deal.II/source/grid/grid_tools.cc +++ b/deal.II/source/grid/grid_tools.cc @@ -1812,6 +1812,138 @@ fix_up_distorted_child_cells (const typename Triangulation::Distor +template +void +GridTools::extract_boundary_mesh (const Triangulation &volume_mesh, + Triangulation &surface_mesh, + std::map::cell_iterator, + typename Triangulation::face_iterator> + &surface_to_volume_mapping, + const std::set &boundary_ids) +{ + + //

Assumptions

+ // We are relying heavily on that + // Triangulation::create_triangulation(...) will keep the order + // pass by CellData and that it will not reorder the vertices. + + //

TODO

+ //
    + //
  • std::set &boundary_ids needs to be set + // to a default value of {0}, no is set to empty + //
  • Add mapping from volume to surface + //
  • Add some assertions + //
  • Create a function (reconcile?) to reconcile the triangulations + // after refinement/coarsening is performed in one of them + //
  • Create a class that manages the triangulations, as direct + // access to refining the triangulations will most likely lead + // to inconsitencies + //
+ + surface_to_volume_mapping.clear(); // Is this reasonable? + + const unsigned int boundary_dim = dim-1; //dimension of the boundary mesh + + // First create surface mesh and mapping from only level(0) cells of volume_mesh + + std::vector::face_iterator> + mapping; // temporary map for level==0 + + + + typename Triangulation::cell_iterator + cell = volume_mesh.begin(0), + endc = volume_mesh.end(0); + typename Triangulation::face_iterator face; + + std::vector< bool > touched (volume_mesh.n_vertices(), false); + std::vector< CellData< boundary_dim > > cells; + std::vector< Point > vertices; + + std::map map_vert_index; //volume vertex indices to surf ones + + unsigned int v_index; + CellData< boundary_dim > c_data; + + for (; cell!=endc; ++cell) + + for (unsigned int i=0; + i < GeometryInfo::faces_per_cell; ++i){ + + face = cell->face(i); + + if ( face->at_boundary() + && + (boundary_ids.empty() || + ( boundary_ids.find(face->boundary_indicator()) != boundary_ids.end())) + ) { + for (unsigned int j=0; + j::vertices_per_cell; ++j) { + + v_index = face->vertex_index(j); + + if ( !touched[v_index] ){ + vertices.push_back(face->vertex(j)); + map_vert_index[v_index] = vertices.size() - 1; + touched[v_index] = true; + } + + c_data.vertices[j] = map_vert_index[v_index]; + c_data.material_id = 0; + + } + + cells.push_back(c_data); + mapping.push_back(face); + + } + } + + // create level 0 surface triangulation + surface_mesh.create_triangulation (vertices, cells, SubCellData()); + + { // Make the actual mapping + typename Triangulation::active_cell_iterator + cell = surface_mesh.begin(0), + endc = surface_mesh.end(0); + for (; cell!=endc; ++cell) + surface_to_volume_mapping[cell] = mapping.at(cell->index()); + } + + do { + bool changed = false; + typename Triangulation::active_cell_iterator + cell = surface_mesh.begin_active(), + endc = surface_mesh.end(); + + for (; cell!=endc; ++cell) + if (surface_to_volume_mapping[cell]->has_children() == true ) { + cell->set_refine_flag (); + changed = true; + } + + + + if (changed){ + + surface_mesh.execute_coarsening_and_refinement (); + + typename Triangulation::cell_iterator + cell = surface_mesh.begin(), + endc = surface_mesh.end(); + for (; cell!=endc; ++cell) + for (unsigned int c=0; cn_children(); c++) + if (surface_to_volume_mapping.find(cell->child(c)) == surface_to_volume_mapping.end()) + surface_to_volume_mapping[cell->child(c)] + = surface_to_volume_mapping[cell]->child(c); + } + else + break; + } + while (true); + +} + // explicit instantiations #include "grid_tools.inst" diff --git a/deal.II/source/grid/grid_tools.inst.in b/deal.II/source/grid/grid_tools.inst.in index 0a1cc0795d..4464e78ed5 100644 --- a/deal.II/source/grid/grid_tools.inst.in +++ b/deal.II/source/grid/grid_tools.inst.in @@ -57,6 +57,15 @@ for (deal_II_dimension : DIMENSIONS) template double GridTools::diameter (const Triangulation &); + template + void + GridTools::extract_boundary_mesh (const Triangulation &volume_mesh, + Triangulation &surface_mesh, + std::map< Triangulation::cell_iterator, + Triangulation::face_iterator> + &surface_to_volume_mapping, + const std::set &boundary_ids + = std::set()); #endif #if deal_II_dimension == 2 @@ -204,6 +213,7 @@ for (deal_II_dimension : DIMENSIONS) get_finest_common_cells (const MGDoFHandler &mesh_1, const MGDoFHandler &mesh_2); + #endif }