Merge lp:~jobh/dolfin/parallel-dirichlet into lp:~fenics-core/dolfin/trunk

[Joachim Haga]

This branch does the following to Dirichlet boundary conditions:

(a) Make DirichletBC pointwise search faster (by not doing unnecessary initialisation of facets, and only visiting vertices once).

(b) Make DirichletBC topological and geometric search safe in parallel for SystemAssembler (by adding a gather() method which exchanges boundary values with neighbour processes).

==

(a) is a simple change, by just delaying initialisation of facets until they are required. The timings for Poisson on a 48^3 cube are:

         before [s] after [s]

Topological search, compiled boundary subdomain
          20.3 19.9

Topological search, python boundary subdomain
          29.2 29.5

Pointwise search, compiled boundary subdomain
          21.4 19.9

Pointwise search, python boundary subdomain
          46.4 20.5 (!)

The recorded times include mesh create, 2x assemble, 2x DirichletBC create, 2x DirichletBC apply. The times do not change very much except for the case of pointwise search with a Python subdomain, which is dramatically improved. Pointwise search is now (much) faster than topological search when using Python subdomains.

(b) is more involved, and includes adding neighbour / dof sharing information to DofMap, and MPI functionality for exchanging data with neighbours (a new distribute() call, and a new send_recv because the current one deadlocks when posting different number of send_recvs on different processes).

It passes the unit and system tests (except KrylovSolver, but that one fails for me in trunk too).

[Garth Wells (garth-wells)]

This crashes the auto-adaptive-poisson demo, which is probably due to a bug in the adaptive solver (see below). I suspect that the solver fails to properly initialise the mesh, but the the initialisation was previously done by chance in the bc code.

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** https://answers.launchpad.net/dolfin
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to initialize mapping of degrees of freedom.
*** Reason: Missing entities of dimension 1. Try calling mesh.init(1).
*** Where: This error was encountered inside DofMap.cpp.
*** -------------------------------------------------------------------------

review: Needs Fixing

[Garth Wells (garth-wells)]

review: Approve

=== modified file 'dolfin/common/MPI.cpp'
--- dolfin/common/MPI.cpp	2011-11-18 12:14:50 +0000
+++ dolfin/common/MPI.cpp	2012-02-29 14:31:44 +0000
@@ -19,9 +19,10 @@
 // Modified by Anders Logg 2007-2011
 // Modified by Ola Skavhaug 2008-2009
 // Modified by Niclas Jansson 2009
+// Modified by Joachim B Haga 2012
 //
 // First added:  2007-11-30
-// Last changed: 2011-08-25
+// Last changed: 2012-02-29
 
 #include <numeric>
 #include <dolfin/log/dolfin_log.h>
@@ -145,6 +146,23 @@
   return r + (index - r * (n + 1)) / n;
 }
 //-----------------------------------------------------------------------------
+dolfin::MPINonblocking::MPINonblocking()
+{}
+//-----------------------------------------------------------------------------
+dolfin::MPINonblocking::~MPINonblocking()
+{
+  wait_all();
+}
+//-----------------------------------------------------------------------------
+void dolfin::MPINonblocking::wait_all()
+{
+  if (!reqs.empty())
+  {
+    boost::mpi::wait_all(reqs.begin(), reqs.end());
+    reqs.clear();
+  }
+}
+//-----------------------------------------------------------------------------
 //-----------------------------------------------------------------------------
 #else
 //-----------------------------------------------------------------------------
@@ -217,4 +235,15 @@
   return 0;
 }
 //-----------------------------------------------------------------------------
+dolfin::MPINonblocking::Nonblocking()
+{}
+//-----------------------------------------------------------------------------
+dolfin::MPINonblocking::~Nonblocking()
+{}
+//-----------------------------------------------------------------------------
+void dolfin::MPINonblocking::wait_all()
+{
+  error_no_mpi("call MPINonblocking::wait_all");
+}
+
 #endif
 
=== modified file 'dolfin/common/MPI.h'
--- dolfin/common/MPI.h	2011-12-04 11:33:42 +0000
+++ dolfin/common/MPI.h	2012-02-29 14:31:44 +0000
@@ -18,9 +18,10 @@
 // Modified by Ola Skavhaug 2008-2009
 // Modified by Anders Logg 2008-2011
 // Modified by Niclas Jansson 2009
+// Modified by Joachim B Haga 2012
 //
 // First added:  2007-11-30
-// Last changed: 2011-09-28
+// Last changed: 2012-02-29
 
 #ifndef __MPI_DOLFIN_WRAPPER_H
 #define __MPI_DOLFIN_WRAPPER_H
@@ -31,6 +32,7 @@
 #include <dolfin/log/dolfin_log.h>
 
 #ifdef HAS_MPI
