diff --git a/tests/test_integral_operators.py b/tests/test_integral_operators.py
index 2a0979e..d4c7ebf 100644
--- a/tests/test_integral_operators.py
+++ b/tests/test_integral_operators.py
@@ -1,43 +1,70 @@
 #!/usr/bin/env python
 # coding: utf-8
 # -----------------------------------------------------------------------------
 # @author Lucas Frérot <lucas.frerot@epfl.ch>
 #
 # @section LICENSE
 #
 # Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
 # Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
 # Solides)
 #
 # Tamaas is free  software: you can redistribute it and/or  modify it under the
 # terms  of the  GNU Lesser  General Public  License as  published by  the Free
 # Software Foundation, either version 3 of the License, or (at your option) any
 # later version.
 #
 # Tamaas 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  Lesser General  Public License  for more
 # details.
 #
 # You should  have received  a copy  of the GNU  Lesser General  Public License
 # along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
 # -----------------------------------------------------------------------------
 import tamaas as tm
 import numpy as np
 import matplotlib.pyplot as plt
+from pyevtk.hl import gridToVTK
+
+def dumpGrid(name, coords, field, func=lambda x: x):
+    field_vtk = [np.zeros([coords[i].shape[0] for i in range(3)]) for j in range(3)]
+    for i in range(3):
+        field_vtk[i][:, :, :] = func(field[:, :, :, i])
+
+    gridToVTK(name,
+              *coords,
+              pointData={
+                  'X':field_vtk[0],
+                  'Y':field_vtk[1],
+                  'Z':field_vtk[2],
+              })
+
+N = 16
+
+domain = [1.] * 3
+discretization = [N] * 3
 
 model = tm.ModelFactory.createModel(tm.model_type.volume_2d,
-                                    [1.]*3,
-                                    [8]*3)
+                                    domain,
+                                    discretization)
+
 displacement = model.getDisplacement()
 
 engine = tm._tamaas._test_features.Kelvin(model)
 
 source = np.zeros_like(displacement)
-source[4, 0, 0, 2] = 1
+source[N//2, 0, 0, 2] = 1
 
+print("Starting computation")
 engine.applyVolumeForcePotential(source, displacement)
+print("End computation")
+
+coords = [np.linspace(0, domain[i], discretization[i]) for i in range(3)]
+
+dumpGrid('kelvin', coords, displacement, lambda x: np.einsum('ijk->jki', x))
+dumpGrid('source', coords, source, lambda x: np.einsum('ijk->jki', x))
 
-plt.imshow(displacement[4, :, :, 2])
+plt.imshow(displacement[:, 4, :, 2])
 plt.colorbar()
-plt.show()
+#plt.show()