Index: ThirdParty/vmg/src/units/particle/bspline.hpp
===================================================================
--- ThirdParty/vmg/src/units/particle/bspline.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
+++ ThirdParty/vmg/src/units/particle/bspline.hpp	(revision 7faa5c793164ff799487bb1193ab2a18a3a7e2fb)
@@ -0,0 +1,313 @@
+/*
+ *    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   bspline.hpp
+ * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
+ * @date   Mon Nov 21 13:27:22 2011
+ *
+ * @brief  B-Splines for molecular dynamics.
+ *
+ */
+
+#ifndef BSPLINE_HPP_
+#define BSPLINE_HPP_
+
+#include <cmath>
+
+#include "base/helper.hpp"
+#include "base/index.hpp"
+#include "base/polynomial.hpp"
+#include "base/vector.hpp"
+#include "grid/grid.hpp"
+#include "units/particle/particle.hpp"
+
+namespace VMG
+{
+
+namespace Particle
+{
+
+class BSpline
+{
+public:
+  BSpline(const int& near_field_cells, const vmg_float& h);
+
+  vmg_float EvaluateSpline(const vmg_float& val) const
+  {
+    for (unsigned int i=0; i<intervals.size(); ++i)
+      if (val < intervals[i])
+        return spline_nom[i](val) / spline_denom[i](val);
+    return 0.0;
+  }
+
+  void SetSpline(Grid& grid, const Particle& p) const
+  {
+    assert(p.Pos().X() >= grid.Extent().Begin().X() && p.Pos().X() < grid.Extent().End().X());
+    assert(p.Pos().Y() >= grid.Extent().Begin().Y() && p.Pos().Y() < grid.Extent().End().Y());
+    assert(p.Pos().Z() >= grid.Extent().Begin().Z() && p.Pos().Z() < grid.Extent().End().Z());
+
+    vmg_float *vals = new vmg_float[Helper::intpow(2*near_field_cells+1,3)];
+
+    vmg_float temp_val;
+    vmg_float int_val = 0.0;
+    int c = 0;
+
+    const int index_global_x = grid.Global().GlobalBegin().X() + std::floor((p.Pos().X() - grid.Extent().Begin().X()) / grid.Extent().MeshWidth().X());
+    const int index_global_y = grid.Global().GlobalBegin().Y() + std::floor((p.Pos().Y() - grid.Extent().Begin().Y()) / grid.Extent().MeshWidth().Y());
+    const int index_global_z = grid.Global().GlobalBegin().Z() + std::floor((p.Pos().Z() - grid.Extent().Begin().Z()) / grid.Extent().MeshWidth().Z());
+
+    assert(index_global_x >= grid.Global().LocalBegin().X() && index_global_x < grid.Global().LocalEnd().X());
+    assert(index_global_y >= grid.Global().LocalBegin().Y() && index_global_y < grid.Global().LocalEnd().Y());
+    assert(index_global_z >= grid.Global().LocalBegin().Z() && index_global_z < grid.Global().LocalEnd().Z());
+
+    const int index_local_x = index_global_x - grid.Global().LocalBegin().X() + grid.Local().HaloSize1().X();
+    const int index_local_y = index_global_y - grid.Global().LocalBegin().Y() + grid.Local().HaloSize1().Y();
+    const int index_local_z = index_global_z - grid.Global().LocalBegin().Z() + grid.Local().HaloSize1().Z();
+
+    assert(index_local_x >= grid.Local().Begin().X() && index_local_x < grid.Local().End().X());
+    assert(index_local_y >= grid.Local().Begin().Y() && index_local_y < grid.Local().End().Y());
+    assert(index_local_z >= grid.Local().Begin().Z() && index_local_z < grid.Local().End().Z());
+
+    const vmg_float pos_beg_x =  grid.Extent().Begin().X() + grid.Extent().MeshWidth().X() * (index_global_x - grid.Global().GlobalBegin().X() - near_field_cells) - p.Pos().X();
+    const vmg_float pos_beg_y =  grid.Extent().Begin().Y() + grid.Extent().MeshWidth().Y() * (index_global_y - grid.Global().GlobalBegin().Y() - near_field_cells) - p.Pos().Y();
+    const vmg_float pos_beg_z =  grid.Extent().Begin().Z() + grid.Extent().MeshWidth().Z() * (index_global_z - grid.Global().GlobalBegin().Z() - near_field_cells) - p.Pos().Z();
+
+    const vmg_float& h_x = grid.Extent().MeshWidth().X();
+    const vmg_float& h_y = grid.Extent().MeshWidth().Y();
+    const vmg_float& h_z = grid.Extent().MeshWidth().Z();
+
+    // Iterate over all grid points which lie in the support of the interpolating B-Spline
+    vmg_float dir_x = pos_beg_x;
+    for (int i=-1*near_field_cells; i<=near_field_cells; ++i) {
+      vmg_float dir_y = pos_beg_y;
+      for (int j=-1*near_field_cells; j<=near_field_cells; ++j) {
+        vmg_float dir_z = pos_beg_z;
+        for (int k=-1*near_field_cells; k<=near_field_cells; ++k) {
+
+          // Compute distance from grid point to particle
+          temp_val = EvaluateSpline(std::sqrt(dir_x*dir_x+dir_y*dir_y+dir_z*dir_z));
+          vals[c++] = temp_val * p.Charge();
+          int_val += temp_val;
+
+          dir_z += h_z;
+        }
+        dir_y += h_y;
+      }
+      dir_x += h_x;
+    }
+
+    // Reciprocal value of the numerically integrated spline
+    int_val = 1.0 / (int_val * h_x * h_y * h_z);
+
+	c = 0;
+    for (int i=-1*near_field_cells; i<=near_field_cells; ++i)
+      for (int j=-1*near_field_cells; j<=near_field_cells; ++j)
+        for (int k=-1*near_field_cells; k<=near_field_cells; ++k)
+          grid(index_local_x + i,
+              index_local_y + j,
+              index_local_z + k) += vals[c++] * int_val;
+    assert( c == Helper::intpow(2*near_field_cells+1,3) );
+
+    delete [] vals;
+  }
+
+  void changeGridBySelfInducedPotential(Grid& grid, Particle& p, const vmg_float &sign) const
+  {
+	assert(p.Pos().X() >= grid.Extent().Begin().X() && p.Pos().X() < grid.Extent().End().X());
+	assert(p.Pos().Y() >= grid.Extent().Begin().Y() && p.Pos().Y() < grid.Extent().End().Y());
+	assert(p.Pos().Z() >= grid.Extent().Begin().Z() && p.Pos().Z() < grid.Extent().End().Z());
+
+	const int index_global_x = grid.Global().GlobalBegin().X() + std::floor((p.Pos().X() - grid.Extent().Begin().X()) / grid.Extent().MeshWidth().X());
+	const int index_global_y = grid.Global().GlobalBegin().Y() + std::floor((p.Pos().Y() - grid.Extent().Begin().Y()) / grid.Extent().MeshWidth().Y());
+	const int index_global_z = grid.Global().GlobalBegin().Z() + std::floor((p.Pos().Z() - grid.Extent().Begin().Z()) / grid.Extent().MeshWidth().Z());
+
+	assert(index_global_x >= grid.Global().LocalBegin().X() && index_global_x < grid.Global().LocalEnd().X());
+	assert(index_global_y >= grid.Global().LocalBegin().Y() && index_global_y < grid.Global().LocalEnd().Y());
+	assert(index_global_z >= grid.Global().LocalBegin().Z() && index_global_z < grid.Global().LocalEnd().Z());
+
+	const int index_local_x = index_global_x - grid.Global().LocalBegin().X() + grid.Local().HaloSize1().X();
+	const int index_local_y = index_global_y - grid.Global().LocalBegin().Y() + grid.Local().HaloSize1().Y();
+	const int index_local_z = index_global_z - grid.Global().LocalBegin().Z() + grid.Local().HaloSize1().Z();
+
+	assert(index_local_x >= grid.Local().Begin().X() && index_local_x < grid.Local().End().X());
+	assert(index_local_y >= grid.Local().Begin().Y() && index_local_y < grid.Local().End().Y());
+	assert(index_local_z >= grid.Local().Begin().Z() && index_local_z < grid.Local().End().Z());
+
+	const vmg_float pos_beg_x =  grid.Extent().Begin().X() + grid.Extent().MeshWidth().X() * (index_global_x - grid.Global().GlobalBegin().X() - near_field_cells) - p.Pos().X();
+	const vmg_float pos_beg_y =  grid.Extent().Begin().Y() + grid.Extent().MeshWidth().Y() * (index_global_y - grid.Global().GlobalBegin().Y() - near_field_cells) - p.Pos().Y();
+	const vmg_float pos_beg_z =  grid.Extent().Begin().Z() + grid.Extent().MeshWidth().Z() * (index_global_z - grid.Global().GlobalBegin().Z() - near_field_cells) - p.Pos().Z();
+
+	const vmg_float& h_x = grid.Extent().MeshWidth().X();
+	const vmg_float& h_y = grid.Extent().MeshWidth().Y();
+	const vmg_float& h_z = grid.Extent().MeshWidth().Z();
+
+	vmg_float dir_x = pos_beg_x;
+    for (int i=-1*near_field_cells; i<=near_field_cells; ++i) {
+  	  vmg_float dir_y = pos_beg_y;
+      for (int j=-1*near_field_cells; j<=near_field_cells; ++j) {
+  		vmg_float dir_z = pos_beg_z;
+        for (int k=-1*near_field_cells; k<=near_field_cells; ++k) {
+  		  const double length_sq = dir_x*dir_x+dir_y*dir_y+dir_z*dir_z;
+		  const double length = std::sqrt(length_sq);
+          grid(index_local_x + i,
+              index_local_y + j,
+              index_local_z + k) += sign * p.Charge() * EvaluatePotential(length);
+
+          dir_z += h_z;
+		}
+	    dir_y += h_y;
+	  }
+	  dir_x += h_x;
+	}
+  }
+
+  void SubtractSelfInducedForces(const Grid& grid, Particle& p) const
+  {
+	assert(p.Pos().X() >= grid.Extent().Begin().X() && p.Pos().X() < grid.Extent().End().X());
+	assert(p.Pos().Y() >= grid.Extent().Begin().Y() && p.Pos().Y() < grid.Extent().End().Y());
+	assert(p.Pos().Z() >= grid.Extent().Begin().Z() && p.Pos().Z() < grid.Extent().End().Z());
+
+	vmg_float temp_val = 0.;
+
+	const int index_global_x = grid.Global().GlobalBegin().X() + std::floor((p.Pos().X() - grid.Extent().Begin().X()) / grid.Extent().MeshWidth().X());
+	const int index_global_y = grid.Global().GlobalBegin().Y() + std::floor((p.Pos().Y() - grid.Extent().Begin().Y()) / grid.Extent().MeshWidth().Y());
+	const int index_global_z = grid.Global().GlobalBegin().Z() + std::floor((p.Pos().Z() - grid.Extent().Begin().Z()) / grid.Extent().MeshWidth().Z());
+
+	assert(index_global_x >= grid.Global().LocalBegin().X() && index_global_x < grid.Global().LocalEnd().X());
+	assert(index_global_y >= grid.Global().LocalBegin().Y() && index_global_y < grid.Global().LocalEnd().Y());
+	assert(index_global_z >= grid.Global().LocalBegin().Z() && index_global_z < grid.Global().LocalEnd().Z());
+
+	const int index_local_x = index_global_x - grid.Global().LocalBegin().X() + grid.Local().HaloSize1().X();
+	const int index_local_y = index_global_y - grid.Global().LocalBegin().Y() + grid.Local().HaloSize1().Y();
+	const int index_local_z = index_global_z - grid.Global().LocalBegin().Z() + grid.Local().HaloSize1().Z();
+
+	assert(index_local_x >= grid.Local().Begin().X() && index_local_x < grid.Local().End().X());
+	assert(index_local_y >= grid.Local().Begin().Y() && index_local_y < grid.Local().End().Y());
+	assert(index_local_z >= grid.Local().Begin().Z() && index_local_z < grid.Local().End().Z());
+
+	const vmg_float pos_beg_x =  grid.Extent().Begin().X() + grid.Extent().MeshWidth().X() * (index_global_x - grid.Global().GlobalBegin().X() - near_field_cells) - p.Pos().X();
+	const vmg_float pos_beg_y =  grid.Extent().Begin().Y() + grid.Extent().MeshWidth().Y() * (index_global_y - grid.Global().GlobalBegin().Y() - near_field_cells) - p.Pos().Y();
+	const vmg_float pos_beg_z =  grid.Extent().Begin().Z() + grid.Extent().MeshWidth().Z() * (index_global_z - grid.Global().GlobalBegin().Z() - near_field_cells) - p.Pos().Z();
+
+	const vmg_float& h_x = grid.Extent().MeshWidth().X();
+	const vmg_float& h_y = grid.Extent().MeshWidth().Y();
+	const vmg_float& h_z = grid.Extent().MeshWidth().Z();
+
+	vmg_float int_val = 0.;
+    vmg_float dir_x = pos_beg_x;
+    for (int i=-1*near_field_cells; i<=near_field_cells; ++i) {
+      vmg_float dir_y = pos_beg_y;
+      for (int j=-1*near_field_cells; j<=near_field_cells; ++j) {
+        vmg_float dir_z = pos_beg_z;
+        for (int k=-1*near_field_cells; k<=near_field_cells; ++k) {
+
+  		  const double length_sq = dir_x*dir_x+dir_y*dir_y+dir_z*dir_z;
+  		  if (fabs(length_sq) > std::numeric_limits<double>::epsilon()) {
+			  // Compute distance from grid point to particle
+			  int_val += EvaluateSpline(std::sqrt(length_sq));
+  		  }
+
+          dir_z += h_z;
+        }
+        dir_y += h_y;
+      }
+      dir_x += h_x;
+    }
+    // Reciprocal value of the numerically integrated spline
+    if (fabs(int_val) > std::numeric_limits<double>::epsilon())
+    	int_val = 1. / (int_val * h_x * h_y * h_z);
+    else
+    	int_val = 1. / (h_x * h_y * h_z);
+
+    // Iterate over all grid points which lie in the support of the interpolating B-Spline
+    vmg_float test_int_val = 0.;
+	dir_x = pos_beg_x;
+	for (int i=-1*near_field_cells; i<=near_field_cells; ++i) {
+	  vmg_float dir_y = pos_beg_y;
+	  for (int j=-1*near_field_cells; j<=near_field_cells; ++j) {
+		vmg_float dir_z = pos_beg_z;
+		for (int k=-1*near_field_cells; k<=near_field_cells; ++k) {
+
+		  // Compute distance from grid point to particle
+		  const double length_sq = dir_x*dir_x+dir_y*dir_y+dir_z*dir_z;
+		  if (fabs(length_sq) > std::numeric_limits<double>::epsilon()) {
+		    const double length = std::sqrt(length_sq);
+
+		    temp_val = h_x * h_y * h_z * int_val * EvaluateSpline(length);
+		    test_int_val += temp_val;
+		    p.Field()[0] -= p.Charge()* temp_val * dir_x / (length_sq*length);
+		    p.Field()[1] -= p.Charge()* temp_val * dir_y / (length_sq*length);
+		    p.Field()[2] -= p.Charge()* temp_val * dir_z / (length_sq*length);
+		  } else {
+			  std::cerr << "Value very small " << length_sq << "=("
+					  << dir_x << ")^2+(" << dir_y << ")^2+(" << dir_z << ")^2 for particle at ("
+					  << p.Pos().X() << ","
+					  << p.Pos().Y() << ","
+					  << p.Pos().Z() << ")"
+					  << "\n";
+		  }
+		  dir_z += h_z;
+		}
+		dir_y += h_y;
+	  }
+	  dir_x += h_x;
+	}
+	if ( fabs(test_int_val -1.) > std::numeric_limits<double>::epsilon()*1e4 ) {
+		std::cerr << "Integrated spline value should be 1 but is " << test_int_val << std::endl;
+		assert(0);
+	}
+  }
+
+  vmg_float EvaluatePotential(const vmg_float& val) const
+  {
+    for (unsigned int i=0; i<intervals.size(); ++i)
+      if (val < intervals[i])
+        return potential_nom[i](val) / potential_denom[i](val);
+    return potential_nom.back()(val) / potential_denom.back()(val);
+  }
+
+  vmg_float EvaluateField(const vmg_float& val) const
+  {
+    for (unsigned int i=0; i<intervals.size(); ++i)
+      if (val < intervals[i])
+        return field_nom[i](val) / field_denom[i](val);
+    return 0.0;
+  }
+
+  const vmg_float& GetAntiDerivativeAtZero() const
+  {
+    return antid;
+  }
+
+private:
+  std::vector<Polynomial> spline_nom, spline_denom;
+  std::vector<Polynomial> potential_nom, potential_denom;
+  std::vector<Polynomial> field_nom, field_denom;
+  vmg_float antid;
+  std::vector<vmg_float> intervals;
+
+  const vmg_float R;
+  const int near_field_cells;
+};
+
+}
+
+}
+
+#endif /* BSPLINE_HPP_ */