+#include <boost/serialization/utility.hpp>
 #include <boost/mpi.hpp>
 #include <mpi.h>
 #endif
@@ -58,6 +60,32 @@
   };
   #endif
 
+  /// This class provides stateful (single communicator) non-blocking MPI functionality.
+
+  class MPINonblocking
+  {
+  public:
+    /// Create instance with associated MPI_WORLD communicator
+    MPINonblocking();
+
+    /// Destroy instance (waits for outstanding requests)
+    ~MPINonblocking();
+
+    /// Non-blocking send and receive
+    template<typename T>
+    void send_recv(const T& send_value, uint dest,
+                   T& recv_value, uint source);
+
+    /// Wait for all requests to finish
+    void wait_all();
+
+  #ifdef HAS_MPI
+  private:
+    MPICommunicator mpi_comm;
+    std::vector<boost::mpi::request> reqs;
+  #endif
+  };
+
   /// This class provides utility functions for easy communcation with MPI.
 
   class MPI
@@ -100,13 +128,22 @@
       distribute(in_values, destinations, out_values, sources);
     }
 
+    /// Distribute local arrays on a group of processes (typically neighbours
+    /// from GenericDofMap::neighbours()). It is important that each process'
+    /// group includes exactly the processes that has it in their groups,
+    /// otherwise it will deadlock.
+    template<typename T>
+    static void distribute(const std::set<uint> group,
+                           const std::map<uint,T>& in_values_per_dest,
+                           std::map<uint,T>& out_values_per_src);
+
     // Broadcast value from broadcaster process to all processes
     template<typename T>
     static void broadcast(T& value, uint broadcaster=0)
     {
       #ifdef HAS_MPI
       MPICommunicator mpi_comm;
-      boost::mpi::communicator comm(*mpi_comm, boost::mpi::comm_duplicate);
+      boost::mpi::communicator comm(*mpi_comm, boost::mpi::comm_attach);
       boost::mpi::broadcast(comm, value, broadcaster);
       #endif
     }
