From daeed038f37bc442bd403e2b31c9f2428a17e8cd Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Mon, 30 Nov 2015 17:21:06 -0600
Subject: [PATCH] Also copy edge boundary ids to manifold ids in
 copy_boundary_to_manifold_id().

---
 doc/news/changes.h                |  7 ++++++
 include/deal.II/grid/grid_tools.h | 14 +++++++-----
 source/grid/grid_tools.cc         | 38 ++++++++++++++++++++++++++-----
 3 files changed, 47 insertions(+), 12 deletions(-)

diff --git a/doc/news/changes.h b/doc/news/changes.h
index a6215cb618..a33ee863d3 100644
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -531,6 +531,13 @@ inconvenience this causes.
   (Martin Kronbichler, 2015/12/05)
   </li>
 
+  <li> Fixed: The GridTools::copy_boundary_to_manifold_id() function
+  only copied boundary indicators from faces, but in 3d forgot about
+  edges. This is now fixed.
+  <br>
+  (Wolfgang Bangerth, 2015/11/30)
+  </li>
+
   <li> Fixed: The constructor of SymmetricTensor that takes an array
   of initializing elements led to a compiler error. This is now
   fixed.
diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h
index c8909f8d22..69c7c2fc15 100644
--- a/include/deal.II/grid/grid_tools.h
+++ b/include/deal.II/grid/grid_tools.h
@@ -1253,17 +1253,19 @@ namespace GridTools
   /*@{*/
 
   /**
-   * Copy boundary ids to manifold ids. The default manifold_id for new
-   * Triangulation objects is numbers::invalid_manifold_id. This function
-   * copies the boundary_ids of the boundary faces to the manifold_ids of the
-   * same faces, allowing the user to change the boundary_ids and use them for
+   * Copy boundary ids to manifold ids on faces and edges at the boundary.
+   * The default manifold_id for new Triangulation objects is
+   * numbers::invalid_manifold_id. This function copies the boundary_ids of
+   * the boundary faces and edges to the manifold_ids of the same faces
+   * and edges, allowing the user to change the boundary_ids and use them for
    * boundary conditions regardless of the geometry, which will use
    * manifold_ids to create new points. Only active cells will be iterated
    * over. This is a function you'd typically call when there is only one
-   * active level on your Triangulation.
+   * active level on your Triangulation. Mesh refinement will then inherit
+   * these indicators to child cells, faces, and edges.
    *
    * The optional parameter @p reset_boundary_ids, indicates whether this
-   * function should reset the boundary_ids of the Triangulation to its
+   * function should reset the boundary_ids of boundary faces and edges to its
    * default value 0 after copying its value to the manifold_id. By default,
    * boundary_ids are left untouched.
    *
diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index 8f569d8588..a9481555f8 100644
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -3252,25 +3252,51 @@ next_cell:
 #endif
   }
 
+
+
   template <int dim, int spacedim>
   void copy_boundary_to_manifold_id(Triangulation<dim, spacedim> &tria,
                                     const bool reset_boundary_ids)
   {
+    // in 3d, we not only have to copy boundary ids of faces, but also of edges
+    // because we see them twice (once from each adjacent boundary face),
+    // we cannot immediately reset their boundary ids. thus, copy first
+    // and reset later
+    if (dim >= 3)
+      for (typename Triangulation<dim,spacedim>::active_cell_iterator
+           cell=tria.begin_active();
+           cell != tria.end(); ++cell)
+        for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+          if (cell->face(f)->at_boundary())
+            for (unsigned int e=0; e<GeometryInfo<dim>::lines_per_face; ++e)
+              cell->face(f)->line(e)->set_manifold_id
+              (static_cast<types::manifold_id>(cell->face(f)->line(e)->boundary_id()));
 
-    typename Triangulation<dim,spacedim>::active_cell_iterator
-    cell=tria.begin_active(), endc=tria.end();
-
-    for (; cell != endc; ++cell)
+    // now do cells
+    for (typename Triangulation<dim,spacedim>::active_cell_iterator
+         cell=tria.begin_active();
+         cell != tria.end(); ++cell)
       for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
         if (cell->face(f)->at_boundary())
           {
+            // copy boundary to manifold ids
             cell->face(f)->set_manifold_id
             (static_cast<types::manifold_id>(cell->face(f)->boundary_id()));
+
+            // then reset boundary ids if so desired, and in 3d also that
+            // of edges
             if (reset_boundary_ids == true)
-              cell->face(f)->set_boundary_id(0);
+              {
+                cell->face(f)->set_boundary_id(0);
+                if (dim >= 3)
+                  for (unsigned int e=0; e<GeometryInfo<dim>::lines_per_face; ++e)
+                    cell->face(f)->line(e)->set_boundary_id(0);
+              }
           }
   }
 
+
+
   template <int dim, int spacedim>
   void copy_material_to_manifold_id(Triangulation<dim, spacedim> &tria,
                                     const bool compute_face_ids)
@@ -3285,7 +3311,7 @@ next_cell:
           {
             for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
               {
-                if (cell->neighbor(f) != endc)
+                if (cell->at_boundary(f) == false)
                   cell->face(f)->set_manifold_id
                   (std::min(cell->material_id(),
                             cell->neighbor(f)->material_id()));
-- 
2.39.5