diff src/Function.cc @ 256:8fe68d94ab76

Add parallel assembly of matrices and vectors
author Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
date Wed, 30 Jul 2014 20:59:37 +0200
parents 5e9b5bbdc56b
children 61830a4f9ab9
line wrap: on
line diff
--- a/src/Function.cc	Wed Jul 30 18:09:13 2014 +0200
+++ b/src/Function.cc	Wed Jul 30 20:59:37 2014 +0200
@@ -19,6 +19,13 @@
 #include "functionspace.h"
 #include "function.h"
 #include "dolfin_compat.h"
+#ifdef LATEST_DOLFIN
+#include <boost/mpi.hpp>
+#include <boost/serialization/vector.hpp>
+#include <boost/serialization/utility.hpp>
+#endif
+
+void copy_vector (Array <double> const&, SHARED_PTR <dolfin::GenericVector>);
 
 DEFUN_DLD (Function, args, , "-*- texinfo -*-\n\
 @deftypefn {Function File} {[@var{func}]} = \
@@ -90,20 +97,11 @@
                        and of the vector must agree");
               else
                 {
-#ifdef LATEST_DOLFIN
-                  dolfin::Vector du (MPI_COMM_WORLD, myu.length ());
-#else
-                  dolfin::Vector du(myu.length ());
-#endif
-                  for (std::size_t i = 0; i < myu.length (); ++i)
-                    du.setitem (i, myu(i));
+                  dolfin::Function aux (V);
+                  copy_vector (myu, aux.vector ());
 
-                  SHARED_PTR <dolfin::Vector> 
-                    uu (new dolfin::Vector(du));
-
-                  SHARED_PTR <const dolfin::Function>
-                    u (new dolfin::Function(V, uu));
-
+                  SHARED_PTR <dolfin::Function const>
+                    u (new dolfin::Function (aux));
                   retval = new function (str, u);
                }
             }
@@ -140,3 +138,64 @@
     }
   return retval;
 }
+
+void
+copy_vector (Array <double> const& arr, SHARED_PTR <dolfin::GenericVector> vec)
+{
+  unsigned const size =
+#ifdef LATEST_DOLFIN
+    dolfin::MPI::size (MPI_COMM_WORLD);
+#else
+    dolfin::MPI::num_processes ();
+#endif
+
+  if (size > 1)
+    {
+      std::vector <double> locvals;
+      boost::mpi::communicator world;
+      unsigned const rank = world.rank ();
+
+      std::pair <std::size_t, std::size_t> const loc_range =
+        vec->local_range ();
+
+      for (unsigned p = 1; p < size; ++p)
+        if (rank == p)
+          world.send (0, p, loc_range);
+
+      if (rank == 0)
+        {
+          for (unsigned p = 0; p < size; ++p)
+            {
+              std::pair <std::size_t, std::size_t> range;
+              if (p == 0)
+                range = loc_range;
+              else
+                world.recv (p, p, range);
+
+              std::vector <double> entries;
+              for (std::size_t i = range.first; i < range.second; ++i)
+                entries.push_back (arr(i));
+
+              if (p == 0)
+                locvals = entries;
+              else
+                world.send (p, p, entries);
+            }
+        }
+
+      for (unsigned p = 1; p < size; ++p)
+        {
+          if (rank == p)
+            world.recv (0, p, locvals);
+        }
+
+      vec->set_local (locvals);
+    }
+  else
+    for (std::size_t i = 0; i < arr.length (); ++i)
+      vec->setitem (i, arr(i));
+
+  vec->apply ("insert");
+
+  return;
+}