diff --git a/python/damask/_geom.py b/python/damask/_geom.py
index b40bedb04..5d89d0750 100644
--- a/python/damask/_geom.py
+++ b/python/damask/_geom.py
@@ -849,37 +849,30 @@ class Geom:
         Parameters
         ----------
         periodic : Boolean, optional
-            Show boundaries at periodic nodes too. Defalts to False.
+            Show boundaries at periodic nodes too. Defaults to False.
 
         """
         v = VTK.from_polyData(grid_filters.node_coord0(self.grid,self.size,self.origin).reshape(-1,3,order='F'))
         cells = vtk.vtkCellArray()
 
-        #Grid to nodes
-        a=self.grid[0]+1
-        b=self.grid[1]+1
-        c=self.grid[2]+1
-        
         if periodic:
-            FaceYZ_per=np.empty()
-            FaceXY_per=np.empty()
-            FaceXZ_per=np.empty()
+            FaceYZ_per=np.empty([])
+            FaceXY_per=np.empty([])
+            FaceXZ_per=np.empty([])
         else:          
-            # get nodes of periodic boundaries in X
-            FaceYZ_per=np.concatenate((  np.arange(0,a*b*c,a) , np.arange(self.grid[0],   a*b*c,a)  ))  
-            # get nodes of periodic boundaries in Z
-            FaceXY_per=np.concatenate((  np.arange(0,a*b)     , np.arange(a*b*c -  a*b, a*b*c )  ))
-            # get nodes of periodic boundaries in Y
-            FaceXZ_0=np.arange(0, a) # face Y=0
-            for i in range(c):
-                add=np.arange( (i+1)*(a)*(b), (i+1)*(a)*(b)+ (a))
-                FaceXZ_0=np.append(FaceXZ_0, add)
-            FaceXZ_end=np.arange(a*b-b, a*b) # face Y=end
-            for i in range(c):
-                add= np.arange( (i+1)*a*b-b, (i+1)*a*b)
-                FaceXZ_end=np.append(FaceXZ_end, add)
-            FaceXZ_per=np.concatenate (( FaceXZ_0,  FaceXZ_end ))
-        
+            FaceYZ_per=np.concatenate((
+                        np.arange(0                                        , np.prod(self.grid+1)     , (self.grid[0]+1) ),
+                        np.arange(self.grid[0]                             , np.prod(self.grid+1)     , (self.grid[0]+1) )  ))
+            FaceXY_per=np.concatenate((  
+                        np.arange(0                                        , np.prod(self.grid[:2]+1) , 1                ), 
+                        np.arange(np.prod(self.grid+1) - np.prod(self.grid[:2]+1), \
+                                                    np.prod(self.grid+1)     , 1                )  ))
+            FaceXZ_per=np.concatenate((   
+                        np.concatenate([np.arange((i)*np.prod(self.grid[:2]+1), \
+                                       (i)*np.prod(self.grid[:2]+1)+ (self.grid[0]+1)) for i in range(self.grid[2]+1)]) ,
+                        np.concatenate([np.arange( (i+1)*np.prod(self.grid[:2]+1)-(self.grid[0]+1),\
+                                                (i+1)*np.prod(self.grid[:2]+1))  for i in range(self.grid[2]+1)])  ))
+          
         for d_s in [0,1,2]:
             base_nodes = np.where(self.material==np.roll(self.material,1,d_s),False,True)
             for d in [0,1,2]:
@@ -898,17 +891,17 @@ class Geom:
                     q = vtk.vtkQuad()
                     q.GetPointIds().SetId(0, p[0])
                     if d_s == 0:
-                        q.GetPointIds().SetId(1, p[0]+a)
-                        q.GetPointIds().SetId(2, p[0]+a*b+a)
-                        q.GetPointIds().SetId(3, p[0]+a*b)
+                        q.GetPointIds().SetId(1, p[0]+(self.grid[0]+1))
+                        q.GetPointIds().SetId(2, p[0]+np.prod(self.grid[:2]+1)+(self.grid[0]+1))
+                        q.GetPointIds().SetId(3, p[0]+np.prod(self.grid[:2]+1))
                     if d_s == 1:
-                        q.GetPointIds().SetId(1, p[0]+a*b)
-                        q.GetPointIds().SetId(2, p[0]+a*b+1)
+                        q.GetPointIds().SetId(1, p[0]+np.prod(self.grid[:2]+1))
+                        q.GetPointIds().SetId(2, p[0]+np.prod(self.grid[:2]+1)+1)
                         q.GetPointIds().SetId(3, p[0]+1)
                     if d_s == 2:
                         q.GetPointIds().SetId(1, p[0]+1)
-                        q.GetPointIds().SetId(2, p[0]+a+1)
-                        q.GetPointIds().SetId(3, p[0]+a)
+                        q.GetPointIds().SetId(2, p[0]+(self.grid[0]+1)+1)
+                        q.GetPointIds().SetId(3, p[0]+(self.grid[0]+1))
                     cells.InsertNextCell(q)
         v.vtk_data.SetPolys(cells)
         v.save('GrainBoundaries')