From 0b8acf158186abfeefcc21bca1f26c75912e203f Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Mon, 11 Jan 2016 15:49:59 -0600
Subject: [PATCH] Make sure DataOutFaces also works in parallel.

---
 doc/news/changes.h                |  6 ++++
 source/numerics/data_out_faces.cc | 47 +++++++++++++++++--------------
 2 files changed, 32 insertions(+), 21 deletions(-)

diff --git a/doc/news/changes.h b/doc/news/changes.h
index ab1f6ea546..73b5e5bfa2 100644
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -572,6 +572,12 @@ inconvenience this causes.
 
 
 <ol>
+  <li> Fixed: The DataOutFaces class should now also work with triangulations
+  of type parallel::distributed::Triangulation.
+  <br>
+  (Heikki Virtanen, Wolfgang Bangerth, 2016/01/11)
+  </li>
+
   <li> Fixed: AlignedVector<T>::fill() (and thus, Table<N,T>::reinit) did not
   correctly call the destructor of T() and could leak memory for complicated
   class types that depend on their constructor to free memory.
diff --git a/source/numerics/data_out_faces.cc b/source/numerics/data_out_faces.cc
index 170656c2d9..9a77385892 100644
--- a/source/numerics/data_out_faces.cc
+++ b/source/numerics/data_out_faces.cc
@@ -355,18 +355,17 @@ template <int dim, typename DoFHandlerType>
 typename DataOutFaces<dim,DoFHandlerType>::FaceDescriptor
 DataOutFaces<dim,DoFHandlerType>::first_face ()
 {
-  // simply find first active cell
-  // with a face on the boundary
+  // simply find first active cell with a face on the boundary
   typename Triangulation<dimension,space_dimension>::active_cell_iterator cell = this->triangulation->begin_active();
   for (; cell != this->triangulation->end(); ++cell)
-    for (unsigned int f=0; f<GeometryInfo<dimension>::faces_per_cell; ++f)
-      if (!surface_only || cell->face(f)->at_boundary())
-        return FaceDescriptor(cell, f);
-
-  // ups, triangulation has no
-  // boundary? impossible!
-  Assert (false, ExcInternalError());
+    if (cell->is_locally_owned())
+      for (unsigned int f=0; f<GeometryInfo<dimension>::faces_per_cell; ++f)
+        if (!surface_only || cell->face(f)->at_boundary())
+          return FaceDescriptor(cell, f);
 
+  // just return an invalid descriptor if we haven't found a locally
+  // owned face. this can happen in parallel where all boundary
+  // faces are owned by other processors
   return FaceDescriptor();
 }
 
@@ -378,14 +377,16 @@ DataOutFaces<dim,DoFHandlerType>::next_face (const FaceDescriptor &old_face)
 {
   FaceDescriptor face = old_face;
 
-  // first check whether the present cell has more faces on the boundary
+  // first check whether the present cell has more faces on the boundary. since
+  // we started with this face, its cell must clearly be locally owned
+  Assert (face.first->is_locally_owned(), ExcInternalError());
   for (unsigned int f=face.second+1; f<GeometryInfo<dimension>::faces_per_cell; ++f)
     if (!surface_only || face.first->face(f)->at_boundary())
       // yup, that is so, so return it
       {
         face.second = f;
         return face;
-      };
+      }
 
   // otherwise find the next active cell that has a face on the boundary
 
@@ -399,17 +400,21 @@ DataOutFaces<dim,DoFHandlerType>::next_face (const FaceDescriptor &old_face)
   // while there are active cells
   while (active_cell != this->triangulation->end())
     {
-      // check all the faces of this active cell
-      for (unsigned int f=0; f<GeometryInfo<dimension>::faces_per_cell; ++f)
-        if (!surface_only || active_cell->face(f)->at_boundary())
-          {
-            face.first  = active_cell;
-            face.second = f;
-            return face;
-          };
-      // the present cell had no faces on the boundary, so check next cell
+      // check all the faces of this active cell. but skip it altogether
+      // if it isn't locally owned
+      if (active_cell->is_locally_owned())
+        for (unsigned int f=0; f<GeometryInfo<dimension>::faces_per_cell; ++f)
+          if (!surface_only || active_cell->face(f)->at_boundary())
+            {
+              face.first  = active_cell;
+              face.second = f;
+              return face;
+            }
+
+      // the present cell had no faces on the boundary (or was not locally
+      // owned), so check next cell
       ++active_cell;
-    };
+    }
 
   // we fell off the edge, so return with invalid pointer
   face.first  = this->triangulation->end();
-- 
2.39.5