@@ -118,7 +155,7 @@
     {
       #ifdef HAS_MPI
       MPICommunicator mpi_comm;
-      boost::mpi::communicator comm(*mpi_comm, boost::mpi::comm_duplicate);
+      boost::mpi::communicator comm(*mpi_comm, boost::mpi::comm_attach);
       boost::mpi::scatter(comm, in_values, out_value, sending_process);
       #else
       dolfin_assert(sending_process == 0);
@@ -134,7 +171,7 @@
     {
       #ifdef HAS_MPI
       MPICommunicator mpi_comm;
-      boost::mpi::communicator comm(*mpi_comm, boost::mpi::comm_duplicate);
+      boost::mpi::communicator comm(*mpi_comm, boost::mpi::comm_attach);
       boost::mpi::gather(comm, in_value, out_values, receiving_process);
       #else
       out_values.clear();
@@ -148,7 +185,7 @@
     {
       #ifdef HAS_MPI
       MPICommunicator mpi_comm;
-      boost::mpi::communicator comm(*mpi_comm, boost::mpi::comm_duplicate);
+      boost::mpi::communicator comm(*mpi_comm, boost::mpi::comm_attach);
       boost::mpi::all_gather(comm, in_value, out_values);
       #else
       out_values.clear();
@@ -191,7 +228,7 @@
     {
       #ifdef HAS_MPI
       MPICommunicator mpi_comm;
-      boost::mpi::communicator comm(*mpi_comm, boost::mpi::comm_duplicate);
+      boost::mpi::communicator comm(*mpi_comm, boost::mpi::comm_attach);
       T out;
       boost::mpi::all_reduce(comm, value, out, op);
       return out;
@@ -207,20 +244,17 @@
     /// reduction op)
     static uint global_offset(uint range, bool exclusive);
 
-    /// Send-receive and data
+    /// Send-receive data. Note that if the number of posted send-receives may
+    /// differ between processes, another interface (such as MPINonblocking::send_recv)
+    /// must be used since duplicating the communicator requires participation
+    /// from all processes.
     template<typename T>
     static void send_recv(const T& send_value, uint dest,
                           T& recv_value, uint source)
     {
       #ifdef HAS_MPI
-      MPICommunicator mpi_comm;
-      boost::mpi::communicator comm(*mpi_comm, boost::mpi::comm_duplicate);
-
-      // Non-blocking send-receive
-      boost::mpi::request reqs[2];
-      reqs[0] = comm.isend(dest, 0, send_value);
-      reqs[1] = comm.irecv(source, 0, recv_value);
-      boost::mpi::wait_all(reqs, reqs + 2);
+      MPINonblocking mpi;
+      mpi.send_recv(send_value, dest, recv_value, source);
       #else
       dolfin_error("MPI.h",
                    "call MPI::send_recv",
@@ -246,6 +280,13 @@
 
   private:
 
+    #ifndef HAS_MPI
+    static void error_no_mpi(const char *where)
+    {
+      dolfin_error("MPI.h", msg, "DOLFIN has been configured without MPI support");
+    }
+    #endif
+
     #ifdef HAS_MPI
     // Return MPI data type
     template<typename T> static MPI_Datatype mpi_type()
@@ -259,6 +300,7 @@
 
   };
 
+
   #ifdef HAS_MPI
   // Specialisations for MPI_Datatypes
   template<> inline MPI_Datatype MPI::mpi_type<double>() { return MPI_DOUBLE; }
@@ -337,7 +379,63 @@
                  "DOLFIN has been configured without MPI support");
   }
   #endif
-  //---------------------------------------------------------------------------
+
+  //-----------------------------------------------------------------------------
+  template<typename T>
+  void dolfin::MPI::distribute(const std::set<uint> group,
+                               const std::map<uint,T>& in_values_per_dest,
+                               std::map<uint,T>& out_values_per_src)
+  {
+  #ifdef HAS_MPI
+    typedef typename std::map<uint,T>::const_iterator map_const_iterator;
+    typedef typename std::map<uint,T>::iterator map_iterator;
+    dolfin::MPINonblocking mpi;
+    const T no_data;
+
+    // Send and receive values to all processes in groups (non-blocking). If a
+    // given process is not found in in_values_per_dest, send empty data.
+
+    out_values_per_src.clear();
+    for (std::set<uint>::const_iterator dest = group.begin(); dest != group.end(); ++dest)
+    {
+      map_const_iterator values = in_values_per_dest.find(*dest);
+      if (values != in_values_per_dest.end())
+        mpi.send_recv(values->second, *dest, out_values_per_src[*dest], *dest);
+      else
+        mpi.send_recv(no_data, *dest, out_values_per_src[*dest], *dest);
+    }
+
+    // Wait for all MPI calls before modifying out_values_per_src
+
+    mpi.wait_all();
+
+    // Remove received no_data entries.
+
+    map_iterator it = out_values_per_src.begin();
+    while (it != out_values_per_src.end())
+    {
+      map_iterator tmp = it++;
+      if (tmp->second.empty())
+        out_values_per_src.erase(tmp); // map::erase only invalidates current iterator
+    }
+  #else
+    error_no_mpi("call MPI::distribute");
+  #endif
+  }
+
+  //-----------------------------------------------------------------------------
+  template<typename T>
+  void dolfin::MPINonblocking::send_recv(const T& send_value, uint dest,
+                                         T& recv_value, uint source)
+  {
+  #ifdef HAS_MPI
+    boost::mpi::communicator comm(*mpi_comm, boost::mpi::comm_attach);
+    reqs.push_back(comm.isend(dest, 0, send_value));
+    reqs.push_back(comm.irecv(source, 0, recv_value));
+  #else
+    error_no_mpi("call MPINonblocking::send_recv");
+  #endif
+  }
 
 }
 
 
=== modified file 'dolfin/fem/DirichletBC.cpp'
--- dolfin/fem/DirichletBC.cpp	2011-12-15 00:03:23 +0000
+++ dolfin/fem/DirichletBC.cpp	2012-02-29 14:31:44 +0000
@@ -18,13 +18,15 @@
 // Modified by Kristian Oelgaard, 2008
 // Modified by Martin Sandve Alnes, 2008
 // Modified by Johan Hake, 2009
+// Modified by Joachim B Haga, 2009
 //
 // First added:  2007-04-10
-// Last changed: 2011-09-19
+// Last changed: 2012-02-29
 
 #include <map>
 #include <utility>
 #include <boost/assign/list_of.hpp>
+#include <boost/serialization/utility.hpp>
 
 #include <dolfin/common/constants.h>
 #include <dolfin/common/Array.h>
@@ -67,7 +69,6 @@
 {
   check();
   parameters = default_parameters();
-  init_from_sub_domain(_user_sub_domain);
 }
 //-----------------------------------------------------------------------------
 DirichletBC::DirichletBC(boost::shared_ptr<const FunctionSpace> V,
@@ -80,7 +81,6 @@
 {
   check();
   parameters = default_parameters();
-  init_from_sub_domain(_user_sub_domain);
 }
 //-----------------------------------------------------------------------------
 DirichletBC::DirichletBC(const FunctionSpace& V, const GenericFunction& g,
@@ -89,11 +89,12 @@
   : BoundaryCondition(V),
     Hierarchical<DirichletBC>(*this),
     g(reference_to_no_delete_pointer(g)),
-    _method(method)
+    _method(method),
+    _user_mesh_function(reference_to_no_delete_pointer(sub_domains)),
+    _user_sub_domain_marker(sub_domain)
 {
   check();
   parameters = default_parameters();
-  init_from_mesh_function(sub_domains, sub_domain);
 }
 //-----------------------------------------------------------------------------
 DirichletBC::DirichletBC(boost::shared_ptr<const FunctionSpace> V,
@@ -103,22 +104,23 @@
                          std::string method)
   : BoundaryCondition(V),
     Hierarchical<DirichletBC>(*this),
-    g(g), _method(method)
+    g(g), _method(method),
+    _user_mesh_function(sub_domains),
+    _user_sub_domain_marker(sub_domain)
 {
   check();
   parameters = default_parameters();
-  init_from_mesh_function(*sub_domains, sub_domain);
 }
 //-----------------------------------------------------------------------------
 DirichletBC::DirichletBC(const FunctionSpace& V, const GenericFunction& g,
                          uint sub_domain, std::string method)
   : BoundaryCondition(V),
     Hierarchical<DirichletBC>(*this),
-    g(reference_to_no_delete_pointer(g)), _method(method)
+    g(reference_to_no_delete_pointer(g)), _method(method),
+    _user_sub_domain_marker(sub_domain)
 {
   check();
   parameters = default_parameters();
-  init_from_mesh(sub_domain);
 }
 //-----------------------------------------------------------------------------
 DirichletBC::DirichletBC(boost::shared_ptr<const FunctionSpace> V,
@@ -126,11 +128,11 @@
                          uint sub_domain, std::string method)
   : BoundaryCondition(V),
     Hierarchical<DirichletBC>(*this),
-    g(g), _method(method)
+    g(g), _method(method),
+    _user_sub_domain_marker(sub_domain)
 {
   check();
   parameters = default_parameters();
-  init_from_mesh(sub_domain);
 }
 //-----------------------------------------------------------------------------
 DirichletBC::DirichletBC(boost::shared_ptr<const FunctionSpace> V,
@@ -201,6 +203,43 @@
   apply(&A, &b, &x);
 }
 //-----------------------------------------------------------------------------
+void DirichletBC::gather(Map& boundary_values) const
+{
+  typedef std::vector<std::pair<uint, double> > bv_vec_type;
+  typedef std::map<uint, bv_vec_type> map_type;
+
+  typedef boost::unordered_map<uint, std::vector<uint> > shared_dof_type;
+  typedef shared_dof_type::const_iterator shared_dof_iterator;
+  typedef std::vector<uint>::const_iterator proc_iterator;
+
+  dolfin_assert(_function_space->dofmap());
+  const GenericDofMap& dofmap = *_function_space->dofmap();
+  const shared_dof_type& shared_dofs = dofmap.shared_dofs();
+
+  // Create list of boundary values to send to each processor
+
+  map_type proc_map;
+  for (Map::const_iterator bv = boundary_values.begin(); bv != boundary_values.end(); ++bv)
+  {
+    // If the boundary value is attached to a shared dof, add it to the list of
+    // boundary values for each of the processors that share it
+    shared_dof_iterator shared_dof = shared_dofs.find(bv->first);
+    if (shared_dof != shared_dofs.end())
+      for (proc_iterator proc = shared_dof->second.begin(); proc != shared_dof->second.end(); ++proc)
+        proc_map[*proc].push_back(*bv);
+  }
+
+  // Distribute the lists between neighbours
+
+  map_type received_bvs;
+  MPI::distribute(dofmap.neighbours(), proc_map, received_bvs);
+
+  // Add the received boundary values to the local boundary values
+
+  for (map_type::const_iterator it = received_bvs.begin(); it != received_bvs.end(); ++it)
+    boundary_values.insert(it->second.begin(), it->second.end());
+}
+//-----------------------------------------------------------------------------
 void DirichletBC::get_boundary_values(Map& boundary_values,
                          std::string method) const
 {
@@ -564,7 +603,20 @@
   }
 }
 //-----------------------------------------------------------------------------
-void DirichletBC::init_from_sub_domain(boost::shared_ptr<const SubDomain> sub_domain)
+void DirichletBC::init_facets() const
+{
+  if (facets.size() > 0)
+    return;
+
+  if (_user_sub_domain)
+    init_from_sub_domain(_user_sub_domain);
+  else if (_user_mesh_function)
+    init_from_mesh_function(*_user_mesh_function, _user_sub_domain_marker);
+  else
+    init_from_mesh(_user_sub_domain_marker);
+}
+//-----------------------------------------------------------------------------
+void DirichletBC::init_from_sub_domain(boost::shared_ptr<const SubDomain> sub_domain) const
 {
   dolfin_assert(facets.size() == 0);
 
@@ -595,7 +647,7 @@
 }
 //-----------------------------------------------------------------------------
 void DirichletBC::init_from_mesh_function(const MeshFunction<uint>& sub_domains,
-                                          uint sub_domain)
+                                          uint sub_domain) const
 {
   dolfin_assert(facets.size() == 0);
 
@@ -624,7 +676,7 @@
   }
 }
 //-----------------------------------------------------------------------------
-void DirichletBC::init_from_mesh(uint sub_domain)
+void DirichletBC::init_from_mesh(uint sub_domain) const
 {
   dolfin_assert(facets.size() == 0);
 
@@ -678,6 +730,8 @@
   dolfin_assert(_function_space);
   dolfin_assert(g);
 
+  init_facets();
+
   // Special case
   if (facets.size() == 0)
   {
@@ -724,18 +778,19 @@
       const double value = data.w[data.facet_dofs[i]];
       boundary_values[global_dof] = value;
     }
-
     p++;
   }
 }
 //-----------------------------------------------------------------------------
 void DirichletBC::compute_bc_geometric(Map& boundary_values,
-                                      BoundaryCondition::LocalData& data) const
+                                       BoundaryCondition::LocalData& data) const
 {
   dolfin_assert(_function_space);
   dolfin_assert(_function_space->element());
   dolfin_assert(g);
 
+  init_facets();
+
   // Special case
   if (facets.size() == 0)
   {
@@ -814,7 +869,11 @@
   dolfin_assert(_function_space);
   dolfin_assert(_function_space->element());
   dolfin_assert(g);
-  dolfin_assert(_user_sub_domain);
+
+  if (!_user_sub_domain)
+    dolfin_error("DirichletBC.cpp",
+                 "computing Dirichlet boundary values, pointwise search",
+                 "A SubDomain is required for pointwise search");
 
   // Get mesh and dofmap
   dolfin_assert(_function_space->mesh());
@@ -828,9 +887,13 @@
   // Create UFC cell object
   UFCCell ufc_cell(mesh);
 
+  // Speeder-upper
+  std::pair<uint,uint> local_range = dofmap.ownership_range();
+  std::vector<bool> already_visited(local_range.second - local_range.first);
+  std::fill(already_visited.begin(), already_visited.end(), false);
+
   // Iterate over cells
   Progress p("Computing Dirichlet boundary values, pointwise search", mesh.num_cells());
-  Array<double> x(gdim);
   for (CellIterator cell(mesh); !cell.end(); ++cell)
   {
     // Update UFC cell
@@ -848,9 +911,17 @@
     // Loop all dofs on cell
     for (uint i = 0; i < dofmap.cell_dimension(cell->index()); ++i)
     {
-      // Check if the coordinates are part of the sub domain (calls user-defined 'indside' function)
-      for (uint j = 0; j < gdim; ++j)
-        x[j] = data.coordinates[i][j];
+      const uint global_dof = cell_dofs[i];
+      if (global_dof >= local_range.first && global_dof < local_range.second)
+      {
+        const uint dof_index = global_dof - local_range.first;
+        if (already_visited[dof_index])
+          continue;
+        already_visited[dof_index] = true;
+      }
+
+      // Check if the coordinates are part of the sub domain (calls user-defined 'inside' function)
+      Array<double> x(gdim, &data.coordinates[i][0]);
       if (!_user_sub_domain->inside(x, false))
         continue;
 
@@ -863,7 +934,6 @@
       }
 
       // Set boundary value
-      const uint global_dof = cell_dofs[i];
       const double value = data.w[i];
       boundary_values[global_dof] = value;
     }
 
=== modified file 'dolfin/fem/DirichletBC.h'
--- dolfin/fem/DirichletBC.h	2011-10-26 08:59:01 +0000
+++ dolfin/fem/DirichletBC.h	2012-02-29 14:31:44 +0000
@@ -17,9 +17,10 @@
 //
 // Modified by Kristian Oelgaard, 2007
 // Modified by Johan Hake, 2009
+// Modified by Joachim B Haga, 2012
 //
 // First added:  2007-04-10
-// Last changed: 2011-04-13
+// Last changed: 2012-02-29
 //
 // FIXME: This class needs some cleanup, in particular collecting
 // FIXME: all data from different representations into a common
@@ -94,8 +95,6 @@
   /// three possibilties are "topological", "geometric" and
   /// "pointwise".
 
-  /// This class specifies the interface for setting (strong)
-
   class DirichletBC : public BoundaryCondition, public Hierarchical<DirichletBC>
   {
 
@@ -301,6 +300,16 @@
     void get_boundary_values(Map& boundary_values,
                              std::string method="default") const;
 
+
+    /// Get boundary values from neighbour processes. If a method other than
+    /// "pointwise" is used, this is necessary to ensure all boundary dofs are
+    /// marked on all processes.
+    ///
+    /// *Arguments*
+    ///     boundary_values (boost::unordered_map<uint, double>)
+    ///         Map from dof to boundary value.
+    void gather(Map& boundary_values) const;
+
     /// Make rows of matrix associated with boundary condition zero,
     /// useful for non-diagonal matrices in a block matrix.
     ///
@@ -401,15 +410,18 @@
     // Check input data to constructor
     void check() const;
 
+    // Initialize facets (from sub domain, mesh, etc)
+    void init_facets() const;
+
     // Initialize sub domain markers from sub domain
-    void init_from_sub_domain(boost::shared_ptr<const SubDomain> sub_domain);
+    void init_from_sub_domain(boost::shared_ptr<const SubDomain> sub_domain) const;
 
     // Initialize sub domain markers from MeshFunction
     void init_from_mesh_function(const MeshFunction<uint>& sub_domains,
-                                 uint sub_domain);
+                                 uint sub_domain) const;
 
     // Initialize sub domain markers from mesh
-    void init_from_mesh(uint sub_domain);
+    void init_from_mesh(uint sub_domain) const;
 
     // Compute dofs and values for application of boundary conditions using
     // given method
@@ -444,8 +456,13 @@
     boost::shared_ptr<const SubDomain> _user_sub_domain;
 
     // Boundary facets, stored as pairs (cell, local facet number)
-    std::vector<std::pair<uint, uint> > facets;
-
+    mutable std::vector<std::pair<uint, uint> > facets;
+
+    // User defined mesh function
+    boost::shared_ptr<const MeshFunction<uint> > _user_mesh_function;
+
+    // User defined sub domain marker for mesh or mesh function
+    uint _user_sub_domain_marker;
   };
 
 }
 
