diff --git a/processing/post/DADF5_vtk_points.py b/processing/post/DADF5_vtk_points.py
index ee71f787d..87c1ad93e 100755
--- a/processing/post/DADF5_vtk_points.py
+++ b/processing/post/DADF5_vtk_points.py
@@ -40,19 +40,14 @@ if options.con is None: options.con=[]
 for filename in options.filenames:
   results = damask.DADF5(filename)
   
-  Polydata = vtk.vtkPolyData()
   Points   = vtk.vtkPoints()
-  Vertices = vtk.vtkCellArray()
-  
-  if results.structured:                                                                            # for grid solvers calculate points
-    delta = results.size/results.grid*0.5
-    for z in np.linspace(delta[2],results.size[2]-delta[2],results.grid[2]):
-      for y in np.linspace(delta[1],results.size[1]-delta[1],results.grid[1]):
-        for x in np.linspace(delta[0],results.size[0]-delta[0],results.grid[0]):
-          pointID = Points.InsertNextPoint([x,y,z])
-          Vertices.InsertNextCell(1)
-          Vertices.InsertCellPoint(pointID)
+  Vertices = vtk.vtkCellArray()  
+  for c in results.cell_coordinates():
+    pointID = Points.InsertNextPoint(c)
+    Vertices.InsertNextCell(1)
+    Vertices.InsertCellPoint(pointID)
 
+  Polydata = vtk.vtkPolyData()
   Polydata.SetPoints(Points)
   Polydata.SetVerts(Vertices)
   Polydata.Modified()
diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py
index 8444fa259..955765fcc 100644
--- a/python/damask/dadf5.py
+++ b/python/damask/dadf5.py
@@ -335,6 +335,20 @@ class DADF5():
     return dataset
 
 
+  def cell_coordinates(self):
+    """Initial coordinates of the cell centers."""
+    if self.structured:
+      delta = self.size/self.grid*0.5
+      z, y, x = np.meshgrid(np.linspace(delta[2],self.size[2]-delta[2],self.grid[2]),
+                            np.linspace(delta[1],self.size[1]-delta[1],self.grid[1]),
+                            np.linspace(delta[0],self.size[0]-delta[0],self.grid[0]),
+                           )
+      return np.concatenate((x[:,:,:,None],y[:,:,:,None],y[:,:,:,None]),axis = 3).reshape([np.product(self.grid),3])
+    else:
+      with h5py.File(self.filename,'r') as f:
+        return f['geometry/x_c'][()]
+
+
   def add_Cauchy(self,P='P',F='F'):
     """
     Adds Cauchy stress calculated from 1st Piola-Kirchhoff stress and deformation gradient.