Adding missing files

src/interpolation/interpolation.hpp

+/*
+ * mp4_interpolation.hpp
+ *
+ *  Created on: May 4, 2017
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_HPP_
+#define OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_HPP_
+
+#include "NN/CellList/MemFast.hpp"
+#include "NN/CellList/CellList.hpp"
+
+#define INTERPOLATION_ERROR_OBJECT std::runtime_error("Runtime interpolation error");
+
+constexpr int inte_m2p = 0;
+constexpr int inte_p2m = 1;
+
+template<typename T>
+struct mul_inte
+{
+	inline static void value(T & result, const T & coeff, const T & src)
+	{
+		result += coeff * src;
+	}
+};
+
+template<typename T, unsigned int N1>
+struct mul_inte<T[N1]>
+{
+	inline static void value(T (& result)[N1], const T & coeff, const T (& src)[N1])
+	{
+		for (size_t i = 0 ; i < N1 ; i++)
+			result[i] += coeff * src[i];
+	}
+};
+
+template<typename T, unsigned int N1, unsigned int N2>
+struct mul_inte<T[N1][N2]>
+{
+	inline static void value(T (& result)[N1][N2], const T & coeff, const T (& src)[N1][N2])
+	{
+		for (size_t i = 0 ; i < N1 ; i++)
+			for (size_t j = 0 ; j < N2 ; j++)
+				result[i][j] += coeff * src[i][j];
+	}
+};
+
+template<unsigned int prp_g, unsigned int prp_v,unsigned int m2p_or_p2m>
+struct inte_template
+{
+	template<typename grid, typename vector, typename iterator, typename key_type> inline static void value(grid & gd,
+			                                                          vector & vd,
+																	  const grid_dist_key_dx<vector::dims> & k_dist,
+																	  iterator & key_p,
+																	  const grid_cpu<vector::dims,aggregate<typename vector::stype>> & a_int,
+																	  const key_type & key)
+	{
+		mul_inte<typename std::remove_reference<decltype(gd.template get<prp_g>(k_dist))>::type>::value(gd.template get<prp_g>(k_dist),a_int.template get<0>(key),vd.template getProp<prp_v>(key_p));
+	}
+};
+
+template<unsigned int prp_g, unsigned int prp_v>
+struct inte_template<prp_g,prp_v,inte_m2p>
+{
+	template<typename grid, typename vector, typename iterator, typename key_type> inline static void value(grid & gd,
+			                                                          vector & vd,
+																	  const grid_dist_key_dx<vector::dims> & k_dist,
+																	  iterator & key_p,
+																	  const grid_cpu<vector::dims,aggregate<typename vector::stype>> & a_int,
+																	  const key_type & key)
+	{
+		mul_inte<typename std::remove_reference<decltype(gd.template get<prp_g>(k_dist))>::type>::value(vd.template getProp<prp_v>(key_p),a_int.template get<0>(key),gd.template get<prp_g>(k_dist));
+	}
+};
+
+
+/*! \brief return the sub-domain where this particle must be interpolated
+ *
+ * \param p particle position
+ *
+ * \return the sub-domain id
+ *
+ */
+template<typename vector, typename grid>
+inline size_t getSub(Point<vector::dims,typename vector::stype> & p,
+		             const CellList<vector::dims,typename vector::stype,Mem_fast<vector::dims,typename vector::stype>,shift<vector::dims,typename vector::stype>> & geo_cell,
+					 grid & gd)
+{
+	size_t cell = geo_cell.getCell(p);
+
+	for (size_t i = 0 ; i < geo_cell.getNelements(cell) ; i++)
+	{
+		size_t ns = geo_cell.get(cell,i);
+
+		if (gd.getDecomposition().getSubDomain(ns).isInside(p))
+			return ns;
+	}
+
+#ifdef SE_CLASS1
+	std::cerr << __FILE__ << ":" << __LINE__ << " Error: " << std::endl;
+	ACTION_ON_ERROR(INTERPOLATION_ERROR_OBJECT)
+#endif
+
+	return (size_t)-1;
+}
+
+
+template<unsigned int prp_g, unsigned int prp_v, unsigned int m2p_or_p2m, typename kernel, typename iterator, typename vector, typename grid, typename grid_inte>
+																 inline void inte_calc(iterator & it,
+		                     	 	 	 	 	 	 	 	 	 const vector & vd,
+																 const Box<vector::dims,typename vector::stype> & domain,
+																 int (& ip)[vector::dims][kernel::np],
+																 grid & gd,
+																 const typename vector::stype (& dx)[vector::dims],
+																 typename vector::stype (& xp)[vector::dims],
+																 const grid_inte & a_int,
+																 typename vector::stype (& a)[vector::dims][kernel::np],
+																 typename vector::stype (& x)[vector::dims][kernel::np],
+																 size_t (& sz)[vector::dims],
+																 const CellList<vector::dims,typename vector::stype,Mem_fast<vector::dims,typename vector::stype>,shift<vector::dims,typename vector::stype>> & geo_cell)
+{
+	auto key_p = it.get();
+
+	Point<vector::dims,typename vector::stype> p = vd.getPos(key_p);
+
+	// On which sub-domain we interpolate the particle
+	size_t sub = getSub<vector>(p,geo_cell,gd);
+
+	typename vector::stype x0[vector::dims];
+
+	// calculate the position of the particle in the grid
+	// coordinated
+	for (size_t i = 0 ; i < vector::dims ; i++)
+		x0[i] = (p.get(i)-domain.getLow(i))*dx[i];
+
+	// convert into integer
+	for (size_t i = 0 ; i < vector::dims ; i++)
+		ip[i][0] = (int)x0[i];
+
+	// convert the global grid position into local grid position
+	grid_key_dx<vector::dims> base;
+
+	for (size_t i = 0 ; i < vector::dims ; i++)
+		base.set_d(i,ip[i][0] - gd.getLocalGridsInfo().get(sub).origin.get(i) - (long int)kernel::np/2 + 1);
+
+	// convenient grid of points
+
+	for (size_t j = 0 ; j < kernel::np-1 ; j++)
+	{
+		for (size_t i = 0 ; i < vector::dims ; i++)
+			ip[i][j+1] = (int)ip[i][j]+1;
+	}
+
+	for (size_t i = 0 ; i < vector::dims ; i++)
+		xp[i] = x0[i] - ip[i][0];
+
+	for (long int j = 0 ; j < kernel::np ; j++)
+	{
+		for (size_t i = 0 ; i < vector::dims ; i++)
+			x[i][j] = - xp[i] + typename vector::stype((long int)j - (long int)kernel::np/2 + 1);
+	}
+
+	for (size_t j = 0 ; j < kernel::np ; j++)
+	{
+		for (size_t i = 0 ; i < vector::dims ; i++)
+			a[i][j] = kernel::value(x[i][j],j);
+	}
+
+	grid_sm<vector::dims,void> gs(sz);
+	grid_key_dx_iterator<vector::dims> kit(gs);
+
+	while (kit.isNext())
+	{
+		auto key = kit.get();
+
+		a_int.template get<0>(key) = 1;
+
+		for (size_t i = 0 ; i < vector::dims ; i++)
+			a_int.template get<0>(key) *= a[i][key.get(i)];
+
+		++kit;
+	}
+
+	grid_key_dx_iterator<vector::dims> kit2(gs);
+	grid_dist_key_dx<vector::dims> k_dist;
+	k_dist.setSub(sub);
+
+	while (kit2.isNext())
+	{
+		auto key = kit2.get();
+
+		for (size_t i = 0 ; i < vector::dims ; i++)
+			k_dist.getKeyRef().set_d(i,key.get(i) + base.get(i));
+
+		inte_template<prp_g,prp_v,m2p_or_p2m>::value(gd,vd,k_dist,key_p,a_int,key);
+
+		++kit2;
+	}
+}
+
+
+template<typename vector,typename grid, typename kernel>
+class interpolate
+{
+	CellList<vector::dims,typename vector::stype,Mem_fast<vector::dims,typename vector::stype>,shift<vector::dims,typename vector::stype>> geo_cell;
+
+	struct Box_vol
+	{
+		Box<vector::dims,size_t> bv;
+
+		size_t vol;
+
+		void operator<(const Box_vol & bv)
+		{
+			return vol < bv.vol;
+		}
+	};
+
+	//! particles
+	vector & vd;
+
+	//! grid
+	grid & gd;
+
+
+public:
+
+	interpolate(vector & vd, grid & gd)
+	:vd(vd),gd(gd)
+	{
+		// get the processor bounding box in grid units
+		Box<vector::dims,typename vector::stype> bb = gd.getDecomposition().getProcessorBounds();
+		Box<vector::dims,typename vector::stype> bunit = gd.getDecomposition().getCellDecomposer().getCellBox();
+
+		size_t div[vector::dims];
+
+		for (size_t i = 0 ; i < vector::dims ; i++)
+			div[i] = (bb.getHigh(i) - bb.getLow(i)) / bunit.getHigh(i);
+
+		geo_cell.Initialize(bb,div);
+
+		// Now draw the domain into the cell list
+
+		auto & dec = gd.getDecomposition();
+
+		for (size_t i = 0 ; i < dec.getNSubDomain() ; i++)
+		{
+			const Box<vector::dims,typename vector::stype> & bx = dec.getSubDomain(i);
+
+			// get the cells this box span
+			const grid_key_dx<vector::dims> p1 = geo_cell.getCellGrid(bx.getP1());
+			const grid_key_dx<vector::dims> p2 = geo_cell.getCellGrid(bx.getP2());
+
+			// Get the grid and the sub-iterator
+			auto & gi = geo_cell.getGrid();
+			grid_key_dx_iterator_sub<vector::dims> g_sub(gi,p1,p2);
+
+			// add the box-id to the cell list
+			while (g_sub.isNext())
+			{
+				auto key = g_sub.get();
+				geo_cell.addCell(gi.LinId(key),i);
+				++g_sub;
+			}
+		}
+	};
+
+	template<unsigned int prp_v, unsigned int prp_g> void p2m(vector & vd, grid & gd)
+	{
+#ifdef SE_CLASS1
+
+		if (!vd.getDecomposition().is_equal_ng(gd.getDecomposition()) )
+		{
+			std::cerr << __FILE__ << ":" << __LINE__ << " Error: the distribution of the vector of particles" <<
+					" and the grid is different. In order to interpolate the two data structure must have the" <<
+					" same decomposition" << std::endl;
+
+			ACTION_ON_ERROR(INTERPOLATION_ERROR_OBJECT)
+		}
+
+#endif
+
+		Box<vector::dims,typename vector::stype> domain = vd.getDecomposition().getDomain();
+
+		// grid spacing
+		typename vector::stype dx[vector::dims];
+
+		for (size_t i = 0 ; i < vector::dims ; i++)
+			dx[i] = 1.0/gd.spacing(i);
+
+		// point position
+		typename vector::stype xp[vector::dims];
+
+		int ip[vector::dims][kernel::np];
+		typename vector::stype x[vector::dims][kernel::np];
+		typename vector::stype a[vector::dims][kernel::np];
+
+		size_t sz[vector::dims];
+
+		for (size_t i = 0 ; i < vector::dims ; i++)
+			sz[i] = kernel::np;
+
+		grid_cpu<vector::dims,aggregate<typename vector::stype>> a_int(sz);
+		a_int.setMemory();
+
+		auto it = vd.getDomainIterator();
+
+		while (it.isNext() == true)
+		{
+			inte_calc<prp_g,prp_v,inte_p2m,kernel>(it,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell);
+
+			++it;
+		}
+	}
+
+	template<unsigned int prp_v, unsigned int prp_g> void m2p(grid & gd, vector & vd)
+	{
+#ifdef SE_CLASS1
+
+		if (!vd.getDecomposition().is_equal_ng(gd.getDecomposition()) )
+		{
+			std::cerr << __FILE__ << ":" << __LINE__ << " Error: the distribution of the vector of particles" <<
+					" and the grid is different. In order to interpolate the two data structure must have the" <<
+					" same decomposition" << std::endl;
+
+			ACTION_ON_ERROR(INTERPOLATION_ERROR_OBJECT)
+		}
+
+#endif
+
+		Box<vector::dims,typename vector::stype> domain = vd.getDecomposition().getDomain();
+
+		// grid spacing
+		typename vector::stype dx[vector::dims];
+
+		for (size_t i = 0 ; i < vector::dims ; i++)
+			dx[i] = 1.0/gd.spacing(i);
+
+		// point position
+		typename vector::stype xp[vector::dims];
+
+		int ip[vector::dims][kernel::np];
+		typename vector::stype x[vector::dims][kernel::np];
+		typename vector::stype a[vector::dims][kernel::np];
+
+		size_t sz[vector::dims];
+
+		for (size_t i = 0 ; i < vector::dims ; i++)
+			sz[i] = kernel::np;
+
+		grid_cpu<vector::dims,aggregate<typename vector::stype>> a_int(sz);
+		a_int.setMemory();
+
+		auto it = vd.getDomainIterator();
+
+		while (it.isNext() == true)
+		{
+			inte_calc<prp_g,prp_v,inte_m2p,kernel>(it,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell);
+
+			++it;
+		}
+	}
+
+};
+
+#endif /* OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_HPP_ */

