Generalize the extract_boundary_mesh function to arbitrary container types.

git-svn-id: https://svn.dealii.org/trunk@22801 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 ea4db5ce47ff83f8901664d31c5afa92985ec93d..459d10f5cf25e30329fd2f87e4f37e78a1800351 100644 (file)
--- a/deal.II/include/deal.II/grid/grid_tools.h
+++ b/deal.II/include/deal.II/grid/grid_tools.h
@@ -795,7 +795,7 @@ 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
@@ -805,7 +805,7 @@ class GridTools
                                      * 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.
@@ -832,17 +832,42 @@ class GridTools
                                      * rather than any curved boundary object
                                      * you may want to use to determine the
                                      * location of new vertices.
+                                     *
+                                     * @note Oftentimes, the
+                                     * <code>Container</code>
+                                     * template type will be of kind
+                                     * Triangulation; in that case,
+                                     * the map that is returned will
+                                     * be between Triangulation cell
+                                     * iterators of the surface mesh
+                                     * and Triangulation face
+                                     * iterators of the volume
+                                     * mesh. However, one often needs
+                                     * to have this mapping between
+                                     * DoFHandler (or hp::DoFHandler)
+                                     * iterators. In that case, you
+                                     * can pass DoFHandler arguments
+                                     * as first and second parameter;
+                                     * the function will in that case
+                                     * re-build the triangulation
+                                     * underlying the second argument
+                                     * and return a map between
+                                     * DoFHandler iterators. However,
+                                     * the function will not actually
+                                     * distribute degrees of freedom
+                                     * on this newly created surface
+                                     * mesh.
                                      */
-    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> 
+    template <template <int,int> class Container, int dim, int spacedim>
+    static void
+    extract_boundary_mesh (const Container<dim,spacedim> &volume_mesh,
+                          Container<dim-1,spacedim>     &surface_mesh,
+                          std::map<typename Container<dim-1,spacedim>::cell_iterator,
+                          typename Container<dim,spacedim>::face_iterator>
     &surface_to_volume_mapping,
                           const std::set<unsigned char> &boundary_ids
                           = std::set<unsigned char>());
-  
+
                                      /**
                                       * Exception
                                       */

diff --git a/deal.II/source/grid/grid_tools.cc b/deal.II/source/grid/grid_tools.cc
index 47d29ea58e271301a21cfa420ca38ae8aae421a8..646769dabbc7a4b08515bdc24d838738452eb568 100644 (file)
--- a/deal.II/source/grid/grid_tools.cc
+++ b/deal.II/source/grid/grid_tools.cc
@@ -43,17 +43,34 @@ DEAL_II_NAMESPACE_OPEN
 // and the like
 namespace
 {
-    template<int dim, int spacedim>
-    const Triangulation<dim, spacedim> &get_tria(const Triangulation<dim, spacedim> &tria)
-   {
-      return tria;
-   }
-
-    template<int dim, template<int, int> class Container, int spacedim>
-    const Triangulation<dim> &get_tria(const Container<dim,spacedim> &container)
-   {
-      return container.get_tria();
-   }
+  template<int dim, int spacedim>
+  const Triangulation<dim, spacedim> &
+  get_tria(const Triangulation<dim, spacedim> &tria)
+  {
+    return tria;
+  }
+
+  template<int dim, template<int, int> class Container, int spacedim>
+  const Triangulation<dim,spacedim> &
+  get_tria(const Container<dim,spacedim> &container)
+  {
+    return container.get_tria();
+  }
+
+
+  template<int dim, int spacedim>
+  Triangulation<dim, spacedim> &
+  get_tria(Triangulation<dim, spacedim> &tria)
+  {
+    return tria;
+  }
+
+  template<int dim, template<int, int> class Container, int spacedim>
+  const Triangulation<dim,spacedim> &
+  get_tria(Container<dim,spacedim> &container)
+  {
+    return container.get_tria();
+  }
 }
 
 
