Index: ThirdParty/vmg/src/discretization/boundary_value.hpp
===================================================================
--- ThirdParty/vmg/src/discretization/boundary_value.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/discretization/boundary_value.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,75 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012-2013 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   boundary_value.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Sun Apr 14 17:23:40 2013
+ *
+ * @brief  Holds a boundary value for open boundary
+ *         conditions.
+ *
+ */
+
+#ifndef BOUNDARY_VALUE_HPP_
+#define BOUNDARY_VALUE_HPP_
+
+#include "base/index.hpp"
+#include "base/vector.hpp"
+#include "grid/grid.hpp"
+
+namespace VMG
+{
+
+class BoundaryValue
+{
+public:
+  BoundaryValue() :
+    _index(0),
+    _pos(0.0),
+    _val(0.0)
+  {}
+
+  BoundaryValue(const Grid& grid, const Index& index) :
+    _index(index),
+    _pos(grid.GetSpatialPosGlobal(index)),
+    _val(0.0)
+  {}
+
+  vmg_float& Val() {return _val;}
+  const vmg_float& Val() const {return _val;}
+
+  const Index& GetIndex() const {return _index;}
+  void SetIndex(const Grid& grid, const Index& index)
+  {
+    _index = index;
+    _pos = grid.GetSpatialPosGlobal(index);
+  }
+
+  const Vector&  GetPos() const {return _pos;}
+
+private:
+  Index _index;
+  Vector _pos;
+  vmg_float _val;
+};
+
+}
+
+
+#endif /* BOUNDARY_VALUE_HPP_ */
Index: ThirdParty/vmg/src/discretization/boundary_value_setter.cpp
===================================================================
--- ThirdParty/vmg/src/discretization/boundary_value_setter.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/discretization/boundary_value_setter.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,126 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012-2013 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   boundary_value_setter.cpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Apr 15 16:21:40 2013
+ *
+ * @brief  Class for setting boundary values with open
+ *         boundary conditions.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include "libvmg_config.h"
+#endif
+
+#include "comm/comm.hpp"
+#include "discretization/boundary_value_setter.hpp"
+#include "grid/grid.hpp"
+#include "grid/grid_index_translations.hpp"
+#include "grid/multigrid.hpp"
+#include "mg.hpp"
+
+using namespace VMG;
+
+void BoundaryValueSetter::ComputeBoundaryValues()
+{
+  const Grid& grid_ml = (*MG::GetRhs())(MG::GetRhs()->MaxLevel());
+  Grid& rhs_gml = (*MG::GetRhs())(MG::GetRhs()->GlobalMaxLevel());
+  Grid& sol_gml = (*MG::GetSol())(MG::GetSol()->GlobalMaxLevel());
+
+  for (std::vector<BoundaryValue>::iterator iter_b = boundary_values.begin(); iter_b != boundary_values.end(); ++iter_b)
+    iter_b->Val() = 0.0;
+
+  for (GridIterator iter_g = grid_ml.Iterators().Local().Begin(); iter_g != grid_ml.Iterators().Local().End(); ++iter_g)
+    for (std::vector<BoundaryValue>::iterator iter_b = boundary_values.begin(); iter_b != boundary_values.end(); ++iter_b)
+      iter_b->Val() += grid_ml.GetVal(*iter_g) * greens_function((iter_b->GetPos() - grid_ml.GetSpatialPos(*iter_g)).Length());
+
+  std::vector<BoundaryValue> list_bv = MG::GetComm()->CommBoundaryValues();
+
+  for (std::vector<BoundaryValue>::iterator iter = list_bv.begin(); iter != list_bv.end(); ++iter) {
+    const Index i = GridIndexTranslations::GlobalToLocal(rhs_gml, iter->GetIndex());
+    const vmg_float val = grid_ml.Extent().MeshWidth().Product() * iter->Val();
+    rhs_gml(i) = val;
+    sol_gml(i) = val;
+  }
+
+  MG::GetComm()->PrintGrid(rhs_gml, "DEBUG GREEN BOUNDARY");
+}
+
+void BoundaryValueSetter::InitBoundaryValueSetter()
+{
+  Comm& comm = *MG::GetComm();
+
+  const Grid& grid = (*MG::GetRhs())(MG::GetRhs()->GlobalMaxLevel());
+
+  Index begin_local, end_local, i;
+  for (int i=0; i<3; ++i) {
+    begin_local[i] = grid.Global().GlobalBegin()[i] + (comm.BoundaryConditions()[i] == Open ? 1 : 0);
+    end_local[i] = grid.Global().GlobalEnd()[i] - (comm.BoundaryConditions()[i] == Open ? 1 : 0);
+  }
+
+  if (comm.BoundaryConditions().X() == Open) {
+
+    i.X() = grid.Global().GlobalBegin().X();
+    for (i.Y() = begin_local.Y(); i.Y() < end_local.Y(); ++i.Y())
+      for (i.Z() = begin_local.Z(); i.Z() < end_local.Z(); ++i.Z())
+        boundary_values.push_back(BoundaryValue(grid, i));
+
+    i.X() = grid.Global().GlobalEnd().X()-1;
+    for (i.Y() = begin_local.Y(); i.Y() < end_local.Y(); ++i.Y())
+      for (i.Z() = begin_local.Z(); i.Z() < end_local.Z(); ++i.Z())
+        boundary_values.push_back(BoundaryValue(grid, i));
+
+  }
+
+  if (comm.BoundaryConditions().Y() == Open) {
+
+    i.Y() = grid.Global().GlobalBegin().Y();
+    for (i.X() = begin_local.X(); i.X() < end_local.X(); ++i.X())
+      for (i.Z() = begin_local.Z(); i.Z() < end_local.Z(); ++i.Z())
+        boundary_values.push_back(BoundaryValue(grid, i));
+
+    i.Y() = grid.Global().GlobalEnd().Y() - 1;
+    for (i.X() = begin_local.X(); i.X() < end_local.X(); ++i.X())
+      for (i.Z() = begin_local.Z(); i.Z() < end_local.Z(); ++i.Z())
+        boundary_values.push_back(BoundaryValue(grid, i));
+
+  }
+
+  if (comm.BoundaryConditions().Z() == Open) {
+
+    i.Z() = grid.Global().GlobalBegin().Z();
+    for (i.X() = begin_local.X(); i.X() < end_local.X(); ++i.X())
+      for (i.Y() = begin_local.Y(); i.Y() < end_local.Y(); ++i.Y())
+        boundary_values.push_back(BoundaryValue(grid, i));
+
+    i.Z() = grid.Global().GlobalEnd().Z() - 1;
+    for (i.X() = begin_local.X(); i.X() < end_local.X(); ++i.X())
+      for (i.Y() = begin_local.Y(); i.Y() < end_local.Y(); ++i.Y())
+        boundary_values.push_back(BoundaryValue(grid, i));
+
+  }
+
+}
+
+BoundaryValueSetter::~BoundaryValueSetter()
+{
+
+}
Index: ThirdParty/vmg/src/discretization/boundary_value_setter.hpp
===================================================================
--- ThirdParty/vmg/src/discretization/boundary_value_setter.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/discretization/boundary_value_setter.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,66 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012-2013 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   boundary_value_setter.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Apr 15 16:21:40 2013
+ *
+ * @brief  Class for setting boundary values with open
+ *         boundary conditions.
+ *
+ */
+
+#ifndef BOUNDARY_VALUE_SETTER_HPP_
+#define BOUNDARY_VALUE_SETTER_HPP_
+
+#include <vector>
+
+#include "base/object.hpp"
+#include "discretization/boundary_value.hpp"
+
+namespace VMG
+{
+
+class Grid;
+
+class BoundaryValueSetter : public Object
+{
+public:
+  BoundaryValueSetter(bool _register = true) :
+    Object("BOUNDARY_VALUE_SETTER", _register)
+  {InitBoundaryValueSetter();}
+
+  virtual ~BoundaryValueSetter();
+
+  void ComputeBoundaryValues();
+
+  const std::vector<BoundaryValue>& BoundaryValues() const {return boundary_values;}
+
+protected:
+  virtual vmg_float greens_function(vmg_float arg) = 0;
+
+private:
+  void InitBoundaryValueSetter();
+
+  std::vector<BoundaryValue> boundary_values;
+};
+
+}
+
+#endif /* BOUNDARY_VALUE_SETTER_HPP_ */
Index: ThirdParty/vmg/src/discretization/boundary_value_setter_open.hpp
===================================================================
--- ThirdParty/vmg/src/discretization/boundary_value_setter_open.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/discretization/boundary_value_setter_open.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,53 @@
+/*
+ *    vmg - a versatile multigrid solver
+ *    Copyright (C) 2012-2013 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   boundary_value_setter_open.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Apr 15 16:21:40 2013
+ *
+ * @brief  For open boundary conditions.
+ *
+ */
+
+#ifndef BOUNDARY_VALUE_SETTER_OPEN_HPP_
+#define BOUNDARY_VALUE_SETTER_OPEN_HPP_
+
+#include "base/math.hpp"
+#include "discretization/boundary_value_setter.hpp"
+
+namespace VMG
+{
+
+class Grid;
+
+class BoundaryValueSetterOpen : public BoundaryValueSetter
+{
+public:
+  BoundaryValueSetterOpen(bool _register = true) :
+    BoundaryValueSetter(_register)
+  {}
+  virtual ~BoundaryValueSetterOpen() {}
+
+protected:
+  vmg_float greens_function(vmg_float arg) {return 1.0 / (4.0 * Math::pi * arg);}
+};
+
+}
+
+#endif /* BOUNDARY_VALUE_SETTER_OPEN_HPP_ */
Index: ThirdParty/vmg/src/discretization/discretization_poisson_fd.cpp
===================================================================
--- ThirdParty/vmg/src/discretization/discretization_poisson_fd.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/discretization/discretization_poisson_fd.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,96 @@
+/*
+ *    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   discretization_poisson_fd.cpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Apr 18 13:03:47 2011
+ *
+ * @brief  Discretization of the poisson equation
+ *         using the standard 7-point stencil or
+ *         the Collatz Mehrstellen Ansatz.
+ *         Discretization error: O(h^2) or O(h^4)
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include "discretization/discretization_poisson_fd.hpp"
+#include "grid/grid.hpp"
+#include "mg.hpp"
+
+using namespace VMG;
+
+void DiscretizationPoissonFD::InitDiscretizationPoissonFD()
+{
+  switch (order)
+    {
+    case 2:
+      stencil.SetDiag(6.0);
+      stencil.push_back(-1,  0,  0, -1.0);
+      stencil.push_back( 1,  0,  0, -1.0);
+      stencil.push_back( 0, -1,  0, -1.0);
+      stencil.push_back( 0,  1,  0, -1.0);
+      stencil.push_back( 0,  0, -1, -1.0);
+      stencil.push_back( 0,  0,  1, -1.0);
+      break;
+    case 4:
+      stencil.SetDiag(24.0/6.0);
+      stencil.push_back(-1,  0,  0, -2.0/6.0);
+      stencil.push_back( 1,  0,  0, -2.0/6.0);
+      stencil.push_back( 0, -1,  0, -2.0/6.0);
+      stencil.push_back( 0,  1,  0, -2.0/6.0);
+      stencil.push_back( 0,  0, -1, -2.0/6.0);
+      stencil.push_back( 0,  0,  1, -2.0/6.0);
+      stencil.push_back(-1, -1,  0, -1.0/6.0);
+      stencil.push_back(-1,  1,  0, -1.0/6.0);
+      stencil.push_back( 1, -1,  0, -1.0/6.0);
+      stencil.push_back( 1,  1,  0, -1.0/6.0);
+      stencil.push_back(-1,  0, -1, -1.0/6.0);
+      stencil.push_back(-1,  0,  1, -1.0/6.0);
+      stencil.push_back( 1,  0, -1, -1.0/6.0);
+      stencil.push_back( 1,  0,  1, -1.0/6.0);
+      stencil.push_back( 0, -1, -1, -1.0/6.0);
+      stencil.push_back( 0, -1,  1, -1.0/6.0);
+      stencil.push_back( 0,  1, -1, -1.0/6.0);
+      stencil.push_back( 0,  1,  1, -1.0/6.0);
+      break;
+    default:
+      assert(0 != "vmg choose discretization order 2 or 4");
+      break;
+    }
+}
+
+void DiscretizationPoissonFD::ModifyRightHandSide()
+{
+  if (order == 4) {
+
+    Stencil stencil(6.0/12.0);
+    stencil.push_back(-1,  0,  0, 1.0/12.0);
+    stencil.push_back( 1,  0,  0, 1.0/12.0);
+    stencil.push_back( 0, -1,  0, 1.0/12.0);
+    stencil.push_back( 0,  1,  0, 1.0/12.0);
+    stencil.push_back( 0,  0, -1, 1.0/12.0);
+    stencil.push_back( 0,  0,  1, 1.0/12.0);
+
+    stencil.Apply(MG::GetRhsMaxLevel());
+
+  }
+}
Index: ThirdParty/vmg/src/discretization/discretization_poisson_fd.hpp
===================================================================
--- ThirdParty/vmg/src/discretization/discretization_poisson_fd.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/discretization/discretization_poisson_fd.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,60 @@
+/*
+ *    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   discretization_poisson_fd.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Apr 18 13:01:40 2011
+ *
+ * @brief  Finite difference discretization of the
+ *         Poisson equation.
+ *
+ */
+
+#ifndef DISCRETIZATION_POISSON_FD_HPP_
+#define DISCRETIZATION_POISSON_FD_HPP_
+
+#include "base/discretization.hpp"
+#include "base/helper.hpp"
+
+namespace VMG
+{
+
+class DiscretizationPoissonFD : public Discretization
+{
+public:
+  DiscretizationPoissonFD(const int& order) :
+    Discretization(order)
+  {
+    InitDiscretizationPoissonFD();
+  }
+
+  vmg_float OperatorPrefactor(const Grid& grid) const
+  {
+    return 1.0 / Helper::pow_2(grid.Extent().MeshWidth().Max());
+  }
+
+  void ModifyRightHandSide();
+
+private:
+  void InitDiscretizationPoissonFD();
+};
+
+}
+
+#endif /* DISCRETIZATION_POISSON_FD_HPP_ */
Index: ThirdParty/vmg/src/discretization/discretization_poisson_fv.cpp
===================================================================
--- ThirdParty/vmg/src/discretization/discretization_poisson_fv.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/discretization/discretization_poisson_fv.cpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,330 @@
+/*
+ *    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   discretization_poisson_fv.cpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Apr 18 13:03:47 2011
+ *
+ * @brief  Finite volume discretization for the Poisson
+ *         equation. Absolutely equivalent to the finite
+ *         difference discretization unless you use
+ *         hierarchically coarsened grids.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <libvmg_config.h>
+#endif
+
+#include "comm/comm.hpp"
+#include "discretization/discretization_poisson_fv.hpp"
+#include "grid/grid_index_translations.hpp"
+
+using namespace VMG;
+
+void DiscretizationPoissonFV::InitDiscretizationPoissonFV()
+{
+  switch (order)
+  {
+  case 2:
+    stencil.SetDiag(6.0);
+    stencil.push_back(-1,  0,  0, -1.0);
+    stencil.push_back( 1,  0,  0, -1.0);
+    stencil.push_back( 0, -1,  0, -1.0);
+    stencil.push_back( 0,  1,  0, -1.0);
+    stencil.push_back( 0,  0, -1, -1.0);
+    stencil.push_back( 0,  0,  1, -1.0);
+    break;
+  case 4:
+    stencil.SetDiag(24.0/6.0);
+    stencil.push_back(-1,  0,  0, -2.0/6.0);
+    stencil.push_back( 1,  0,  0, -2.0/6.0);
+    stencil.push_back( 0, -1,  0, -2.0/6.0);
+    stencil.push_back( 0,  1,  0, -2.0/6.0);
+    stencil.push_back( 0,  0, -1, -2.0/6.0);
+    stencil.push_back( 0,  0,  1, -2.0/6.0);
+    stencil.push_back(-1, -1,  0, -1.0/6.0);
+    stencil.push_back(-1,  1,  0, -1.0/6.0);
+    stencil.push_back( 1, -1,  0, -1.0/6.0);
+    stencil.push_back( 1,  1,  0, -1.0/6.0);
+    stencil.push_back(-1,  0, -1, -1.0/6.0);
+    stencil.push_back(-1,  0,  1, -1.0/6.0);
+    stencil.push_back( 1,  0, -1, -1.0/6.0);
+    stencil.push_back( 1,  0,  1, -1.0/6.0);
+    stencil.push_back( 0, -1, -1, -1.0/6.0);
+    stencil.push_back( 0, -1,  1, -1.0/6.0);
+    stencil.push_back( 0,  1, -1, -1.0/6.0);
+    stencil.push_back( 0,  1,  1, -1.0/6.0);
+    break;
+  default:
+    assert(0 != "vmg choose discretization order 2 or 4");
+    break;
+  }
+}
+
+void DiscretizationPoissonFV::ModifyRightHandSide()
+{
+  if (order == 4) {
+
+    Grid& rhs = MG::GetRhsMaxLevel();
+
+    Stencil stencil(6.0/12.0);
+    stencil.push_back(-1,  0,  0, 1.0/12.0);
+    stencil.push_back( 1,  0,  0, 1.0/12.0);
+    stencil.push_back( 0, -1,  0, 1.0/12.0);
+    stencil.push_back( 0,  1,  0, 1.0/12.0);
+    stencil.push_back( 0,  0, -1, 1.0/12.0);
+    stencil.push_back( 0,  0,  1, 1.0/12.0);
+
+    stencil.Apply(rhs);
+
+  }
+}
+
+void DiscretizationPoissonFV::SetInnerBoundaryCompute(Grid& sol_f, Grid& rhs_f, Grid& sol_c) const
+{
+  Index i_c, i_f;
+
+  Comm& comm = *MG::GetComm();
+
+  const Boundary& bc = comm.BoundaryConditions();
+  const Index off((GridIndexTranslations::LocalToGlobal(sol_f, sol_f.Local().Begin())[0] % 2 == 0 ? 0 : 1),
+                  (GridIndexTranslations::LocalToGlobal(sol_f, sol_f.Local().Begin())[1] % 2 == 0 ? 0 : 1),
+                  (GridIndexTranslations::LocalToGlobal(sol_f, sol_f.Local().Begin())[2] % 2 == 0 ? 0 : 1));
+
+  const Index begin_f = sol_f.Local().Begin() - off;
+  const Index end_f = sol_f.Local().End();
+
+  const Index begin_c = GridIndexTranslations::GlobalFinestToLocal(sol_c, GridIndexTranslations::LocalToGlobalFinest(sol_f, begin_f));
+
+  const Index b1_f = sol_f.Local().BoundaryBegin1();
+  const Index b2_f = sol_f.Local().BoundaryBegin2();
+
+  const Index b1_finest = GridIndexTranslations::LocalToGlobalFinest(sol_f, sol_f.Local().Begin() + off);
+  const Index b2_finest = GridIndexTranslations::LocalToGlobalFinest(sol_f, sol_f.Local().End()-1);
+
+  const Index b1_c = GridIndexTranslations::GlobalFinestToLocal(sol_c, b1_finest) - 1;
+  const Index b2_c = GridIndexTranslations::GlobalFinestToLocal(sol_c, b2_finest) + 1;
+
+  const vmg_float c_1_3 = 1.0 / 3.0;
+  const vmg_float c_2_3 = 2.0 / 3.0;
+  const vmg_float c_4_3 = 4.0 / 3.0;
+
+  comm.CommToGhosts(sol_f);
+  comm.CommToGhosts(sol_c);
+
+  //
+  // X-direction
+  //
+  if (bc.X() == Open) {
+
+    if (sol_f.Local().BoundarySize1().X() > 0) {
+
+      for (i_f.Y()=begin_f.Y(), i_c.Y()=begin_c.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2, ++i_c.Y())
+        for (i_f.Z()=begin_f.Z(), i_c.Z()=begin_c.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2, ++i_c.Z())
+          rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (sol_c(b1_c.X(),i_c.Y(),i_c.Z()) - sol_f(b1_f.X()+1,i_f.Y(),i_f.Z()));
+
+      for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
+        for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2)
+          rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (rhs_f(b1_f.X(),i_f.Y()-1,i_f.Z()) + rhs_f(b1_f.X(),i_f.Y()+1,i_f.Z()));
+
+      for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2)
+        for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
+          rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (rhs_f(b1_f.X(),i_f.Y(),i_f.Z()-1) + rhs_f(b1_f.X(),i_f.Y(),i_f.Z()+1));
+
+      for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
+        for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
+          rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = 0.25 * (rhs_f(b1_f.X(),i_f.Y()-1,i_f.Z()-1) +
+                                                    rhs_f(b1_f.X(),i_f.Y()+1,i_f.Z()-1) +
+                                                    rhs_f(b1_f.X(),i_f.Y()-1,i_f.Z()+1) +
+                                                    rhs_f(b1_f.X(),i_f.Y()+1,i_f.Z()+1));
+
+      for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y())
+        for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z())
+          rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = sol_f(b1_f.X(),i_f.Y(),i_f.Z()) = c_2_3 * rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) +
+                                                                              c_4_3 * sol_f(b1_f.X()+1,i_f.Y(),i_f.Z()) -
+                                                                              c_1_3 * sol_f(b1_f.X()+2,i_f.Y(),i_f.Z());
+
+    }
+
+    if (sol_f.Local().BoundarySize2().X() > 0) {
+
+      for (i_f.Y()=begin_f.Y(), i_c.Y()=begin_c.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2, ++i_c.Y())
+        for (i_f.Z()=begin_f.Z(), i_c.Z()=begin_c.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2, ++i_c.Z())
+          rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (sol_c(b2_c.X(),i_c.Y(),i_c.Z()) - sol_f(b2_f.X()-1,i_f.Y(),i_f.Z()));
+
+      for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
+        for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2)
+          rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (rhs_f(b2_f.X(),i_f.Y()-1,i_f.Z()) + rhs_f(b2_f.X(),i_f.Y()+1,i_f.Z()));
+
+      for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2)
+        for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
+          rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (rhs_f(b2_f.X(),i_f.Y(),i_f.Z()-1) + rhs_f(b2_f.X(),i_f.Y(),i_f.Z()+1));
+
+      for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
+        for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
+          rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = 0.25 * (rhs_f(b2_f.X(),i_f.Y()-1,i_f.Z()-1) +
+                                                    rhs_f(b2_f.X(),i_f.Y()+1,i_f.Z()-1) +
+                                                    rhs_f(b2_f.X(),i_f.Y()-1,i_f.Z()+1) +
+                                                    rhs_f(b2_f.X(),i_f.Y()+1,i_f.Z()+1));
+
+      for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y())
+        for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z())
+          rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = sol_f(b2_f.X(),i_f.Y(),i_f.Z()) = c_2_3 * rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) +
+                                                                              c_4_3 * sol_f(b2_f.X()-1,i_f.Y(),i_f.Z()) -
+                                                                              c_1_3 * sol_f(b2_f.X()-2,i_f.Y(),i_f.Z());
+
+    }
+  }
+
+  //
+  // Y-direction
+  //
+  if (bc.Y() == Open) {
+
+    if (sol_f.Local().BoundarySize1().Y() > 0) {
+
+      for (i_f.X()=begin_f.X(), i_c.X()=begin_c.X(); i_f.X()<end_f.X(); i_f.X()+=2, ++i_c.X())
+        for (i_f.Z()=begin_f.Z(), i_c.Z()=begin_c.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2, ++i_c.Z())
+          rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = 0.5 * (sol_c(i_c.X(),b1_c.Y(),i_c.Z()) - sol_f(i_f.X(),b1_f.Y()+1,i_f.Z()));
+
+      for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
+        for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2)
+          rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = 0.5 * (rhs_f(i_f.X()-1,b1_f.Y(),i_f.Z()) + rhs_f(i_f.X()+1,b1_f.Y(),i_f.Z()));
+
+      for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); i_f.X()+=2)
+        for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
+          rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = 0.5 * (rhs_f(i_f.X(),b1_f.Y(),i_f.Z()-1) + rhs_f(i_f.X(),b1_f.Y(),i_f.Z()+1));
+
+      for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
+        for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
+          rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = 0.25 * (rhs_f(i_f.X()-1,b1_f.Y(),i_f.Z()-1) +
+                                                    rhs_f(i_f.X()+1,b1_f.Y(),i_f.Z()-1) +
+                                                    rhs_f(i_f.X()-1,b1_f.Y(),i_f.Z()+1) +
+                                                    rhs_f(i_f.X()+1,b1_f.Y(),i_f.Z()+1));
+
+      for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
+        for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z())
+          rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = sol_f(i_f.X(),b1_f.Y(),i_f.Z()) = c_2_3 * rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) +
+                                                                              c_4_3 * sol_f(i_f.X(),b1_f.Y()+1,i_f.Z()) -
+                                                                              c_1_3 * sol_f(i_f.X(),b1_f.Y()+2,i_f.Z());
+
+    }
+
+    if (sol_f.Local().BoundarySize2().Y() > 0) {
+
+      for (i_f.X()=begin_f.X(), i_c.X()=begin_c.X(); i_f.X()<end_f.X(); i_f.X()+=2, ++i_c.X())
+        for (i_f.Z()=begin_f.Z(), i_c.Z()=begin_c.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2, ++i_c.Z())
+          rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = 0.5 * (sol_c(i_c.X(),b2_c.Y(),i_c.Z()) - sol_f(i_f.X(),b2_f.Y()-1,i_f.Z()));
+
+      for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
+        for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2)
+          rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = 0.5 * (rhs_f(i_f.X()-1,b2_f.Y(),i_f.Z()) + rhs_f(i_f.X()+1,b2_f.Z(),i_f.Z()));
+
+      for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); i_f.X()+=2)
+        for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
+          rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = 0.5 * (rhs_f(i_f.X(),b2_f.Y(),i_f.Z()-1) + rhs_f(i_f.X(),b2_f.Y(),i_f.Z()+1));
+
+      for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
+        for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
+          rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = 0.25 * (rhs_f(i_f.X()-1,b2_f.Y(),i_f.Z()-1) +
+                                                    rhs_f(i_f.X()+1,b2_f.Y(),i_f.Z()-1) +
+                                                    rhs_f(i_f.X()-1,b2_f.Y(),i_f.Z()+1) +
+                                                    rhs_f(i_f.X()+1,b2_f.Y(),i_f.Z()+1));
+
+      for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
+        for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z())
+          rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = sol_f(i_f.X(),b2_f.Y(),i_f.Z()) = c_2_3 * rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) +
+                                                                              c_4_3 * sol_f(i_f.X(),b2_f.Y()-1,i_f.Z()) -
+                                                                              c_1_3 * sol_f(i_f.X(),b2_f.Y()-2,i_f.Z());
+
+    }
+
+  }
+
+  //
+  // Z-direction
+  //
+  if (bc.Z() == Open) {
+
+    if (sol_f.Local().BoundarySize1().Z() > 0) {
+
+      for (i_f.X()=begin_f.X(), i_c.X()=begin_c.X(); i_f.X()<end_f.X(); i_f.X()+=2, ++i_c.X())
+        for (i_f.Y()=begin_f.Y(), i_c.Y()=begin_c.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2, ++i_c.Y())
+          rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = 0.5 * (sol_c(i_c.X(),i_c.Y(),b1_c.Z()) - sol_f(i_f.X(),i_f.Y(),b1_f.Z()+1));
+
+      for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
+        for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2)
+          rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = 0.5 * (rhs_f(i_f.X()-1,i_f.Y(),b1_f.Z()) + rhs_f(i_f.X()+1,i_f.Y(),b1_f.Z()));
+
+      for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); i_f.X()+=2)
+        for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
+          rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = 0.5 * (rhs_f(i_f.X(),i_f.Y()-1,b1_f.Z()) + rhs_f(i_f.X(),i_f.Y()+1,b1_f.Z()));
+
+      for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
+        for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
+          rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = 0.25 * (rhs_f(i_f.X()-1,i_f.Y()-1,b1_f.Z()) +
+                                                    rhs_f(i_f.X()+1,i_f.Y()-1,b1_f.Z()) +
+                                                    rhs_f(i_f.X()-1,i_f.Y()+1,b1_f.Z()) +
+                                                    rhs_f(i_f.X()+1,i_f.Y()+1,b1_f.Z()));
+
+      for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
+        for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y())
+          rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = sol_f(i_f.X(),i_f.Y(),b1_f.Z()) = c_2_3 * rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) +
+                                                                              c_4_3 * sol_f(i_f.X(),i_f.Y(),b1_f.Z()+1) -
+                                                                              c_1_3 * sol_f(i_f.X(),i_f.Y(),b1_f.Z()+2);
+
+    }
+
+    if (sol_f.Local().BoundarySize2().Z() > 0) {
+
+      for (i_f.X()=begin_f.X(), i_c.X()=begin_c.X(); i_f.X()<end_f.X(); i_f.X()+=2, ++i_c.X())
+        for (i_f.Y()=begin_f.Y(), i_c.Y()=begin_c.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2, ++i_c.Y())
+          rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = 0.5 * (sol_c(i_c.X(),i_c.Y(),b2_c.Z()) - sol_f(i_f.X(),i_f.Y(),b2_f.Z()-1));
+
+      for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
+        for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2)
+          rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = 0.5 * (rhs_f(i_f.X()-1,i_f.Y(),b2_f.Z()) + rhs_f(i_f.X()+1,i_f.Y(),b2_f.Z()));
+
+      for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); i_f.X()+=2)
+        for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
+          rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = 0.5 * (rhs_f(i_f.X(),i_f.Y()-1,b2_f.Z()) + rhs_f(i_f.X(),i_f.Y()+1,b2_f.Z()));
+
+      for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
+        for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
+          rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = 0.25 * (rhs_f(i_f.X()-1,i_f.Y()-1,b2_f.Z()) +
+                                                    rhs_f(i_f.X()+1,i_f.Y()-1,b2_f.Z()) +
+                                                    rhs_f(i_f.X()-1,i_f.Y()+1,b2_f.Z()) +
+                                                    rhs_f(i_f.X()+1,i_f.Y()+1,b2_f.Z()));
+
+      for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
+        for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y())
+          rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = sol_f(i_f.X(),i_f.Y(),b2_f.Z()) = c_2_3 * rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) +
+                                                                              c_4_3 * sol_f(i_f.X(),i_f.Y(),b2_f.Z()-1) -
+                                                                              c_1_3 * sol_f(i_f.X(),i_f.Y(),b2_f.Z()-2);
+
+    }
+
+  }
+
+#ifdef DEBUG_MATRIX_CHECKS
+  rhs_f.IsConsistent();
+  sol_f.IsConsistent();
+#endif
+}
Index: ThirdParty/vmg/src/discretization/discretization_poisson_fv.hpp
===================================================================
--- ThirdParty/vmg/src/discretization/discretization_poisson_fv.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/discretization/discretization_poisson_fv.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,63 @@
+/*
+ *    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   discretization_poisson_fv.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Apr 18 13:03:47 2011
+ *
+ * @brief  Finite volume discretization for the Poisson
+ *         equation. Absolutely equivalent to the finite
+ *         difference discretization unless you use
+ *         hierarchically coarsened grids.
+ *
+ */
+
+#ifndef DISCRETIZATION_POISSON_FV_HPP_
+#define DISCRETIZATION_POISSON_FV_HPP_
+
+#include "base/discretization.hpp"
+#include "base/helper.hpp"
+
+namespace VMG
+{
+
+class DiscretizationPoissonFV : public Discretization
+{
+public:
+  DiscretizationPoissonFV(const int& order) :
+    Discretization(order)
+  {
+    InitDiscretizationPoissonFV();
+  }
+
+  vmg_float OperatorPrefactor(const Grid& grid) const
+  {
+    return 1.0  / Helper::pow_2(grid.Extent().MeshWidth().Max());
+  }
+
+  void ModifyRightHandSide();
+
+private:
+  void InitDiscretizationPoissonFV();
+  void SetInnerBoundaryCompute(Grid& sol_f, Grid& rhs_f, Grid& sol_c) const;
+};
+
+}
+
+#endif /* DISCRETIZATION_POISSON_FV_HPP_ */
