From 219f53f09cae9c41c8390a02e9d3bcf55efcfaa2 Mon Sep 17 00:00:00 2001
From: Pietro Incardona <incardon@mpi-cbg.de>
Date: Tue, 11 Aug 2015 03:00:14 +0200
Subject: [PATCH] Compilable VTKWriter.hpp

---
 src/CSVWriter.hpp            |   2 +-
 src/GraphMLWriter.hpp        |   4 +-
 src/VTKWriter.hpp            | 146 +------------
 src/VTKWriter_graph.hpp      | 148 ++++++++++++-
 src/VTKWriter_grids.hpp      | 406 +++++++++++++++++++++++++++++++++++
 src/VTKWriter_unit_tests.hpp |  81 +++++++
 6 files changed, 638 insertions(+), 149 deletions(-)
 create mode 100644 src/VTKWriter_grids.hpp

diff --git a/src/CSVWriter.hpp b/src/CSVWriter.hpp
index 9085c5c..7edaa41 100644
--- a/src/CSVWriter.hpp
+++ b/src/CSVWriter.hpp
@@ -158,7 +158,7 @@ class CSVWriter
 		std::stringstream str;
 
 		// write positional columns
-		for (int i = 0 ; i < v_pos::value_type::dims ; i++)
+		for (size_t i = 0 ; i < v_pos::value_type::dims ; i++)
 		{
 			if (i == 0)
 				str << "x[" << i << "]";
diff --git a/src/GraphMLWriter.hpp b/src/GraphMLWriter.hpp
index 5c8eb18..e6a1728 100644
--- a/src/GraphMLWriter.hpp
+++ b/src/GraphMLWriter.hpp
@@ -23,7 +23,7 @@ void create_prop(std::string * str)
 	if (has_attributes<T>::value )
 	{
 		// Create properties names based on the attributes name defined
-		for (int i = 0 ; i < T::max_prop ; i++)
+		for (size_t i = 0 ; i < T::max_prop ; i++)
 		{
 			str[i] = std::string(T::attributes::name[i]);
 		}
@@ -31,7 +31,7 @@ void create_prop(std::string * str)
 	else
 	{
 		// Create default properties name
-		for (int i = 0 ; i < T::max_prop ; i++)
+		for (size_t i = 0 ; i < T::max_prop ; i++)
 		{
 			str[i] = "attr" + std::to_string(i);
 		}
diff --git a/src/VTKWriter.hpp b/src/VTKWriter.hpp
index a14386e..8406141 100644
--- a/src/VTKWriter.hpp
+++ b/src/VTKWriter.hpp
@@ -93,152 +93,9 @@ enum file_type
 	ASCII
 };
 
-/*! \brief this class is a functor for "for_each" algorithm
- *
- * This class is a functor for "for_each" algorithm. For each
- * element of the boost::vector the operator() is called.
- * Is mainly used to create a string containing all the vertex
- * properties
- *
- * \tparam G graph type
- * \tparam attr has the vertex attributes
- *
- */
-
-template<typename G, bool attr>
-struct vtk_vertex_node
-{
-	// Vertex spatial type information
-	typedef typename G::V_type::s_type s_type;
-
-	s_type (& x)[3];
-
-	// Vertex object container
-	typename G::V_container & vo;
-
-	// vertex node string
-	std::string & v_node;
-
-	/*! \brief Constructor
-	 *
-	 * Create a vertex properties list
-	 *
-	 * \param v_node std::string that is filled with the graph properties in the GraphML format
-	 * \param n_obj object container to access its properties for example encapc<...>
-	 *
-	 */
-	vtk_vertex_node(std::string & v_node, typename G::V_container & n_obj, s_type (&x)[3])
-	:x(x),vo(n_obj),v_node(v_node)
-	{
-	};
-
-	//! \brief Write collected information
-	void write()
-	{
-		v_node += std::to_string(x[0]) + " " + std::to_string(x[1]) + " " + std::to_string(x[2]) + "\n";
-	}
-
-	//! It call the functor for each member
-    template<typename T>
-    void operator()(T& t)
-    {
-    	// if the attribute name is x y or z, create a string with the value of the properties and store it
-		if (G::V_type::attributes::name[T::value] == "x"){x[0] = convert<typename boost::remove_reference<decltype(vo.template get<T::value>())>::type>::template to<s_type>(vo.template get<T::value>());}
-		else if (G::V_type::attributes::name[T::value] == "y"){x[1] = convert<typename boost::remove_reference<decltype(vo.template get<T::value>())>::type>::template to<s_type>(vo.template get<T::value>());}
-		else if (G::V_type::attributes::name[T::value] == "z"){x[2] = convert<typename boost::remove_reference<decltype(vo.template get<T::value>())>::type>::template to<s_type>(vo.template get<T::value>());}
-    }
-};
-
-/*! \brief this class is a functor for "for_each" algorithm
- *
- * This class is a functor for "for_each" algorithm. For each
- * element of the boost::vector the operator() is called.
- * Is mainly used to create a string containing all the vertex
- * properties
- *
- * Specialization when we do not have vertex attributes
- *
- * \tparam G graph type
- *
- */
-
-template<typename G>
-struct vtk_vertex_node<G,false>
-{
-	// Vertex object container
-	typename G::V_container & vo;
-
-	// vertex node string
-	std::string & v_node;
-
-	/*! \brief Constructor
-	 *
-	 * Create a vertex properties list
-	 *
-	 * \param v_node std::string that is filled with the graph properties in the GraphML format
-	 * \param n_obj object container to access its properties for example encapc<...>
-	 *
-	 */
-	vtk_vertex_node(std::string & v_node, typename G::V_container & n_obj)
-	:vo(n_obj),v_node(v_node)
-	{
-	};
-
-	//! It call the functor for each member
-    template<typename T>
-    void operator()(T& t)
-    {
-    	v_node += "0 0 0\n";
-    }
-};
-
-
-/*! \brief this class is a functor for "for_each" algorithm
- *
- * This class is a functor for "for_each" algorithm. For each
- * element of the boost::vector the operator() is called.
- * Is mainly used to create a string containing all the edge
- * properties
- *
- */
-
-template<typename G>
-struct vtk_edge_node
-{
-	// Vertex object container
-	typename G::E_container & vo;
-
-	// edge node string
-	std::string & e_node;
-
-	/*! \brief Constructor
-	 *
-	 * Create an edge node
-	 *
-	 * \param e_node std::string that is filled with the graph properties in the GraphML format
-	 * \param n_obj object container to access the object properties for example encapc<...>
-	 * \param n_prop number of properties
-	 *
-	 */
-	vtk_edge_node(std::string & e_node, typename G::E_container & n_obj)
-	:vo(n_obj),e_node(e_node)
-	{
-	};
-
-	/*! \brief Create a new node
-	 *
-	 * \param vc node number
-	 *
-	 */
-	void new_node(size_t v_c, size_t s, size_t d)
-	{
-		// start a new node
-		e_node += "2 " + std::to_string(s) + " " + std::to_string(d) + "\n";
-	}
-};
-
 #define GRAPH 1
 #define VECTOR_BOX 2
+#define VECTOR_GRIDS 3
 
 template <typename Graph, unsigned int imp>
 class VTKWriter
@@ -248,5 +105,6 @@ class VTKWriter
 
 #include "VTKWriter_graph.hpp"
 #include "VTKWriter_vector_box.hpp"
+#include "VTKWriter_grids.hpp"
 
 #endif /* VTKWRITER_HPP_ */
diff --git a/src/VTKWriter_graph.hpp b/src/VTKWriter_graph.hpp
index caa4f5d..fa56981 100644
--- a/src/VTKWriter_graph.hpp
+++ b/src/VTKWriter_graph.hpp
@@ -2,12 +2,156 @@
  * VTKWriter_graph.hpp
  *
  *  Created on: May 5, 2015
- *      Author: i-bird
+ *      Author: Pietro Incardona
  */
 
 #ifndef VTKWRITER_GRAPH_HPP_
 #define VTKWRITER_GRAPH_HPP_
 
+/*! \brief this class is a functor for "for_each" algorithm
+ *
+ * This class is a functor for "for_each" algorithm. For each
+ * element of the boost::vector the operator() is called.
+ * Is mainly used to create a string containing all the vertex
+ * properties
+ *
+ * \tparam G graph type
+ * \tparam attr has the vertex attributes
+ *
+ */
+
+template<typename G, bool attr>
+struct vtk_vertex_node
+{
+	// Vertex spatial type information
+	typedef typename G::V_type::s_type s_type;
+
+	s_type (& x)[3];
+
+	// Vertex object container
+	typename G::V_container & vo;
+
+	// vertex node string
+	std::string & v_node;
+
+	/*! \brief Constructor
+	 *
+	 * Create a vertex properties list
+	 *
+	 * \param v_node std::string that is filled with the graph properties in the GraphML format
+	 * \param n_obj object container to access its properties for example encapc<...>
+	 *
+	 */
+	vtk_vertex_node(std::string & v_node, typename G::V_container & n_obj, s_type (&x)[3])
+	:x(x),vo(n_obj),v_node(v_node)
+	{
+	};
+
+	//! \brief Write collected information
+	void write()
+	{
+		v_node += std::to_string(x[0]) + " " + std::to_string(x[1]) + " " + std::to_string(x[2]) + "\n";
+	}
+
+	//! It call the functor for each member
+    template<typename T>
+    void operator()(T& t)
+    {
+    	// if the attribute name is x y or z, create a string with the value of the properties and store it
+		if (G::V_type::attributes::name[T::value] == "x"){x[0] = convert<typename boost::remove_reference<decltype(vo.template get<T::value>())>::type>::template to<s_type>(vo.template get<T::value>());}
+		else if (G::V_type::attributes::name[T::value] == "y"){x[1] = convert<typename boost::remove_reference<decltype(vo.template get<T::value>())>::type>::template to<s_type>(vo.template get<T::value>());}
+		else if (G::V_type::attributes::name[T::value] == "z"){x[2] = convert<typename boost::remove_reference<decltype(vo.template get<T::value>())>::type>::template to<s_type>(vo.template get<T::value>());}
+    }
+};
+
+/*! \brief this class is a functor for "for_each" algorithm
+ *
+ * This class is a functor for "for_each" algorithm. For each
+ * element of the boost::vector the operator() is called.
+ * Is mainly used to create a string containing all the vertex
+ * properties
+ *
+ * Specialization when we do not have vertex attributes
+ *
+ * \tparam G graph type
+ *
+ */
+
+template<typename G>
+struct vtk_vertex_node<G,false>
+{
+	// Vertex object container
+	typename G::V_container & vo;
+
+	// vertex node string
+	std::string & v_node;
+
+	/*! \brief Constructor
+	 *
+	 * Create a vertex properties list
+	 *
+	 * \param v_node std::string that is filled with the graph properties in the GraphML format
+	 * \param n_obj object container to access its properties for example encapc<...>
+	 *
+	 */
+	vtk_vertex_node(std::string & v_node, typename G::V_container & n_obj)
+	:vo(n_obj),v_node(v_node)
+	{
+	};
+
+	//! It call the functor for each member
+    template<typename T>
+    void operator()(T& t)
+    {
+    	v_node += "0 0 0\n";
+    }
+};
+
+
+/*! \brief this class is a functor for "for_each" algorithm
+ *
+ * This class is a functor for "for_each" algorithm. For each
+ * element of the boost::vector the operator() is called.
+ * Is mainly used to create a string containing all the edge
+ * properties
+ *
+ */
+
+template<typename G>
+struct vtk_edge_node
+{
+	// Vertex object container
+	typename G::E_container & vo;
+
+	// edge node string
+	std::string & e_node;
+
+	/*! \brief Constructor
+	 *
+	 * Create an edge node
+	 *
+	 * \param e_node std::string that is filled with the graph properties in the GraphML format
+	 * \param n_obj object container to access the object properties for example encapc<...>
+	 * \param n_prop number of properties
+	 *
+	 */
+	vtk_edge_node(std::string & e_node, typename G::E_container & n_obj)
+	:vo(n_obj),e_node(e_node)
+	{
+	};
+
+	/*! \brief Create a new node
+	 *
+	 * \param vc node number
+	 *
+	 */
+	void new_node(size_t v_c, size_t s, size_t d)
+	{
+		// start a new node
+		e_node += "2 " + std::to_string(s) + " " + std::to_string(d) + "\n";
+	}
+};
+
 /*! \brief This class specialize functions in the case the type T
  * has or not defined attributes
  *
@@ -18,7 +162,7 @@
  *         define or not attributes name
  *
  * \tparam Graph type of graph we are processing
- * \tparam p the property we are going to write
+ * \tparam i the property we are going to write
  *
  */
 
diff --git a/src/VTKWriter_grids.hpp b/src/VTKWriter_grids.hpp
new file mode 100644
index 0000000..2c6f3b5
--- /dev/null
+++ b/src/VTKWriter_grids.hpp
@@ -0,0 +1,406 @@
+/*
+ * VTKWriter_grids.hpp
+ *
+ *  Created on: May 5, 2015
+ *      Author: Pietro Incardona
+ */
+
+#ifndef VTKWRITER_GRIDS_HPP_
+#define VTKWRITER_GRIDS_HPP_
+
+#include <boost/mpl/pair.hpp>
+#include "VTKWriter_grids_util.hpp"
+
+template <typename Grid, typename St>
+class ele_g
+{
+public:
+
+	typedef Grid value_type;
+
+	ele_g(const Grid & g, Point<Grid::dims,St> & offset, Point<Grid::dims,St> & spacing, Box<Grid::dims,St> dom)
+	:g(g),offset(offset),spacing(spacing),dom(dom)
+	{}
+
+	std::string dataset;
+	//! Grid
+	const Grid & g;
+	//! offset where it start
+	Point<Grid::dims,St> offset;
+	// spacing of the grid
+	Point<Grid::dims,St> spacing;
+	// Part of the grid that is real domain
+	Box<Grid::dims,St> dom;
+};
+
+
+
+/*! \brief This class specialize functions in the case the type T
+ * has not defined attributes
+ *
+ * In C++ partial specialization of a function is not allowed so we have to
+ * encapsulate this function in a class
+ *
+ * \tparam has_attributes parameter that specialize the function in case the vertex
+ *         define or not attributes name
+ *
+ * \tparam i id of the property we are going to write
+ *
+ */
+
+template<typename ele_g, typename St, unsigned int i>
+class prop_output_g<false,St,ele_g,i>
+{
+public:
+
+	/*! \brief Given a Graph return the point data header for a typename T
+	 *
+	 * \tparam T type to write
+	 * \param n_node number of the node
+	 *
+	 */
+
+	static std::string get_point_property_header()
+	{
+		//! vertex node output string
+		std::string v_out;
+
+		// Check if T is a supported format
+		// for now we support only scalar of native type
+
+		std::string type = getType<typename boost::fusion::result_of::at<typename ele_g::value_type::value_type::type,boost::mpl::int_<i>>::type>();
+
+		// if the type is not supported return
+		if (type.size() == 0)
+		{return v_out;}
+
+		// Create point data properties
+		v_out += "SCALARS " + get_attributes() + " " + type + "\n";
+
+		// Default lookup table
+		v_out += "LOOKUP_TABLE default\n";
+
+		// return the vertex list
+		return v_out;
+	}
+
+	/*! \brief Get the attributes name
+	 *
+	 */
+	static std::string get_attributes()
+	{
+		return std::string("attr" + std::to_string(i));
+	}
+};
+
+/*! \brief this class is a functor for "for_each" algorithm
+ *
+ * This class is a functor for "for_each" algorithm. For each
+ * element of the boost::vector the operator() is called.
+ * Is mainly used to produce at output for each property
+ *
+ * \tparam Graph graph we are processing
+ *
+ * \param dim Dimensionality
+ * \param S type of grid
+ *
+ */
+
+template<typename ele_g, typename St>
+struct prop_out_g
+{
+	// property output string
+	std::string & v_out;
+
+	// grid that we are processing
+	const openfpm::vector_std< ele_g > & vg;
+
+	/*! \brief constructor
+	 *
+	 * \param v_out string to fill with the vertex properties
+	 *
+	 */
+	prop_out_g(std::string & v_out, const openfpm::vector_std< ele_g > & vg)
+	:v_out(v_out),vg(vg)
+	{};
+
+	//! It produce an output for each property
+    template<typename T>
+    void operator()(T& t) const
+    {
+    	typedef typename boost::mpl::at<typename ele_g::value_type::value_type::type,boost::mpl::int_<T::value>>::type ptype;
+
+    	meta_prop<boost::mpl::int_<T::value> ,ele_g,St, ptype > m(vg,v_out);
+    }
+};
+
+/*!
+ *
+ * It write a VTK format file in case of grids defined on a space
+ *
+ * \tparam boost::mpl::pair<G,S>
+ *
+ * where G is the type of grid S is the type of space, float, double ...
+ *
+ */
+template <typename pair>
+class VTKWriter<pair,VECTOR_GRIDS>
+{
+	//! Vector of grids
+	openfpm::vector< ele_g<typename pair::first,typename pair::second> > vg;
+
+	/*! \brief Get the total number of points
+	 *
+	 * \return the total number
+	 *
+	 */
+	size_t get_total()
+	{
+		size_t tot = 0;
+
+		//! Calculate the full number of vertices
+		for (size_t i = 0 ; i < vg.size() ; i++)
+		{
+			tot += vg.get(i).g.size();
+		}
+		return tot;
+	}
+
+	/*! \brief It get the vertex properties list
+	 *
+	 * It get the vertex properties list of the vertex defined as VTK header
+	 *
+	 * \return a string that define the vertex properties in graphML format
+	 *
+	 */
+
+	std::string get_vertex_properties_list()
+	{
+		//! vertex property output string
+		std::string v_out;
+
+		// write the number of vertex
+		v_out += "VERTICES " + std::to_string(get_total()) + " " + std::to_string(vg.size() * 2) + "\n";
+
+		// return the vertex properties string
+		return v_out;
+	}
+
+	/*! \brief It get the vertex properties list
+	 *
+	 * It get the vertex properties list of the vertex defined as a VTK header
+	 *
+	 * \return a string that define the vertex properties in graphML format
+	 *
+	 */
+
+	std::string get_point_properties_list()
+	{
+		//! vertex property output string
+		std::string v_out;
+
+		// write the number of vertex
+		v_out += "POINTS " + std::to_string(get_total()) + " float" + "\n";
+
+		// return the vertex properties string
+		return v_out;
+	}
+
+	/*! \brief Create the VTK point definition
+	 *
+	 * \tparam s_type spatial type of the data
+	 * \tparam attr false x,y,z are set to 0 for each vertex
+	 *
+	 */
+
+	template <bool attr> std::string get_point_list()
+	{
+		//! vertex node output string
+		std::stringstream v_out;
+
+		//! For each defined grid
+
+		for (size_t i = 0 ; i < vg.size() ; i++)
+		{
+			//! Get the iterator
+			auto it = vg.get(i).g.getIterator();
+
+			//! Where the grid is defined
+			Box<pair::first::dims,typename pair::second> dom;
+
+			// if there is the next element
+			while (it.isNext())
+			{
+				Point<pair::first::dims,typename pair::second> p;
+				p = it.get().toPoint();
+				p = p * vg.get(i).spacing + vg.get(i).offset;
+
+				v_out << p.toString() << "\n";
+
+				// increment the iterator and counter
+				++it;
+			}
+		}
+
+		// return the vertex list
+		return v_out.str();
+	}
+
+	/*! \brief Create the VTK vertex definition
+	 *
+	 * \tparam s_type spatial type of the data
+	 * \tparam attr false x,y,z are set to 0 for each vertex
+	 *
+	 */
+
+	std::string get_vertex_list()
+	{
+		//! vertex node output string
+		std::string v_out;
+
+		size_t k = 0;
+
+		for (size_t i = 0 ; i < vg.size() ; i++)
+		{
+			//! For each grid point create a vertex
+			auto it = vg.get(i).g.getIterator();
+
+			while (it.isNext())
+			{
+				v_out += "1 " + std::to_string(k) + "\n";
+
+				++k;
+				++it;
+			}
+		}
+		// return the vertex list
+		return v_out;
+	}
+
+	/*! \brief Get the point data header
+	 *
+	 * \return a string with the point data header for VTK format
+	 *
+	 */
+
+	std::string get_point_data_header()
+	{
+		std::string v_out;
+
+		v_out += "POINT_DATA " + std::to_string(get_total()) + "\n";
+
+		return v_out;
+	}
+
+public:
+
+	/*!
+	 *
+	 * VTKWriter constructor
+	 *
+	 */
+	VTKWriter()
+	{}
+
+	/*! \brief Add grid dataset
+	 *
+	 * \param g Grid to add
+	 * \param offset grid offset
+	 * \param spacing spacing of the grid
+	 * \param dom part of the spacethat is the domain
+	 *
+	 */
+	void add(const typename pair::first & g, Point<pair::first::dims,typename pair::second> & offset, Point<pair::first::dims,typename pair::second> & spacing, Box<pair::first::dims,typename pair::second> dom)
+	{
+		ele_g<typename pair::first,typename pair::second> t(g,offset,spacing,dom);
+
+		vg.add(t);
+	}
+
+	/*! \brief It write a VTK file from a graph
+	 *
+	 * \tparam prp_out which properties to output [default = -1 (all)]
+	 *
+	 * \param file path where to write
+	 * \param name of the graph
+	 * \param file_type specify if it is a VTK BINARY or ASCII file [default = ASCII]
+	 *
+	 */
+
+	template<int prp = -1> bool write(std::string file, std::string f_name = "grids" , file_type ft = file_type::ASCII)
+	{
+		// Header for the vtk
+		std::string vtk_header;
+		// Point list of the VTK
+		std::string point_list;
+		// Vertex list of the VTK
+		std::string vertex_list;
+		// Graph header
+		std::string vtk_binary_or_ascii;
+		// vertex properties header
+		std::string point_prop_header;
+		// edge properties header
+		std::string vertex_prop_header;
+		// Data point header
+		std::string point_data_header;
+		// Data point
+		std::string point_data;
+
+		// VTK header
+		vtk_header = "# vtk DataFile Version 3.0\n"
+				     + f_name + "\n";
+
+		// Choose if binary or ASCII
+		if (ft == file_type::ASCII)
+		{vtk_header += "ASCII\n";}
+		else
+		{vtk_header += "BINARY\n";}
+
+		// Data type for graph is DATASET POLYDATA
+		vtk_header += "DATASET POLYDATA\n";
+
+		// point properties header
+		point_prop_header = get_point_properties_list();
+
+		// Get point list
+		point_list = get_point_list<has_attributes<typename pair::first::value_type>::value>();
+
+		// vertex properties header
+		vertex_prop_header = get_vertex_properties_list();
+
+		// Get vertex list
+		vertex_list = get_vertex_list();
+
+		// Get the point data header
+		point_data_header = get_point_data_header();
+
+		// For each property in the vertex type produce a point data
+
+		prop_out_g< ele_g<typename pair::first,typename pair::second>, typename pair::second > pp(point_data, vg);
+
+		if (prp == -1)
+			boost::mpl::for_each< boost::mpl::range_c<int,0, pair::first::value_type::max_prop> >(pp);
+		else
+			boost::mpl::for_each< boost::mpl::range_c<int,prp, prp> >(pp);
+
+		// write the file
+		std::ofstream ofs(file);
+
+		// Check if the file is open
+		if (ofs.is_open() == false)
+		{std::cerr << "Error cannot create the VTK file: " + file;}
+
+		ofs << vtk_header << point_prop_header << point_list <<
+				vertex_prop_header << vertex_list << point_data_header << point_data;
+
+		// Close the file
+
+		ofs.close();
+
+		// Completed succefully
+		return true;
+	}
+};
+
+
+#endif /* VTKWRITER_GRAPH_HPP_ */
diff --git a/src/VTKWriter_unit_tests.hpp b/src/VTKWriter_unit_tests.hpp
index eb7ed81..f0accab 100644
--- a/src/VTKWriter_unit_tests.hpp
+++ b/src/VTKWriter_unit_tests.hpp
@@ -221,6 +221,87 @@ BOOST_AUTO_TEST_CASE( vtk_writer_use_vector_box)
 	BOOST_REQUIRE_EQUAL(test,true);
 }
 
+/*! \fill the CPU with some random data
+ *
+ *
+ */
+void fill_grid_some_data(grid_cpu<2,Point_test<float>> & g)
+{
+	typedef Point_test<float> p;
+
+	auto it = g.getIterator();
+
+	while (it.isNext())
+	{
+		g.template get<p::x>(it.get()) = it.get().get(0);
+		g.template get<p::y>(it.get()) = it.get().get(0);
+		g.template get<p::z>(it.get()) = 0;
+		g.template get<p::s>(it.get()) = 1.0;
+		g.template get<p::v>(it.get())[0] = g.getGrid().LinId(it.get());
+		g.template get<p::v>(it.get())[1] = g.getGrid().LinId(it.get());
+		g.template get<p::v>(it.get())[2] = g.getGrid().LinId(it.get());
+
+		g.template get<p::t>(it.get())[0][0] = g.getGrid().LinId(it.get());
+		g.template get<p::t>(it.get())[0][1] = g.getGrid().LinId(it.get());
+		g.template get<p::t>(it.get())[0][2] = g.getGrid().LinId(it.get());
+		g.template get<p::t>(it.get())[1][0] = g.getGrid().LinId(it.get());
+		g.template get<p::t>(it.get())[1][1] = g.getGrid().LinId(it.get());
+		g.template get<p::t>(it.get())[1][2] = g.getGrid().LinId(it.get());
+		g.template get<p::t>(it.get())[2][0] = g.getGrid().LinId(it.get());
+		g.template get<p::t>(it.get())[2][1] = g.getGrid().LinId(it.get());
+		g.template get<p::t>(it.get())[2][2] = g.getGrid().LinId(it.get());
+
+		++it;
+	}
+}
+
+
+BOOST_AUTO_TEST_CASE( vtk_writer_use_grids)
+{
+	// Create box grids
+	Point<2,float> offset1({0.0,0.0});
+	Point<2,float> spacing1({0.1,0.2});
+	Box<2,size_t> d1({1,2},{14,15});
+
+	// Create box grids
+	Point<2,float> offset2({5.0,7.0});
+	Point<2,float> spacing2({0.2,0.1});
+	Box<2,size_t> d2({2,1},{13,15});
+
+	// Create box grids
+	Point<2,float> offset3({0.0,7.0});
+	Point<2,float> spacing3({0.05,0.07});
+	Box<2,size_t> d3({3,2},{11,10});
+
+	// Create box grids
+	Point<2,float> offset4({5.0,0.0});
+	Point<2,float> spacing4({0.1,0.1});
+	Box<2,size_t> d4({1,1},{7,7});
+
+	size_t sz[] = {16,16};
+	grid_cpu<2,Point_test<float>> g1(sz);
+	fill_grid_some_data(g1);
+	grid_cpu<2,Point_test<float>> g2(sz);
+	fill_grid_some_data(g2);
+	grid_cpu<2,Point_test<float>> g3(sz);
+	fill_grid_some_data(g3);
+	grid_cpu<2,Point_test<float>> g4(sz);
+	fill_grid_some_data(g4);
+
+	// Create a writer and write
+	VTKWriter<boost::mpl::pair<grid_cpu<2,Point_test<float>>,float>,VECTOR_GRIDS> vtk_g;
+	vtk_g.add(g1,offset1,spacing1,d1);
+	vtk_g.add(g2,offset2,spacing2,d2);
+	vtk_g.add(g3,offset3,spacing3,d3);
+	vtk_g.add(g4,offset4,spacing4,d4);
+
+	vtk_g.write("vtk_grids.vtk");
+
+	// Check that match
+	bool test = compare("vtk_grids.vtk","vtk_grids_test.vtk");
+	BOOST_REQUIRE_EQUAL(test,true);
+}
+
 BOOST_AUTO_TEST_SUITE_END()
 
 #endif /* VTKWRITER_UNIT_TESTS_HPP_ */
-- 
GitLab