diff --git a/python/damask/_geom.py b/python/damask/_geom.py
index d3b375bf6..03179aac9 100644
--- a/python/damask/_geom.py
+++ b/python/damask/_geom.py
@@ -785,7 +785,7 @@ class Geom:
         """
         def mp(entry,mapper):
             return mapper[entry] if entry in mapper else entry
-       
+
         mp = np.vectorize(mp)
         mapper = dict(zip(from_material,to_material))
 
@@ -796,6 +796,20 @@ class Geom:
                    )
 
 
+    def sort(self):
+        """Sort material indices such that min(material) is located at (0,0,0)."""
+        a = self.material.flatten(order='F')
+        from_ma = pd.unique(a)
+        sort_idx = np.argsort(from_ma)
+        ma = np.unique(a)[sort_idx][np.searchsorted(from_ma,a,sorter = sort_idx)]
+
+        return Geom(material = ma.reshape(self.grid,order='F'),
+                    size     = self.size,
+                    origin   = self.origin,
+                    comments = self.comments+[util.execution_stamp('Geom','sort')],
+                   )
+
+
     def vicinity_offset(self,vicinity=1,offset=None,trigger=[],periodic=True):
         """
         Offset material index of points in the vicinity of xxx.
@@ -822,13 +836,9 @@ class Geom:
         def tainted_neighborhood(stencil,trigger):
 
             me = stencil[stencil.shape[0]//2]
-            if len(trigger) == 0:
-                return np.any(stencil != me)
-            if me in trigger:
-                trigger = set(trigger)
-                trigger.remove(me)
-                trigger = list(trigger)
-            return np.any(np.in1d(stencil,np.array(trigger)))
+            return np.any(stencil != me
+                          if len(trigger) == 0 else
+                          np.in1d(stencil,np.array(list(set(trigger) - {me}))))
 
         offset_ = np.nanmax(self.material) if offset is None else offset
         mask = ndimage.filters.generic_filter(self.material,
diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py
index 0d393a4c4..5c2777da7 100644
--- a/python/tests/test_Geom.py
+++ b/python/tests/test_Geom.py
@@ -6,6 +6,7 @@ from damask import Geom
 from damask import Table
 from damask import Rotation
 from damask import util
+from damask import seeds
 from damask import grid_filters
 
 
@@ -204,6 +205,11 @@ class TestGeom:
         assert np.array_equiv(t,f) or (not geom_equal(modified,default))
         assert geom_equal(default, modified.substitute(t,f))
 
+    def test_sort(self):
+        grid = np.random.randint(5,20,3)
+        m = Geom(np.random.randint(1,20,grid)*3,np.ones(3)).sort().material.flatten(order='F')
+        for i,v in enumerate(m):
+            assert i==0 or v > m[:i].max() or v in m[:i]
 
     @pytest.mark.parametrize('axis_angle',[np.array([1,0,0,86.7]), np.array([0,1,0,90.4]), np.array([0,0,1,90]),
                                            np.array([1,0,0,175]),np.array([0,-1,0,178]),np.array([0,0,1,180])])
@@ -384,3 +390,13 @@ class TestGeom:
         t = Table(np.column_stack((coords,z)),{'coords':3,'z':1})
         g = Geom.from_table(t,'coords',['1_coords','z'])
         assert g.N_materials == g.grid[0]*2 and (g.material[:,:,-1]-g.material[:,:,0] == grid[0]).all()
+
+
+    def test_from_table_recover(self,tmp_path):
+        grid = np.random.randint(60,100,3)
+        size = np.ones(3)+np.random.rand(3)
+        s = seeds.from_random(size,np.random.randint(60,100))
+        geom = Geom.from_Voronoi_tessellation(grid,size,s)
+        coords = grid_filters.cell_coord0(grid,size)
+        t = Table(np.column_stack((coords.reshape(-1,3,order='F'),geom.material.flatten(order='F'))),{'c':3,'m':1})
+        assert geom_equal(geom.sort().renumber(),Geom.from_table(t,'c',['m']))