Index: ThirdParty/vmg/src/comm/comm.cpp
===================================================================
--- ThirdParty/vmg/src/comm/comm.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/comm.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,169 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   comm.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Wed Jun 16 13:21:06 2010
+ *
+ * @brief  Base class for communication.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#ifdef HAVE_BOOST_FILESYSTEM_PATH_HPP
+#include <boost/filesystem.hpp>
+namespace fs = boost::filesystem;
+#endif
+
+#include "base/discretization.hpp"
+#include "base/helper.hpp"
+#include "base/stencil.hpp"
+#include "comm/comm.hpp"
+#include "comm/domain_decomposition.hpp"
+#include "grid/grid.hpp"
+#include "mg.hpp"
+
+using namespace VMG;
+
+vmg_float Comm::ComputeResidualNorm(Multigrid& sol_mg, Multigrid& rhs_mg)
+{
+  Stencil::iterator stencil_iter;
+  vmg_float norm = 0.0;
+
+  const vmg_float prefactor = MG::GetDiscretization()->OperatorPrefactor(sol_mg());
+  const Stencil& A = MG::GetDiscretization()->GetStencil();
+
+  this->CommToGhosts(sol_mg());
+
+  if (sol_mg().Global().BoundaryType() == LocallyRefined)
+    MG::GetDiscretization()->SetInnerBoundary(sol_mg(), rhs_mg(), sol_mg(sol_mg.Level()-1));
+
+  const Grid& sol = sol_mg();
+  const Grid& rhs = rhs_mg();
+
+#ifdef DEBUG_MATRIX_CHECKS
+  sol.IsCompatible(rhs);
+  sol.IsConsistent();
+  rhs.IsConsistent();
+#endif
+
+  for (int i=rhs.Local().Begin().X(); i<rhs.Local().End().X(); ++i)
+    for (int j=rhs.Local().Begin().Y(); j<rhs.Local().End().Y(); ++j)
+      for (int k=rhs.Local().Begin().Z(); k<rhs.Local().End().Z(); ++k) {
+	vmg_float val = rhs.GetVal(i,j,k) - prefactor * A.GetDiag() * sol.GetVal(i,j,k);
+	for (stencil_iter=A.begin(); stencil_iter!=A.end(); ++stencil_iter)
+	  val -= prefactor * stencil_iter->Val() * sol.GetVal(i + stencil_iter->Disp().X(),
+							      j + stencil_iter->Disp().Y(),
+							      k + stencil_iter->Disp().Z());
+	norm += val*val;
+      }
+
+  norm = GlobalSum(norm);
+  norm = std::sqrt(sol.Extent().MeshWidth().Product() * norm);
+
+  return norm;
+}
+
+Grid& Comm::GetParticleGrid()
+{
+  if (particle_grid != NULL)
+    return *particle_grid;
+
+  const Grid& grid = (*MG::GetRhs())(MG::GetRhs()->MaxLevel());
+
+  LocalIndices local = grid.Local();
+  GlobalIndices global = grid.Global();
+  SpatialExtent extent = grid.Extent();
+
+  const int& near_field_cells = MG::GetFactory().GetObjectStorageVal<int>("PARTICLE_NEAR_FIELD_CELLS");
+
+  for (int i=0; i<3; ++i) {
+
+    if (BoundaryConditions()[i] == Open || BoundaryConditions()[i] == Dirichlet) {
+
+      if (global.LocalBegin()[i] == global.GlobalBegin()[i])
+	global.LocalBegin()[i] += 1;
+
+      if (global.LocalEnd()[i] == global.GlobalEnd()[i])
+	global.LocalEnd()[i] -= 1;
+
+      global.LocalSize()[i] = global.LocalEnd()[i] - global.LocalBegin()[i];
+
+      global.GlobalBegin()[i] += 1;
+      global.GlobalEnd()[i] -= 1;
+      global.GlobalSize()[i] = global.GlobalEnd()[i] - global.GlobalBegin()[i];
+
+      extent.Begin()[i] += extent.MeshWidth()[i];
+      extent.End()[i] -= extent.MeshWidth()[i];
+      extent.Size()[i] = extent.End()[i] - extent.Begin()[i];
+
+    }
+
+  }
+
+  local.BoundaryBegin1() = local.BoundaryEnd1() = local.BoundarySize1() = 0;
+  local.BoundaryBegin2() = local.BoundaryEnd2() = local.BoundarySize2() = 0;
+
+  local.End() -= local.Begin();
+  local.Begin() = 0;
+
+  // Set grid size of intermediate temporary grid
+  for (int i=0; i<3; ++i) {
+
+    if (local.HaloSize1()[i] > 0) {
+      local.HaloBegin1()[i] = 0;
+      local.HaloEnd1()[i] = near_field_cells;
+      local.HaloSize1()[i] = near_field_cells;
+      local.Begin()[i] = near_field_cells;
+      local.End()[i] = local.Begin()[i] + local.Size()[i];
+    }
+
+    if (local.HaloSize2()[i] > 0) {
+      local.HaloBegin2()[i] = local.End()[i];
+      local.HaloEnd2()[i] = local.HaloBegin2()[i] + near_field_cells;
+      local.HaloSize2()[i] = near_field_cells;
+    }
+
+  }
+
+  local.SizeTotal() = local.Size() + local.HaloSize1() + local.HaloSize2();
+
+  particle_grid = new Grid(global, local, extent);
+
+  return *particle_grid;
+}
+
+Comm::~Comm()
+{
+  delete dd;
+  delete particle_grid;
+}
+
+const std::string& Comm::OutputPath()
+{
+  if (!output_directory_is_created) {
+    output_path_str = CreateOutputDirectory();
+    output_directory_is_created = true;
+  }
+
+  return output_path_str;
+}
Index: ThirdParty/vmg/src/comm/comm.hpp
===================================================================
--- ThirdParty/vmg/src/comm/comm.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/comm.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,163 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   comm.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Wed Jun 16 13:21:06 2010
+ *
+ * @brief  Base class for communication.
+ *
+ */
+
+#ifndef COMM_HPP_
+#define COMM_HPP_
+
+#include <map>
+#include <vector>
+
+#include "base/defs.hpp"
+#include "base/index.hpp"
+#include "base/object.hpp"
+#include "comm/domain_decomposition.hpp"
+#include "thirdparty/pugixml/pugixml.hpp"
+
+namespace VMG
+{
+
+class BoundaryValue;
+class DomainDecomposition;
+class GlobalIndices;
+class Grid;
+class Interface;
+class Multigrid;
+class Stencil;
+
+class Comm : public Object
+{
+public:
+  Comm(const Boundary& boundary, DomainDecomposition* domain_dec, bool register_ = true) :
+    Object("COMM", register_),
+    bound(boundary),
+    dd(domain_dec),
+    particle_grid(NULL),
+    output_directory_is_created(false),
+    output_count(0)
+  {}
+
+  virtual ~Comm();
+
+  virtual Grid& GetCoarserGrid(Multigrid& multigrid) = 0;
+  virtual Grid& GetFinerGrid(Multigrid& multigrid) = 0;
+  Grid& GetParticleGrid();
+
+  virtual void CommFromGhosts(Grid& grid) = 0;
+  virtual void CommToGhosts(Grid& grid) = 0;
+  virtual void CommSubgrid(Grid& grid_old, Grid& grid_new, const int& direction) = 0;
+  virtual void CommAddSubgrid(Grid& grid_old, Grid& grid_new, const int& direction) = 0;
+  virtual std::vector<BoundaryValue> CommBoundaryValues() = 0;
+
+  virtual void CommToGhostsAsyncStart(Grid& grid) = 0;
+  virtual void CommToGhostsAsyncFinish(Grid& grid) = 0;
+  virtual void CommFromGhostsAsyncStart(Grid& grid) = 0;
+  virtual void CommFromGhostsAsyncFinish(Grid& grid) = 0;
+
+  virtual void Barrier() {}
+
+  virtual vmg_float GlobalSum(vmg_float value) {return value;}
+  virtual vmg_float GlobalSumRoot(vmg_float value) {return value;}
+  virtual void GlobalSumArray(vmg_float* array, const vmg_int& size) {}
+  virtual vmg_float GlobalMax(vmg_float value) {return value;}
+  virtual vmg_float GlobalMaxRoot(vmg_float value) {return value;}
+  virtual void GlobalMaxArray(vmg_float* array, const vmg_int& size) {}
+  virtual void GlobalBroadcast(vmg_float& value) {}
+  virtual void GlobalGather(vmg_float& value, vmg_float* array) {array[0] = value;}
+
+  virtual vmg_int GlobalSum(vmg_int value) {return value;}
+  virtual vmg_int GlobalSumRoot(vmg_int value) {return value;}
+  virtual void GlobalSumArray(vmg_int* array, const vmg_int& size) {}
+  virtual vmg_int GlobalMax(vmg_int value) {return value;}
+  virtual vmg_int GlobalMaxRoot(vmg_int value) {return value;}
+  virtual void GlobalMaxArray(vmg_int* array, const vmg_int& size) {}
+  virtual void GlobalBroadcast(vmg_int& value) {}
+  virtual void GlobalGather(vmg_int& value, vmg_int* array) {array[0] = value;}
+
+  virtual void GlobalBroadcast(char* str) {}
+
+  virtual vmg_float LevelSum(const Grid& grid, vmg_float value) {return value;}
+  virtual vmg_float LevelSumRoot(const Grid& grid, vmg_float value) {return value;}
+  virtual void LevelSumArray(const Grid& grid,  vmg_float* array, const vmg_int& size) {}
+
+  virtual vmg_int LevelSum(const Grid& grid, vmg_int value) {return value;}
+  virtual vmg_int LevelSumRoot(const Grid& grid, vmg_int value) {return value;}
+  virtual void LevelSumArray(const Grid& grid, vmg_int* array, const vmg_int& size) {}
+
+  virtual void Print(const OutputLevel level, const char* format, ...) = 0;
+  virtual void PrintOnce(const OutputLevel level, const char* format, ...) = 0;
+
+  virtual void PrintXML(const std::string& filename, const std::string& xml_data) = 0;
+  virtual void PrintXMLAll(const std::string& filename, const std::string& xml_data) = 0;
+  virtual void PrintAllSettings() = 0;
+  virtual void PrintGrid(Grid& grid, const char* information) = 0;
+  virtual void PrintDefect(Grid& sol, Grid& rhs, const char* information) = 0;
+
+  virtual void DebugPrintError(const Grid& sol, const char* information) {}
+  virtual void DebugPrintErrorNorm(Grid& sol) {}
+  virtual void DebugPrintGridStructure(Multigrid& multigrid) {}
+
+  virtual int GlobalRank() const {return 0;}
+  virtual int GlobalSize() const {return 1;}
+  virtual Index GlobalPos() const {return Index(0);}
+  virtual Index GlobalProcs() const {return Index(1);}
+
+  virtual int Rank(const Grid& grid) const {return 0;}
+  virtual int Size(const Grid& grid) const {return 1;}
+  virtual Index Pos(const Grid& grid) const {return Index(0);}
+  virtual Index Procs(const Grid& grid) const {return Index(1);}
+
+  virtual void PostInit(Multigrid& sol, Multigrid& rhs) = 0;
+
+  const Boundary& BoundaryConditions() const {return bound;}
+
+  vmg_float ComputeResidualNorm(Multigrid& sol, Multigrid& rhs);
+
+  void ComputeDomainDecomposition(const Interface& interface) {dd->Compute(*this, interface, decomposed_global);}
+  const std::vector<GlobalIndices>& DecomposedGlobalMe() const {return decomposed_global.find(GlobalPos())->second;}
+  const std::map<Index, std::vector<GlobalIndices> >& DecomposedGlobal() const {return decomposed_global;}
+
+protected:
+  const std::string& OutputPath();
+  int OutputCount() {return ++output_count;}
+
+  std::map<Index, std::vector<GlobalIndices> > decomposed_global;
+
+  Boundary bound;
+  DomainDecomposition* dd;
+
+private:
+  virtual std::string CreateOutputDirectory() = 0;
+
+  Grid* particle_grid;
+  std::string output_path_str;
+  bool output_directory_is_created;
+  int output_count;
+};
+
+}
+
+#endif /* COMM_HPP_ */
Index: ThirdParty/vmg/src/comm/comm_mpi.cpp
===================================================================
--- ThirdParty/vmg/src/comm/comm_mpi.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/comm_mpi.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,1093 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   comm_mpi.cpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Wed Jun 16 13:21:06 2010
+ *
+ * @brief  Class for MPI-based communication.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#ifdef HAVE_MPI
+
+#include <mpi.h>
+#ifdef HAVE_MARMOT
+#include <enhancempicalls.h>
+#include <sourceinfompicalls.h>
+#endif
+
+#ifdef HAVE_BOOST_FILESYSTEM_PATH_HPP
+#include <boost/filesystem.hpp>
+namespace fs = boost::filesystem;
+#endif
+
+#ifdef HAVE_VTK
+#include <vtkAbstractArray.h>
+#include <vtkImageData.h>
+#include <vtkPointData.h>
+#include <vtkSmartPointer.h>
+#include <vtkXMLImageDataWriter.h>
+#endif
+
+#include <cstring>
+#include <sstream>
+
+#include "base/helper.hpp"
+#include "base/tuple.hpp"
+#include "comm/comm_mpi.hpp"
+#include "comm/mpi/datatypes_local.hpp"
+#include "discretization/boundary_value_setter.hpp"
+#include "grid/grid.hpp"
+#include "grid/multigrid.hpp"
+#include "grid/tempgrid.hpp"
+#include "mg.hpp"
+#include "base/timer.hpp"
+
+static char print_buffer[512];
+
+using namespace VMG;
+
+void CommMPI::IsendAll(Grid& grid, std::vector<VMG::MPI::Datatype>& types, const MPI_Comm& comm, const int& tag_start)
+{
+  for (std::vector<VMG::MPI::Datatype>::const_iterator iter=types.begin(); iter!=types.end(); ++iter)
+    iter->Isend(grid, tag_start, comm, Request());
+}
+
+void CommMPI::IrecvAll(Grid& grid, std::vector<VMG::MPI::Datatype>& types, const MPI_Comm& comm, const int& tag_start)
+{
+  for (std::vector<VMG::MPI::Datatype>::const_iterator iter=types.begin(); iter!=types.end(); ++iter)
+    iter->Irecv(grid, tag_start, comm, Request());
+}
+
+void CommMPI::IsendAllBuffered(const Grid& grid, std::vector<VMG::MPI::Datatype>& types, const MPI_Comm& comm, const int& tag_start)
+{
+  for (std::vector<VMG::MPI::Datatype>::iterator iter=types.begin(); iter!=types.end(); ++iter)
+    iter->IsendBuffered(grid, tag_start, comm, Request());
+}
+
+void CommMPI::IrecvAllBuffered(std::vector<VMG::MPI::Datatype>& types, const MPI_Comm& comm, const int& tag_start)
+{
+  for (std::vector<VMG::MPI::Datatype>::iterator iter=types.begin(); iter!=types.end(); ++iter)
+    iter->IrecvBuffered(tag_start, comm, Request());
+}
+
+void CommMPI::ReplaceBufferAll(Grid& grid, const std::vector<VMG::MPI::Datatype>& types)
+{
+  for (std::vector<VMG::MPI::Datatype>::const_iterator iter=types.begin(); iter!= types.end(); ++iter)
+    iter->GridReplace(grid);
+}
+
+void CommMPI::AddBufferAll(Grid& grid, const std::vector<VMG::MPI::Datatype>& types)
+{
+  for (std::vector<VMG::MPI::Datatype>::const_iterator iter=types.begin(); iter!= types.end(); ++iter)
+    iter->GridSum(grid);
+}
+
+void CommMPI::CommSubgrid(Grid& grid_old, Grid& grid_new, const int& direction)
+{
+  MPI_Comm comm = settings.CommunicatorGlobal(grid_old);
+  if (comm != MPI_COMM_NULL) {
+    VMG::MPI::DatatypesGlobal& datatypes = settings.DatatypesGlobal(grid_old, grid_new, direction);
+    IrecvAllBuffered(datatypes.Receive(), comm, 0411);
+    IsendAllBuffered(grid_old, datatypes.Send(), comm, 0411);
+    WaitAll();
+    ReplaceBufferAll(grid_new, datatypes.Receive());
+  }
+}
+
+void CommMPI::CommAddSubgrid(Grid& grid_old, Grid& grid_new, const int& direction)
+{
+  MPI_Comm comm = settings.CommunicatorGlobal(grid_old);
+  if (comm != MPI_COMM_NULL) {
+    VMG::MPI::DatatypesGlobal& datatypes = settings.DatatypesGlobal(grid_old, grid_new, direction);
+    IrecvAllBuffered(datatypes.Receive(), comm, 1806);
+    IsendAllBuffered(grid_old, datatypes.Send(), comm, 1806);
+    WaitAll();
+    AddBufferAll(grid_new, datatypes.Receive());
+  }
+}
+
+std::vector<BoundaryValue> CommMPI::CommBoundaryValues()
+{
+  MPI_Comm comm = settings.CommunicatorLocal((*MG::GetRhs())(MG::GetRhs()->MaxLevel()));
+  if (comm != MPI_COMM_NULL) {
+
+    const std::vector<BoundaryValue>& bvs = MG::GetBoundaryValueSetter()->BoundaryValues();
+
+    std::vector<vmg_float> val_buffer; val_buffer.reserve(bvs.size());
+
+    int comm_size, comm_rank;
+    MPI_Comm_rank(comm, &comm_rank);
+    MPI_Comm_size(comm, &comm_size);
+
+    std::vector< std::vector<BoundaryValue> > bvs_distributed(comm_size);
+
+    for (unsigned int i=0; i<bvs.size(); ++i)
+      bvs_distributed[settings.BoundaryRanks()[i]].push_back(bvs[i]);
+
+    int bvs_count[comm_size];
+    for (int i=0; i<comm_size; ++i) {
+      bvs_count[i] = bvs_distributed[i].size();
+      for (unsigned int j=0; j<bvs_distributed[i].size(); ++j)
+        val_buffer.push_back(bvs_distributed[i][j].Val());
+    }
+
+    MPI_Reduce_scatter(MPI_IN_PLACE, &val_buffer.front(), bvs_count, MPI_DOUBLE, MPI_SUM, comm);
+
+    int begin = 0;
+    for (int i=0; i<comm_rank; ++i)
+      begin += bvs_distributed[i].size();
+
+    for (unsigned int i=0; i<bvs_distributed[comm_rank].size(); ++i)
+      bvs_distributed[comm_rank][i].Val() = val_buffer[begin+i];
+
+    return bvs_distributed[comm_rank];
+  }
+
+  return std::vector<BoundaryValue>();
+}
+
+void CommMPI::CommToGhosts(Grid& grid)
+{
+  MPI_Comm comm = settings.CommunicatorLocal(grid);
+  if (comm != MPI_COMM_NULL) {
+    VMG::MPI::DatatypesLocal& types = settings.DatatypesLocal(grid);
+    IrecvAllBuffered(types.Halo(), comm, 2310);
+    IsendAllBuffered(grid, types.NB(), comm, 2310);
+    WaitAll();
+    ReplaceBufferAll(grid, types.Halo());
+  }
+}
+
+void CommMPI::CommToGhostsAsyncStart(Grid& grid)
+{
+  MPI_Comm comm = settings.CommunicatorLocal(grid);
+  if (comm != MPI_COMM_NULL) {
+    VMG::MPI::DatatypesLocal& types = settings.DatatypesLocal(grid);
+    IrecvAllBuffered(types.Halo(), comm, 2412);
+    IsendAllBuffered(grid, types.NB(), comm, 2412);
+    TestAll();
+  }
+}
+
+void CommMPI::CommToGhostsAsyncFinish(Grid& grid)
+{
+  WaitAll();
+  ReplaceBufferAll(grid, settings.DatatypesLocal(grid).Halo());
+}
+
+void CommMPI::CommFromGhosts(Grid& grid)
+{
+  MPI_Comm comm = settings.CommunicatorLocal(grid);
+  if (comm != MPI_COMM_NULL) {
+    VMG::MPI::DatatypesLocal& types = settings.DatatypesLocal(grid);
+    IrecvAllBuffered(types.NB(), comm, 1337);
+    IsendAllBuffered(grid, types.Halo(), comm, 1337);
+    WaitAll();
+    AddBufferAll(grid, types.NB());
+  }
+}
+
+void CommMPI::CommFromGhostsAsyncStart(Grid& grid)
+{
+  MPI_Comm comm = settings.CommunicatorLocal(grid);
+  if (comm != MPI_COMM_NULL) {
+    VMG::MPI::DatatypesLocal& types = settings.DatatypesLocal(grid);
+    IrecvAllBuffered(types.NB(), comm, 0xc0ffee);
+    IsendAllBuffered(grid, types.Halo(), comm, 0xc0ffee);
+    TestAll();
+  }
+}
+
+void CommMPI::CommFromGhostsAsyncFinish(Grid& grid)
+{
+  WaitAll();
+  AddBufferAll(grid, settings.DatatypesLocal(grid).NB());
+}
+
+void CommMPI::Barrier()
+{
+  MPI_Barrier(comm_global);
+}
+
+vmg_float CommMPI::GlobalSum(vmg_float value)
+{
+  vmg_float result = 0;
+  MPI_Allreduce(&value, &result, 1, MPI_DOUBLE, MPI_SUM, comm_global);
+  return result;
+}
+
+vmg_float CommMPI::GlobalSumRoot(vmg_float value)
+{
+  vmg_float result = 0;
+  MPI_Reduce(&value, &result, 1, MPI_DOUBLE, MPI_SUM, 0, comm_global);
+  return result;
+}
+
+void CommMPI::GlobalSumArray(vmg_float* array, const vmg_int& size)
+{
+  MPI_Allreduce(MPI_IN_PLACE, array, size, MPI_DOUBLE, MPI_SUM, comm_global);
+}
+
+void CommMPI::GlobalBroadcast(vmg_float& value)
+{
+  MPI_Bcast(&value, 1, MPI_DOUBLE, 0, comm_global);
+}
+
+void CommMPI::GlobalGather(vmg_float& value, vmg_float* array)
+{
+  MPI_Gather(&value, 1, MPI_DOUBLE, array, 1, MPI_DOUBLE, 0, comm_global);
+}
+
+vmg_int CommMPI::GlobalSum(vmg_int value)
+{
+  vmg_int result = 0;
+  MPI_Allreduce(&value, &result, 1, MPI_INT, MPI_SUM, comm_global);
+  return result;
+}
+
+vmg_int CommMPI::GlobalSumRoot(vmg_int value)
+{
+  vmg_int result = 0;
+  MPI_Reduce(&value, &result, 1, MPI_INT, MPI_SUM, 0, comm_global);
+  return result;
+}
+
+void CommMPI::GlobalSumArray(vmg_int* array, const vmg_int& size)
+{
+  MPI_Allreduce(MPI_IN_PLACE, array, size, MPI_INT, MPI_SUM, comm_global);
+}
+
+vmg_float CommMPI::GlobalMax(vmg_float value)
+{
+  vmg_float result = 0.0;
+  MPI_Allreduce(&value, &result, 1, MPI_DOUBLE, MPI_MAX, comm_global);
+  return result;
+}
+
+vmg_float CommMPI::GlobalMaxRoot(vmg_float value)
+{
+  vmg_float result = 0.0;
+  MPI_Reduce(&value, &result, 1, MPI_DOUBLE, MPI_MAX, 0, comm_global);
+  return result;
+}
+
+void CommMPI::GlobalMaxArray(vmg_float* array, const vmg_int& size)
+{
+  MPI_Allreduce(MPI_IN_PLACE, array, size, MPI_DOUBLE, MPI_MAX, comm_global);
+}
+
+vmg_int CommMPI::GlobalMax(vmg_int value)
+{
+  vmg_int result = 0;
+  MPI_Allreduce(&value, &result, 1, MPI_INT, MPI_MAX, comm_global);
+  return result;
+}
+
+vmg_int CommMPI::GlobalMaxRoot(vmg_int value)
+{
+  vmg_int result = 0;
+  MPI_Reduce(&value, &result, 1, MPI_INT, MPI_MAX, 0, comm_global);
+  return result;
+}
+
+void CommMPI::GlobalMaxArray(vmg_int* array, const vmg_int& size)
+{
+  MPI_Allreduce(MPI_IN_PLACE, array, size, MPI_INT, MPI_MAX, comm_global);
+}
+
+void CommMPI::GlobalBroadcast(vmg_int& value)
+{
+  MPI_Bcast(&value, 1, MPI_INT, 0, comm_global);
+}
+
+void CommMPI::GlobalGather(vmg_int& value, vmg_int* array)
+{
+  MPI_Gather(&value, 1, MPI_INT, array, 1, MPI_INT, 0, comm_global);
+}
+
+void CommMPI::GlobalBroadcast(char* str)
+{
+  int size = std::strlen(str) + 1;
+  MPI_Bcast(&size, 1, MPI_INT, 0, comm_global);
+  MPI_Bcast(str, size, MPI_CHAR, 0, comm_global);
+}
+
+vmg_float CommMPI::LevelSum(const Grid& grid, vmg_float value)
+{
+  vmg_float result = 0.0;
+  MPI_Comm comm = settings.CommunicatorLocal(grid);
+  assert(comm != MPI_COMM_NULL);
+  MPI_Allreduce(&value, &result, 1, MPI_DOUBLE, MPI_SUM, comm);
+  return result;
+}
+
+vmg_float CommMPI::LevelSumRoot(const Grid& grid, vmg_float value)
+{
+  vmg_float result = 0.0;
+  MPI_Comm comm = settings.CommunicatorLocal(grid);
+  assert(comm != MPI_COMM_NULL);
+  MPI_Reduce(&value, &result, 1, MPI_DOUBLE, MPI_SUM, 0, comm);
+  return result;
+}
+
+void CommMPI::LevelSumArray(const Grid& grid, vmg_float* array, const vmg_int& size)
+{
+  MPI_Comm comm = settings.CommunicatorLocal(grid);
+  assert(comm != MPI_COMM_NULL);
+  MPI_Allreduce(MPI_IN_PLACE, array, size, MPI_DOUBLE, MPI_SUM, comm);
+}
+
+vmg_int CommMPI::LevelSum(const Grid& grid, vmg_int value)
+{
+  vmg_int result = 0;
+  MPI_Comm comm = settings.CommunicatorLocal(grid);
+  assert(comm != MPI_COMM_NULL);
+  MPI_Allreduce(&value, &result, 1, MPI_INT, MPI_SUM, comm);
+  return result;
+}
+
+vmg_int CommMPI::LevelSumRoot(const Grid& grid, vmg_int value)
+{
+  vmg_int result = 0;
+  MPI_Comm comm = settings.CommunicatorLocal(grid);
+  assert(comm != MPI_COMM_NULL);
+  MPI_Reduce(&value, &result, 1, MPI_INT, MPI_SUM, 0, comm);
+  return result;
+}
+
+void CommMPI::LevelSumArray(const Grid& grid, vmg_int* array, const vmg_int& size)
+{
+  MPI_Comm comm = settings.CommunicatorLocal(grid);
+  assert(comm != MPI_COMM_NULL);
+  MPI_Allreduce(MPI_IN_PLACE, array, size, MPI_INT, MPI_SUM, comm);
+}
+
+//TODO: Variable length buffer
+void CommMPI::Print(const OutputLevel level, const char* format, ...)
+{
+  bool print = (level == Output);
+
+#ifdef OUTPUT_INFO
+  print |= (level == Info);
+#endif
+#ifdef OUTPUT_DEBUG
+  print |= (level == Debug);
+#endif
+#ifdef OUTPUT_TIMING
+  print |= (level == Timing);
+#endif
+
+  if (print) {
+    va_list args;
+    va_start(args, format);
+    vsprintf(print_buffer, format, args);
+    printf("vmg: Rank %d: %s\n", GlobalRank(), print_buffer);
+    va_end(args);
+  }
+}
+
+void CommMPI::PrintOnce(const OutputLevel level, const char* format, ...)
+{
+  bool print = (level == Output);
+
+#ifdef OUTPUT_INFO
+  print |= (level == Info);
+#endif
+#ifdef OUTPUT_DEBUG
+  print |= (level == Debug);
+#endif
+#ifdef OUTPUT_TIMING
+  print |= (level == Timing);
+#endif
+
+  if (GlobalRank() == 0 && print) {
+    va_list args;
+    va_start(args, format);
+    vsprintf(print_buffer, format, args);
+    printf("vmg: Rank %d: %s\n", GlobalRank(), print_buffer);
+    va_end(args);
+  }
+}
+
+void CommMPI::PrintXML(const std::string& filename, const std::string& xml_data)
+{
+  MPI_File file;
+  std::stringstream path, xml_header;
+
+  pugi::xml_document doc;
+  pugi::xml_node node_data = doc.append_child("Global").append_child("NumProcesses").append_child(pugi::node_pcdata);
+  node_data.set_value(Helper::ToString(GlobalProcs().Product()).c_str());
+  doc.save(xml_header);
+
+  path << OutputPath() << filename;
+
+  char* filename_array = Helper::GetCharArray(path.str());
+  char* xml_header_array = Helper::GetCharArray(xml_header.str());
+  char* str_array = Helper::GetCharArray(xml_data);
+
+  MPI_File_open(MPI_COMM_SELF, filename_array, MPI_MODE_WRONLY|MPI_MODE_CREATE, MPI_INFO_NULL, &file);
+  MPI_File_set_size(file, 0);
+  MPI_File_write(file, xml_header_array, xml_header.str().size(), MPI_CHAR, MPI_STATUS_IGNORE);
+  MPI_File_write(file, str_array, xml_data.size(), MPI_CHAR, MPI_STATUS_IGNORE);
+  MPI_File_close(&file);
+
+  delete [] filename_array;
+  delete [] xml_header_array;
+  delete [] str_array;
+}
+
+void CommMPI::PrintXMLAll(const std::string& filename, const std::string& xml_data)
+{
+  MPI_File file;
+  std::stringstream path;
+
+  path << OutputPath() << filename;
+
+  char* filename_array = Helper::GetCharArray(path.str());
+  char* str_array = Helper::GetCharArray(xml_data);
+
+  MPI_File_open(comm_global, filename_array, MPI_MODE_WRONLY|MPI_MODE_CREATE, MPI_INFO_NULL, &file);
+  MPI_File_set_size(file, 0);
+
+  if (GlobalRank() == 0) {
+    std::stringstream xml_header;
+    pugi::xml_document doc;
+    pugi::xml_node node_data = doc.append_child("Global").append_child("NumProcesses").append_child(pugi::node_pcdata);
+    node_data.set_value(Helper::ToString(GlobalProcs().Product()).c_str());
+    doc.save(xml_header);
+
+    char* xml_header_array = Helper::GetCharArray(xml_header.str());
+
+    MPI_File_write_shared(file, xml_header_array, xml_header.str().size(), MPI_CHAR, MPI_STATUS_IGNORE);
+
+    delete [] xml_header_array;
+  }
+
+  MPI_File_write_ordered(file, str_array, xml_data.size(), MPI_CHAR, MPI_STATUS_IGNORE);
+  MPI_File_close(&file);
+
+  delete [] filename_array;
+  delete [] str_array;
+}
+
+void CommMPI::PrintGridInformation(const Grid& grid, char* filename, const std::string& name)
+{
+  std::stringstream buf;
+  MPI_File file;
+  int rank, size;
+  int size_local, size_local_max;
+
+  MPI_Comm comm = settings.CommunicatorGlobal(grid);
+  MPI_Comm comm_local = settings.CommunicatorLocal(grid);
+
+  if (comm_local != MPI_COMM_NULL)
+    MPI_Comm_size(comm_local, &size_local);
+  else
+    size_local = 0;
+
+  if (comm != MPI_COMM_NULL) {
+
+    MPI_Reduce(&size_local, &size_local_max, 1, MPI_INT, MPI_MAX, 0, comm);
+
+    MPI_File_open(comm, filename, MPI_MODE_WRONLY|MPI_MODE_CREATE|MPI_MODE_APPEND, MPI_INFO_NULL, &file);
+
+    MPI_Comm_rank(comm, &rank);
+    MPI_Comm_size(comm, &size);
+
+    if (rank == 0) {
+
+      buf << "###########################################################" << std::endl
+          << "GLOBAL INFORMATION:" << std::endl
+          << "  NAME:                 " << name << std::endl
+          << "  LEVEL:                " << grid.Level() << std::endl
+          << "  COMM_SIZE_GLOBAL:     " << size << std::endl
+          << "  COMM_SIZE_LOCAL:      " << size_local_max << std::endl
+          << "  GLOBAL:" << std::endl
+          << "    GLOBAL_BEGIN:            " << grid.Global().GlobalBegin() << std::endl
+          << "    GLOBAL_END:              " << grid.Global().GlobalEnd() << std::endl
+          << "    GLOBAL_SIZE:             " << grid.Global().GlobalSize() << std::endl
+          << "    GLOBAL_BEGIN_FINEST:     " << grid.Global().GlobalBeginFinest() << std::endl
+          << "    GLOBAL_END_FINEST:       " << grid.Global().GlobalEndFinest() << std::endl
+          << "    GLOBAL_SIZE_FINEST:      " << grid.Global().GlobalSizeFinest() << std::endl
+          << "  EXTENT:" << std::endl
+          << "    BEGIN:              " << grid.Extent().Begin() << std::endl
+          << "    END:                " << grid.Extent().End() << std::endl
+          << "    SIZE:               " << grid.Extent().Size() << std::endl
+          << "    MESH_WIDTH:         " << grid.Extent().MeshWidth() << std::endl
+          << std::endl
+          << "LOCAL INFORMATION:" << std::endl;
+    }
+
+    buf << "RANK " << rank << ":" << std::endl
+        << "  GLOBAL:" << std::endl
+        << "    LOCAL_BEGIN:          " << grid.Global().LocalBegin() << std::endl
+        << "    LOCAL_END:            " << grid.Global().LocalEnd() << std::endl
+        << "    LOCAL_SIZE:           " << grid.Global().LocalSize() << std::endl
+        << "    BOUNDARY_TYPE:        " << grid.Global().BoundaryType() << std::endl
+        << "  LOCAL:" << std::endl
+        << "    BEGIN:                " << grid.Local().Begin() << std::endl
+        << "    END:                  " << grid.Local().End() << std::endl
+        << "    SIZE:                 " << grid.Local().Size() << std::endl
+        << "    SIZE_TOTAL:           " << grid.Local().SizeTotal() << std::endl
+        << "    HALO_BEGIN_1:         " << grid.Local().HaloBegin1() << std::endl
+        << "    HALO_END_1:           " << grid.Local().HaloEnd1() << std::endl
+        << "    HALO_SIZE_1:          " << grid.Local().HaloSize1() << std::endl
+        << "    HALO_BEGIN_2:         " << grid.Local().HaloBegin2() << std::endl
+        << "    HALO_END_2:           " << grid.Local().HaloEnd2() << std::endl
+        << "    HALO_SIZE_2:          " << grid.Local().HaloSize2() << std::endl
+        << "    BOUNDARY_BEGIN_1:     " << grid.Local().BoundaryBegin1() << std::endl
+        << "    BOUNDARY_END_1:       " << grid.Local().BoundaryEnd1() << std::endl
+        << "    BOUNDARY_SIZE_1:      " << grid.Local().BoundarySize1() << std::endl
+        << "    BOUNDARY_BEGIN_2:     " << grid.Local().BoundaryBegin2() << std::endl
+        << "    BOUNDARY_END_2:       " << grid.Local().BoundaryEnd2() << std::endl
+        << "    BOUNDARY_SIZE_2:      " << grid.Local().BoundarySize2() << std::endl;
+
+    if (rank == size-1)
+      buf << "###########################################################" << std::endl;
+
+    char* char_buf = Helper::GetCharArray(buf.str());
+    MPI_File_write_ordered(file, char_buf, buf.str().size(), MPI_CHAR, MPI_STATUS_IGNORE);
+    delete [] char_buf;
+
+    MPI_File_close(&file);
+
+  }
+}
+
+void CommMPI::PrintDatatypes(char* filename)
+{
+  std::stringstream buf;
+  MPI_File file;
+
+  if (comm_global != MPI_COMM_NULL) {
+
+    MPI_File_open(comm_global, filename, MPI_MODE_WRONLY|MPI_MODE_CREATE|MPI_MODE_APPEND, MPI_INFO_NULL, &file);
+
+    buf << settings << std::endl;
+
+    char* char_buf = Helper::GetCharArray(buf.str());
+    MPI_File_write_ordered(file, char_buf, buf.str().size(), MPI_CHAR, MPI_STATUS_IGNORE);
+    delete [] char_buf;
+
+    MPI_File_close(&file);
+
+  }
+}
+
+void CommMPI::PrintAllSettings()
+{
+  std::stringstream buf;
+  MPI_File file;
+
+  const Multigrid& mg = *MG::GetSol();
+
+  buf << OutputPath() << "settings.txt";
+  char *filename = Helper::GetCharArray(buf.str());
+
+  MPI_File_open(comm_global, filename, MPI_MODE_WRONLY|MPI_MODE_CREATE, MPI_INFO_NULL, &file);
+  MPI_File_set_size(file, 0);
+  MPI_File_close(&file);
+
+  for (int i=mg.MinLevel(); i<=mg.MaxLevel(); ++i)
+    PrintGridInformation(mg(i), filename, "MULTIGRID");
+
+  for (int i=mg.MinLevel()+1; i<=mg.MaxLevel(); ++i)
+    PrintGridInformation(settings.CoarserGrid(mg(i)), filename, "COARSER_GRID");
+
+  for (int i=mg.MinLevel(); i<mg.MaxLevel(); ++i)
+    PrintGridInformation(settings.FinerGrid(mg(i)), filename, "FINER_GRID");
+
+  PrintDatatypes(filename);
+
+  delete [] filename;
+
+}
+
+void CommMPI::PrintGrid(Grid& grid, const char* information)
+{
+  int output_count = OutputCount();
+
+#ifdef HAVE_VTK
+
+  if (settings.CommunicatorLocal(grid) != MPI_COMM_NULL) {
+
+    vtkSmartPointer<vtkImageData> image = vtkSmartPointer<vtkImageData>::New();
+    image->SetExtent(grid.Global().LocalBegin().X(), grid.Global().LocalEnd().X()-1,
+                     grid.Global().LocalBegin().Y(), grid.Global().LocalEnd().Y()-1,
+                     grid.Global().LocalBegin().Z(), grid.Global().LocalEnd().Z()-1);
+
+    image->SetSpacing(grid.Extent().MeshWidth().vec());
+    image->SetOrigin(grid.Extent().Begin().vec());
+    image->SetScalarTypeToDouble();
+    image->SetNumberOfScalarComponents(1);
+    image->AllocateScalars();
+    image->GetPointData()->GetScalars()->SetName(information);
+
+    Index i;
+    Index begin, end, offset;
+
+    for (int j=0; j<3; ++j) {
+      if (grid.Local().BoundarySize1()[j] > 0) {
+        begin[j] = grid.Local().BoundaryBegin1()[j];
+        offset[j] = grid.Global().LocalBegin()[j];
+      } else {
+        begin[j] = grid.Local().Begin()[j];
+        offset[j] =  grid.Global().LocalBegin()[j] - grid.Local().Begin()[j];
+      }
+      end[j] = (grid.Local().BoundarySize2()[j] > 0 ? grid.Local().BoundaryEnd2()[j] : grid.Local().End()[j]);
+    }
+
+    for (i.X()=begin.X(); i.X()<end.X(); ++i.X())
+      for (i.Y()=begin.Y(); i.Y()<end.Y(); ++i.Y())
+        for (i.Z()=begin.Z(); i.Z()<end.Z(); ++i.Z())
+          image->SetScalarComponentFromDouble(i.X() + offset.X(),
+                                              i.Y() + offset.Y(),
+                                              i.Z() + offset.Z(),
+                                              0, grid.GetVal(i));
+
+    image->Update();
+
+    int rank, size;
+    MPI_Comm_rank(comm_global, &rank);
+    MPI_Comm_size(comm_global, &size);
+
+    std::stringstream filename;
+    filename << OutputPath() << std::setw(4) << std::setfill('0') << output_count << "_" << rank << ".vti";
+
+    vtkSmartPointer<vtkXMLImageDataWriter> writer = vtkSmartPointer<vtkXMLImageDataWriter>::New();
+    writer->SetFileName(filename.str().c_str());
+    writer->SetCompressorTypeToNone();
+    writer->SetDataModeToAscii();
+    writer->SetInput(image);
+    writer->Update();
+    writer->Write();
+
+  }
+
+//FIXME
+#else /* HAVE_VTK */
+  Index i;
+  std::stringstream buf;
+
+  Index begin, end;
+  Index begin_local, end_local, begin_global, end_global;
+
+  CommToGhosts(grid);
+
+  for (int i=0; i<3; ++i) {
+    end[i] = grid.Local().End()[i] + (grid.Global().LocalEnd()[i] == grid.Global().GlobalSize()[i] ? 0 : grid.Local().HaloSize1()[i]);
+    end_local[i] = grid.Global().LocalEnd()[i] - (grid.Global().LocalEnd()[i] == grid.Global().GlobalSize()[i] ? 1 : 0);
+  }
+
+  begin = grid.Local().Begin();
+  begin_local = grid.Global().LocalBegin();
+  begin_global = 0;
+  end_global = grid.Global().GlobalSize()-1;
+
+  for (i.Z()=begin.Z(); i.Z()<end.Z(); ++i.Z())
+    for (i.Y()=begin.Y(); i.Y()<end.Y(); ++i.Y())
+      for (i.X()=begin.X(); i.X()<end.X(); ++i.X())
+        buf << std::scientific << grid.GetVal(i) << " ";
+
+  CreateOutputFiles(grid, buf, information,
+      begin_global, end_global,
+      begin_local, end_local,
+      output_count);
+#endif /* HAVE_VTK */
+}
+
+void CommMPI::PrintDefect(Grid& sol, Grid& rhs, const char* information)
+{
+  TempGrid *temp = MG::GetTempGrid();
+  temp->SetProperties(sol);
+  temp->ImportFromResidual(sol, rhs);
+  PrintGrid(*temp, information);
+}
+
+void CommMPI::CreateOutputFiles(const Grid& grid, const std::stringstream& serial_data, const char* information,
+    const Index& begin_global, const Index& end_global,
+    const Index& begin_local, const Index& end_local,
+    const int& output_count)
+{
+  MPI_Comm comm = settings.CommunicatorGlobal(grid);
+
+  if (comm != MPI_COMM_NULL) {
+
+    MPI_File file;
+    std::string conv_information = Helper::ReplaceWhitespaces(information, "_");
+
+    CreateParallelOutputFile(grid, comm, output_count, conv_information.c_str(),
+        begin_global, end_global, begin_local, end_local);
+
+    file = CreateSerialOutputFile(grid, comm, output_count, conv_information.c_str(),
+        begin_global, end_global, begin_local, end_local);
+
+    char *char_buf = Helper::GetCharArray(serial_data.str());
+    MPI_File_write(file, char_buf, serial_data.str().size(), MPI_CHAR, MPI_STATUS_IGNORE);
+    delete [] char_buf;
+
+    FinalizeSerialOutputFile(file);
+
+  }
+}
+
+void CommMPI::CreateParallelOutputFile(const Grid& grid, MPI_Comm& comm,
+    const int& output_count, const char* information,
+    const Index& begin_global, const Index& end_global,
+    const Index& begin_local, const Index& end_local)
+{
+  int rank;
+  MPI_File file;
+  char parallel_filename[513], serial_filename[513];
+  std::stringstream buf;
+
+  MPI_Comm_rank(comm, &rank);
+
+  sprintf(parallel_filename, "%s%04d.pvti", OutputPath().c_str(), output_count);
+  sprintf(serial_filename, "%04d_%d.vti", output_count, rank);
+
+  MPI_File_open(comm, parallel_filename, MPI_MODE_WRONLY|MPI_MODE_CREATE, MPI_INFO_NULL, &file);
+  MPI_File_set_size(file, 0);
+
+  if (rank == 0) {
+
+    buf << "<?xml version=\"1.0\"?>" << std::endl
+        << "<VTKFile type=\"PImageData\" version=\"0.1\" byte_order=\"LittleEndian\">" << std::endl
+        << "  <PImageData WholeExtent=\"";
+
+    for (int i=0; i<3; ++i)
+      buf << begin_global[i] << " " << end_global[i] << " ";
+
+    buf << "\"" << std::endl
+        << "              GhostLevel=\"0\" Origin=\"0 0 0\" Spacing=\"";
+
+    for (int i=0; i<3; ++i)
+      buf << grid.Extent().MeshWidth()[i] << " ";
+
+    buf << "\">" << std::endl
+        << "    <PPointData Scalars=\"" << information << "\">" << std::endl
+        << "      <PDataArray type=\"Float32\" Name=\"" << information << "\"/>" << std::endl
+        << "    </PPointData>" << std::endl;
+
+    char* char_buf = Helper::GetCharArray(buf.str());
+
+    MPI_File_write_shared(file, char_buf, buf.str().size(), MPI_CHAR, MPI_STATUS_IGNORE);
+
+    delete [] char_buf;
+  }
+
+  buf.str("");
+
+  if ((end_local-begin_local).Product() > 0) {
+    buf << "    <Piece Extent=\"";
+
+    for (int i=0; i<3; ++i)
+      buf << begin_local[i] << " " << end_local[i]  << " ";
+
+    buf << "\" Source=\"" << serial_filename << "\"/>" << std::endl;
+  }
+
+  char* char_buf = Helper::GetCharArray(buf.str());
+
+  MPI_File_write_ordered(file, char_buf, buf.str().size(), MPI_CHAR, MPI_STATUS_IGNORE);
+
+  delete [] char_buf;
+
+  if (rank == 0) {
+
+    buf.str("");
+
+    buf << "  </PImageData>" << std::endl
+        << "</VTKFile>" << std::endl;
+
+    char* char_buf = Helper::GetCharArray(buf.str());
+
+    MPI_File_write_shared(file, char_buf, buf.str().size(), MPI_CHAR, MPI_STATUS_IGNORE);
+
+    delete [] char_buf;
+  }
+
+  MPI_File_close(&file);
+}
+
+MPI_File CommMPI::CreateSerialOutputFile(const Grid& grid, MPI_Comm& comm,
+    const int& output_count, const char* information,
+    const Index& begin_global, const Index& end_global,
+    const Index& begin_local, const Index& end_local)
+{
+  char serial_filename[513];
+  int rank;
+  MPI_File file;
+  std::stringstream buf;
+
+  MPI_Comm_rank(comm_global, &rank);
+
+  sprintf(serial_filename, "%s%04d_%d.vti", OutputPath().c_str(), output_count, rank);
+
+  MPI_File_open(MPI_COMM_SELF, serial_filename, MPI_MODE_WRONLY|MPI_MODE_CREATE, MPI_INFO_NULL, &file);
+
+  buf << "<?xml version=\"1.0\"?>" << std::endl
+      << "<VTKFile type=\"ImageData\" version=\"0.1\" byte_order=\"LittleEndian\">" << std::endl
+      << "  <ImageData WholeExtent=\"";
+
+  for (int i=0; i<3; ++i)
+    buf << begin_global[i] << " " << end_global[i] << " ";
+
+  buf << "\"" << std::endl
+      << "             Origin=\"0 0 0\" Spacing=\"";
+
+  for (int i=0; i<3; ++i)
+    buf << grid.Extent().MeshWidth()[i] << " ";
+
+  buf << "\">" << std::endl
+      << "    <Piece Extent=\"";
+
+  for (int i=0; i<3; ++i)
+    buf << begin_local[i] << " " << end_local[i] << " ";
+
+  buf << "\">" << std::endl
+      << "      <PointData Scalars=\"" << information << "\">" << std::endl
+      << "        <DataArray type=\"Float32\" Name=\"" << information << "\" format=\"ascii\">" << std::endl
+      << "          ";
+
+  char* char_buf = Helper::GetCharArray(buf.str());
+  MPI_File_write(file, char_buf, buf.str().size(), MPI_CHAR, MPI_STATUS_IGNORE);
+  delete [] char_buf;
+
+  return file;
+}
+
+void CommMPI::FinalizeSerialOutputFile(MPI_File& file)
+{
+  std::stringstream buf;
+
+  buf << std::endl
+      << "        </DataArray>" << std::endl
+      << "      </PointData>" << std::endl
+      << "    </Piece>" << std::endl
+      << "  </ImageData>" << std::endl
+      << "</VTKFile>" << std::endl;
+
+  char* char_buf = Helper::GetCharArray(buf.str());
+  MPI_File_write(file, char_buf, buf.str().size(), MPI_CHAR, MPI_STATUS_IGNORE);
+  delete [] char_buf;
+
+  MPI_File_close(&file);
+}
+
+int CommMPI::GlobalRank() const
+{
+  int rank;
+  MPI_Comm_rank(comm_global, &rank);
+  return rank;
+}
+
+int CommMPI::GlobalSize() const
+{
+  int size;
+  MPI_Comm_size(comm_global, &size);
+  return size;
+}
+
+Index CommMPI::GlobalPos() const
+{
+  Index dims, periods, coords;
+  MPI_Cart_get(comm_global, 3, dims.vec(), periods.vec(), coords.vec());
+  return coords;
+}
+
+Index CommMPI::GlobalProcs() const
+{
+  Index dims, periods, coords;
+  MPI_Cart_get(comm_global, 3, dims.vec(), periods.vec(), coords.vec());
+  return dims;
+}
+
+int CommMPI::Rank(const Grid& grid) const
+{
+  int rank;
+  MPI_Comm comm = settings.CommunicatorLocal(grid);
+  assert(comm != MPI_COMM_NULL);
+  MPI_Comm_rank(comm, &rank);
+  return rank;
+}
+
+int CommMPI::Size(const Grid& grid) const
+{
+  int size;
+  MPI_Comm comm = settings.CommunicatorLocal(grid);
+  assert(comm != MPI_COMM_NULL);
+  MPI_Comm_size(comm, &size);
+  return size;
+}
+
+Index CommMPI::Pos(const Grid& grid) const
+{
+  Index dims, periods, coords;
+  MPI_Comm comm = settings.CommunicatorLocal(grid);
+  assert(comm != MPI_COMM_NULL);
+  MPI_Cart_get(comm, 3, dims.vec(), periods.vec(), coords.vec());
+  return coords;
+}
+
+Index CommMPI::Procs(const Grid& grid) const
+{
+  Index dims, periods, coords;
+  MPI_Comm comm = settings.CommunicatorLocal(grid);
+  assert(comm != MPI_COMM_NULL);
+  MPI_Cart_get(comm, 3, dims.vec(), periods.vec(), coords.vec());
+  return dims;
+}
+
+void CommMPI::InitCommMPI(const MPI_Comm& comm)
+{
+  int status, size, rank;
+  int dims[3] = {0, 0, 0};
+  int periods[3];
+
+  for (int i=0; i<3; ++i)
+    periods[i] = (BoundaryConditions()[i] == Periodic ? 1 : 0);
+
+  MPI_Comm_size(comm, &size);
+  MPI_Comm_rank(comm, &rank);
+
+  MPI_Topo_test(comm, &status);
+
+  if (status == MPI_CART) {
+
+    comm_global = comm;
+
+  }else {
+
+    const int log2 = Helper::log_2(size);
+
+    if (Helper::intpow(2, log2) == size) {
+      for (int i=0; i<3; ++i)
+        dims[i] = Helper::intpow(2, log2 / 3 + (log2%3 > i ? 1 : 0));
+    }else {
+      MPI_Dims_create(size, 3, dims);
+    }
+
+#ifdef OUTPUT_DEBUG
+    if (rank == 0)
+      std::printf("Process grid: %d %d %d\n", dims[0], dims[1], dims[2]);
+#endif
+
+    MPI_Cart_create(comm, 3, dims, periods, 1, &comm_global);
+
+  }
+
+  MPI_Info_create(&info);
+  char key[] = "no_locks";
+  char val[] = "true";
+  MPI_Info_set(info, key, val);
+
+}
+
+CommMPI::~CommMPI()
+{
+  MPI_Comm_free(&comm_global);
+#ifdef VMG_ONE_SIDED
+  if (win_created)
+    MPI_Win_free(&win);
+#endif
+  MPI_Info_free(&info);
+}
+
+Grid& CommMPI::GetCoarserGrid(Multigrid& multigrid)
+{
+  return settings.CoarserGrid(multigrid(multigrid.Level()));
+}
+
+Grid& CommMPI::GetFinerGrid(Multigrid& multigrid)
+{
+  return settings.FinerGrid(multigrid(multigrid.Level()));
+}
+
+std::string CommMPI::CreateOutputDirectory()
+{
+#ifdef HAVE_BOOST_FILESYSTEM_PATH_HPP
+  std::string path, unique_path;
+  std::stringstream unique_suffix;
+  int suffix_counter = 0;
+  char buffer[129];
+  time_t rawtime;
+  struct tm *timeinfo;
+  int path_size;
+  char* path_array;
+
+  if (GlobalRank() == 0) {
+
+    time(&rawtime);
+    timeinfo = localtime(&rawtime);
+    strftime(buffer, 128, "./output/%Y_%m_%d_%H_%M_%S/", timeinfo);
+    path = buffer;
+
+    if (!fs::exists("output"))
+      fs::create_directory("output");
+
+    unique_path = path;
+
+    while (fs::exists(unique_path.c_str())) {
+
+      unique_suffix.str("");
+      unique_suffix << "_" << suffix_counter++ << "/";
+
+      unique_path = path;
+      unique_path.replace(unique_path.size()-1, 1, unique_suffix.str());
+
+    }
+
+    fs::create_directory(unique_path.c_str());
+
+    path_size = unique_path.size() + 1;
+    path_array = Helper::GetCharArray(unique_path);
+
+    MPI_Bcast(&path_size, 1, MPI_INT, 0, comm_global);
+
+  }else {
+
+    MPI_Bcast(&path_size, 1, MPI_INT, 0, comm_global);
+    path_array = new char[path_size];
+
+  }
+
+  MPI_Bcast(path_array, path_size, MPI_CHAR, 0, comm_global);
+
+  unique_path = path_array;
+
+  delete [] path_array;
+
+  return unique_path;
+
+#else
+
+  return "./";
+
+#endif
+}
+
+
+#endif /* HAVE_MPI */
Index: ThirdParty/vmg/src/comm/comm_mpi.hpp
===================================================================
--- ThirdParty/vmg/src/comm/comm_mpi.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/comm_mpi.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,197 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   comm_mpi.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Wed Jun 16 13:21:06 2010
+ *
+ * @brief  Class for MPI-based communication.
+ *
+ */
+
+#ifndef COMM_MPI_HPP_
+#define COMM_MPI_HPP_
+
+#ifndef HAVE_MPI
+#error You need MPI in order to compile CommMPI
+#endif
+
+#include <cstdarg>
+#include <cstdio>
+#include <list>
+#include <map>
+#include <vector>
+
+#include "base/has_tempgrids.hpp"
+#include "comm/comm.hpp"
+#include "comm/mpi/settings.hpp"
+#include "comm/mpi/has_request_vec.hpp"
+
+namespace VMG
+{
+
+class DomainDecomposition;
+class GridIteratorSet;
+class TempGrid;
+
+class CommMPI : public Comm, public HasTempGrids, public HasRequestVec
+{
+public:
+  CommMPI(const Boundary& boundary, DomainDecomposition* domain_dec, const MPI_Comm& mpi_comm, bool register_ = true) :
+    Comm(boundary, domain_dec, register_)
+#ifdef VMG_ONE_SIDED
+    ,win_created(false)
+#endif
+  {
+    InitCommMPI(mpi_comm);
+  }
+
+  CommMPI(const Boundary& boundary, DomainDecomposition* domain_dec, bool register_ = true) :
+    Comm(boundary, domain_dec, register_)
+#ifdef VMG_ONE_SIDED
+    ,win_created(false)
+#endif
+  {
+    InitCommMPI(MPI_COMM_WORLD);
+  }
+
+  virtual ~CommMPI();
+
+  Grid& GetCoarserGrid(Multigrid& multigrid);
+  Grid& GetFinerGrid(Multigrid& multigrid);
+
+  void CommFromGhosts(Grid& grid);
+  void CommToGhosts(Grid& grid);
+  void CommSubgrid(Grid& grid_old, Grid& grid_new, const int& direction);
+  void CommAddSubgrid(Grid& grid_old, Grid& grid_new, const int& direction);
+  std::vector<BoundaryValue> CommBoundaryValues();
+
+  void CommToGhostsAsyncStart(Grid& grid);
+  void CommToGhostsAsyncFinish(Grid& grid);
+  void CommFromGhostsAsyncStart(Grid& grid);
+  void CommFromGhostsAsyncFinish(Grid& grid);
+
+  void Barrier();
+
+  vmg_float GlobalSum(vmg_float value);
+  vmg_float GlobalSumRoot(vmg_float value);
+  void GlobalSumArray(vmg_float* array, const vmg_int& size);
+  vmg_float GlobalMax(vmg_float value);
+  vmg_float GlobalMaxRoot(vmg_float value);
+  void GlobalMaxArray(vmg_float* array, const vmg_int& size);
+  void GlobalBroadcast(vmg_float& value);
+  void GlobalGather(vmg_float& value, vmg_float* array);
+
+  vmg_int GlobalSum(vmg_int value);
+  vmg_int GlobalSumRoot(vmg_int value);
+  void GlobalSumArray(vmg_int* array, const vmg_int& size);
+  vmg_int GlobalMax(vmg_int value);
+  vmg_int GlobalMaxRoot(vmg_int value);
+  void GlobalMaxArray(vmg_int* array, const vmg_int& size);
+  void GlobalBroadcast(vmg_int& value);
+  void GlobalGather(vmg_int& value, vmg_int* array);
+
+  void GlobalBroadcast(char* str);
+
+  vmg_float LevelSum(const Grid& grid, vmg_float value);
+  vmg_float LevelSumRoot(const Grid& grid, vmg_float value);
+  void LevelSumArray(const Grid& grid,  vmg_float* array, const vmg_int& size);
+
+  vmg_int LevelSum(const Grid& grid, vmg_int value);
+  vmg_int LevelSumRoot(const Grid& grid, vmg_int value);
+  void LevelSumArray(const Grid& grid,  vmg_int* array, const vmg_int& size);
+
+  void Print(const OutputLevel level, const char* format, ...);
+  void PrintOnce(const OutputLevel level, const char* format, ...);
+
+  void PrintXML(const std::string& filename, const std::string& xml_data);
+  void PrintXMLAll(const std::string& filename, const std::string& xml_data);
+  void PrintAllSettings();
+  void PrintGrid(Grid& grid, const char* information);
+  void PrintDefect(Grid& sol, Grid& rhs, const char* information);
+
+  virtual int GlobalRank() const;
+  virtual int GlobalSize() const;
+  virtual Index GlobalPos() const;
+  virtual Index GlobalProcs() const;
+
+  virtual int Rank(const Grid& grid) const;
+  virtual int Size(const Grid& grid) const;
+  virtual Index Pos(const Grid& grid) const;
+  virtual Index Procs(const Grid& grid) const;
+
+  void PostInit(Multigrid& sol, Multigrid& rhs)
+  {
+    settings.ComputeSettings(sol, rhs, comm_global);
+  }
+
+  virtual void DebugPrintError(const Grid& sol, const char* information) {}
+  virtual void DebugPrintErrorNorm(Grid& sol) {}
+  virtual void DebugPrintGridStructure(Multigrid& multigrid) {}
+
+private:
+  void InitCommMPI(const MPI_Comm& comm);
+
+  void CreateOutputFiles(const Grid& grid, const std::stringstream& serial_data, const char* information,
+			 const Index& begin_global, const Index& end_global,
+			 const Index& begin_local, const Index& end_local,
+			 const int& output_count);
+
+  void CreateParallelOutputFile(const Grid& grid, MPI_Comm& comm,
+				const int& output_count, const char* information,
+				const Index& begin_global, const Index& end_global,
+				const Index& begin_local, const Index& end_local);
+
+  MPI_File CreateSerialOutputFile(const Grid& grid, MPI_Comm& comm,
+				  const int& output_count, const char* information,
+				  const Index& begin_global, const Index& end_global,
+				  const Index& begin_local, const Index& end_local);
+
+  void FinalizeSerialOutputFile(MPI_File& file);
+
+  std::string CreateOutputDirectory();
+
+protected:
+  void IsendAll(Grid& grid, std::vector<VMG::MPI::Datatype>& types, const MPI_Comm& comm, const int& tag_start);
+  void IrecvAll(Grid& grid, std::vector<VMG::MPI::Datatype>& types, const MPI_Comm& comm, const int& tag_start);
+
+  void IsendAllBuffered(const Grid& grid, std::vector<VMG::MPI::Datatype>& types, const MPI_Comm& comm, const int& tag_start);
+  void IrecvAllBuffered(std::vector<VMG::MPI::Datatype>& types, const MPI_Comm& comm, const int& tag_start);
+
+  void ReplaceBufferAll(Grid& grid, const std::vector<VMG::MPI::Datatype>& types);
+  void AddBufferAll(Grid& grid, const std::vector<VMG::MPI::Datatype>& types);
+
+  void PrintGridInformation(const Grid& grid, char* filename, const std::string& name);
+  void PrintDatatypes(char* filename);
+
+  VMG::MPI::Settings settings;
+
+  MPI_Comm comm_global;
+  MPI_Info info;
+
+#ifdef VMG_ONE_SIDED
+  bool win_created;
+  MPI_Win win;
+#endif
+
+};
+
+}
+
+#endif /* COMM_MPI_HPP_ */
Index: ThirdParty/vmg/src/comm/comm_serial.cpp
===================================================================
--- ThirdParty/vmg/src/comm/comm_serial.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/comm_serial.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,574 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   comm_serial.cpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Apr 18 12:28:12 2011
+ *
+ * @brief  VMG::CommSerial
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#ifdef HAVE_BOOST_FILESYSTEM_PATH_HPP
+#include <boost/filesystem.hpp>
+namespace fs = boost::filesystem;
+#endif
+
+#include <cstdarg>
+#include <cstdio>
+#include <sstream>
+
+#include "base/discretization.hpp"
+#include "base/helper.hpp"
+#include "base/index.hpp"
+#include "base/stencil.hpp"
+#include "base/vector.hpp"
+#include "comm/comm_serial.hpp"
+#include "discretization/boundary_value_setter.hpp"
+#include "grid/multigrid.hpp"
+#include "grid/tempgrid.hpp"
+#include "mg.hpp"
+
+using namespace VMG;
+
+static char print_buffer[512];
+
+void CommSerial::InitCommSerial()
+{
+  error_norm_count = 0;
+  residual_count = 0;
+}
+
+void CommSerial::CommSubgrid(Grid& grid_old, Grid& grid_new, const int& direction)
+{
+  Grid::iterator iter;
+
+  if (&grid_old == &grid_new)
+    return;
+
+  for (iter = grid_old.Iterators().CompleteGrid().Begin(); iter != grid_old.Iterators().CompleteGrid().End(); ++iter)
+    grid_new(*iter) = grid_old.GetVal(*iter);
+}
+
+void CommSerial::CommAddSubgrid(Grid& grid_old, Grid& grid_new, const int& direction)
+{
+  Grid::iterator iter;
+
+  if (&grid_old == &grid_new)
+    return;
+
+  for (iter = grid_old.Iterators().CompleteGrid().Begin(); iter != grid_old.Iterators().CompleteGrid().End(); ++iter)
+    grid_new(*iter) += grid_old.GetVal(*iter);
+}
+
+std::vector<BoundaryValue> CommSerial::CommBoundaryValues()
+{
+  return MG::GetBoundaryValueSetter()->BoundaryValues();
+}
+
+Grid& CommSerial::GetFinerGrid(Multigrid& multigrid)
+{
+  return multigrid(multigrid.Level()+1);
+}
+
+Grid& CommSerial::GetCoarserGrid(Multigrid& multigrid)
+{
+  return multigrid(multigrid.Level()-1);
+}
+
+void CommSerial::CommFromGhosts(Grid& grid)
+{
+  Grid::iterator iter;
+
+  const GridIteratorSuite& iterators = grid.Iterators();
+  const Index& size = grid.Local().Size();
+
+  if (BoundaryConditions().X() == Periodic) {
+
+    for (iter = iterators.Halo1().X().Begin(); iter != iterators.Halo1().X().End(); ++iter)
+      grid((*iter).X()+size.X(), (*iter).Y(), (*iter).Z()) += grid.GetVal(*iter);
+
+    for (iter = iterators.Halo2().X().Begin(); iter != iterators.Halo2().X().End(); ++iter)
+      grid((*iter).X()-size.X(), (*iter).Y(), (*iter).Z()) += grid.GetVal(*iter);
+  }
+
+  if (BoundaryConditions().Y() == Periodic) {
+
+    for (iter = iterators.Halo1().Y().Begin(); iter != iterators.Halo1().Y().End(); ++iter)
+      grid((*iter).X(), (*iter).Y()+size.Y(), (*iter).Z()) += grid.GetVal(*iter);
+
+    for (iter = iterators.Halo2().Y().Begin(); iter != iterators.Halo2().Y().End(); ++iter)
+      grid((*iter).X(), (*iter).Y()-size.Y(), (*iter).Z()) += grid.GetVal(*iter);
+
+  }
+
+  if (BoundaryConditions().Z() == Periodic) {
+
+    for (iter = iterators.Halo1().Z().Begin(); iter != iterators.Halo1().Z().End(); ++iter)
+      grid((*iter).X(), (*iter).Y(), (*iter).Z()+size.Z()) += grid.GetVal(*iter);
+
+    for (iter = iterators.Halo2().Z().Begin(); iter != iterators.Halo2().Z().End(); ++iter)
+      grid((*iter).X(), (*iter).Y(), (*iter).Z()-size.Z()) += grid.GetVal(*iter);
+
+  }
+}
+
+void CommSerial::CommToGhosts(Grid& grid)
+{
+  Grid::iterator iter;
+
+  const GridIteratorSuite& iterators = grid.Iterators();
+  const Index& size = grid.Local().Size();
+
+  if (BoundaryConditions().Z() == Periodic) {
+
+    for (iter = iterators.Halo1().Z().Begin(); iter != iterators.Halo1().Z().End(); ++iter)
+      grid(*iter) = grid.GetVal((*iter).X(), (*iter).Y(), (*iter).Z()+size.Z());
+
+    for (iter = iterators.Halo2().Z().Begin(); iter != iterators.Halo2().Z().End(); ++iter)
+      grid(*iter) = grid.GetVal((*iter).X(), (*iter).Y(), (*iter).Z()-size.Z());
+
+  }
+
+  if (BoundaryConditions().Y() == Periodic) {
+
+    for (iter = iterators.Halo1().Y().Begin(); iter != iterators.Halo1().Y().End(); ++iter)
+      grid(*iter) = grid.GetVal((*iter).X(), (*iter).Y()+size.Y(), (*iter).Z());
+
+    for (iter = iterators.Halo2().Y().Begin(); iter != iterators.Halo2().Y().End(); ++iter)
+      grid(*iter) = grid.GetVal((*iter).X(), (*iter).Y()-size.Y(), (*iter).Z());
+
+  }
+
+  if (BoundaryConditions().X() == Periodic) {
+
+    for (iter = iterators.Halo1().X().Begin(); iter != iterators.Halo1().X().End(); ++iter)
+      grid(*iter) = grid.GetVal((*iter).X()+size.X(), (*iter).Y(), (*iter).Z());
+
+    for (iter = iterators.Halo2().X().Begin(); iter != iterators.Halo2().X().End(); ++iter)
+      grid(*iter) = grid.GetVal((*iter).X()-size.X(), (*iter).Y(), (*iter).Z());
+
+  }
+}
+
+void CommSerial::CommToGhostsAsyncStart(Grid& grid)
+{
+  CommToGhosts(grid);
+}
+
+void CommSerial::CommToGhostsAsyncFinish(Grid& grid)
+{
+}
+
+void CommSerial::CommFromGhostsAsyncStart(Grid& grid)
+{
+  CommFromGhosts(grid);
+}
+
+void CommSerial::CommFromGhostsAsyncFinish(Grid& grid)
+{
+}
+
+void CommSerial::Print(const OutputLevel level, const char* format, ...)
+{
+  bool print = (level == Output);
+
+#ifdef OUTPUT_INFO
+  print |= (level == Info);
+#endif
+#ifdef OUTPUT_DEBUG
+  print |= (level == Debug);
+#endif
+#ifdef OUTPUT_TIMING
+  print |= (level == Timing);
+#endif
+
+  if (print) {
+    va_list args;
+    va_start(args, format);
+    vsprintf(print_buffer, format, args);
+    printf("VMG: %s\n", print_buffer);
+    va_end(args);
+  }
+}
+
+void CommSerial::PrintOnce(const OutputLevel level, const char* format, ...)
+{
+  bool print = (level == Output);
+
+#ifdef OUTPUT_INFO
+  print |= (level == Info);
+#endif
+#ifdef OUTPUT_DEBUG
+  print |= (level == Debug);
+#endif
+#ifdef OUTPUT_TIMING
+  print |= (level == Timing);
+#endif
+
+  if (print) {
+    va_list args;
+    va_start(args, format);
+    vsprintf(print_buffer, format, args);
+    printf("VMG: %s\n", print_buffer);
+    va_end(args);
+  }
+}
+
+void CommSerial::PrintXML(const std::string& filename, const std::string& xml_data)
+{
+  pugi::xml_document doc;
+  doc.load(xml_data.c_str());
+
+  pugi::xml_node node_global = doc.prepend_child("Global");
+  pugi::xml_node node_num_processes = node_global.append_child("NumProcesses");
+  pugi::xml_node node_num_processes_data = node_num_processes.append_child(pugi::node_pcdata);
+  node_num_processes_data.set_value(Helper::ToString(GlobalProcs().Product()).c_str());
+
+  doc.save_file(filename.c_str());
+}
+
+void CommSerial::PrintXMLAll(const std::string& filename, const std::string& xml_data)
+{
+  PrintXML(filename, xml_data);
+}
+
+void CommSerial::PrintGrid(Grid& grid, const char* information)
+{
+  Index i;
+  std::stringstream out;
+  std::ofstream out_file;
+
+  OpenFileAndPrintHeader(out_file, grid, information);
+
+  for (i.Z()=grid.Local().Begin().Z(); i.Z()<grid.Local().End().Z(); ++i.Z())
+    for (i.Y()=grid.Local().Begin().Y(); i.Y()<grid.Local().End().Y(); ++i.Y())
+      for (i.X()=grid.Local().Begin().X(); i.X()<grid.Local().End().X(); ++i.X())
+        out << std::scientific << grid.GetVal(i) << std::endl;
+
+  out_file << out.str();
+  out_file.close();
+}
+
+void CommSerial::PrintDefect(Grid& sol, Grid& rhs, const char* information)
+{
+  Grid::iterator iter;
+  std::ofstream out_file;
+  std::stringstream out;
+
+  TempGrid *temp = MG::GetTempGrid();
+  temp->SetProperties(sol);
+  temp->ImportFromResidual(sol, rhs);
+
+  OpenFileAndPrintHeader(out_file, *temp, information);
+
+  for (iter=temp->Iterators().Local().Begin(); iter!=temp->Iterators().Local().End(); ++iter)
+    out << std::scientific << temp->GetVal(*iter) << std::endl;
+
+  out_file << out.str();
+  out_file.close();
+}
+
+void CommSerial::OpenFileAndPrintHeader(std::ofstream& out, const Grid& grid, const char* information)
+{
+  char path_str[129];
+  int count = OutputCount();
+
+  sprintf(path_str, "%s%04d.dat", OutputPath().c_str(), count);
+
+  out.open(path_str, std::ios::trunc);
+
+  out << "# vtk DataFile Version 2.0" << std::endl
+      << count << ": " << information << std::endl
+      << "ASCII" << std::endl
+      << "DATASET STRUCTURED_POINTS" << std::endl
+      << grid.Local().Size().X() << " "
+      << grid.Local().Size().Y() << " "
+      << grid.Local().Size().Z() << std::endl
+      << "ORIGIN 0 0 0" << std::endl
+      << "SPACING " << grid.Extent().MeshWidth().X() << " "
+      << grid.Extent().MeshWidth() << " "
+      << grid.Extent().MeshWidth() << std::endl
+      << "POINT_DATA " << grid.Local().Size().Product() << std::endl
+      << "SCALARS residual double 1" << std::endl
+      << "LOOKUP_TABLE default" << std::endl;
+}
+
+void CommSerial::PrintAllSettings()
+{
+  Multigrid* mg = MG::GetFactory().Get("SOL")->Cast<Multigrid>();
+  std::stringstream buf;
+  std::ofstream out;
+
+  buf << OutputPath() << "settings.txt";
+
+  out.open(buf.str().c_str(), std::ios::trunc);
+
+  for (int i=mg->MinLevel(); i<=mg->MaxLevel(); ++i) {
+
+    out << "###########################################################" << std::endl
+        << "LEVEL:                " << i << std::endl
+        << "GLOBAL:" << std::endl
+        << "  LOCAL_BEGIN:         " << (*mg)(i).Global().LocalBegin() << std::endl
+        << "  LOCAL_END:           " << (*mg)(i).Global().LocalEnd() << std::endl
+        << "  LOCAL_SIZE:          " << (*mg)(i).Global().LocalSize() << std::endl
+        << "  GLOBAL_BEGIN:        " << (*mg)(i).Global().GlobalBegin() << std::endl
+        << "  GLOBAL_END:          " << (*mg)(i).Global().GlobalEnd() << std::endl
+        << "  GLOBAL_SIZE:         " << (*mg)(i).Global().GlobalSize() << std::endl
+        << "  GLOBAL_BEGIN_FINEST: " << (*mg)(i).Global().GlobalBeginFinest() << std::endl
+        << "  GLOBAL_END_FINEST:   " << (*mg)(i).Global().GlobalEndFinest() << std::endl
+        << "  GLOBAL_SIZE_FINEST:  " << (*mg)(i).Global().GlobalSizeFinest() << std::endl
+        << "  BOUNDARY_TYPE:       " << (*mg)(i).Global().BoundaryType() << std::endl
+        << "LOCAL:" << std::endl
+        << "  BEGIN:              " << (*mg)(i).Local().Begin() << std::endl
+        << "  END:                " << (*mg)(i).Local().End() << std::endl
+        << "  SIZE:               " << (*mg)(i).Local().Size() << std::endl
+        << "  SIZE_TOTAL:         " << (*mg)(i).Local().SizeTotal() << std::endl
+        << "  HALO_BEGIN_1:       " << (*mg)(i).Local().HaloBegin1() << std::endl
+        << "  HALO_END_1:         " << (*mg)(i).Local().HaloEnd1() << std::endl
+        << "  HALO_SIZE_1:        " << (*mg)(i).Local().HaloSize1() << std::endl
+        << "  HALO_BEGIN_2:       " << (*mg)(i).Local().HaloBegin2() << std::endl
+        << "  HALO_END_2:         " << (*mg)(i).Local().HaloEnd2() << std::endl
+        << "  HALO_SIZE_2:        " << (*mg)(i).Local().HaloSize2() << std::endl
+        << "  BOUNDARY_BEGIN_1:   " << (*mg)(i).Local().BoundaryBegin1() << std::endl
+        << "  BOUNDARY_END_1:     " << (*mg)(i).Local().BoundaryEnd1() << std::endl
+        << "  BOUNDARY_SIZE_1:    " << (*mg)(i).Local().BoundarySize1() << std::endl
+        << "  BOUNDARY_BEGIN_2:   " << (*mg)(i).Local().BoundaryBegin2() << std::endl
+        << "  BOUNDARY_END_2:     " << (*mg)(i).Local().BoundaryEnd2() << std::endl
+        << "  BOUNDARY_SIZE_2:    " << (*mg)(i).Local().BoundarySize2() << std::endl
+        << "EXTENT:" << std::endl
+        << "  BEGIN:              " << (*mg)(i).Extent().Begin() << std::endl
+        << "  END:                " << (*mg)(i).Extent().End() << std::endl
+        << "  SIZE:               " << (*mg)(i).Extent().Size() << std::endl
+        << "  MESH_WIDTH:         " << (*mg)(i).Extent().MeshWidth() << std::endl
+        << "###########################################################" << std::endl;
+
+  }
+
+  assert(out.good());
+  out.close();
+}
+
+inline int GetIndex(const Grid& grid, int i, int j, int k)
+{
+  if (grid.Global().BoundaryType() == LocallyRefined)
+    return k + grid.Local().Size().Z() * (j + grid.Local().Size().Y() * i);
+  else
+    return k + grid.Local().SizeTotal().Z() * (j + grid.Local().SizeTotal().Y() * i);
+}
+
+void CommSerial::PrintGridStructureLevel(Grid& grid, std::ofstream& out)
+{
+  const Vector& sp = grid.Extent().MeshWidth();
+  int numLines;
+
+  if (grid.Global().BoundaryType() == LocallyRefined)
+    numLines = grid.Local().Size().X() * grid.Local().Size().Y() +
+    grid.Local().Size().Y() * grid.Local().Size().Z() +
+    grid.Local().Size().X() * grid.Local().Size().Z();
+  else
+    numLines = grid.Local().SizeTotal().X() * grid.Local().SizeTotal().Y() +
+    grid.Local().SizeTotal().Y() * grid.Local().SizeTotal().Z() +
+    grid.Local().SizeTotal().X() * grid.Local().SizeTotal().Z();
+
+  out << "    <Piece";
+
+  if (grid.Global().BoundaryType() == LocallyRefined) {
+    out << " NumberOfPoints=\"" << grid.Local().Size().Product()  << "\""
+        << " NumberOfVerts=\"" << grid.Local().Size().Product()  << "\"";
+  }else {
+    out << " NumberOfPoints=\"" << grid.Local().SizeTotal().Product()  << "\""
+        << " NumberOfVerts=\"" << grid.Local().SizeTotal().Product()  << "\"";
+  }
+
+  out << " NumberOfLines=\"" << numLines << "\""
+      << " NumberOfStrips=\"0\""
+      << " NumberOfPolys=\"0\"" << ">" << std::endl
+      << "      <Points>" << std::endl
+      << "        <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl
+      << "          ";
+
+  if (grid.Global().BoundaryType() == LocallyRefined) {
+    for (int i=0; i<grid.Local().Size().X(); i++)
+      for (int j=0; j<grid.Local().Size().Y(); j++)
+        for (int k=0; k<grid.Local().Size().Z(); k++)
+          out << grid.Extent().Begin().X() + i * sp.X() << " "
+          << grid.Extent().Begin().Y() + j * sp.Y() << " "
+          << grid.Extent().Begin().Z() + k * sp.Z() << " ";
+  }else {
+    for (int i=0; i<grid.Local().SizeTotal().X(); i++)
+      for (int j=0; j<grid.Local().SizeTotal().Y(); j++)
+        for (int k=0; k<grid.Local().SizeTotal().Z(); k++)
+          out << grid.Extent().Begin().X() + i * sp.X() << " "
+          << grid.Extent().Begin().Y() + j * sp.Y() << " "
+          << grid.Extent().Begin().Z() + k * sp.Z() << " ";
+  }
+
+  out << std::endl
+      << "        </DataArray>" << std::endl
+      << "      </Points>" << std::endl
+      << "      <Verts>" << std::endl
+      << "        <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << std::endl
+      << "          ";
+
+  if (grid.Global().BoundaryType() == LocallyRefined) {
+    for (int i=0; i<grid.Local().Size().Product(); i++)
+      out << i << " ";
+  }else {
+    for (int i=0; i<grid.Local().SizeTotal().Product(); i++)
+      out << i << " ";
+  }
+
+  out << std::endl
+      << "        </DataArray>" << std::endl
+      << "        <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << std::endl
+      << "          ";
+
+  if (grid.Global().BoundaryType() == LocallyRefined) {
+    for (int i=1; i<=grid.Local().Size().Product(); i++)
+      out  << i << " ";
+  }else {
+    for (int i=1; i<=grid.Local().SizeTotal().Product(); i++)
+      out  << i << " ";
+  }
+
+  out << std::endl
+      << "        </DataArray>" << std::endl
+      << "      </Verts>" << std::endl
+      << "      <Lines>" << std::endl
+      << "        <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << std::endl
+      << "          ";
+
+  if (grid.Global().BoundaryType() == LocallyRefined) {
+    for (int i=0; i<grid.Local().Size().X(); i++)
+      for (int j=0; j<grid.Local().Size().Y(); j++)
+        out << GetIndex(grid,i,j,0) << " "
+        << GetIndex(grid,i,j,grid.Local().Size().Z()-1) << " ";
+
+    for (int j=0; j<grid.Local().Size().Y(); j++)
+      for (int k=0; k<grid.Local().Size().Z(); k++)
+        out << GetIndex(grid,0,j,k) << " "
+        << GetIndex(grid,grid.Local().Size().X()-1,j,k) << " ";
+
+    for (int i=0; i<grid.Local().Size().X(); i++)
+      for (int k=0; k<grid.Local().Size().Z(); k++)
+        out << GetIndex(grid,i,0,k) << " "
+        << GetIndex(grid,i,grid.Local().Size().Y()-1,k) << " ";
+  }else {
+    for (int i=0; i<grid.Local().SizeTotal().X(); i++)
+      for (int j=0; j<grid.Local().SizeTotal().Y(); j++)
+        out << GetIndex(grid,i,j,0) << " "
+        << GetIndex(grid,i,j,grid.Local().SizeTotal().Z()-1) << " ";
+
+    for (int j=0; j<grid.Local().SizeTotal().Y(); j++)
+      for (int k=0; k<grid.Local().SizeTotal().Z(); k++)
+        out << GetIndex(grid,0,j,k) << " "
+        << GetIndex(grid,grid.Local().SizeTotal().X()-1,j,k) << " ";
+
+    for (int i=0; i<grid.Local().SizeTotal().X(); i++)
+      for (int k=0; k<grid.Local().SizeTotal().Z(); k++)
+        out << GetIndex(grid,i,0,k) << " "
+        << GetIndex(grid,i,grid.Local().SizeTotal().Y()-1,k) << " ";
+  }
+
+  out << std::endl
+      << "        </DataArray>" << std::endl
+      << "        <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << std::endl
+      << "          ";
+
+  for (int i=1; i<=numLines; i++)
+    out << 2*i << " ";
+
+  out << std::endl
+      << "        </DataArray>" << std::endl
+      << "      </Lines>" << std::endl
+      << "    </Piece>" << std::endl;
+}
+
+void CommSerial::DebugPrintGridStructure(Multigrid& grid)
+{
+  std::ofstream out;
+  char path_str[129];
+
+  sprintf(path_str, "%sgrid.vtp", OutputPath().c_str());
+
+  out.open(path_str, std::ios::trunc);
+
+  if (!out.good()) {
+    Print(Info, "File %s not accessible.", path_str);
+    return;
+  }
+
+  out << "<?xml version=\"1.0\"?>" << std::endl
+      << "<VTKFile type=\"PolyData\">" << std::endl
+      << "  <PolyData>" << std::endl;
+
+  for (int i=grid.MinLevel(); i<=grid.MaxLevel(); i++)
+    PrintGridStructureLevel(grid(i), out);
+
+  out << "  </PolyData>" << std::endl
+      << "</VTKFile>" << std::endl;
+
+  out.close();
+}
+
+std::string CommSerial::CreateOutputDirectory()
+{
+#ifdef HAVE_BOOST_FILESYSTEM_PATH_HPP
+  std::string path, unique_path;
+  std::stringstream unique_suffix;
+  int suffix_counter = 0;
+  char buffer[129];
+  time_t rawtime;
+  struct tm *timeinfo;
+
+  time(&rawtime);
+  timeinfo = localtime(&rawtime);
+  strftime(buffer, 128, "./output/%Y_%m_%d_%H_%M_%S/", timeinfo);
+  path = buffer;
+
+  if (!fs::exists("output"))
+    fs::create_directory("output");
+
+  unique_path = path;
+
+  while (fs::exists(unique_path.c_str())) {
+
+    unique_suffix.str("");
+    unique_suffix << "_" << suffix_counter++ << "/";
+
+    unique_path = path;
+    unique_path.replace(unique_path.size()-1, 1, unique_suffix.str());
+
+  }
+
+  fs::create_directory(unique_path.c_str());
+
+  return unique_path;
+
+#else
+
+  return "./";
+
+#endif
+}
Index: ThirdParty/vmg/src/comm/comm_serial.hpp
===================================================================
--- ThirdParty/vmg/src/comm/comm_serial.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/comm_serial.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,94 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   comm_serial.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Apr 18 12:28:29 2011
+ *
+ * @brief  VMG::CommSerial
+ *
+ */
+
+#ifndef COMM_SERIAL_HPP_
+#define COMM_SERIAL_HPP_
+
+#include <fstream>
+#include <map>
+
+#include "comm/comm.hpp"
+#include "comm/domain_decomposition_serial.hpp"
+
+namespace VMG
+{
+
+class DomainDecomposition;
+class Grid;
+class Multigrid;
+
+class CommSerial : public Comm
+{
+public:
+  CommSerial(const Boundary& boundary, bool register_ = true) :
+    Comm(boundary, new DomainDecompositionSerial(), register_)
+  {
+    InitCommSerial();
+  }
+
+  Grid& GetCoarserGrid(Multigrid& multigrid);
+  Grid& GetFinerGrid(Multigrid& multigrid);
+
+  void CommFromGhosts(Grid& grid);
+  void CommToGhosts(Grid& grid);
+  void CommSubgrid(Grid& grid_old, Grid& grid_new, const int& direction);
+  void CommAddSubgrid(Grid& grid_old, Grid& grid_new, const int& direction);
+  std::vector<BoundaryValue> CommBoundaryValues();
+
+  void CommToGhostsAsyncStart(Grid& grid);
+  void CommToGhostsAsyncFinish(Grid& grid);
+  void CommFromGhostsAsyncStart(Grid& grid);
+  void CommFromGhostsAsyncFinish(Grid& grid);
+
+  void Print(const OutputLevel level, const char* format, ...);
+  void PrintOnce(const OutputLevel level, const char* format, ...);
+  void PrintXML(const std::string& filename, const std::string& xml_data);
+  void PrintXMLAll(const std::string& filename, const std::string& xml_data);
+  void PrintAllSettings();
+  void PrintGrid(Grid& grid, const char* information);
+  void PrintDefect(Grid& sol, Grid& rhs, const char* information);
+
+  void DebugPrintGridStructure(Multigrid& multigrid);
+
+  void PostInit(Multigrid& sol, Multigrid& rhs) {}
+
+private:
+  void InitCommSerial();
+
+  void OpenFileAndPrintHeader(std::ofstream& out, const Grid& mesh, const char* information);
+  void PrintGridStructureLevel(Grid& grid, std::ofstream& out);
+  std::string CreateOutputDirectory();
+
+  int error_norm_count, residual_count;
+
+  std::map< const Grid*, Grid*> finer, coarser;
+
+};
+
+}
+
+#endif /* COMM_SERIAL_HPP_ */
Index: ThirdParty/vmg/src/comm/domain_decomposition.hpp
===================================================================
--- ThirdParty/vmg/src/comm/domain_decomposition.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/domain_decomposition.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,56 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   domain_decomposition.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Jun 27 12:22:57 2011
+ *
+ * @brief  Computes a domain decomposition on a global
+ *         multigrid
+ *
+ */
+
+#ifndef DOMAIN_DECOMPOSITION_HPP_
+#define DOMAIN_DECOMPOSITION_HPP_
+
+#include <map>
+#include <vector>
+
+#include "base/index.hpp"
+#include "grid/grid_properties.hpp"
+
+namespace VMG
+{
+
+class Comm;
+class Interface;
+class Multigrid;
+
+class DomainDecomposition
+{
+public:
+  DomainDecomposition() {}
+  virtual ~DomainDecomposition() {}
+
+  virtual void Compute(Comm& comm, const Interface& interface, std::map<Index, std::vector<GlobalIndices> >& global) = 0;
+};
+
+}
+
+#endif /* DOMAIN_DECOMPOSITION_HPP_ */
Index: ThirdParty/vmg/src/comm/domain_decomposition_mpi.cpp
===================================================================
--- ThirdParty/vmg/src/comm/domain_decomposition_mpi.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/domain_decomposition_mpi.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,160 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   domain_decomposition_mpi.cpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Jun 27 12:53:50 2011
+ *
+ * @brief  Computes a domain decomposition which separates
+ *         the finest grid equally for all processes.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include "base/interface.hpp"
+#include "comm/comm.hpp"
+#include "comm/domain_decomposition_mpi.hpp"
+#include "grid/grid.hpp"
+#include "grid/multigrid.hpp"
+
+using namespace VMG;
+
+void DomainDecompositionMPI::Compute(Comm& comm, const Interface& interface, std::map<Index, std::vector<GlobalIndices> >& global)
+{
+
+  GlobalIndices global_l;
+  Index pos, remainder, procs;
+  Index last_procs = comm.GlobalProcs();
+
+  global.clear();
+
+  for (unsigned int i=0; i<interface.Global().size(); ++i) {
+
+    for (pos.X() = 0; pos.X() < comm.GlobalProcs().X(); ++pos.X())
+      for (pos.Y() = 0; pos.Y() < comm.GlobalProcs().Y(); ++pos.Y())
+        for (pos.Z() = 0; pos.Z() < comm.GlobalProcs().Z(); ++pos.Z()) {
+
+          global_l = interface.Global()[i];
+
+          if (IsActive(global_l.GlobalSize(), pos, procs, comm.GlobalProcs())) {
+
+            if (i == 0) {
+
+              remainder = global_l.GlobalSize() % procs;
+
+              global_l.LocalSize() = global_l.GlobalSize() / procs;
+              for (int j=0; j<3; ++j)
+                if (pos[j] < remainder[j])
+                  ++(global_l.LocalSize()[j]);
+
+              global_l.LocalBegin() = global_l.GlobalBegin() + pos * global_l.LocalSize();
+              for (int j=0; j<3; ++j)
+                if (pos[j] >= remainder[j])
+                  global_l.LocalBegin()[j] += remainder[j];
+
+              global_l.LocalEnd() = global_l.LocalBegin() + global_l.LocalSize();
+
+            } else {
+
+              for (int j=0; j<3; ++j) {
+
+                if (procs[j] == last_procs[j]) {
+
+                  if (global[pos].back().LocalBegin()[j] == global[pos].back().GlobalBegin()[j])
+                    global_l.LocalBegin()[j] = global_l.GlobalBegin()[j];
+                  else
+                    global_l.LocalBegin()[j] = global[pos].back().LocalBegin()[j] / 2;
+
+                  if (global[pos].back().LocalEnd()[j] == global[pos].back().GlobalEnd()[j])
+                    global_l.LocalEnd()[j] = global_l.GlobalEnd()[j];
+                  else
+                    global_l.LocalEnd()[j] = global[pos].back().LocalEnd()[j] / 2;
+
+                  global_l.LocalSize()[j] = global_l.LocalEnd()[j] - global_l.LocalBegin()[j];
+
+                } else {
+
+                  remainder[j] = global_l.GlobalSize()[j] % procs[j];
+
+                  global_l.LocalSize()[j] = global_l.GlobalSize()[j] / procs[j];
+                  if (pos[j] < remainder[j])
+                    ++(global_l.LocalSize()[j]);
+
+                  global_l.LocalBegin()[j] = global_l.GlobalBegin()[j] + pos[j] * global_l.LocalSize()[j];
+                  if (pos[j] >= remainder[j])
+                    global_l.LocalBegin()[j] += remainder[j];
+
+                  global_l.LocalEnd()[j] = global_l.LocalBegin()[j] + global_l.LocalSize()[j];
+
+                }
+              }
+            }
+          }else {
+            global_l.LocalBegin() = 0;
+            global_l.LocalEnd() = 0;
+            global_l.LocalSize() = 0;
+          }
+
+          global[pos].push_back(global_l);
+
+        }
+
+    last_procs = procs;
+
+  }
+}
+
+bool DomainDecompositionMPI::IsActive(const Index& size_global, const Index& pos, Index& procs, const Index& max_procs)
+{
+  bool is_active = true;
+  const int points_min = 5;
+
+  procs = size_global / points_min + 1;
+
+  for (int i=0; i<3; ++i) {
+    procs[i] = std::min(procs[i], max_procs[i]);
+    is_active &= pos[i] < procs[i];
+  }
+
+  return is_active;
+}
+
+void DomainDecompositionMPI::FineToCoarse(Comm& comm, int& begin, int& end, int levels)
+{
+  int last_point = end - 1;
+
+  for (int i=0; i<levels; ++i) {
+
+    if (begin % 2 == 0)
+      begin /= 2;
+    else
+      begin = (begin+1) / 2;
+
+    if (last_point % 2 == 0)
+      last_point /= 2;
+    else
+      last_point = (last_point-1) / 2;
+
+  }
+
+  end = last_point + 1;
+}
Index: ThirdParty/vmg/src/comm/domain_decomposition_mpi.hpp
===================================================================
--- ThirdParty/vmg/src/comm/domain_decomposition_mpi.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/domain_decomposition_mpi.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,51 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   domain_decomposition_mpi.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Jun 27 12:53:50 2011
+ *
+ * @brief  Computes a domain decomposition which separates
+ *         the finest grid equally for all processes.
+ *
+ */
+
+#ifndef DOMAIN_DECOMPOSITION_MPI_HPP_
+#define DOMAIN_DECOMPOSITION_MPI_HPP_
+
+#include "comm/domain_decomposition.hpp"
+
+namespace VMG
+{
+
+class Index;
+
+class DomainDecompositionMPI : public DomainDecomposition
+{
+public:
+  void Compute(Comm& comm, const Interface& interface, std::map<Index, std::vector<GlobalIndices> >& global);
+
+private:
+  bool IsActive(const Index& size_global, const Index& pos, Index& procs, const Index& max_procs);
+  void FineToCoarse(Comm& comm, int& begin, int& end, int levels);
+};
+
+}
+
+#endif /* DOMAIN_DECOMPOSITION_MPI_HPP_ */
Index: ThirdParty/vmg/src/comm/domain_decomposition_serial.cpp
===================================================================
--- ThirdParty/vmg/src/comm/domain_decomposition_serial.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/domain_decomposition_serial.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,43 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   domain_decomposition_serial.cpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Jun 27 12:33:10 2011
+ *
+ * @brief  Domain decomposition which gives all grid points
+ *         to one process
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include "base/interface.hpp"
+#include "comm/domain_decomposition_serial.hpp"
+#include "grid/grid.hpp"
+#include "grid/multigrid.hpp"
+
+using namespace VMG;
+
+void DomainDecompositionSerial::Compute(Comm& comm, const Interface& interface, std::map<Index, std::vector<GlobalIndices> >& global)
+{
+  global[0] = interface.Global();
+}
Index: ThirdParty/vmg/src/comm/domain_decomposition_serial.hpp
===================================================================
--- ThirdParty/vmg/src/comm/domain_decomposition_serial.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/domain_decomposition_serial.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,45 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   domain_decomposition_serial.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Jun 27 12:33:10 2011
+ *
+ * @brief  Domain decomposition which gives all grid points
+ *         to one process
+ *
+ */
+
+#ifndef DOMAIN_DECOMPOSITION_SERIAL_HPP_
+#define DOMAIN_DECOMPOSITION_SERIAL_HPP_
+
+#include "comm/domain_decomposition.hpp"
+
+namespace VMG
+{
+
+class DomainDecompositionSerial : public DomainDecomposition
+{
+public:
+  void Compute(Comm& comm, const Interface& interface, std::map<Index, std::vector<GlobalIndices> >& global);
+};
+
+}
+
+#endif /* DOMAIN_DECOMPOSITION_SERIAL_HPP_ */
Index: ThirdParty/vmg/src/comm/mpi/datatype.cpp
===================================================================
--- ThirdParty/vmg/src/comm/mpi/datatype.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/mpi/datatype.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,219 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   datatype.cpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Tue Apr 09 12:19:58 2013
+ *
+ * @brief  Helper class to handle MPI datatypes.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#ifndef HAVE_MPI
+#error MPI is needed to compile VMG::MPI::Datatype
+#endif
+
+#include <mpi.h>
+#ifdef HAVE_MARMOT
+#include <enhancempicalls.h>
+#include <sourceinfompicalls.h>
+#endif
+
+#include <cstring>
+#include <sstream>
+
+#include "comm/mpi/datatype.hpp"
+
+using namespace VMG;
+
+void VMG::MPI::Datatype::Send(Grid& grid, const int& tag, const MPI_Comm& comm) const
+{
+  if (Feasible())
+    MPI_Send(&grid(0), 1, _type, _rank, _tag_send+tag, comm);
+}
+
+void VMG::MPI::Datatype::Isend(Grid& grid, const int& tag, const MPI_Comm& comm, MPI_Request& request) const
+{
+  if (Feasible())
+    MPI_Isend(&grid(0), 1, _type, _rank, _tag_send+tag, comm, &request);
+}
+
+void VMG::MPI::Datatype::Recv(Grid& grid, const int& tag, const MPI_Comm& comm) const
+{
+  if (Feasible())
+    MPI_Recv(&grid(0), 1 ,_type, _rank, _tag_recv+tag, comm, MPI_STATUS_IGNORE);
+}
+
+void VMG::MPI::Datatype::Irecv(Grid& grid, const int& tag, const MPI_Comm& comm, MPI_Request& request) const
+{
+  if (Feasible())
+    MPI_Irecv(&grid(0), 1, _type, _rank, _tag_recv+tag, comm, &request);
+}
+
+void VMG::MPI::Datatype::SendBuffered(const Grid& grid, const int& tag, const MPI_Comm& comm)
+{
+  if (Feasible()) {
+
+    Index i;
+    int c = 0;
+    const Index end = _starts + _subsizes;
+    const size_t memcpy_size = _subsizes.Z() * sizeof(vmg_float);
+
+    for (i.X()=_starts.X(); i.X()<end.X();  ++i.X())
+      for (i.Y()=_starts.Y(); i.Y()<end.Y(); ++i.Y()) {
+        std::memcpy(&_buffer[c], &grid.GetVal(i.X(), i.Y(), _starts.Z()), memcpy_size);
+        c += _subsizes.Z();
+      }
+
+    MPI_Send(&_buffer.front(), _buffer.size(), MPI_DOUBLE, _rank, _tag_send+tag, comm);
+  }
+}
+
+void VMG::MPI::Datatype::IsendBuffered(const Grid& grid, const int& tag, const MPI_Comm& comm, MPI_Request& request)
+{
+  if (Feasible()) {
+
+    Index i;
+    unsigned int c = 0;
+    const Index end = _starts + _subsizes;
+    const size_t memcpy_size = _subsizes.Z() * sizeof(vmg_float);
+
+    for (i.X()=_starts.X(); i.X()<end.X();  ++i.X())
+      for (i.Y()=_starts.Y(); i.Y()<end.Y(); ++i.Y()) {
+        std::memcpy(&_buffer[c], &grid.GetVal(i.X(), i.Y(), _starts.Z()), memcpy_size);
+        c += _subsizes.Z();
+      }
+
+    assert(c == _buffer.size());
+
+    MPI_Isend(&_buffer.front(), _buffer.size(), MPI_DOUBLE, _rank, _tag_send+tag, comm, &request);
+  }
+}
+
+void VMG::MPI::Datatype::RecvBuffered(const int& tag, const MPI_Comm& comm)
+{
+  if (Feasible())
+    MPI_Recv(&_buffer.front(), _buffer.size(), MPI_DOUBLE, _rank, _tag_recv+tag, comm, MPI_STATUS_IGNORE);
+}
+
+void VMG::MPI::Datatype::IrecvBuffered(const int& tag, const MPI_Comm& comm, MPI_Request& request)
+{
+  if (Feasible())
+    MPI_Irecv(&_buffer.front(), _buffer.size(), MPI_DOUBLE, _rank, _tag_recv+tag, comm, &request);
+}
+
+void VMG::MPI::Datatype::GridReplace(Grid& grid) const
+{
+  if (Feasible()) {
+
+    Index i;
+    unsigned int c = 0;
+    const Index end = _starts + _subsizes;
+    const size_t memcpy_size = _subsizes.Z() * sizeof(vmg_float);
+
+    for (i.X()=_starts.X(); i.X()<end.X(); ++i.X())
+      for (i.Y()=_starts.Y(); i.Y()<end.Y(); ++i.Y()) {
+        std::memcpy(&grid(i.X(), i.Y(), _starts.Z()), &_buffer[c], memcpy_size);
+        c += _subsizes.Z();
+      }
+
+    assert(c == _buffer.size());
+  }
+}
+
+void VMG::MPI::Datatype::GridSum(Grid& grid) const
+{
+  if (Feasible()) {
+
+    Index i;
+    const Index end = _starts + _subsizes;
+    std::vector<vmg_float>::const_iterator iter = _buffer.begin();
+
+    for (i.X()=_starts.X(); i.X()<end.X(); ++i.X())
+      for (i.Y()=_starts.Y(); i.Y()<end.Y(); ++i.Y())
+        for (i.Z()=_starts.Z(); i.Z()<end.Z(); ++i.Z())
+          grid(i) += *iter++;
+
+    assert(iter == _buffer.end());
+  }
+}
+
+void VMG::MPI::Datatype::Set(const GridIteratorSet& bounds, const Grid& grid, const int& rank,
+    const int& tag_send, const int& tag_receive)
+{
+  _sizes = grid.Local().SizeTotal();
+  _subsizes = bounds.Begin().GetEnd() - bounds.Begin().GetBegin();
+  _starts = bounds.Begin().GetBegin();
+  _rank = rank;
+  _tag_send = tag_send;
+  _tag_recv = tag_receive;
+
+  if (_type != MPI_DATATYPE_NULL)
+    MPI_Type_free(&_type);
+
+  InitDatatype();
+}
+
+void VMG::MPI::Datatype::Set(const Index& sizes, const Index& subsizes, const Index& starts, const int& rank,
+    const int& tag_send, const int& tag_receive)
+{
+  _sizes = sizes;
+  _subsizes = subsizes;
+  _starts = starts;
+  _rank = rank;
+  _tag_send = tag_send;
+  _tag_recv = tag_receive;
+
+  if (_type != MPI_DATATYPE_NULL)
+    MPI_Type_free(&_type);
+
+  InitDatatype();
+}
+
+void VMG::MPI::Datatype::InitDatatype()
+{
+  if (Feasible()) {
+    MPI_Type_create_subarray(3, _sizes.vec(), _subsizes.vec(), _starts.vec(), MPI_ORDER_C, MPI_DOUBLE, &_type);
+    MPI_Type_commit(&_type);
+    if (_alloc_buffer)
+      _buffer.resize(_subsizes.Product());
+  }else  {
+    _type = MPI_DATATYPE_NULL;
+  }
+}
+
+std::string VMG::MPI::Datatype::ToString() const
+{
+  std::stringstream str;
+  str << "VMG::MPI::Datatype {" << std::endl
+      << "  Rank:        " << _rank << std::endl
+      << "  TagSend:     " << _tag_send << std::endl
+      << "  TagReceive:  " << _tag_recv << std::endl
+      << "  Sizes:       " << _sizes << std::endl
+      << "  Subsizes:    " << _subsizes << std::endl
+      << "  Starts:      " << _starts << std::endl
+      << "  AllocBuffer: " << _alloc_buffer << std::endl
+      << "}";
+
+  return str.str();
+}
Index: ThirdParty/vmg/src/comm/mpi/datatype.hpp
===================================================================
--- ThirdParty/vmg/src/comm/mpi/datatype.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/mpi/datatype.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,156 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   datatype.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Tue Apr 09 12:19:58 2013
+ *
+ * @brief  Helper class to handle MPI datatypes.
+ *
+ */
+
+#ifndef DATATYPE_HPP_
+#define DATATYPE_HPP_
+
+#include <string>
+#include <vector>
+
+#include "base/index.hpp"
+#include "grid/grid.hpp"
+
+namespace VMG
+{
+
+namespace MPI
+{
+
+class Datatype
+{
+public:
+  Datatype() :
+    _sizes(0),
+    _subsizes(0),
+    _starts(0),
+    _rank(-1),
+    _tag_send(0),
+    _tag_recv(0),
+    _type(MPI_DATATYPE_NULL),
+    _alloc_buffer(false)
+  {}
+
+  Datatype(Index sizes, Index subsizes, Index starts, const int& rank,
+	   const int& tag_send, const int& tag_receive,
+           const bool& alloc_buffer) :
+    _sizes(sizes),
+    _subsizes(subsizes),
+    _starts(starts),
+    _rank(rank),
+    _tag_send(tag_send),
+    _tag_recv(tag_receive),
+    _alloc_buffer(alloc_buffer)
+  {
+    InitDatatype();
+  }
+
+  Datatype(const GridIteratorSet& bounds, const Grid& grid, const int& rank,
+	   const int& tag_send, const int& tag_receive,
+           const bool& alloc_buffer) :
+    _sizes(grid.Local().SizeTotal()),
+    _subsizes(bounds.Begin().GetEnd() - bounds.Begin().GetBegin()),
+    _starts(bounds.Begin().GetBegin()),
+    _rank(rank),
+    _tag_send(tag_send),
+    _tag_recv(tag_receive),
+    _alloc_buffer(alloc_buffer)
+  {
+    InitDatatype();
+  }
+
+  Datatype(const Datatype& other) :
+    _sizes(other._sizes),
+    _subsizes(other._subsizes),
+    _starts(other._starts),
+    _rank(other._rank),
+    _tag_send(other._tag_send),
+    _tag_recv(other._tag_recv),
+    _alloc_buffer(other._alloc_buffer)
+  {
+    InitDatatype();
+  }
+
+  ~Datatype()
+  {
+    if (_type != MPI_DATATYPE_NULL)
+      MPI_Type_free(&_type);
+  }
+
+  void Set(const GridIteratorSet& bounds, const Grid& grid, const int& rank,
+	   const int& tag_send, const int& tag_receive);
+
+  void Set(const Index& sizes, const Index& subsizes, const Index& starts, const int& rank,
+	   const int& tag_send, const int& tag_receive);
+
+  const Index& Sizes() const {return _sizes;}
+  const Index& Subsizes() const {return _subsizes;}
+  const Index& Starts() const {return _starts;}
+
+  const int& Rank() const {return _rank;}
+  const int& TagSend() const {return _tag_send;}
+  const int& TagReceive() const {return _tag_recv;}
+
+  const MPI_Datatype& Type() const {return _type;}
+
+  std::vector<vmg_float>& Buffer() {return _buffer;}
+  const std::vector<vmg_float>& Buffer() const {return _buffer;}
+
+  void Send(Grid& grid, const int& tag, const MPI_Comm& comm) const;
+  void Isend(Grid& grid, const int& tag, const MPI_Comm& comm, MPI_Request& request) const;
+  void Recv(Grid& grid, const int& tag, const MPI_Comm& comm) const;
+  void Irecv(Grid& grid, const int& tag, const MPI_Comm& comm, MPI_Request& request) const;
+
+  void SendBuffered(const Grid& grid, const int& tag, const MPI_Comm& comm);
+  void IsendBuffered(const Grid& grid, const int& tag, const MPI_Comm& comm, MPI_Request& request);
+  void RecvBuffered(const int& tag, const MPI_Comm& comm);
+  void IrecvBuffered(const int& tag, const MPI_Comm& comm, MPI_Request& request);
+
+  void GridReplace(Grid& grid) const;
+  void GridSum(Grid& grid) const;
+
+  bool Feasible() const
+  {
+    return _sizes.Product() > 0 && _subsizes.Product() > 0;
+  }
+
+  std::string ToString() const;
+
+private:
+  void InitDatatype();
+
+  Index _sizes, _subsizes, _starts;
+  int _rank, _tag_send, _tag_recv;
+  MPI_Datatype _type;
+  std::vector<vmg_float> _buffer;
+  bool _alloc_buffer;
+};
+
+}
+
+}
+
+#endif /* DATATYPE_HPP_ */
Index: ThirdParty/vmg/src/comm/mpi/datatypes_global.cpp
===================================================================
--- ThirdParty/vmg/src/comm/mpi/datatypes_global.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/mpi/datatypes_global.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,69 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   datatypes_global.cpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Tue Apr 09 12:19:58 2013
+ *
+ * @brief  Saves MPI datatypes for global grid communication.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include "libvmg_config.h"
+#endif
+
+#ifdef HAVE_MPI
+#include <mpi.h>
+#endif
+
+#include <sstream>
+#include <string>
+
+#include "comm/mpi/datatypes_global.hpp"
+
+using namespace VMG;
+
+std::string VMG::MPI::DatatypesGlobal::ToString() const
+{
+  std::stringstream str;
+  std::vector<Datatype>::const_iterator iter;
+
+  str << "==== DATATYPES_GLOBAL ====" << std::endl
+      << "****      SEND        ****" << std::endl;
+
+  for (iter = send.begin(); iter != send.end(); ++iter)
+    str << iter->ToString() << std::endl;
+
+  str << "**************************" << std::endl
+      << "****     RECEIVE      ****" << std::endl;
+
+  for (iter = recv.begin(); iter != recv.end(); ++iter)
+    str << iter->ToString() << std::endl;
+
+  str << "**************************" << std::endl
+      << "==========================";
+
+  return str.str();
+}
+
+std::ostream& VMG::MPI::operator<<(std::ostream& out, const VMG::MPI::DatatypesGlobal& dt_global)
+{
+  return out << dt_global.ToString();
+}
Index: ThirdParty/vmg/src/comm/mpi/datatypes_global.hpp
===================================================================
--- ThirdParty/vmg/src/comm/mpi/datatypes_global.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/mpi/datatypes_global.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,70 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   datatypes_global.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Jan 2 18:45:22 2012
+ *
+ * @brief  Saves MPI datatypes for global grid communication.
+ *
+ */
+
+#ifndef DATATYPES_GLOBAL_HPP_
+#define DATATYPES_GLOBAL_HPP_
+
+#include <vector>
+
+#include "comm/mpi/datatype.hpp"
+
+namespace VMG
+{
+
+namespace MPI
+{
+
+class DatatypesGlobal
+{
+public:
+  DatatypesGlobal() {}
+  ~DatatypesGlobal() {}
+
+  typedef std::vector<Datatype>::iterator iterator;
+  typedef std::vector<Datatype>::const_iterator const_iterator;
+  typedef std::vector<Datatype>::reverse_iterator reverse_iterator;
+  typedef std::vector<Datatype>::const_reverse_iterator const_reverse_iterator;
+
+  std::vector<Datatype>& Send() {return send;}
+  const std::vector<Datatype>& Send() const {return send;}
+
+  std::vector<Datatype>& Receive() {return recv;}
+  const std::vector<Datatype>& Receive() const {return recv;}
+
+  std::string ToString() const;
+
+private:
+  std::vector<Datatype> send, recv;
+};
+
+std::ostream& operator<<(std::ostream& out, const DatatypesGlobal& dt_global);
+
+}
+
+}
+
+#endif /* DATATYPES_GLOBAL_HPP_ */
Index: ThirdParty/vmg/src/comm/mpi/datatypes_local.cpp
===================================================================
--- ThirdParty/vmg/src/comm/mpi/datatypes_local.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/mpi/datatypes_local.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,76 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   datatypes_local.cpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Jan 2 18:45:22 2012
+ *
+ * @brief  Saves MPI datatypes for ghost cell communication.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include "libvmg_config.h"
+#endif
+
+#ifdef HAVE_MPI
+#include <mpi.h>
+#endif
+
+#include <sstream>
+#include <string>
+
+#include "comm/mpi/datatypes_local.hpp"
+
+using namespace VMG;
+
+std::string VMG::MPI::DatatypesLocal::ToString() const
+{
+  std::stringstream str;
+  std::vector<Datatype>::const_iterator iter_dt;
+  std::vector<Index>::const_iterator iter_ind;
+
+  str << "==== DATATYPES_LOCAL  ====" << std::endl
+      << "****      HALO        ****" << std::endl;
+
+  for (iter_dt =_halo.begin(); iter_dt != _halo.end(); ++iter_dt)
+    str << iter_dt->ToString() << std::endl;
+
+  str << "**************************" << std::endl
+      << "****  NEAR_BOUNDARY   ****" << std::endl;
+
+  for (iter_dt = _nb.begin(); iter_dt != _nb.end(); ++iter_dt)
+    str << iter_dt->ToString() << std::endl;
+
+  str << "**************************" << std::endl
+      << "****     OFFSET       ****" << std::endl;
+
+  for (iter_ind = _offset.begin(); iter_ind != _offset.end(); ++iter_ind)
+    str << *iter_ind << std::endl;
+
+  str << "**************************" << std::endl
+      << "==========================";
+
+  return str.str();
+}
+
+std::ostream& VMG::MPI::operator<<(std::ostream& out, const VMG::MPI::DatatypesLocal& dt_local)
+{
+  return out << dt_local.ToString();
+}
Index: ThirdParty/vmg/src/comm/mpi/datatypes_local.hpp
===================================================================
--- ThirdParty/vmg/src/comm/mpi/datatypes_local.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/mpi/datatypes_local.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,435 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   datatypes_local.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Jan 2 18:45:22 2012
+ *
+ * @brief  Saves MPI datatypes for ghost cell communication.
+ *
+ */
+
+#ifndef DATATYPES_LOCAL_HPP_
+#define DATATYPES_LOCAL_HPP_
+
+#ifdef HAVE_MPI
+#include <mpi.h>
+
+#include <vector>
+
+#include "base/index.hpp"
+#include "comm/mpi/datatype.hpp"
+
+namespace VMG
+{
+
+namespace MPI
+{
+
+class DatatypesLocal
+{
+public:
+  template <class T>
+  DatatypesLocal(const T& grid, const MPI_Comm& comm, const bool& alloc)
+  {
+  InitDatatypesLocal(grid, comm, alloc);
+  }
+
+  ~DatatypesLocal() {}
+
+  std::vector<Datatype>& Halo() {return _halo;}
+  std::vector<Datatype>& NB() {return _nb;}
+
+  const std::vector<Datatype>& Halo() const {return _halo;}
+  const std::vector<Datatype>& NB() const {return _nb;}
+
+  const std::vector<Index>& Offset() const {return _offset;}
+
+  std::string ToString() const;
+
+private:
+  template <class T>
+  void InitDatatypesLocal(const T& grid, const MPI_Comm& comm, const bool& alloc_buffer);
+
+  std::vector<Datatype> _halo, _nb;
+  std::vector<Index> _offset;
+};
+
+std::ostream& operator<<(std::ostream& out, const DatatypesLocal& dt_local);
+
+namespace
+{
+
+inline int to_1d(const Index& i)
+{
+  return i.Z()+3*(i.Y()+3*i.X());
+}
+
+inline bool _is_valid(const Index& coord, const Index& dims, const Index& periods)
+{
+  return (periods[0] || (coord[0] >= 0 && coord[0] < dims[0])) &&
+      (periods[1] || (coord[1] >= 0 && coord[1] < dims[1])) &&
+      (periods[2] || (coord[2] >= 0 && coord[2] < dims[2]));
+}
+
+}
+
+template <class T>
+void DatatypesLocal::InitDatatypesLocal(const T& grid, const MPI_Comm& comm, const bool& alloc_buffer)
+{
+  if (comm != MPI_COMM_NULL) {
+
+    int index, ranks[27];
+    Index dims, periods, coords;
+    Index sizes, subsizes, starts;
+    Index offset, i;
+    const LocalIndices& l = grid.Local();
+
+    MPI_Cart_get(comm, 3, dims.vec(), periods.vec(), coords.vec());
+
+    for (i.X()=-1; i.X()<=1; ++i.X())
+      for (i.Y()=-1; i.Y()<=1; ++i.Y())
+        for (i.Z()=-1; i.Z()<=1; ++i.Z())
+          if (_is_valid(coords + i, dims, periods))
+            MPI_Cart_rank(comm, (coords + i).vec(), &ranks[to_1d(i+1)]);
+
+    sizes = l.SizeTotal();
+
+    /* -1 0 0 */
+    offset = Index(-1,0,0);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.HaloSize1().X(), l.Size().Y(), l.Size().Z());
+      starts = Index(0, l.Begin().Y(), l.Begin().Z());
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 0, 1, alloc_buffer));
+      starts = l.Begin();
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 0, 1, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* 1 0 0 */
+    offset = Index(1,0,0);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.HaloSize2().X(), l.Size().Y(), l.Size().Z());
+      starts = Index(l.End().X(), l.Begin().Y(), l.Begin().Z());
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 1, 0, alloc_buffer));
+      starts = Index(l.End().X()-l.HaloSize2().X(), l.Begin().Y(), l.Begin().Z());
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 1, 0, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* 0 -1 0 */
+    offset = Index(0,-1,0);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.Size().X(), l.HaloSize1().Y(), l.Size().Z());
+      starts = Index(l.Begin().X(), 0, l.Begin().Z());
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 2, 3, alloc_buffer));
+      starts = l.Begin();
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 2, 3, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* 0 1 0 */
+    offset = Index(0,1,0);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.Size().X(), l.HaloSize2().Y(), l.Size().Z());
+      starts = Index(l.Begin().X(), l.End().Y(), l.Begin().Z());
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 3, 2, alloc_buffer));
+      starts = Index(l.Begin().X(), l.End().Y()-l.HaloSize2().Y(), l.Begin().Z());
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 3, 2, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* 0 0 -1 */
+    offset = Index(0,0,-1);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.Size().X(), l.Size().Y(), l.HaloSize1().Z());
+      starts = Index(l.Begin().X(), l.Begin().Y(), 0);
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 4, 5, alloc_buffer));
+      starts = l.Begin();
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 4, 5, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* 0 0 1 */
+    offset = Index(0,0,1);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.Size().X(), l.Size().Y(), l.HaloSize2().Z());
+      starts = Index(l.Begin().X(), l.Begin().Y(), l.End().Z());
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 5, 4, alloc_buffer));
+      starts = Index(l.Begin().X(), l.Begin().Y(), l.End().Z()-l.HaloSize2().Z());
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 5, 4, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* -1 -1 0 */
+    offset = Index(-1,-1,0);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.HaloSize1().X(), l.HaloSize1().Y(), l.Size().Z());
+      starts = Index(0, 0, l.Begin().Z());
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 6, 7, alloc_buffer));
+      starts = l.Begin();
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 6, 7, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* -1 1 0 */
+    offset = Index(-1,1,0);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.HaloSize1().X(), l.HaloSize2().Y(), l.Size().Z());
+      starts = Index(0, l.End().Y(), l.Begin().Z());
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 8, 9, alloc_buffer));
+      starts = Index(l.Begin().X(), l.End().Y()-l.HaloSize2().Y(), l.Begin().Z());
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 8, 9, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* 1 -1 0 */
+    offset = Index(1,-1,0);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.HaloSize2().X(), l.HaloSize1().Y(), l.Size().Z());
+      starts = Index(l.End().X(), 0, l.Begin().Z());
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 9, 8, alloc_buffer));
+      starts = Index(l.End().X()-l.HaloSize2().X(), l.Begin().Y(), l.Begin().Z());
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 9, 8, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* 1 1 0 */
+    offset = Index(1,1,0);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.HaloSize2().X(), l.HaloSize2().Y(), l.Size().Z());
+      starts = Index(l.End().X(), l.End().Y(), l.Begin().Z());
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 7, 6, alloc_buffer));
+      starts = Index(l.End().X()-l.HaloSize2().X(), l.End().Y()-l.HaloSize2().Y(), l.Begin().Z());
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 7, 6, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* -1 0 -1 */
+    offset = Index(-1,0,-1);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.HaloSize1().X(), l.Size().Y(), l.HaloSize1().Z());
+      starts = Index(0, l.Begin().Y(), 0);
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 10, 11, alloc_buffer));
+      starts = l.Begin();
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 10, 11, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* -1 0 1 */
+    offset = Index(-1,0,1);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.HaloSize1().X(), l.Size().Y(), l.HaloSize2().Z());
+      starts = Index(0, l.Begin().Y(), l.End().Z());
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 12, 13, alloc_buffer));
+      starts = Index(l.Begin().X(), l.Begin().Y(), l.End().Z()-l.HaloSize2().Z());
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 12, 13, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* 1 0 -1 */
+    offset = Index(1,0,-1);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.HaloSize2().X(), l.Size().Y(), l.HaloSize1().Z());
+      starts = Index(l.End().X(), l.Begin().Y(), 0);
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 13, 12, alloc_buffer));
+      starts = Index(l.End().X()-l.HaloSize2().X(), l.Begin().Y(), l.Begin().Z());
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 13, 12, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* 1 0 1 */
+    offset = Index(1,0,1);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.HaloSize2().X(), l.Size().Y(), l.HaloSize2().Z());
+      starts = Index(l.End().X(), l.Begin().Y(), l.End().Z());
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 11, 10, alloc_buffer));
+      starts = Index(l.End().X()-l.HaloSize2().X(), l.Begin().Y(), l.End().Z()-l.HaloSize2().Z());
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 11, 10, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* 0 -1 -1 */
+    offset = Index(0,-1,-1);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.Size().X(), l.HaloSize1().Y(), l.HaloSize1().Z());
+      starts = Index(l.Begin().X(), 0, 0);
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 14, 15, alloc_buffer));
+      starts = l.Begin();
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 14, 15, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* 0 -1 1 */
+    offset = Index(0,-1,1);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.Size().X(), l.HaloSize1().Y(), l.HaloSize2().Z());
+      starts = Index(l.Begin().X(), 0, l.End().Z());
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 16, 17, alloc_buffer));
+      starts = Index(l.Begin().X(), l.Begin().Y(), l.End().Z()-l.HaloSize2().Z());
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 16, 17, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* 0 1 -1 */
+    offset = Index(0,1,-1);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.Size().X(), l.HaloSize2().Y(), l.HaloSize1().Z());
+      starts = Index(l.Begin().X(), l.End().Y(), 0);
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 17, 16, alloc_buffer));
+      starts = Index(l.Begin().X(), l.End().Y()-l.HaloSize2().Y(), l.Begin().Z());
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 17, 16, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* 0 1 1 */
+    offset = Index(0,1,1);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.Size().X(), l.HaloSize2().Y(), l.HaloSize2().Z());
+      starts = Index(l.Begin().X(), l.End().Y(), l.End().Z());
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 15, 14, alloc_buffer));
+      starts = Index(l.Begin().X(), l.End().Y()-l.HaloSize2().Y(), l.End().Z()-l.HaloSize2().Z());
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 15, 14, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* -1 -1 -1 */
+    offset = Index(-1,-1,-1);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = l.HaloSize1();
+      starts = 0;
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 18, 19, alloc_buffer));
+      starts = l.Begin();
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 18, 19, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* -1 -1 1 */
+    offset = Index(-1,-1,1);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.HaloSize1().X(), l.HaloSize1().Y(), l.HaloSize2().Z());
+      starts = Index(0, 0, l.End().Z());;
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 20, 21, alloc_buffer));
+      starts = Index(l.Begin().X(), l.Begin().Y(), l.End().Z()-l.HaloSize2().Z());
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 20, 21, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* -1 1 -1 */
+    offset = Index(-1,1,-1);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.HaloSize1().X(), l.HaloSize2().Y(), l.HaloSize1().Z());
+      starts = Index(0, l.End().Y(), 0);
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 22, 23, alloc_buffer));
+      starts = Index(l.Begin().X(), l.End().Y()-l.HaloSize2().Y(), l.Begin().Z());
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 22, 23, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* 1 -1 -1 */
+    offset = Index(1,-1,-1);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.HaloSize2().X(), l.HaloSize1().Y(), l.HaloSize1().Z());
+      starts = Index(l.End().X(), 0, 0);
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 24, 25, alloc_buffer));
+      starts = Index(l.End().X()-l.HaloSize2().X(), l.Begin().Y(), l.Begin().Z());
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 24, 25, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* -1 1 1 */
+    offset = Index(-1,1,1);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.HaloSize1().X(), l.HaloSize2().Y(), l.HaloSize2().Z());
+      starts = Index(0, l.End().Y(), l.End().Z());
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 25, 24, alloc_buffer));
+      starts = Index(l.Begin().X(), l.End().Y()-l.HaloSize2().Y(), l.End().Z()-l.HaloSize2().Z());
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 25, 24, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* 1 -1 1 */
+    offset = Index(1,-1,1);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.HaloSize2().X(), l.HaloSize1().Y(), l.HaloSize2().Z());
+      starts = Index(l.End().X(), 0, l.End().Z());
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 23, 22, alloc_buffer));
+      starts = Index(l.End().X()-l.HaloSize2().X(), l.Begin().Y(), l.End().Z()-l.HaloSize2().Z());
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 23, 22, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* 1 1 -1 */
+    offset = Index(1,1,-1);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = Index(l.HaloSize2().X(), l.HaloSize2().Y(), l.HaloSize1().Z());
+      starts = Index(l.End().X(), l.End().Y(), 0);
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 21, 20, alloc_buffer));
+      starts = Index(l.End().X()-l.HaloSize2().X(), l.End().Y()-l.HaloSize2().Y(), l.Begin().Z());
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 21, 20, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+    /* 1 1 1 */
+    offset = Index(1,1,1);
+    if (_is_valid(coords + offset, dims, periods)) {
+      index = to_1d(offset+1);
+      subsizes = l.HaloSize2();
+      starts = l.End();
+      _halo.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 19, 18, alloc_buffer));
+      starts = l.End()-l.HaloSize2();
+      _nb.push_back(VMG::MPI::Datatype(sizes, subsizes, starts, ranks[index], 19, 18, alloc_buffer));
+      _offset.push_back(offset);
+    }
+
+  }
+}
+
+}
+
+}
+
+#endif /* HAVE_MPI */
+
+#endif /* DATATYPES_LOCAL_HPP_ */
Index: ThirdParty/vmg/src/comm/mpi/error_handler.cpp
===================================================================
--- ThirdParty/vmg/src/comm/mpi/error_handler.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/mpi/error_handler.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,66 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   error_handler.cpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Nov 21 13:27:22 2011
+ *
+ * @brief  Convert MPI errors to C++ exceptions. Used to be
+ *         able to call a debugger when MPI crashes
+ *         internally.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#ifdef DEBUG
+
+#ifdef HAVE_MPI
+
+#include <mpi.h>
+#ifdef HAVE_MARMOT
+#include <enhancempicalls.h>
+#include <sourceinfompicalls.h>
+#endif
+
+#include <sstream>
+
+#include "comm/mpi/error_handler.hpp"
+
+void VMG::MPI::ConvertToException(MPI_Comm* comm, int* err, ...)
+{
+  int rank;
+  char error_string[MPI_MAX_ERROR_STRING+1], comm_name[MPI_MAX_OBJECT_NAME+1];
+  int error_len, comm_name_len;
+
+  MPI_Comm_rank(*comm, &rank);
+  MPI_Error_string(*err, error_string, &error_len);
+  MPI_Comm_get_name(*comm, comm_name, &comm_name_len);
+
+  std::stringstream error_concat;
+  error_concat << "vmg: Rank " << rank << " on communicator \"" << comm_name << "\": " << error_string;
+
+  throw Exception(error_concat.str().c_str());
+}
+
+#endif /* DEBUG */
+
+#endif /* HAVE_MPI */
Index: ThirdParty/vmg/src/comm/mpi/error_handler.hpp
===================================================================
--- ThirdParty/vmg/src/comm/mpi/error_handler.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/mpi/error_handler.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,76 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   error_handler.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Nov 21 13:27:22 2011
+ *
+ * @brief  Convert MPI errors to C++ exceptions. Used to be
+ *         able to call a debugger when MPI crashes
+ *         internally.
+ *
+ */
+
+#ifndef ERROR_HANDLER_HPP_
+#define ERROR_HANDLER_HPP_
+
+#ifdef DEBUG
+
+#ifdef HAVE_MPI
+
+#include <stdexcept>
+
+namespace VMG
+{
+
+namespace MPI
+{
+
+class Exception : public std::exception
+{
+public:
+  Exception(const std::string& what) :
+    std::exception(),
+    m_what(what)
+  {}
+
+  virtual ~Exception() throw()
+  {}
+
+  const char* what() const throw()
+  {
+    return m_what.c_str();
+  }
+
+protected:
+  std::string m_what;
+
+};
+
+  void ConvertToException(MPI_Comm* comm, int* err, ...);
+
+}
+
+}
+
+#endif /* HAVE_MPI */
+
+#endif /* DEBUG */
+
+#endif /* ERROR_HANDLER_HPP_ */
Index: ThirdParty/vmg/src/comm/mpi/has_request_vec.hpp
===================================================================
--- ThirdParty/vmg/src/comm/mpi/has_request_vec.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/mpi/has_request_vec.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,83 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   has_request_vec.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Nov 21 13:27:22 2011
+ *
+ * @brief  A convenience base class for classes using MPI_Requests.
+ *
+ */
+
+
+#ifndef HAS_REQUEST_VEC_HPP_
+#define HAS_REQUEST_VEC_HPP_
+
+#include <vector>
+
+#ifdef HAVE_MPI
+
+namespace VMG
+{
+
+class HasRequestVec
+{
+protected:
+
+MPI_Request& Request()
+{
+  request_vec.push_back(MPI_Request());
+  return request_vec.back();
+}
+
+void WaitAll()
+{
+  if (!request_vec.empty()) {
+#ifndef NDEBUG
+    int rval = MPI_Waitall(static_cast<int>(request_vec.size()), &request_vec.front(), MPI_STATUSES_IGNORE);
+    assert(rval == MPI_SUCCESS);
+#else
+    MPI_Waitall(static_cast<int>(request_vec.size()), &request_vec.front(), MPI_STATUSES_IGNORE);
+#endif
+    request_vec.clear();
+  }
+}
+
+int TestAll()
+{
+  int flag = 1;
+  if (!request_vec.empty())
+    MPI_Testall(request_vec.size(), &request_vec.front(), &flag, MPI_STATUSES_IGNORE);
+  return flag;
+}
+
+size_t RequestsPending()
+{
+  return request_vec.size();
+}
+
+private:
+  std::vector<MPI_Request> request_vec;
+};
+
+}
+
+#endif /* HAVE_MPI */
+
+#endif /* HAS_REQUEST_VEC_HPP_ */
Index: ThirdParty/vmg/src/comm/mpi/key.hpp
===================================================================
--- ThirdParty/vmg/src/comm/mpi/key.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/mpi/key.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,231 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   key.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Nov 21 13:27:22 2011
+ *
+ * @brief  Class to distinguish the different entities in maps.
+ *
+ */
+
+#ifndef KEY_HPP_
+#define KEY_HPP_
+
+#ifndef HAVE_MPI
+#error MPI is needed to compile this class
+#endif /* HAVE_MPI */
+
+#include <iostream>
+#include <list>
+
+#include "grid/grid.hpp"
+
+namespace VMG
+{
+
+class Grid;
+class Index;
+
+namespace MPI
+{
+
+class KeyStorage
+{
+public:
+  KeyStorage(const Grid& grid) :
+    begin(grid.Global().LocalBegin()),
+    end(grid.Global().LocalEnd()),
+    size_local(grid.Local().SizeTotal()),
+    size_global(grid.Global().GlobalSizeFinest()),
+    level(grid.Level())
+  {}
+
+  KeyStorage(const Index& begin, const Index& end,
+	     const Index& size_local, const Index& size_global,
+	     const int& level) :
+    begin(begin),
+    end(end),
+    size_local(size_local),
+    size_global(size_global),
+    level(level)
+  {}
+
+  KeyStorage(const KeyStorage& other) :
+    begin(other.begin),
+    end(other.end),
+    size_local(other.size_local),
+    size_global(other.size_global),
+    level(other.level)
+  {}
+
+  ~KeyStorage() {}
+
+  bool operator==(const KeyStorage& other) const
+  {
+    return this->begin == other.begin &&
+           this->end == other.end &&
+           this->size_local == other.size_local &&
+           this->size_global == other.size_global &&
+           this->level == other.level;
+  }
+
+  bool operator!=(const KeyStorage& other) const
+  {
+    return !(*this==other);
+  }
+
+  bool operator<(const KeyStorage& other) const
+  {
+    if (this->begin < other.begin) return true;
+    if (this->begin != other.begin) return false;
+    if (this->end < other.end) return true;
+    if (this->end != other.end) return false;
+    if (this->size_local < other.size_local) return true;
+    if (this->size_local != other.size_local) return false;
+    if (this->size_global < other.size_global) return true;
+    if (this->size_global != other.size_global) return false;
+    if (this->level < other.level) return true;
+    return false;
+  }
+
+private:
+  const Index begin, end, size_local, size_global;
+  const int level;
+};
+
+class KeySorted {
+public:
+  KeySorted(const KeySorted& other)
+  {
+    std::list<KeyStorage>::const_iterator i;
+    this->keys.clear();
+    for (i=other.keys.begin(); i!=other.keys.end(); ++i)
+      this->keys.push_back(*i);
+  }
+
+  KeySorted(const Grid& grid)
+  {
+    keys.push_back(KeyStorage(grid));
+  }
+
+  KeySorted(const Grid& grid_1, const Grid& grid_2)
+  {
+    keys.push_back(KeyStorage(grid_1));
+    keys.push_back(KeyStorage(grid_2));
+    keys.sort();
+  }
+
+  KeySorted(const Index& begin, const Index& end,
+	    const Index& size_local, const Index& size_global,
+	    const int& level)
+  {
+    keys.push_back(KeyStorage(begin, end, size_local, size_global, level));
+  }
+
+  ~KeySorted() {}
+
+  bool operator<(const KeySorted& other) const
+  {
+    if (this->keys.size() < other.keys.size()) return true;
+    if (this->keys.size() > other.keys.size()) return false;
+
+    std::list<KeyStorage>::const_iterator i1, i2;
+    for (i1=this->keys.begin(), i2=other.keys.begin();
+	 i1!=this->keys.end(), i2!=other.keys.end();
+	 ++i1, ++i2) {
+      if (*i1 < *i2) return true;
+      if (*i1 != *i2) return false;
+    }
+
+    return false;
+  }
+
+private:
+  std::list<KeyStorage> keys;
+};
+
+
+  /*
+   * direction: 0 - from multigrid to temporary grid
+   *            1 - from temporary grid to multigrid
+   *            for single grid datatypes always 0
+   */
+class KeyUnsorted {
+public:
+  KeyUnsorted(const KeyUnsorted& other)
+  {
+    std::list<KeyStorage>::const_iterator i;
+    this->keys.clear();
+    for (i=other.keys.begin(); i!=other.keys.end(); ++i)
+      this->keys.push_back(*i);
+    this->direction = other.direction;
+  }
+
+  KeyUnsorted(const Grid& grid, const int& direction)
+  {
+    keys.push_back(KeyStorage(grid));
+    this->direction = direction;
+  }
+
+  KeyUnsorted(const Grid& grid_1, const Grid& grid_2, const int& direction)
+  {
+    keys.push_back(KeyStorage(grid_1));
+    keys.push_back(KeyStorage(grid_2));
+    this->direction = direction;
+  }
+
+  KeyUnsorted(const Index& begin, const Index& end,
+	      const Index& size_local, const Index& size_global,
+	      const int& level, const int& direction)
+  {
+    keys.push_back(KeyStorage(begin, end, size_local, size_global, level));
+    this->direction = direction;
+  }
+
+  ~KeyUnsorted() {}
+
+  bool operator<(const KeyUnsorted& other) const
+  {
+    if (this->keys.size() < other.keys.size()) return true;
+    if (this->keys.size() > other.keys.size()) return false;
+
+    std::list<KeyStorage>::const_iterator i1, i2;
+    for (i1=this->keys.begin(), i2=other.keys.begin();
+	 i1!=this->keys.end(), i2!=other.keys.end();
+	 ++i1, ++i2) {
+      if (*i1 < *i2) return true;
+      if (*i1 != *i2) return false;
+    }
+
+    if (this->direction < other.direction) return true;
+
+    return false;
+  }
+
+private:
+  std::list<KeyStorage> keys;
+  int direction;
+};
+
+}
+
+}
+
+#endif /* KEY_HPP_ */
Index: ThirdParty/vmg/src/comm/mpi/settings.cpp
===================================================================
--- ThirdParty/vmg/src/comm/mpi/settings.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/mpi/settings.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,525 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   settings.cpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Jan 2 18:45:22 2012
+ *
+ * @brief  Computes several MPI-related settings.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#ifdef HAVE_MPI
+#include <mpi.h>
+#ifdef HAVE_MARMOT
+#include <enhancempicalls.h>
+#include <sourceinfompicalls.h>
+#endif
+#else
+#error MPI is needed to compile CommMPISettings
+#endif
+
+#include <cassert>
+#include <sstream>
+#include <string>
+
+#include "base/index.hpp"
+#include "comm/comm.hpp"
+#include "comm/mpi/settings.hpp"
+#include "discretization/boundary_value_setter.hpp"
+#include "grid/multigrid.hpp"
+#include "grid/tempgrid.hpp"
+#include "mg.hpp"
+
+using namespace VMG;
+
+Grid& VMG::MPI::Settings::FinerGrid(const Grid& grid)
+{
+  assert(finer_grids.find(&grid) != finer_grids.end());
+  return *finer_grids.find(&grid)->second;
+}
+
+Grid& VMG::MPI::Settings::CoarserGrid(const Grid& grid)
+{
+  assert(coarser_grids.find(&grid) != coarser_grids.end());
+  return *coarser_grids.find(&grid)->second;
+}
+
+MPI_Comm VMG::MPI::Settings::CommunicatorGlobal(const Grid& grid) const
+{
+  assert(communicators_global.find(grid.Level()) != communicators_global.end());
+  return communicators_global.find(grid.Level())->second;
+}
+
+MPI_Comm VMG::MPI::Settings::CommunicatorLocal(const Grid& grid) const
+{
+  KeyUnsorted k(grid, 0);
+  assert(communicators_local.find(k) != communicators_local.end());
+  return communicators_local.find(k)->second;
+}
+
+VMG::MPI::DatatypesGlobal& VMG::MPI::Settings::DatatypesGlobal(const Grid& grid_old, const Grid& grid_new, const int& direction)
+{
+  KeyUnsorted k(grid_old, grid_new, direction);
+  assert(datatypes_global.find(k) != datatypes_global.end());
+  return datatypes_global.find(k)->second;
+}
+
+VMG::MPI::DatatypesLocal& VMG::MPI::Settings::DatatypesLocal(const Grid& grid)
+{
+  KeyUnsorted k(grid, 0);
+  assert(datatypes_local.find(k) != datatypes_local.end());
+  return datatypes_local.find(k)->second;
+}
+
+void VMG::MPI::Settings::ComputeSettings(Multigrid& sol, Multigrid& rhs, MPI_Comm& comm_global)
+{
+
+  std::map<const Grid*, Grid*>::const_iterator grid_iter;
+
+  Comm& comm = *MG::GetComm();
+
+  /*
+   * Create coarser grids
+   */
+  for (int i=sol.MaxLevel(); i>sol.MinLevel(); --i) {
+
+    TempGrid* temp_grid = new TempGrid();
+    temp_grid->SetPropertiesToCoarser(sol(i), comm.BoundaryConditions());
+
+    if (temp_grid->Global().LocalBegin().IsComponentwiseGreaterOrEqual(sol(i-1).Global().LocalBegin()) &&
+        temp_grid->Global().LocalEnd().IsComponentwiseLessOrEqual(sol(i-1).Global().LocalEnd())) {
+      delete temp_grid;
+      coarser_grids.insert(std::make_pair(&sol(i), &sol(i-1)));
+    }else {
+      coarser_grids.insert(std::make_pair(&sol(i), temp_grid));
+    }
+
+  }
+
+  //FIXME
+  for (int i=rhs.MaxLevel(); i>rhs.MinLevel(); --i) {
+
+    TempGrid* temp_grid = new TempGrid();
+    temp_grid->SetPropertiesToCoarser(rhs(i), comm.BoundaryConditions());
+
+    if (temp_grid->Global().LocalBegin().IsComponentwiseGreaterOrEqual(rhs(i-1).Global().LocalBegin()) &&
+        temp_grid->Global().LocalEnd().IsComponentwiseLessOrEqual(rhs(i-1).Global().LocalEnd())) {
+      delete temp_grid;
+      coarser_grids.insert(std::make_pair(&rhs(i), &rhs(i-1)));
+    }else {
+      coarser_grids.insert(std::make_pair(&rhs(i), temp_grid));
+    }
+
+  }
+
+  /*
+   * Create finer grids
+   */
+  for (int i=sol.MinLevel(); i<sol.MaxLevel(); ++i) {
+
+    TempGrid* temp_grid = new TempGrid();
+    temp_grid->SetPropertiesToFiner(sol(i), comm.BoundaryConditions());
+    if (temp_grid->Global().LocalBegin() == sol(i+1).Global().LocalBegin() &&
+        temp_grid->Global().LocalEnd() == sol(i+1).Global().LocalEnd()) {
+      delete temp_grid;
+      finer_grids.insert(std::make_pair(&sol(i), &sol(i+1)));
+    }else {
+      finer_grids.insert(std::make_pair(&sol(i), temp_grid));
+    }
+
+  }
+
+  for (int i=rhs.MinLevel(); i<rhs.MaxLevel(); ++i) {
+
+    TempGrid* temp_grid = new TempGrid();
+    temp_grid->SetPropertiesToFiner(rhs(i), comm.BoundaryConditions());
+    if (temp_grid->Global().LocalBegin() == rhs(i+1).Global().LocalBegin() &&
+        temp_grid->Global().LocalEnd() == rhs(i+1).Global().LocalEnd()) {
+      delete temp_grid;
+      finer_grids.insert(std::make_pair(&rhs(i), &rhs(i+1)));
+    }else {
+      finer_grids.insert(std::make_pair(&rhs(i), temp_grid));
+    }
+
+  }
+
+  /*
+   * Create global communicators
+   */
+  for (int i=sol.MinLevel()+1; i<sol.MaxLevel(); ++i)
+    CreateGlobalCommunicator(comm_global, &sol(i), &CoarserGrid(sol(i+1)), &FinerGrid(sol(i-1)));
+
+  CreateGlobalCommunicator(comm_global, &sol(sol.MinLevel()), &CoarserGrid(sol(sol.MinLevel()+1)));
+  CreateGlobalCommunicator(comm_global, &sol(sol.MaxLevel()), &FinerGrid(sol(sol.MaxLevel()-1)));
+
+  MPI_Comm my_comm_global;
+  MPI_Comm_dup(comm_global, &my_comm_global);
+  communicators_global.insert(std::make_pair(0, my_comm_global));
+
+  /*
+   * Create local communicators
+   */
+  for (int i=sol.MinLevel(); i<=sol.MaxLevel(); ++i)
+    CreateLocalCommunicator(comm_global, sol(i));
+
+  for (int i=sol.MinLevel(); i<sol.MaxLevel(); ++i)
+    CreateLocalCommunicator(comm_global, FinerGrid(sol(i)));
+
+  for (int i=sol.MaxLevel(); i>sol.MinLevel(); --i)
+    CreateLocalCommunicator(comm_global, CoarserGrid(sol(i)));
+
+  if (MG::GetFactory().TestObject("PARTICLE_NEAR_FIELD_CELLS"))
+    CreateLocalCommunicator(comm_global, comm.GetParticleGrid());
+
+  /*
+   * Create single grid datatypes
+   */
+  for (int i=sol.MinLevel(); i<=sol.MaxLevel(); ++i)
+    datatypes_local.insert(std::make_pair(KeyUnsorted(sol(i), 0), VMG::MPI::DatatypesLocal(sol(i), CommunicatorLocal(sol(i)), true)));
+
+  for (grid_iter=finer_grids.begin(); grid_iter!=finer_grids.end(); ++grid_iter)
+    datatypes_local.insert(std::make_pair(KeyUnsorted(*grid_iter->second, 0), VMG::MPI::DatatypesLocal(*grid_iter->second, CommunicatorLocal(*grid_iter->second), true)));
+
+  for (grid_iter=coarser_grids.begin(); grid_iter!=coarser_grids.end(); ++grid_iter)
+    datatypes_local.insert(std::make_pair(KeyUnsorted(*grid_iter->second, 0), VMG::MPI::DatatypesLocal(*grid_iter->second, CommunicatorLocal(*grid_iter->second), true)));
+
+  if (MG::GetFactory().TestObject("PARTICLE_NEAR_FIELD_CELLS"))
+    datatypes_local.insert(std::make_pair(KeyUnsorted(comm.GetParticleGrid(), 0), VMG::MPI::DatatypesLocal(comm.GetParticleGrid(), CommunicatorLocal(comm.GetParticleGrid()), true)));
+
+  /*
+   * Create two grid datatypes
+   */
+  for (int i=sol.MinLevel(); i<sol.MaxLevel(); ++i) {
+    AddDatatypeGlobal(sol(i), CoarserGrid(sol(i+1)), 0);
+    AddDatatypeGlobal(CoarserGrid(sol(i+1)), sol(i), 1);
+  }
+
+  for (int i=sol.MaxLevel(); i>sol.MinLevel(); --i) {
+    AddDatatypeGlobal(sol(i), FinerGrid(sol(i-1)), 0);
+    AddDatatypeGlobal(FinerGrid(sol(i-1)), sol(i), 1);
+  }
+
+  InitializeBoundaryValues();
+}
+
+void VMG::MPI::Settings::CreateGlobalCommunicator(MPI_Comm& comm_global, const Grid* grid_1, const Grid* grid_2, const Grid* grid_3)
+{
+  int rank;
+  MPI_Comm comm_new;
+
+  const bool in_communicator = (grid_1->Global().LocalSize().Product() > 0) ||
+      (grid_2 && grid_2->Global().LocalSize().Product() > 0) ||
+      (grid_3 && grid_3->Global().LocalSize().Product() > 0);
+
+  MPI_Comm_rank(comm_global, &rank);
+
+  if (in_communicator) {
+    Index dims, periods, coords;
+    MPI_Comm comm_temp;
+    MPI_Cart_get(comm_global, 3, dims.vec(), periods.vec(), coords.vec());
+    MPI_Comm_split(comm_global, 1, rank, &comm_temp);
+    dims = GlobalDims(comm_temp, coords);
+    MPI_Cart_create(comm_temp, 3, dims.vec(), periods.vec(), 0, &comm_new);
+    MPI_Comm_free(&comm_temp);
+  }else {
+    MPI_Comm_split(comm_global, MPI_UNDEFINED, rank, &comm_new);
+  }
+
+  communicators_global.insert(std::make_pair(grid_1->Level(), comm_new));
+}
+
+void VMG::MPI::Settings::CreateLocalCommunicator(MPI_Comm& comm_global, const Grid& grid)
+{
+  int rank, comm_equal;
+  MPI_Comm comm_new;
+  std::set<MPI_Comm>::iterator iter;
+
+  MPI_Comm_rank(comm_global, &rank);
+
+  if (grid.Global().LocalSize().Product() > 0) {
+    Index dims, periods, coords;
+    MPI_Comm comm_temp;
+    MPI_Cart_get(comm_global, 3, dims.vec(), periods.vec(), coords.vec());
+    MPI_Comm_split(comm_global, 1, rank, &comm_temp);
+    dims = GlobalDims(comm_temp, coords);
+    MPI_Cart_create(comm_temp, 3, dims.vec(), periods.vec(), 0, &comm_new);
+    MPI_Comm_free(&comm_temp);
+  }else {
+    MPI_Comm_split(comm_global, MPI_UNDEFINED, rank, &comm_new);
+  }
+
+  if (comm_new != MPI_COMM_NULL) {
+    for (iter=communicators_local_unique.begin(); iter!=communicators_local_unique.end(); ++iter) {
+      if (*iter != MPI_COMM_NULL) {
+        MPI_Comm_compare(comm_new, *iter, &comm_equal);
+        assert(comm_equal != MPI_SIMILAR);
+        if (comm_equal == MPI_IDENT || comm_equal == MPI_CONGRUENT) {
+          MPI_Comm_free(&comm_new);
+          comm_new = *iter;
+          break;
+        }
+      }
+    }
+  }
+
+  std::pair<std::set<MPI_Comm>::iterator, bool> insert_result = communicators_local_unique.insert(comm_new);
+  communicators_local.insert(std::make_pair(KeyUnsorted(grid, 0), *insert_result.first));
+}
+
+void VMG::MPI::Settings::AddDatatypeGlobal(const Grid& grid_old, const Grid& grid_new, const int& direction)
+{
+  MPI_Comm comm = CommunicatorGlobal(grid_old);
+
+  // Insert into map
+  std::pair< std::map<VMG::MPI::KeyUnsorted, VMG::MPI::DatatypesGlobal>::iterator, bool > insert_result =
+      datatypes_global.insert(std::make_pair(VMG::MPI::KeyUnsorted(grid_old, grid_new, direction), VMG::MPI::DatatypesGlobal()));
+  VMG::MPI::DatatypesGlobal& dt_global = insert_result.first->second;
+  bool dt_is_new = insert_result.second;
+
+
+  if (comm != MPI_COMM_NULL) {
+
+    Index begin, end, offset_old, offset_new;
+    int rank, size;
+
+    MPI_Comm_rank(comm, &rank);
+    MPI_Comm_size(comm, &size);
+
+    std::vector<int> buffer;
+    buffer.resize(6*size);
+
+    // Compute local offset for ghost cells
+    for (int i=0; i<3; ++i) {
+      offset_old[i] = (grid_old.Local().HaloSize1()[i] > 0 ? grid_old.Local().Begin()[i] : 0);
+      offset_new[i] = (grid_new.Local().HaloSize1()[i] > 0 ? grid_new.Local().Begin()[i] : 0);
+    }
+
+    // Publish which grid part this process can offer
+    if (&grid_old == &grid_new) {
+      for (int i=0; i<6; ++i)
+        buffer[6*rank+i] = 0;
+    }else {
+      for (int i=0; i<3; ++i) {
+        buffer[6*rank+i] = grid_old.Global().LocalBegin()[i];
+        buffer[6*rank+i+3] = grid_old.Global().LocalEnd()[i];
+      }
+    }
+
+    MPI_Allgather(MPI_IN_PLACE, 6, MPI_INT, &buffer.front(), 6, MPI_INT, comm);
+
+    if (dt_is_new) {
+
+      // Decide who offers a useful grid part
+      for (int i=0; i<size; ++i) {
+        for (int j=0; j<3; ++j) {
+          begin[j] = buffer[6*i+j];
+          end[j] = buffer[6*i+j+3];
+        }
+
+        begin = begin.Clamp(grid_new.Global().LocalBegin(), grid_new.Global().LocalEnd());
+        end = end.Clamp(grid_new.Global().LocalBegin(), grid_new.Global().LocalEnd());
+
+        if ((end-begin).Product() > 0) {
+          // This process has a useful part
+          dt_global.Receive().push_back(VMG::MPI::Datatype(grid_new.Local().SizeTotal(),
+              end - begin,
+              begin - grid_new.Global().LocalBegin() + offset_new,
+              i, 0, 0, true));
+        }
+      }
+    }
+
+    // Publish which grid parts this process needs
+    for (int i=0; i<3; ++i) {
+      buffer[6*rank+i] = grid_new.Global().LocalBegin()[i];
+      buffer[6*rank+i+3] = grid_new.Global().LocalEnd()[i];
+    }
+
+    MPI_Allgather(MPI_IN_PLACE, 6, MPI_INT, &buffer.front(), 6, MPI_INT, comm);
+
+    if (dt_is_new) {
+
+      // Decide who needs a part of my grid
+      for (int i=0; i<size; ++i) {
+
+        if ((i == rank) && (&grid_old == &grid_new))
+          continue;
+
+        for (int j=0; j<3; ++j) {
+          begin[j] = buffer[6*i+j];
+          end[j] = buffer[6*i+j+3];
+        }
+
+        begin = begin.Clamp(grid_old.Global().LocalBegin(), grid_old.Global().LocalEnd());
+        end = end.Clamp(grid_old.Global().LocalBegin(), grid_old.Global().LocalEnd());
+
+        if ((end-begin).Product() > 0) {
+          // This process needs one of my parts
+          dt_global.Send().push_back(VMG::MPI::Datatype(grid_old.Local().SizeTotal(),
+              end - begin,
+              begin - grid_old.Global().LocalBegin() + offset_old,
+              i, 0, 0, true));
+        }
+      }
+    }
+  }
+}
+
+MPI_Datatype& VMG::MPI::Settings::Datatype(const Index& begin, const Index& end,
+    const Index& size_local, const Index& size_global,
+    const int& level)
+{
+  KeyUnsorted k(begin, end, size_local, size_global, level, 0);
+  std::map<KeyUnsorted, MPI_Datatype>::iterator iter = datatypes.find(k);
+
+  if (iter != datatypes.end())
+    return iter->second;
+
+  MPI_Datatype dt;
+  Index sizes =  size_local;
+  Index subsizes = end - begin;
+  Index starts = begin;
+
+  MPI_Type_create_subarray(3, sizes.vec(), subsizes.vec(), starts.vec(), MPI_ORDER_C, MPI_DOUBLE, &dt);
+  MPI_Type_commit(&dt);
+
+  return datatypes.insert(std::make_pair(k, dt)).first->second;
+}
+
+Index VMG::MPI::Settings::GlobalDims(MPI_Comm comm, Index pos)
+{
+  std::set<int> unique_set[3];
+  Index dims;
+
+  int size;
+  MPI_Comm_size(comm, &size);
+
+  int* coordinates = new int[3*size];
+
+  MPI_Allgather(pos.vec(), 3, MPI_INT, coordinates, 3, MPI_INT, comm);
+
+  for (int i=0; i<size; ++i)
+    for (int j=0; j<3; ++j)
+      unique_set[j].insert(coordinates[3*i+j]);
+
+  for (int j=0; j<3; ++j)
+    dims[j] = static_cast<int>(unique_set[j].size());
+
+  delete [] coordinates;
+
+  return dims;
+}
+
+void VMG::MPI::Settings::InitializeBoundaryValues()
+{
+  assert(bv_ranks.size() == 0);
+
+  if (MG::GetFactory().TestObject("BOUNDARY_VALUE_SETTER")) {
+
+    Index coord;
+
+    const int level_index = MG::GetRhs()->MaxLevel() - MG::GetRhs()->GlobalMaxLevel();
+    const std::vector<BoundaryValue>& bvs = MG::GetBoundaryValueSetter()->BoundaryValues();
+    const std::map<Index, std::vector<GlobalIndices> >& global = MG::GetComm()->DecomposedGlobal();
+
+    assert(global.find(0)->second[level_index].BoundaryType() == GlobalMax);
+
+    MPI_Comm comm = CommunicatorLocal((*MG::GetRhs())(MG::GetRhs()->GlobalMaxLevel()));
+
+    bv_ranks.reserve(bvs.size());
+
+    for (std::vector<BoundaryValue>::const_iterator iter_b = bvs.begin(); iter_b != bvs.end(); ++iter_b) {
+      for (std::map<Index, std::vector<GlobalIndices> >::const_iterator iter_g = global.begin(); iter_g != global.end(); ++iter_g) {
+        if (iter_b->GetIndex().IsComponentwiseGreaterOrEqual(iter_g->second[level_index].LocalBegin()) &&
+            iter_b->GetIndex().IsComponentwiseLess(iter_g->second[level_index].LocalEnd())) {
+            bv_ranks.push_back(0);
+            coord = iter_g->first;
+            MPI_Cart_rank(comm, coord.vec(), &bv_ranks.back());
+          break;
+        }
+      }
+    }
+  }
+}
+
+std::string VMG::MPI::Settings::ToString() const
+{
+  std::stringstream str;
+  std::map<KeyUnsorted, VMG::MPI::DatatypesGlobal>::const_iterator iter_global;
+  std::map<KeyUnsorted, VMG::MPI::DatatypesLocal>::const_iterator iter_local;
+
+  str << "#### MPI_SETTINGS ####" << std::endl;
+
+  for (iter_global = datatypes_global.begin(); iter_global != datatypes_global.end(); ++iter_global)
+    str << iter_global->second << std::endl;
+
+  for (iter_local = datatypes_local.begin(); iter_local != datatypes_local.end(); ++iter_local)
+    str << iter_local->second << std::endl;
+
+  str << "######################";
+
+  return str.str();
+}
+
+std::ostream& VMG::MPI::operator<<(std::ostream& out, const VMG::MPI::Settings& settings)
+{
+  return out << settings.ToString();
+}
+
+VMG::MPI::Settings::Settings()
+{
+}
+
+VMG::MPI::Settings::~Settings()
+{
+  std::map<int, MPI_Comm>::iterator iter_comm_global;
+  for (iter_comm_global=communicators_global.begin(); iter_comm_global!=communicators_global.end(); ++iter_comm_global)
+    if (iter_comm_global->second != MPI_COMM_NULL)
+      MPI_Comm_free(&iter_comm_global->second);
+
+  /*
+   * We simply copied some communicators so we have to make sure that we free
+   * each communicator exactly once
+   */
+  std::set<MPI_Comm>::iterator iter_comm_set;
+  for (iter_comm_set=communicators_local_unique.begin(); iter_comm_set!=communicators_local_unique.end(); ++iter_comm_set)
+    if (*iter_comm_set != MPI_COMM_NULL) {
+      MPI_Comm comm_temp = *iter_comm_set;
+      MPI_Comm_free(&comm_temp);
+    }
+
+  std::map<KeyUnsorted, MPI_Datatype>::iterator iter_dt;
+  for (iter_dt=datatypes.begin(); iter_dt!=datatypes.end(); ++iter_dt)
+    if (iter_dt->second != MPI_DATATYPE_NULL)
+      MPI_Type_free(&iter_dt->second);
+
+  std::map<const Grid*, Grid*>::iterator iter_grid;
+  for (iter_grid=finer_grids.begin(); iter_grid!=finer_grids.end(); ++iter_grid)
+    if (iter_grid->second->Father() == NULL)
+      delete iter_grid->second;
+
+  for (iter_grid=coarser_grids.begin(); iter_grid!=coarser_grids.end(); ++iter_grid)
+    if (iter_grid->second->Father() == NULL)
+      delete iter_grid->second;
+}
Index: ThirdParty/vmg/src/comm/mpi/settings.hpp
===================================================================
--- ThirdParty/vmg/src/comm/mpi/settings.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/comm/mpi/settings.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,103 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
+ *
+ *  vmg is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  vmg is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * @file   settings.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Jan 2 18:45:22 2012
+ *
+ * @brief  Computes several MPI-related settings.
+ *
+ */
+
+#ifndef SETTINGS_HPP_
+#define SETTINGS_HPP_
+
+#ifndef HAVE_MPI
+#error You need MPI in order to compile VMG::MPI::Settings
+#endif
+
+#include <map>
+#include <set>
+
+#include "comm/mpi/datatypes_global.hpp"
+#include "comm/mpi/datatypes_local.hpp"
+#include "comm/mpi/key.hpp"
+
+namespace VMG
+{
+
+class TempGrid;
+
+namespace MPI
+{
+
+class DatatypesGlobal;
+
+class Settings
+{
+public:
+  Settings();
+  ~Settings();
+
+  void ComputeSettings(Multigrid& sol, Multigrid& rhs, MPI_Comm& comm);
+
+  Grid& FinerGrid(const Grid& grid);
+  Grid& CoarserGrid(const Grid& grid);
+
+  MPI_Comm CommunicatorGlobal(const Grid& grid) const;
+  MPI_Comm CommunicatorLocal(const Grid& grid) const;
+
+  MPI_Datatype& Datatype(const Index& begin, const Index& end,
+			 const Index& size_local, const Index& size_global,
+			 const int& level);
+  VMG::MPI::DatatypesGlobal& DatatypesGlobal(const Grid& grid_old, const Grid& grid_new, const int& direction);
+  VMG::MPI::DatatypesLocal& DatatypesLocal(const Grid& grid);
+
+  const std::vector<int>& BoundaryRanks() const {return bv_ranks;}
+
+  std::string ToString() const;
+
+private:
+  Index GlobalDims(MPI_Comm comm, Index pos);
+  void AddDatatypeGlobal(const Grid& grid_old, const Grid& grid_new, const int& direction);
+
+  void CreateGlobalCommunicator(MPI_Comm& comm_global, const Grid* grid_1, const Grid* grid_2=NULL, const Grid* grid_3=NULL);
+  void CreateLocalCommunicator(MPI_Comm& comm_global, const Grid& grid);
+
+  void InitializeBoundaryValues();
+
+  std::map<int, MPI_Comm> communicators_global;
+  std::map<KeyUnsorted, MPI_Comm> communicators_local;
+  std::set<MPI_Comm> communicators_local_unique;
+  std::map<const Grid*, Grid*> finer_grids, coarser_grids;
+  std::map<KeyUnsorted, MPI_Datatype> datatypes;
+  std::map<KeyUnsorted, VMG::MPI::DatatypesGlobal> datatypes_global;
+  std::map<KeyUnsorted, VMG::MPI::DatatypesLocal> datatypes_local;
+
+  std::vector<int> bv_ranks;
+
+};
+
+std::ostream& operator<<(std::ostream& out, const VMG::MPI::Settings& settings);
+
+}
+
+}
+
+#endif /* SETTINGS_HPP_ */
