def grid2vtk(list_gx, list_gy, filename):
	f = open(filename,'wt')
	
	# create a dictionary to store the vertices position and its index
	# it does not count twice vertices with the same x,y position
	vertices = {}
	vt_idx = 0
	id_to_vertices = {}
	for k, X in enumerate(list_gx):
		Y = list_gy[k]
		n,m = X.shape
		# iterate over all cells for both matrices and create a new vertice
		# for each cell in the matrices
		for i in range(n):
			for j in range(m):
				if (X[i,j], Y[i,j]) in vertices.keys():
					continue
				vertices.update({(X[i,j], Y[i,j]): vt_idx})
				id_to_vertices.update({vt_idx: (X[i,j], Y[i,j])})
				vt_idx += 1

	f.write('# vtk DataFile Version 2.0\n')
	f.write('file generated by grid2vtk\n')
	f.write('ASCII\n')
	f.write('DATASET UNSTRUCTURED_GRID\n')
	dim = [i.shape for i in list_gx]
	f.write('POINTS '+str(len(vertices))+' float\n')
	
	for v in range(len(vertices)):
		f.write("{:.8f}".format(id_to_vertices[v][0])+' '+"{:.8f}".format(id_to_vertices[v][1])+' 0.0\n')
	f.write('\n')
	
	nc = sum([(i[0]-1)*(i[1]-1) for i in dim])
	f.write('CELLS '+str(nc)+' '+str(nc*5)+'\n')
	
	for k, X in enumerate(list_gx):
		Y = list_gy[k]
		n, m = X.shape
		for i in range(n-1):
			for j in range(m-1):
				f.write('4 %d %d %d %d\n'%(
						vertices[( X[i,j], Y[i,j] )],
						vertices[( X[i,j+1], Y[i,j+1] )],
						vertices[( X[i+1,j+1], Y[i+1,j+1] )],
						vertices[( X[i+1,j], Y[i+1,j] )],
						))

	f.write('\n')
	f.write('CELL_TYPES '+str(nc)+'\n')
	f.write(nc*'7 ')
	f.write('\n')

	f.close()