@@ -1812,14 +1829,14 @@ fix_up_distorted_child_cells (const typename Triangulation<dim,spacedim>::Distor
 
 
 
-template <int dim, int spacedim>
+template <template <int,int> class Container, 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)
+GridTools::extract_boundary_mesh (const Container<dim,spacedim> &volume_mesh,
+                                 Container<dim-1,spacedim>     &surface_mesh,
+                                 std::map<typename Container<dim-1,spacedim>::cell_iterator,
+                                 typename Container<dim,spacedim>::face_iterator>
+                                   &surface_to_volume_mapping,
+                                 const std::set<unsigned char> &boundary_ids)
 {
 // Assumption:
 //    We are relying below on the fact that Triangulation::create_triangulation(...) will keep the order
@@ -1831,11 +1848,11 @@ GridTools::extract_boundary_mesh (const Triangulation<dim,spacedim> &volume_mesh
 
                                   // First create surface mesh and mapping from only level(0) cells of volume_mesh
 
-  std::vector<typename Triangulation<dim,spacedim>::face_iterator>
+  std::vector<typename Container<dim,spacedim>::face_iterator>
     mapping;  // temporary map for level==0
 
 
-  std::vector< bool >                 touched (volume_mesh.n_vertices(), false);
+  std::vector< bool > touched (get_tria(volume_mesh).n_vertices(), false);
   std::vector< CellData< boundary_dim > > cells;
   std::vector< Point<spacedim> >      vertices;
 
@@ -1844,13 +1861,13 @@ GridTools::extract_boundary_mesh (const Triangulation<dim,spacedim> &volume_mesh
   unsigned int v_index;
   CellData< boundary_dim > c_data;
 
-  for (typename Triangulation<dim,spacedim>::cell_iterator
+  for (typename Container<dim,spacedim>::cell_iterator
         cell = volume_mesh.begin(0);
        cell != volume_mesh.end(0);
        ++cell)
     for (unsigned int i=0; i < GeometryInfo<dim>::faces_per_cell; ++i)
       {
-       const typename Triangulation<dim,spacedim>::face_iterator
+       const typename Container<dim,spacedim>::face_iterator
          face = cell->face(i);
 
        if ( face->at_boundary()
@@ -1880,10 +1897,11 @@ GridTools::extract_boundary_mesh (const Triangulation<dim,spacedim> &volume_mesh
       }
 
                                   // create level 0 surface triangulation
-  surface_mesh.create_triangulation (vertices, cells, SubCellData());
+  const_cast<Triangulation<dim-1,spacedim>&>(get_tria(surface_mesh))
+    .create_triangulation (vertices, cells, SubCellData());
 
                                   // Make the actual mapping
-  for (typename Triangulation<dim-1,spacedim>::active_cell_iterator
+  for (typename Container<dim-1,spacedim>::active_cell_iterator
         cell = surface_mesh.begin(0);
        cell!=surface_mesh.end(0); ++cell)
     surface_to_volume_mapping[cell] = mapping.at(cell->index());
@@ -1891,7 +1909,7 @@ GridTools::extract_boundary_mesh (const Triangulation<dim,spacedim> &volume_mesh
   do
     {
       bool changed = false;
-      typename Triangulation<dim-1,spacedim>::active_cell_iterator
+      typename Container<dim-1,spacedim>::active_cell_iterator
        cell = surface_mesh.begin_active(),
        endc = surface_mesh.end();
 
@@ -1904,9 +1922,10 @@ GridTools::extract_boundary_mesh (const Triangulation<dim,spacedim> &volume_mesh
 
       if (changed)
        {
-         surface_mesh.execute_coarsening_and_refinement ();
+         const_cast<Triangulation<dim-1,spacedim>&>(get_tria(surface_mesh))
+           .execute_coarsening_and_refinement();
 
-         typename Triangulation<dim-1,spacedim>::cell_iterator
+         typename Container<dim-1,spacedim>::cell_iterator
            cell = surface_mesh.begin(),
            endc = surface_mesh.end();
          for (; cell!=endc; ++cell)

diff --git a/deal.II/source/grid/grid_tools.inst.in b/deal.II/source/grid/grid_tools.inst.in
index 0c5da4e24046d61fb48001f7269299f10a700820..e4b4f392d98202e6fc1a044e5fd9ebf0b70c0364 100644 (file)
--- a/deal.II/source/grid/grid_tools.inst.in
+++ b/deal.II/source/grid/grid_tools.inst.in
@@ -57,14 +57,6 @@ 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);
 #endif
 
 #if deal_II_dimension == 2
@@ -212,7 +204,36 @@ 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
+  }
+
+// TODO: Merge this and the next block by introducing a TRIA_AND_DOFHANDLER_TEMPLATES list.
+for (deal_II_dimension : DIMENSIONS)
+  {
+#if deal_II_dimension != 1
+    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);
 #endif
   }
 
+for (deal_II_dimension : DIMENSIONS; Container : DOFHANDLER_TEMPLATES)
+  {
+
+#if deal_II_dimension != 1
+    template
+      void
+      GridTools::extract_boundary_mesh (const Container<deal_II_dimension> &volume_mesh,
+                                       Container<deal_II_dimension-1,deal_II_dimension>  &surface_mesh,
+                                       std::map<  Container<deal_II_dimension-1,deal_II_dimension>::cell_iterator,
+                                                  Container<deal_II_dimension>::face_iterator>
+      &surface_to_volume_mapping,
+                                       const std::set<unsigned char> &boundary_ids);
+#endif
+  }
\ No newline at end of file