]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Added the function GridTools::extract_boundary_mesh (...).
authorpauletti <pauletti@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 15 Nov 2010 17:09:16 +0000 (17:09 +0000)
committerpauletti <pauletti@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 15 Nov 2010 17:09:16 +0000 (17:09 +0000)
git-svn-id: https://svn.dealii.org/trunk@22747 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/grid/grid_tools.h
deal.II/source/grid/grid_tools.cc
deal.II/source/grid/grid_tools.inst.in

index 606709fc684aa86344a13c40e0f953295bfd06f6..15237f947587e2a6d91ccbb91bf6c7651176a079 100644 (file)
@@ -795,6 +795,28 @@ class GridTools
     fix_up_distorted_child_cells (const typename Triangulation<dim,spacedim>::DistortedCellList &distorted_cells,
                                  Triangulation<dim,spacedim> &triangulation);
 
+  /** 
+      This function implements a boundary subgrid extraction.
+      Given a <dim,spacedim>-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 <int dim, int spacedim>
+  static void 
+  extract_boundary_mesh (const Triangulation<dim,spacedim> &volume_mesh,
+                        Triangulation<dim-1,spacedim>     &surface_mesh,
+                        std::map<typename Triangulation<dim-1,spacedim>::cell_iterator,
+                                 typename Triangulation<dim,spacedim>::face_iterator> 
+                        &surface_to_volume_mapping,
+                        const std::set<unsigned char> &boundary_ids
+                        = std::set<unsigned char>());
+  
                                      /**
                                       * 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<Point<2> > &all_vertices,
 //                     const unsigned int (&vertex_indices) [GeometryInfo<2>::vertices_per_cell]);
 
 
+
+
+
+
+
 DEAL_II_NAMESPACE_CLOSE
 
 /*----------------------------   grid_tools.h     ---------------------------*/
index 1cd240ce10aeac70b0732f57f2081074fef1dbfc..862fe7a2526ed4e664b06a7134e0f0813b57d35b 100644 (file)
@@ -1812,6 +1812,138 @@ fix_up_distorted_child_cells (const typename Triangulation<dim,spacedim>::Distor
 
 
 
+template <int dim, int spacedim>
+void 
+GridTools::extract_boundary_mesh (const Triangulation<dim,spacedim> &volume_mesh,
+                             Triangulation<dim-1,spacedim>     &surface_mesh,
+                             std::map<typename Triangulation<dim-1,spacedim>::cell_iterator,
+                                      typename Triangulation<dim,spacedim>::face_iterator> 
+                               &surface_to_volume_mapping,
+                             const std::set<unsigned char> &boundary_ids)
+{
+
+  // <h3>Assumptions</h3>
+  //    We are relying heavily on that 
+  //    Triangulation::create_triangulation(...) will keep the order 
+  //    pass by CellData and that it will not reorder the vertices.
+  
+  //    <h3>TODO</h3>
+  //    <ul>
+  //    <li> std::set<unsigned char> &boundary_ids needs to be set 
+  //         to a default value of {0}, no is set to empty
+  //    <li> Add mapping from volume to surface
+  //    <li> Add some assertions
+  //    <li> Create a function (reconcile?) to reconcile the triangulations
+  //         after refinement/coarsening is performed in one of them
+  //    <li> Create a class that manages the triangulations, as direct
+  //         access to refining the triangulations will most likely lead
+  //   to inconsitencies
+  //    </ul>
+
+  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<typename Triangulation<dim,spacedim>::face_iterator> 
+    mapping;  // temporary map for level==0 
+
+
+
+  typename Triangulation<dim,spacedim>::cell_iterator
+    cell = volume_mesh.begin(0),
+    endc = volume_mesh.end(0);
+  typename Triangulation<dim,spacedim>::face_iterator face;
+
+  std::vector< bool >                 touched (volume_mesh.n_vertices(), false);
+  std::vector< CellData< boundary_dim > > cells;
+  std::vector< Point<spacedim> >      vertices;
+
+  std::map<unsigned int,unsigned int> 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<dim>::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<GeometryInfo<boundary_dim>::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<dim-1,spacedim>::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<dim-1,spacedim>::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<dim-1,spacedim>::cell_iterator
+      cell = surface_mesh.begin(),
+      endc = surface_mesh.end();
+      for (; cell!=endc; ++cell)
+       for (unsigned int c=0; c<cell->n_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"
 
index 0a1cc0795d3c8841cbb769843c02604c88c1249a..4464e78ed5bd9098009d9168892ef953fbe674eb 100644 (file)
@@ -57,6 +57,15 @@ for (deal_II_dimension : DIMENSIONS)
     template
       double
       GridTools::diameter<deal_II_dimension> (const Triangulation<deal_II_dimension> &);
+    template 
+      void
+      GridTools::extract_boundary_mesh (const Triangulation<deal_II_dimension> &volume_mesh,
+                                       Triangulation<deal_II_dimension-1,deal_II_dimension>  &surface_mesh,
+                                       std::map<  Triangulation<deal_II_dimension-1,deal_II_dimension>::cell_iterator,
+                                                  Triangulation<deal_II_dimension>::face_iterator> 
+      &surface_to_volume_mapping,
+                                       const std::set<unsigned char> &boundary_ids
+                                       = std::set<unsigned char>());
 #endif
 
 #if deal_II_dimension == 2
@@ -204,6 +213,7 @@ for (deal_II_dimension : DIMENSIONS)
       get_finest_common_cells (const MGDoFHandler<deal_II_dimension,deal_II_dimension+1> &mesh_1,
                               const MGDoFHandler<deal_II_dimension,deal_II_dimension+1> &mesh_2);
 
+    
 #endif
   }
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.