=== modified file 'dolfin/fem/DofMap.cpp'
--- dolfin/fem/DofMap.cpp	2011-11-18 12:14:50 +0000
+++ dolfin/fem/DofMap.cpp	2012-02-29 14:31:44 +0000
@@ -19,9 +19,10 @@
 // Modified by Kent-Andre Mardal, 2009
 // Modified by Ola Skavhaug, 2009
 // Modified by Niclas Jansson, 2009
+// Modified by Joachim B Haga, 2012
 //
 // First added:  2007-03-01
-// Last changed: 2011-10-31
+// Last changed: 2012-02-29
 
 #include <dolfin/common/MPI.h>
 #include <dolfin/common/NoDeleter.h>
@@ -173,6 +174,8 @@
   // Modify dofmap for non-UFC numbering
   ufc_map_to_dofmap.clear();
   _off_process_owner.clear();
+  _shared_dofs.clear();
+  _neighbours.clear();
   if (parent_dofmap.ufc_map_to_dofmap.size() > 0)
   {
     boost::unordered_map<uint, uint>::const_iterator ufc_to_current_dof;
@@ -196,6 +199,14 @@
         boost::unordered_map<uint, uint>::const_iterator parent_off_proc = parent_dofmap._off_process_owner.find(*dof);
         if (parent_off_proc != parent_dofmap._off_process_owner.end())
           _off_process_owner.insert(*parent_off_proc);
+
+        // Add to shared-dof process map, and update the set of neighbours
+        boost::unordered_map<uint, std::vector<uint> >::const_iterator parent_shared = parent_dofmap._shared_dofs.find(*dof);
+        if (parent_shared != parent_dofmap._shared_dofs.end())
+        {
+          _shared_dofs.insert(*parent_shared);
+          _neighbours.insert(parent_shared->second.begin(), parent_shared->second.end());
+        }
       }
     }
   }