src/interpolation/interpolation_unit_tests.hpp

+/*
+ * interpolation_unit_tests.hpp
+ *
+ *  Created on: May 5, 2017
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_UNIT_TESTS_HPP_
+#define OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_UNIT_TESTS_HPP_
+
+#include "interpolation/mp4_kernel.hpp"
+#include "interpolation/interpolation.hpp"
+#include "interpolation/z_spline.hpp"
+
+BOOST_AUTO_TEST_SUITE( interpolation_test )
+
+template<typename grid, unsigned int mom_p> void momenta_grid(grid & gd,typename grid::stype (& mom_tot)[grid::dims])
+{
+	auto it = gd.getDomainGhostIterator();
+
+	for (size_t i = 0 ; i < grid::dims ; i++)
+		mom_tot[i] = 0.0;
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+		auto key_g = gd.getGKey(key);
+
+		for (size_t i = 0 ; i < grid::dims ; i++)
+		{
+			typename grid::stype coord = gd.spacing(i)*key_g.get(i);
+
+			mom_tot[i] += boost::math::pow<mom_p>(coord)*gd.template getProp<0>(key);
+		}
+
+		++it;
+	}
+}
+
+template<typename grid, unsigned int mom_p> void momenta_grid_domain(grid & gd,typename grid::stype (& mom_tot)[grid::dims])
+{
+	auto it = gd.getDomainIterator();
+
+	for (size_t i = 0 ; i < grid::dims ; i++)
+		mom_tot[i] = 0.0;
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+		auto key_g = gd.getGKey(key);
+
+		for (size_t i = 0 ; i < grid::dims ; i++)
+		{
+			typename grid::stype coord = gd.spacing(i)*key_g.get(i);
+
+			mom_tot[i] += boost::math::pow<mom_p>(coord)*gd.template getProp<0>(key);
+		}
+
+		++it;
+	}
+}
+
+template<typename vector, unsigned int mom_p> void momenta_vector(vector & vd,typename vector::stype (& mom_tot)[vector::dims])
+{
+	auto it = vd.getDomainIterator();
+
+	for (size_t i = 0 ; i < vector::dims ; i++)
+		mom_tot[i] = 0.0;
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		for (size_t i = 0 ; i < vector::dims ; i++)
+		{
+			typename vector::stype coord = vd.getPos(key)[i];
+
+			mom_tot[i] += boost::math::pow<mom_p>(coord)*vd.template getProp<0>(key);
+		}
+
+		++it;
+	}
+}
+
+
+BOOST_AUTO_TEST_CASE( interpolation_full_test )
+{
+	Box<2,float> domain({0.0,0.0},{1.0,1.0});
+	size_t sz[2] = {64,64};
+
+	Ghost<2,long int> gg(3);
+	Ghost<2,float> gv(0.01);
+
+	size_t bc_v[2] = {PERIODIC,PERIODIC};
+
+	{
+	vector_dist<2,float,aggregate<float>> vd(4096,domain,bc_v,gv);
+	grid_dist_id<2,float,aggregate<float>> gd(vd.getDecomposition(),sz,gg);
+
+	// set one particle on vd
+
+	auto it = vd.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto p = it.get();
+
+		vd.getPos(p)[0] = (double)rand()/RAND_MAX;
+		vd.getPos(p)[1] = (double)rand()/RAND_MAX;
+		vd.getPos(p)[2] = (double)rand()/RAND_MAX;
+
+		vd.getProp<0>(p) = 5.0/*(double)rand()/RAND_MAX*/;
+
+		++it;
+	}
+
+	vd.map();
+
+	// Reset the grid
+
+	auto it2 = gd.getDomainGhostIterator();
+
+	while (it2.isNext())
+	{
+		auto key = it2.get();
+
+		gd.template get<0>(key) = 0.0;
+
+		++it2;
+	}
+
+	interpolate<decltype(vd),decltype(gd),mp4_kernel<float>> inte(vd,gd);
+
+	inte.p2m<0,0>(vd,gd);
+
+	float mg[2];
+	float mv[2];
+
+	momenta_grid<decltype(gd),0>(gd,mg);
+	momenta_vector<decltype(vd),0>(vd,mv);
+
+	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
+	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+
+	momenta_grid<decltype(gd),1>(gd,mg);
+	momenta_vector<decltype(vd),1>(vd,mv);
+
+	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
+	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+
+	momenta_grid<decltype(gd),2>(gd,mg);
+	momenta_vector<decltype(vd),2>(vd,mv);
+
+	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
+	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+
+	// We have to do a ghost get before interpolating m2p
+	// Before doing mesh to particle particle must be arranged
+	// into a grid like
+
+	vd.clear();
+
+	auto it4 = vd.getGridIterator(sz);
+
+	while(it4.isNext())
+	{
+		auto key = it4.get();
+
+		vd.add();
+
+		vd.getLastPos()[0] = key.get(0) * it4.getSpacing(0) + domain.getLow(0) + 0.1*it4.getSpacing(0);
+		vd.getLastPos()[1] = key.get(1) * it4.getSpacing(1) + domain.getLow(1) + 0.1*it4.getSpacing(1);
+
+		vd.getLastProp<0>() = 0.0;
+
+		++it4;
+	}
+
+	// Reset also the grid
+
+	auto it5 = gd.getDomainGhostIterator();
+
+	while(it5.isNext())
+	{
+		auto key = it5.get();
+
+		gd.get<0>(key) = 0.0;
+
+		++it5;
+	}
+	gd.ghost_get<0>();
+
+	grid_key_dx<2> start({3,3});
+	grid_key_dx<2> stop({(long int)gd.size(0) - 3,(long int)gd.size(1) - 3});
+
+	auto it6 = gd.getSubDomainIterator(start,stop);
+	while(it6.isNext())
+	{
+		auto key = it6.get();
+
+		gd.get<0>(key) = (double)rand()/RAND_MAX;;
+
+		++it6;
+	}
+	gd.ghost_get<0>();
+
+	vd.map();
+	gd.ghost_get<0>();
+	inte.m2p<0,0>(gd,vd);
+
+	momenta_grid_domain<decltype(gd),0>(gd,mg);
+	momenta_vector<decltype(vd),0>(vd,mv);
+
+	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
+	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+
+	momenta_grid_domain<decltype(gd),1>(gd,mg);
+	momenta_vector<decltype(vd),1>(vd,mv);
+
+	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
+	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+
+	momenta_grid_domain<decltype(gd),2>(gd,mg);
+	momenta_vector<decltype(vd),2>(vd,mv);
+
+	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
+	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+
+	}
+}
+
+BOOST_AUTO_TEST_CASE( int_kernel_test )
+{
+		mp4_kernel<float> mp4;
+
+		float tot = 0.0;
+
+		// Check momenta 0
+
+		tot += mp4.value(-1.3f,0);
+		tot += mp4.value(-0.3f,1);
+		tot += mp4.value(0.7f,2);
+		tot += mp4.value(1.7f,3);
+
+		BOOST_REQUIRE_CLOSE(tot,1.0f,0.001);
+
+		// Check momenta 1
+
+		tot = 0.0;
+
+		tot += -1.3f*mp4.value(-1.3f,0);
+		tot += -0.3f*mp4.value(-0.3f,1);
+		tot += 0.7f*mp4.value(0.7f,2);
+		tot += 1.7f*mp4.value(1.7f,3);
+
+		BOOST_REQUIRE_SMALL(tot,0.001f);
+
+		// Check momenta 2
+
+		tot = 0.0;
+
+		tot += (1.3f)*(1.3f)*mp4.value(-1.3f,0);
+		tot += (0.3f)*(0.3f)*mp4.value(-0.3f,1);
+		tot += (0.7f)*(0.7f)*mp4.value(0.7f,2);
+		tot += (1.7f)*(1.7f)*mp4.value(1.7f,3);
+
+		BOOST_REQUIRE_SMALL(tot,0.001f);
+
+
+		//////// Check zeta 1
+
+		tot = 0.0;
+
+		z_kernel<float,1> zk1;
+
+		tot += zk1.value(-0.3f,0);
+		tot += zk1.value(0.7f,1);
+
+		BOOST_REQUIRE_CLOSE(tot,1.0f,0.001);
+
+		//////// zeta 2 is equivalent to mp4 we do not test
+
+		//////// zeta 3
+
+		z_kernel<float,3> zk3;
+
+		tot = 0.0;
+
+		// Check momenta 0
+
+		tot += zk3.value(-2.3f,0);
+		tot += zk3.value(-1.3f,1);
+		tot += zk3.value(-0.3f,2);
+		tot += zk3.value(0.7f,3);
+		tot += zk3.value(1.7f,4);
+		tot += zk3.value(2.7f,5);
+
+		BOOST_REQUIRE_CLOSE(tot,1.0f,0.001);
+
+		// Check momenta 1
+
+		tot = 0.0;
+
+		tot += -2.3*zk3.value(-2.3f,0);
+		tot += -1.3*zk3.value(-1.3f,1);
+		tot += -0.3*zk3.value(-0.3f,2);
+		tot += 0.7*zk3.value(0.7f,3);
+		tot += 1.7*zk3.value(1.7f,4);
+		tot += 2.7*zk3.value(2.7f,5);
+
+		BOOST_REQUIRE_SMALL(tot,0.001f);
+
+		// Check momenta 2
+
+		tot = 0.0;
+
+		tot += 2.3*2.3*zk3.value(-2.3f,0);
+		tot += 1.3*1.3*zk3.value(-1.3f,1);
+		tot += 0.3*0.3*zk3.value(-0.3f,2);
+		tot += 0.7*0.7*zk3.value(0.7f,3);
+		tot += 1.7*1.7*zk3.value(1.7f,4);
+		tot += 2.7*2.7*zk3.value(2.7f,5);
+
+		BOOST_REQUIRE_SMALL(tot,0.001f);
+
+		// Check momenta 3
+
+		tot = 0.0;
+
+		tot += -2.3*-2.3*-2.3*zk3.value(-2.3f,0);
+		tot += -1.3*-1.3*-1.3*zk3.value(-1.3f,1);
+		tot += -0.3*-0.3*-0.3*zk3.value(-0.3f,2);
+		tot += 0.7*0.7*0.7*zk3.value(0.7f,3);
+		tot += 1.7*1.7*1.7*zk3.value(1.7f,4);
+		tot += 2.7*2.7*2.7*zk3.value(2.7f,5);
+
+		BOOST_REQUIRE_SMALL(tot,0.001f);
+
+
+		// z4
+
+		z_kernel<float,4> zk4;
+
+		// Check momenta 0
+
+		tot = 0.0;
+
+		tot += zk4.value(-3.3f,0);
+		tot += zk4.value(-2.3f,1);
+		tot += zk4.value(-1.3f,2);
+		tot += zk4.value(-0.3f,3);
+		tot += zk4.value(0.7f,4);
+		tot += zk4.value(1.7f,5);
+		tot += zk4.value(2.7f,6);
+		tot += zk4.value(3.7f,7);
+
+		BOOST_REQUIRE_CLOSE(tot,1.0f,0.001);
+
+		// Check momenta 1
+
+		tot = 0.0;
+
+		tot += -3.3*zk4.value(-3.3f,0);
+		tot += -2.3*zk4.value(-2.3f,1);
+		tot += -1.3*zk4.value(-1.3f,2);
+		tot += -0.3*zk4.value(-0.3f,3);
+		tot += 0.7*zk4.value(0.7f,4);
+		tot += 1.7*zk4.value(1.7f,5);
+		tot += 2.7*zk4.value(2.7f,6);
+		tot += 3.7*zk4.value(3.7f,7);
+
+		BOOST_REQUIRE_SMALL(tot,0.001f);
+
+		// Check momenta 2
+
+		tot = 0.0;
+
+		tot += 3.3*3.3*zk4.value(-3.3f,0);
+		tot += 2.3*2.3*zk4.value(-2.3f,1);
+		tot += 1.3*1.3*zk4.value(-1.3f,2);
+		tot += 0.3*0.3*zk4.value(-0.3f,3);
+		tot += 0.7*0.7*zk4.value(0.7f,4);
+		tot += 1.7*1.7*zk4.value(1.7f,5);
+		tot += 2.7*2.7*zk4.value(2.7f,6);
+		tot += 3.7*3.7*zk4.value(3.7f,7);
+
+		BOOST_REQUIRE_SMALL(tot,0.001f);
+
+		// Check momenta 3
+
+		tot = 0.0;
+
+		tot += -3.3*-3.3*-3.3*zk4.value(-3.3f,0);
+		tot += -2.3*-2.3*-2.3*zk4.value(-2.3f,1);
+		tot += -1.3*-1.3*-1.3*zk4.value(-1.3f,2);
+		tot += -0.3*-0.3*-0.3*zk4.value(-0.3f,3);
+		tot += 0.7*0.7*0.7*zk4.value(0.7f,4);
+		tot += 1.7*1.7*1.7*zk4.value(1.7f,5);
+		tot += 2.7*2.7*2.7*zk4.value(2.7f,6);
+		tot += 3.7*3.7*3.7*zk4.value(3.7f,7);
+
+		BOOST_REQUIRE_SMALL(tot,0.001f);
+
+		// Check momenta 4
+
+		tot = 0.0;
+
+		tot += -3.3*-3.3*-3.3*-3.3*zk4.value(-3.3f,0);
+		tot += -2.3*-2.3*-2.3*-2.3*zk4.value(-2.3f,1);
+		tot += -1.3*-1.3*-1.3*-1.3*zk4.value(-1.3f,2);
+		tot += -0.3*-0.3*-0.3*-0.3*zk4.value(-0.3f,3);
+		tot += 0.7*0.7*0.7*0.7*zk4.value(0.7f,4);
+		tot += 1.7*1.7*1.7*1.7*zk4.value(1.7f,5);
+		tot += 2.7*2.7*2.7*2.7*zk4.value(2.7f,6);
+		tot += 3.7*3.7*3.7*3.7*zk4.value(3.7f,7);
+
+		BOOST_REQUIRE_SMALL(tot,0.001f);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+
+#endif /* OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_UNIT_TESTS_HPP_ */

src/interpolation/z_spline.hpp

+/*
+ * z_spline.hpp
+ *
+ *  Created on: May 8, 2017
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_NUMERICS_SRC_INTERPOLATION_Z_SPLINE_HPP_
+#define OPENFPM_NUMERICS_SRC_INTERPOLATION_Z_SPLINE_HPP_
+
+#include "mp4_kernel.hpp"
+
+template<typename st, unsigned int ord>
+class z_kernel
+{
+public:
+
+	static const int np = 4;
+
+	static inline st value(st x, size_t i)
+	{
+		std::cerr << __FILE__ << ":" << __LINE__ << ": Error this order has not been implemented" << std::endl;
+		return 0.0;
+	}
+};
+
+template<typename st>
+class z_kernel<st,1>
+{
+public:
+
+	static const int np = 2;
+
+	static inline st value(st x, size_t i)
+	{
+		if (i == 0)
+			return x + 1;
+		else if (i == 1)
+			return 1-x;
+		return 0.0;
+	}
+};
+
+template<typename st>
+class z_kernel<st,2>: public mp4_kernel<st>
+{
+
+};
+
+template<typename st>
+class z_kernel<st,3>
+{
+public:
+
+	static const int np = 6;
+
+	static inline st value(st x, size_t i)
+	{
+		if (i == 0)
+			return 18.0 + ( 38.25 + (31.875 + ( 313.0/24.0 + (2.625 + 5.0/24.0*x)*x)*x)*x)*x;
+		else if (i == 1)
+			return -4.0 + (-18.75 + (-30.625 + (-545.0/24.0 + ( -7.875 - 25.0/24.0*x)*x)*x)*x)*x;
+		else if (i == 2)
+			return 1.0 + (-1.25 +(35.0/12.0 +(5.25 + 25.0/12.0*x)*x)*x)*x*x;
+		else if (i == 3)
+			return 1.0 + (-1.25 +(-35.0/12.0 +(5.25 - 25.0/12.0*x)*x)*x)*x*x;
+		else if (i == 4)
+			return -4.0 + (18.75 + (-30.625 + (545.0/24.0 + ( -7.875 + 25.0/24.0*x)*x)*x)*x)*x;
+		else if (i == 5)
+			return 18.0 + (-38.25 + (31.875 + (-313.0/24.0 + (2.625 - 5.0/24.0*x)*x)*x)*x)*x;
+
+		return 0.0;
+	}
+};
+
+template<typename st>
+class z_kernel<st,4>
+{
+public:
+
+	static const int np = 8;
+
+	static inline st value(st x, size_t i)
+	{
+		if (i == 0)
+			return 726.4 + ( 1491.2 + (58786.0/45.0 + ( 633.0 + (26383.0/144.0 + (22807.0/720.0 + (727.0/240.0 + 89.0/720.0*x)*x)*x)*x)*x)*x)*x;
+		else if (i == 1)
+			return -440 +( -1297.45 + ( -117131/72.0 + ( -1123.5 + ( -66437.0/144.0 + ( -81109.0/720.0 + ( -727.0/48.0 - 623.0/720.0*x)*x)*x)*x)*x)*x)*x;
+		else if (i == 2)
+			return 27.6 + (8617.0/60.0 +(321.825 +(395.5 +( 284.8125 + (119.7875 +(27.2625 + 623.0/240.0*x)*x)*x)*x)*x)*x)*x;
+		else if (i == 3)
+			return 1.0 + (( -49.0/36.0 + (( -959.0/144.0 + ( -2569.0/144.0 + ( -727.0/48.0 - 623.0/144.0*x)*x)*x)*x)*x)*x)*x;
+		else if (i == 4)
+			return 1.0 + (( -49.0/36.0 + (( -959.0/144.0 + ( 2569.0/144.0 + ( -727.0/48.0 + 623.0/144.0*x)*x)*x)*x)*x)*x)*x;
+		else if (i == 5)
+			return 27.6 + (-8617.0/60.0 +(321.825 +(-395.5 +( 284.8125 + (-119.7875 +(27.2625 - 623.0/240.0*x)*x)*x)*x)*x)*x)*x;
+		else if (i == 6)
+			return -440 +( 1297.45 + ( -117131/72.0 + ( 1123.5 + ( -66437.0/144.0 + ( 81109.0/720.0 + ( -727.0/48.0 + 623.0/720.0*x)*x)*x)*x)*x)*x)*x;
+		else if (i == 7)
+			return 726.4 + ( -1491.2 + (58786.0/45.0 + ( -633.0 + (26383.0/144.0 + (-22807.0/720.0 + (727.0/240.0 - 89.0/720.0*x)*x)*x)*x)*x)*x)*x;
+
+		return 0.0;
+	}
+};
+
+#endif /* OPENFPM_NUMERICS_SRC_INTERPOLATION_Z_SPLINE_HPP_ */