@@ -309,6 +320,16 @@
   return _off_process_owner;
 }
 //-----------------------------------------------------------------------------
+const boost::unordered_map<uint, std::vector<uint> >& DofMap::shared_dofs() const
+{
+  return _shared_dofs;
+}
+//-----------------------------------------------------------------------------
+const std::set<uint>& DofMap::neighbours() const
+{
+  return _neighbours;
+}
+//-----------------------------------------------------------------------------
 void DofMap::tabulate_facet_dofs(uint* dofs, uint local_facet) const
 {
   dolfin_assert(_ufc_dofmap);
 
=== modified file 'dolfin/fem/DofMap.h'
--- dolfin/fem/DofMap.h	2011-11-18 12:14:50 +0000
+++ dolfin/fem/DofMap.h	2012-02-29 14:31:44 +0000
@@ -18,9 +18,10 @@
 // Modified by Martin Alnes, 2008
 // Modified by Kent-Andre Mardal, 2009
 // Modified by Ola Skavhaug, 2009
+// Modified by Joachim B Haga, 2012
 //
 // First added:  2007-03-01
-// Last changed: 2011-10-31
+// Last changed: 2012-02-29
 
 #ifndef __DOLFIN_DOF_MAP_H
 #define __DOLFIN_DOF_MAP_H
@@ -178,6 +179,21 @@
     ///         The map from non-local dofs.
     const boost::unordered_map<unsigned int, unsigned int>& off_process_owner() const;
 
+    /// Return map from all shared dofs to the processes (not including the current
+    /// process) that share it.
+    ///
+    /// *Returns*
+    ///     boost::unordered_map<uint, std::vector<uint> >
+    ///         The map from dofs to list of processes
+    const boost::unordered_map<uint, std::vector<uint> >& shared_dofs() const;
+
+    /// Return set of all neighbouring processes.
+    ///
+    /// *Returns*
+    ///     boost::unordered_set<uint>
+    ///         The set of processes
+    const std::set<uint>& neighbours() const;
+
     /// Local-to-global mapping of dofs on a cell
     ///
     /// *Arguments*
@@ -346,6 +362,12 @@
     // this process
     boost::unordered_map<uint, uint> _off_process_owner;
 
+    // List of processes that share a given dof
+    boost::unordered_map<uint, std::vector<uint> > _shared_dofs;
+
+    // Neighbours (processes that we share dofs with)
+    std::set<uint> _neighbours;
+
     // True iff sub-dofmap (a view, i.e. not collapsed)
     bool _is_view;
 
 
=== modified file 'dolfin/fem/DofMapBuilder.cpp'
--- dolfin/fem/DofMapBuilder.cpp	2011-11-18 12:14:50 +0000
+++ dolfin/fem/DofMapBuilder.cpp	2012-02-29 14:31:44 +0000
@@ -17,9 +17,10 @@
 //
 // Modified by Niclas Jansson 2009.
 // Modified by Garth N. Wells 2010.
+// Modified by Joachim B Haga, 2012.
 //
 // First added:  2008-08-12
-// Last changed: 2011-11-14
+// Last changed: 2012-02-29
 
 #include <ufc.h>
 #include <boost/random.hpp>
@@ -87,18 +88,20 @@
 {
   // Create data structures
   DofMapBuilder::set owned_dofs, shared_owned_dofs, shared_unowned_dofs;
+  DofMapBuilder::vec_map shared_dof_processes;
 
   // Computed owned and shared dofs (and owned and un-owned)
   compute_ownership(owned_dofs, shared_owned_dofs, shared_unowned_dofs,
-                    dofmap, global_dofs, mesh);
+                    shared_dof_processes, dofmap, global_dofs, mesh);
 
   // Renumber dofs owned dofs and received new numbering for unowned shared dofs
   parallel_renumber(owned_dofs, shared_owned_dofs, shared_unowned_dofs,
-                    dofmap, mesh);
+                    shared_dof_processes, dofmap, mesh);
 }
 //-----------------------------------------------------------------------------
 void DofMapBuilder::compute_ownership(set& owned_dofs, set& shared_owned_dofs,
                                       set& shared_unowned_dofs,
+                                      vec_map& shared_dof_processes,
                                       const DofMap& dofmap,
                                       const DofMapBuilder::set& global_dofs,
                                       const Mesh& mesh)
@@ -174,6 +177,7 @@
     {
       const uint received_dof  = recv_buffer[i];
       const uint received_vote = recv_buffer[i + 1];
+
       if (shared_owned_dofs.find(received_dof) != shared_owned_dofs.end())
       {
         // Move dofs with higher ownership votes from shared to shared but not owned
@@ -191,6 +195,14 @@
                        "compute mapping of degrees of freedom",
                        "Cannot decide on dof ownership; votes are equal");
         }
+
+        // Remember the sharing of the dof
+        shared_dof_processes[received_dof].push_back(src);
+      }
+      else if (shared_unowned_dofs.find(received_dof) != shared_unowned_dofs.end())
+      {
+        // Remember the sharing of the dof
+        shared_dof_processes[received_dof].push_back(src);
       }
     }
   }
@@ -240,6 +252,7 @@
 void DofMapBuilder::parallel_renumber(const set& owned_dofs,
                              const set& shared_owned_dofs,
                              const set& shared_unowned_dofs,
+                             const vec_map& shared_dof_processes,
                              DofMap& dofmap, const Mesh& mesh)
 {
   log(TRACE, "Renumber dofs for parallel dof map");
@@ -317,6 +330,17 @@
     }
   }
 
+  // Insert the shared-dof-to-process mapping into the dofmap, renumbering as necessary
+  for (vec_map::const_iterator it = shared_dof_processes.begin(); it != shared_dof_processes.end(); ++it)
+  {
+    boost::unordered_map<uint, uint>::const_iterator new_index = old_to_new_dof_index.find(it->first);
+    if (new_index == old_to_new_dof_index.end())
+      dofmap._shared_dofs.insert(*it);
+    else
+      dofmap._shared_dofs.insert(std::make_pair(new_index->second, it->second));
+    dofmap._neighbours.insert(it->second.begin(), it->second.end());
+  }
+
   // Build new dof map
   for (CellIterator cell(mesh); !cell.end(); ++cell)
   {
 
=== modified file 'dolfin/fem/DofMapBuilder.h'
--- dolfin/fem/DofMapBuilder.h	2011-11-24 17:51:44 +0000
+++ dolfin/fem/DofMapBuilder.h	2012-02-29 14:31:44 +0000
@@ -25,6 +25,7 @@
 
 #include <set>
 #include <boost/unordered_set.hpp>
+#include <boost/unordered_map.hpp>
 #include <dolfin/common/Set.h>
 
 namespace ufc
@@ -55,6 +56,8 @@
 
     typedef std::vector<dolfin::uint>::const_iterator vector_it;
 
+    typedef boost::unordered_map<uint, std::vector<uint> > vec_map;
+
   public:
 
     static void build(DofMap& dofmap, const Mesh& dolfin_mesh,
@@ -70,6 +73,7 @@
 
     static void compute_ownership(set& owned_dofs, set& shared_owned_dofs,
                                   set& shared_unowned_dofs,
+                                  vec_map& shared_dof_processes,
                                   const DofMap& dofmap,
                                   const DofMapBuilder::set& global_dofs,
                                   const Mesh& mesh);
@@ -77,6 +81,7 @@
     static void parallel_renumber(const set& owned_dofs,
                                   const set& shared_owned_dofs,
                                   const set& shared_unowned_dofs,
+                                  const vec_map& shared_dof_processes,
                                   DofMap& dofmap, const Mesh& mesh);
 
     /// Compute set of global dofs (e.g. Reals associated with global
 
=== modified file 'dolfin/fem/GenericDofMap.h'
--- dolfin/fem/GenericDofMap.h	2011-09-27 07:14:59 +0000
+++ dolfin/fem/GenericDofMap.h	2012-02-29 14:31:44 +0000
@@ -15,8 +15,10 @@
 // You should have received a copy of the GNU Lesser General Public License
 // along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
 //
+// Modified by Joachim B Haga, 2012
+//
 // First added:  2010-05-26
-// Last changed: 2011-02-23
+// Last changed: 2012-02-29
 
 #ifndef __GENERIC_DOF_MAP_H
 #define __GENERIC_DOF_MAP_H
@@ -108,6 +110,13 @@
     /// Return the set of dof indices
     virtual boost::unordered_set<uint> dofs() const = 0;
 
+    /// Return map from shared dofs to the processes (not including the current
+    /// process) that share it.
+    virtual const boost::unordered_map<uint, std::vector<uint> >& shared_dofs() const = 0;
+
+    /// Return set of all processes that share dofs with the current process.
+    virtual const std::set<uint>& neighbours() const = 0;
+
     /// Re-number based on provided re-numbering map
     virtual void renumber(const std::vector<uint>& renumbering_map) = 0;
 
 
=== modified file 'dolfin/fem/SystemAssembler.cpp'
--- dolfin/fem/SystemAssembler.cpp	2012-02-22 09:40:00 +0000
+++ dolfin/fem/SystemAssembler.cpp	2012-02-29 14:31:44 +0000
@@ -16,9 +16,10 @@
 // along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
 //
 // Modified by Anders Logg 2008-2011
+// Modified by Joachim B Haga 2012
 //
 // First added:  2009-06-22
-// Last changed: 2011-11-14
+// Last changed: 2012-02-29
 
 #include <dolfin/log/dolfin_log.h>
 #include <dolfin/common/Timer.h>
@@ -168,20 +169,9 @@
   DirichletBC::Map boundary_values;
   for (uint i = 0; i < bcs.size(); ++i)
   {
-    // Methods other than 'pointwise' are not robust in parallel since a vertex
-    // can have a bc applied, but the partition might not have a facet on the boundary.
+    bcs[i]->get_boundary_values(boundary_values);
     if (MPI::num_processes() > 1 && bcs[i]->method() != "pointwise")
-    {
-      if (MPI::process_number() == 0)
-      {
-        warning("Dirichlet boundary condition method '%s' is not robust in parallel with symmetric assembly.", bcs[i]->method().c_str());
-        //warning("Caution: 'on_boundary' does not work with 'pointwise' boundary conditions,");
-      }
-      bcs[i]->get_boundary_values(boundary_values);
-      //bcs[i]->get_boundary_values(boundary_values, "pointwise");
-    }
-    else
-      bcs[i]->get_boundary_values(boundary_values);
+      bcs[i]->gather(boundary_values);
   }
 
   // Modify boundary values for incremental (typically nonlinear) problems