diff --git a/PRIVATE b/PRIVATE
index 72d526e57..8dde2a685 160000
--- a/PRIVATE
+++ b/PRIVATE
@@ -1 +1 @@
-Subproject commit 72d526e5750366a9efe4d1fd9d92e0d1ecd2cd38
+Subproject commit 8dde2a68538b7cffbe9d370e2b60be90a31627ab
diff --git a/VERSION b/VERSION
index 3ba40668a..6d4d1897c 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-v2.0.3-2504-gcee9daff
+v2.0.3-2514-g873b9fa8
diff --git a/processing/post/addSchmidfactors.py b/processing/post/addSchmidfactors.py
index b3c4c7500..d1899e9a0 100755
--- a/processing/post/addSchmidfactors.py
+++ b/processing/post/addSchmidfactors.py
@@ -2,6 +2,7 @@
 
 import os
 import sys
+from io import StringIO
 from optparse import OptionParser
 
 import numpy as np
@@ -15,91 +16,82 @@ scriptID   = ' '.join([scriptName,damask.version])
 slipSystems = {
 'fcc':
     np.array([
-    # Slip direction     Plane normal
-     [ 0, 1,-1,     1, 1, 1, ],
-     [-1, 0, 1,     1, 1, 1, ],
-     [ 1,-1, 0,     1, 1, 1, ],
-     [ 0,-1,-1,    -1,-1, 1, ],
-     [ 1, 0, 1,    -1,-1, 1, ],
-     [-1, 1, 0,    -1,-1, 1, ],
-     [ 0,-1, 1,     1,-1,-1, ],
-     [-1, 0,-1,     1,-1,-1, ],
-     [ 1, 1, 0,     1,-1,-1, ],
-     [ 0, 1, 1,    -1, 1,-1, ],
-     [ 1, 0,-1,    -1, 1,-1, ],
-     [-1,-1, 0,    -1, 1,-1, ],
-    ],'f'),
+    [+0,+1,-1 , +1,+1,+1],
+    [-1,+0,+1 , +1,+1,+1],
+    [+1,-1,+0 , +1,+1,+1],
+    [+0,-1,-1 , -1,-1,+1],
+    [+1,+0,+1 , -1,-1,+1],
+    [-1,+1,+0 , -1,-1,+1],
+    [+0,-1,+1 , +1,-1,-1],
+    [-1,+0,-1 , +1,-1,-1],
+    [+1,+1,+0 , +1,-1,-1],
+    [+0,+1,+1 , -1,+1,-1],
+    [+1,+0,-1 , -1,+1,-1],
+    [-1,-1,+0 , -1,+1,-1],
+   ],'d'),
 'bcc':
     np.array([
-    # Slip system <111>{110} 
-     [ 1,-1, 1,     0, 1, 1, ],
-     [-1,-1, 1,     0, 1, 1, ],
-     [ 1, 1, 1,     0,-1, 1, ],
-     [-1, 1, 1,     0,-1, 1, ],
-     [-1, 1, 1,     1, 0, 1, ],
-     [-1,-1, 1,     1, 0, 1, ],
-     [ 1, 1, 1,    -1, 0, 1, ],
-     [ 1,-1, 1,    -1, 0, 1, ],
-     [-1, 1, 1,     1, 1, 0, ],
-     [-1, 1,-1,     1, 1, 0, ],
-     [ 1, 1, 1,    -1, 1, 0, ],
-     [ 1, 1,-1,    -1, 1, 0, ],
-    # Slip system <111>{112}
-     [-1, 1, 1,     2, 1, 1, ],
-     [ 1, 1, 1,    -2, 1, 1, ],
-     [ 1, 1,-1,     2,-1, 1, ],
-     [ 1,-1, 1,     2, 1,-1, ],
-     [ 1,-1, 1,     1, 2, 1, ],
-     [ 1, 1,-1,    -1, 2, 1, ],
-     [ 1, 1, 1,     1,-2, 1, ],
-     [-1, 1, 1,     1, 2,-1, ],
-     [ 1, 1,-1,     1, 1, 2, ],
-     [ 1,-1, 1,    -1, 1, 2, ],
-     [-1, 1, 1,     1,-1, 2, ],
-     [ 1, 1, 1,     1, 1,-2, ],
-    ],'f'),
+    [+1,-1,+1 , +0,+1,+1],
+    [-1,-1,+1 , +0,+1,+1],
+    [+1,+1,+1 , +0,-1,+1],
+    [-1,+1,+1 , +0,-1,+1],
+    [-1,+1,+1 , +1,+0,+1],
+    [-1,-1,+1 , +1,+0,+1],
+    [+1,+1,+1 , -1,+0,+1],
+    [+1,-1,+1 , -1,+0,+1],
+    [-1,+1,+1 , +1,+1,+0],
+    [-1,+1,-1 , +1,+1,+0],
+    [+1,+1,+1 , -1,+1,+0],
+    [+1,+1,-1 , -1,+1,+0],
+    [-1,+1,+1 , +2,+1,+1],
+    [+1,+1,+1 , -2,+1,+1],
+    [+1,+1,-1 , +2,-1,+1],
+    [+1,-1,+1 , +2,+1,-1],
+    [+1,-1,+1 , +1,+2,+1],
+    [+1,+1,-1 , -1,+2,+1],
+    [+1,+1,+1 , +1,-2,+1],
+    [-1,+1,+1 , +1,+2,-1],
+    [+1,+1,-1 , +1,+1,+2],
+    [+1,-1,+1 , -1,+1,+2],
+    [-1,+1,+1 , +1,-1,+2],
+    [+1,+1,+1 , +1,+1,-2],
+   ],'d'),
 'hex':
     np.array([
-    # Basal systems <11.0>{00.1} (independent of c/a-ratio, Bravais notation (4 coordinate base))
-     [ 2, -1, -1,  0,     0,  0,  0,  1, ],
-     [-1,  2, -1,  0,     0,  0,  0,  1, ],
-     [-1, -1,  2,  0,     0,  0,  0,  1, ],
-    # 1st type prismatic systems <11.0>{10.0}  (independent of c/a-ratio)
-     [ 2, -1, -1,  0,     0,  1, -1,  0, ],
-     [-1,  2, -1,  0,    -1,  0,  1,  0, ],
-     [-1, -1,  2,  0,     1, -1,  0,  0, ],
-    # 2nd type prismatic systems <10.0>{11.0} -- a slip; plane normals independent of c/a-ratio
-     [ 0,  1,  -1, 0,     2, -1, -1,  0, ],
-     [-1,  0,  1,  0,    -1,  2, -1,  0, ],
-     [ 1, -1,  0,  0,    -1, -1,  2,  0, ],
-    # 1st type 1st order pyramidal systems <11.0>{-11.1} -- plane normals depend on the c/a-ratio
-     [ 2, -1, -1,  0,     0,  1, -1,  1, ],
-     [-1,  2, -1,  0,    -1,  0,  1,  1, ],
-     [-1, -1,  2,  0,     1, -1,  0,  1, ],
-     [ 1,  1, -2,  0,    -1,  1,  0,  1, ],
-     [-2,  1,  1,  0,     0, -1,  1,  1, ],
-     [ 1, -2,  1,  0,     1,  0, -1,  1, ],
-    # pyramidal system: c+a slip <11.3>{-10.1} -- plane normals depend on the c/a-ratio
-     [ 2, -1, -1,  3,    -1,  1,  0,  1, ],
-     [ 1, -2,  1,  3,    -1,  1,  0,  1, ],
-     [-1, -1,  2,  3,     1,  0, -1,  1, ],
-     [-2,  1,  1,  3,     1,  0, -1,  1, ],
-     [-1,  2, -1,  3,     0, -1,  1,  1, ],
-     [ 1,  1, -2,  3,     0, -1,  1,  1, ],
-     [-2,  1,  1,  3,     1, -1,  0,  1, ],
-     [-1,  2, -1,  3,     1, -1,  0,  1, ],
-     [ 1,  1, -2,  3,    -1,  0,  1,  1, ],
-     [ 2, -1, -1,  3,    -1,  0,  1,  1, ],
-     [ 1, -2,  1,  3,     0,  1, -1,  1, ],
-     [-1, -1,  2,  3,     0,  1, -1,  1, ],
-    # pyramidal system: c+a slip <11.3>{-1-1.2} -- as for hexagonal ice (Castelnau et al. 1996, similar to twin system found below) 
-     [ 2, -1, -1,  3,    -2,  1,  1,  2, ], # sorted according to similar twin system
-     [-1,  2, -1,  3,     1, -2,  1,  2, ], # <11.3>{-1-1.2} shear = 2((c/a)^2-2)/(3 c/a)
-     [-1, -1,  2,  3,     1,  1, -2,  2, ],
-     [-2,  1,  1,  3,     2, -1, -1,  2, ],
-     [ 1, -2,  1,  3,    -1,  2, -1,  2, ],
-     [ 1,  1, -2,  3,    -1, -1,  2,  2, ],
-     ],'f'),
+    [+2,-1,-1,+0 , +0,+0,+0,+1],
+    [-1,+2,-1,+0 , +0,+0,+0,+1],
+    [-1,-1,+2,+0 , +0,+0,+0,+1],
+    [+2,-1,-1,+0 , +0,+1,-1,+0],
+    [-1,+2,-1,+0 , -1,+0,+1,+0],
+    [-1,-1,+2,+0 , +1,-1,+0,+0],
+    [-1,+1,+0,+0 , +1,+1,-2,+0],
+    [+0,-1,+1,+0 , -2,+1,+1,+0],
+    [+1,+0,-1,+0 , +1,-2,+1,+0],
+    [-1,+2,-1,+0 , +1,+0,-1,+1],
+    [-2,+1,+1,+0 , +0,+1,-1,+1],
+    [-1,-1,+2,+0 , -1,+1,+0,+1],
+    [+1,-2,+1,+0 , -1,+0,+1,+1],
+    [+2,-1,-1,+0 , +0,-1,+1,+1],
+    [+1,+1,-2,+0 , +1,-1,+0,+1],
+    [-2,+1,+1,+3 , +1,+0,-1,+1],
+    [-1,-1,+2,+3 , +1,+0,-1,+1],
+    [-1,-1,+2,+3 , +0,+1,-1,+1],
+    [+1,-2,+1,+3 , +0,+1,-1,+1],
+    [+1,-2,+1,+3 , -1,+1,+0,+1],
+    [+2,-1,-1,+3 , -1,+1,+0,+1],
+    [+2,-1,-1,+3 , -1,+0,+1,+1],
+    [+1,+1,-2,+3 , -1,+0,+1,+1],
+    [+1,+1,-2,+3 , +0,-1,+1,+1],
+    [-1,+2,-1,+3 , +0,-1,+1,+1],
+    [-1,+2,-1,+3 , +1,-1,+0,+1],
+    [-2,+1,+1,+3 , +1,-1,+0,+1],
+    [-1,-1,+2,+3 , +1,+1,-2,+2],
+    [+1,-2,+1,+3 , -1,+2,-1,+2],
+    [+2,-1,-1,+3 , -2,+1,+1,+2],
+    [+1,+1,-2,+3 , -1,-1,+2,+2],
+    [-1,+2,-1,+3 , +1,-2,+1,+2],
+    [-2,+1,+1,+3 , +2,-1,-1,+2],
+   ],'d'),
 }
 
 # --------------------------------------------------------------------
@@ -111,11 +103,11 @@ Add columns listing Schmid factors (and optional trace vector of selected system
 
 """, version = scriptID)
 
-latticeChoices = ('fcc','bcc','hex')
+lattice_choices = list(slipSystems.keys())
 parser.add_option('-l',
                   '--lattice',
-                  dest = 'lattice', type = 'choice', choices = latticeChoices, metavar='string',
-                  help = 'type of lattice structure [%default] {}'.format(latticeChoices))
+                  dest = 'lattice', type = 'choice', choices = lattice_choices, metavar='string',
+                  help = 'type of lattice structure [%default] {}'.format(lattice_choices))
 parser.add_option('--covera',
                   dest = 'CoverA', type = 'float', metavar = 'float',
                   help = 'C over A ratio for hexagonal systems [%default]')
@@ -138,88 +130,63 @@ parser.add_option('-o',
 parser.set_defaults(force = (0.0,0.0,1.0),
                     quaternion='orientation',
                     normal = None,
-                    lattice = latticeChoices[0],
+                    lattice = lattice_choices[0],
                     CoverA = np.sqrt(8./3.),
                    )
 
 (options, filenames) = parser.parse_args()
-
-force = np.array(options.force)
-force /= np.linalg.norm(force)
-
-if options.normal is not None:
-  normal = np.array(options.normal)
-  normal /= np.linalg.norm(normal)
-  if abs(np.dot(force,normal)) > 1e-3:
-    parser.error('stress plane normal not orthogonal to force direction')
-else:
-  normal = force
-
-slip_direction = np.zeros((len(slipSystems[options.lattice]),3),'f')
-slip_normal    = np.zeros_like(slip_direction)
-
-
-if options.lattice in latticeChoices[:2]:
-  slip_direction = slipSystems[options.lattice][:,:3]
-  slip_normal    = slipSystems[options.lattice][:,3:]
-elif options.lattice == latticeChoices[2]:
-  # convert 4 Miller index notation of hex to orthogonal 3 Miller index notation
-  for i in range(len(slip_direction)):
-    slip_direction[i] = np.array([slipSystems['hex'][i,0]*1.5,
-                              (slipSystems['hex'][i,0] + 2.*slipSystems['hex'][i,1])*0.5*np.sqrt(3),
-                               slipSystems['hex'][i,3]*options.CoverA,
-                              ])
-    slip_normal[i]    = np.array([slipSystems['hex'][i,4],
-                              (slipSystems['hex'][i,4] + 2.*slipSystems['hex'][i,5])/np.sqrt(3),
-                               slipSystems['hex'][i,7]/options.CoverA,
-                              ])
-
-slip_direction /= np.tile(np.linalg.norm(slip_direction,axis=1),(3,1)).T
-slip_normal    /= np.tile(np.linalg.norm(slip_normal   ,axis=1),(3,1)).T
-
-# --- loop over input files ------------------------------------------------------------------------
-
 if filenames == []: filenames = [None]
 
+force = np.array(options.force)/np.linalg.norm(options.force)
+
+if options.normal is not None:
+    normal = np.array(options.normal)/np.linalg.norm(options.ormal)
+    if abs(np.dot(force,normal)) > 1e-3:
+          parser.error('stress plane normal not orthogonal to force direction')
+else:
+    normal = force
+
+
+if   options.lattice in ['bcc','fcc']:
+    slip_direction = slipSystems[options.lattice][:,:3]
+    slip_normal    = slipSystems[options.lattice][:,3:]
+elif options.lattice == 'hex':
+    slip_direction = np.zeros((len(slipSystems['hex']),3),'d')
+    slip_normal    = np.zeros_like(slip_direction)
+    # convert 4 Miller index notation of hex to orthogonal 3 Miller index notation
+    for i in range(len(slip_direction)):
+        slip_direction[i] = np.array([slipSystems['hex'][i,0]*1.5,
+                                      (slipSystems['hex'][i,0] + 2.*slipSystems['hex'][i,1])*0.5*np.sqrt(3),
+                                      slipSystems['hex'][i,3]*options.CoverA,
+                                     ])
+        slip_normal[i]    = np.array([slipSystems['hex'][i,4],
+                                      (slipSystems['hex'][i,4] + 2.*slipSystems['hex'][i,5])/np.sqrt(3),
+                                      slipSystems['hex'][i,7]/options.CoverA,
+                                     ])
+
+slip_direction /= np.linalg.norm(slip_direction,axis=1,keepdims=True)
+slip_normal    /= np.linalg.norm(slip_normal,   axis=1,keepdims=True)
+
+labels = ['S[{direction[0]:.1g}_{direction[1]:.1g}_{direction[2]:.1g}]'
+           '({normal[0]:.1g}_{normal[1]:.1g}_{normal[2]:.1g})'\
+           .format(normal = theNormal, direction = theDirection,
+           ) for theNormal,theDirection in zip(slip_normal,slip_direction)]
+
 for name in filenames:
-  try:
-    table = damask.ASCIItable(name = name)
-  except IOError:
-    continue
-  damask.util.report(scriptName,name)
+    damask.util.report(scriptName,name)
 
-# ------------------------------------------ read header ------------------------------------------
+    table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
 
-  table.head_read()
+    o = damask.Rotation.from_quaternion(table.get(options.quaternion))
 
-# ------------------------------------------ sanity checks ----------------------------------------
-  if not table.label_dimension(options.quaternion) == 4:
-    damask.util.croak('input {} does not have dimension 4.'.format(options.quaternion))
-    table.close(dismiss = True)                                                                     # close ASCIItable and remove empty file
-    continue
+    force  = np.broadcast_to(force, o.shape+(3,))
+    normal = np.broadcast_to(normal,o.shape+(3,))
+    slip_direction = np.broadcast_to(slip_direction,o.shape+slip_direction.shape)
+    slip_normal    = np.broadcast_to(slip_normal,   o.shape+slip_normal.shape)
+    S = np.abs(np.einsum('ijk,ik->ij',slip_direction,(o@force))*
+               np.einsum('ijk,ik->ij',slip_normal,   (o@normal)))
 
-  column = table.label_index(options.quaternion)
+    for i,label in enumerate(labels):
+        table.add(label,S[:,i],scriptID+' '+' '.join(sys.argv[1:]))
 
-# ------------------------------------------ assemble header ---------------------------------------
-
-  table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
-  table.labels_append(['S[{direction[0]:.1g}_{direction[1]:.1g}_{direction[2]:.1g}]'
-                       '({normal[0]:.1g}_{normal[1]:.1g}_{normal[2]:.1g})'\
-                       .format(normal = theNormal,
-                               direction = theDirection,
-                              ) for theNormal,theDirection in zip(slip_normal,slip_direction)])
-  table.head_write()
-
-# ------------------------------------------ process data ------------------------------------------
-
-  outputAlive = True
-  while outputAlive and table.data_read():                                                          # read next data line of ASCII table
-    o = damask.Rotation(np.array(list(map(float,table.data[column:column+4]))))
-
-    table.data_append(  np.abs(  np.sum(slip_direction * (o * force) ,axis=1) \
-                               * np.sum(slip_normal    * (o * normal),axis=1)))
-    outputAlive = table.data_write()                                                                # output processed line
-
-# ------------------------------------------ output finalization -----------------------------------  
-
-  table.close()                                                                                     # close ASCII tables
+    table.to_ASCII(sys.stdout if name is None else name)
diff --git a/processing/pre/geom_canvas.py b/processing/pre/geom_canvas.py
index edd5fe622..b309364e8 100755
--- a/processing/pre/geom_canvas.py
+++ b/processing/pre/geom_canvas.py
@@ -40,35 +40,22 @@ parser.add_option('-f','--fill',
 parser.set_defaults(offset = (0,0,0))
 
 (options, filenames) = parser.parse_args()
-
-
 if filenames == []: filenames = [None]
 
+
 for name in filenames:
   damask.util.report(scriptName,name)
 
   geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
-  origin = geom.get_origin()
-  size   = geom.get_size()
-  old = new = geom.get_grid()
   offset = np.asarray(options.offset)
 
   if options.grid is not None:
-    new = np.maximum(1,
-                     np.array([int(o*float(n.lower().replace('x',''))) if n.lower().endswith('x') \
-                          else int(n) for o,n in zip(old,options.grid)],dtype=int))
+      grid = np.maximum(1,
+                        np.array([int(o*float(n.lower().replace('x',''))) if n.lower().endswith('x') \
+                                   else int(n) for o,n in zip(geom.grid,options.grid)],dtype=int))
+  else:
+      grid = np.array(options.grid,dtype=int)
 
-  canvas = np.full(new,options.fill if options.fill is not None
-                  else np.nanmax(geom.microstructure)+1,geom.microstructure.dtype)
-
-  l = np.clip( offset,    0,np.minimum(old  +offset,new))                                           # noqa
-  r = np.clip( offset+old,0,np.minimum(old*2+offset,new))
-  L = np.clip(-offset,    0,np.minimum(new  -offset,old))
-  R = np.clip(-offset+new,0,np.minimum(new*2-offset,old))
-  canvas[l[0]:r[0],l[1]:r[1],l[2]:r[2]] = geom.microstructure[L[0]:R[0],L[1]:R[1],L[2]:R[2]]
-
-
-  damask.util.croak(geom.update(canvas,origin=origin+offset*size/old,rescale=True))
+  damask.util.croak(geom.canvas(grid,np.asarray(options.offset),options.fill))
   geom.add_comments(scriptID + ' ' + ' '.join(sys.argv[1:]))
-
   geom.to_file(sys.stdout if name is None else name,pack=False)
diff --git a/processing/pre/geom_renumber.py b/processing/pre/geom_renumber.py
index 6e51062a5..8eab9064a 100755
--- a/processing/pre/geom_renumber.py
+++ b/processing/pre/geom_renumber.py
@@ -22,14 +22,12 @@ Renumber sorted microstructure indices to 1,...,N.
 """, version=scriptID)
 
 (options, filenames) = parser.parse_args()
-
-
 if filenames == []: filenames = [None]
 
 for name in filenames:
-  damask.util.report(scriptName,name)
+    damask.util.report(scriptName,name)
 
-  geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
-  damask.util.croak(geom.renumber())
-  geom.add_comments(scriptID + ' ' + ' '.join(sys.argv[1:]))
-  geom.to_file(sys.stdout if name is None else name,pack=False)
+    geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
+    damask.util.croak(geom.renumber())
+    geom.add_comments(scriptID + ' ' + ' '.join(sys.argv[1:]))
+    geom.to_file(sys.stdout if name is None else name,pack=False)
diff --git a/processing/pre/geom_rotate.py b/processing/pre/geom_rotate.py
index dc2e8e2a3..9a65a0e2b 100755
--- a/processing/pre/geom_rotate.py
+++ b/processing/pre/geom_rotate.py
@@ -5,7 +5,6 @@ import sys
 from io import StringIO
 from optparse import OptionParser
 
-from scipy import ndimage
 import numpy as np
 
 import damask
@@ -52,47 +51,27 @@ parser.add_option('-f', '--fill',
 parser.set_defaults(degrees = False)
 
 (options, filenames) = parser.parse_args()
-
-if [options.rotation,options.eulers,options.matrix,options.quaternion].count(None) < 3:
-  parser.error('more than one rotation specified.')
-if [options.rotation,options.eulers,options.matrix,options.quaternion].count(None) > 3:
-  parser.error('no rotation specified.')
-
-if options.quaternion is not None:
-  rot = damask.Rotation.from_quaternion(np.array(options.quaternion))                               # we might need P=+1 here, too...
-if options.rotation is not None:
-  rot = damask.Rotation.from_axis_angle(np.array(options.rotation),degrees=options.degrees,normalise=True,P=+1)
-if options.matrix is not None:
-  rot = damask.Rotation.from_matrix(np.array(options.Matrix))
-if options.eulers is not None:
-  rot = damask.Rotation.from_Eulers(np.array(options.eulers),degrees=options.degrees)
-
-eulers = rot.as_Eulers(degrees=True)
-
-
 if filenames == []: filenames = [None]
 
+if [options.rotation,options.eulers,options.matrix,options.quaternion].count(None) < 3:
+    parser.error('more than one rotation specified.')
+if [options.rotation,options.eulers,options.matrix,options.quaternion].count(None) > 3:
+    parser.error('no rotation specified.')
+
+if options.quaternion is not None:
+    rot = damask.Rotation.from_quaternion(np.array(options.quaternion))                             # we might need P=+1 here, too...
+if options.rotation is not None:
+    rot = damask.Rotation.from_axis_angle(np.array(options.rotation),degrees=options.degrees,normalise=True,P=+1)
+if options.matrix is not None:
+    rot = damask.Rotation.from_matrix(np.array(options.Matrix))
+if options.eulers is not None:
+    rot = damask.Rotation.from_Eulers(np.array(options.eulers),degrees=options.degrees)
+
+
 for name in filenames:
-  damask.util.report(scriptName,name)
+    damask.util.report(scriptName,name)
 
-  geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
-  size   = geom.get_size()
-  grid   = geom.get_grid()
-  origin = geom.get_origin()
-  microstructure = geom.get_microstructure()
-  fill = np.nanmax(microstructure)+1 if options.fill is None else options.fill
-  dtype = float if np.isnan(fill) or int(fill) != fill or microstructure.dtype==np.float else int
-
-  # These rotations are always applied in the reference coordinate system, i.e. (z,x,z) not (z,x',z'')
-  # this seems to be ok, see https://www.cs.utexas.edu/~theshark/courses/cs354/lectures/cs354-14.pdf
-  microstructure = ndimage.rotate(microstructure,eulers[2],(0,1),order=0,
-                                  prefilter=False,output=dtype,cval=fill)            # rotation around z
-  microstructure = ndimage.rotate(microstructure,eulers[1],(1,2),order=0,
-                                  prefilter=False,output=dtype,cval=fill)            # rotation around x
-  microstructure = ndimage.rotate(microstructure,eulers[0],(0,1),order=0,
-                                  prefilter=False,output=dtype,cval=fill)            # rotation around z
-
-  damask.util.croak(geom.update(microstructure,origin=origin-(np.asarray(microstructure.shape)-grid)/2*size/grid,rescale=True))
-  geom.add_comments(scriptID + ' ' + ' '.join(sys.argv[1:]))
-
-  geom.to_file(sys.stdout if name is None else name,pack=False)
+    geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
+    damask.util.croak(geom.rotate(rot,options.fill))
+    geom.add_comments(scriptID + ' ' + ' '.join(sys.argv[1:]))
+    geom.to_file(sys.stdout if name is None else name,pack=False)
diff --git a/processing/pre/geom_translate.py b/processing/pre/geom_translate.py
index 1d08c23dc..44c7e8718 100755
--- a/processing/pre/geom_translate.py
+++ b/processing/pre/geom_translate.py
@@ -40,22 +40,16 @@ parser.set_defaults(origin = (0.0,0.0,0.0),
                    )
 
 (options, filenames) = parser.parse_args()
+if filenames == []: filenames = [None]
 
 sub = list(map(int,options.substitute))
 
-
-if filenames == []: filenames = [None]
-
 for name in filenames:
-  damask.util.report(scriptName,name)
+    damask.util.report(scriptName,name)
 
-  geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
-
-  substituted = geom.get_microstructure()
-  for old,new in zip(sub[0::2],sub[1::2]): substituted[geom.microstructure==old] = new              # substitute microstructure indices
-  substituted += options.microstructure                                                             # constant shift
-
-  damask.util.croak(geom.update(substituted,origin=geom.get_origin()+options.origin))
-  geom.add_comments(scriptID + ' ' + ' '.join(sys.argv[1:]))
-
-  geom.to_file(sys.stdout if name is None else name,pack=False)
+    geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
+    geom.renumber(sub[0::2],sub[1::2],origin=geom.origin+options.origin)
+    geom.microstructure+= options.microstructure
+    damask.util.croak(geom)
+    geom.add_comments(scriptID + ' ' + ' '.join(sys.argv[1:]))
+    geom.to_file(sys.stdout if name is None else name,pack=False)
diff --git a/python/damask/__init__.py b/python/damask/__init__.py
index 9a01e8e62..0f343a581 100644
--- a/python/damask/__init__.py
+++ b/python/damask/__init__.py
@@ -1,4 +1,4 @@
-"""Main aggregator."""
+"""Tools for pre and post processing of DAMASK simulations."""
 import os as _os
 import re as _re
 
diff --git a/python/damask/_geom.py b/python/damask/_geom.py
index 7c2fbd82e..fe4e9a17f 100644
--- a/python/damask/_geom.py
+++ b/python/damask/_geom.py
@@ -581,3 +581,58 @@ class Geom:
 
         #self.add_comments('geom.py:renumber v{}'.format(version)
         return self.update(renumbered)
+
+
+    def rotate(self,R,fill=None):
+        """Rotate microstructure (pad if required)."""
+        if fill is None: fill = np.nanmax(self.microstructure) + 1
+        dtype = float if np.isnan(fill) or int(fill) != fill or self.microstructure.dtype==np.float else int
+
+        Eulers = R.as_Eulers(degrees=True)
+        microstructure_in = self.get_microstructure()
+
+        # These rotations are always applied in the reference coordinate system, i.e. (z,x,z) not (z,x',z'')
+        # see https://www.cs.utexas.edu/~theshark/courses/cs354/lectures/cs354-14.pdf
+        for angle,axes in zip(Eulers[::-1], [(0,1),(1,2),(0,1)]):
+            microstructure_out = ndimage.rotate(microstructure_in,angle,axes,order=0,
+                                                prefilter=False,output=dtype,cval=fill)
+            if np.prod(microstructure_in.shape) == np.prod(microstructure_out.shape):
+                # avoid scipy interpolation errors for rotations close to multiples of 90°
+                microstructure_in = np.rot90(microstructure_in,k=np.rint(angle/90.).astype(int),axes=axes)
+            else:
+                microstructure_in = microstructure_out
+
+        origin = self.origin-(np.asarray(microstructure_in.shape)-self.grid)*.5 * self.size/self.grid
+
+        #self.add_comments('geom.py:rotate v{}'.format(version)
+        return self.update(microstructure_in,origin=origin,rescale=True)
+
+
+    def canvas(self,grid=None,offset=None,fill=None):
+        """Crop or enlarge/pad microstructure."""
+        if fill is None: fill = np.nanmax(self.microstructure) + 1
+        if offset is None: offset = 0
+        dtype = float if int(fill) != fill or self.microstructure.dtype==np.float else int
+
+        canvas = np.full(self.grid if grid is None else grid,
+                         fill if fill is not None else np.nanmax(self.microstructure)+1,dtype)
+
+        l = np.clip( offset,          0,np.minimum(self.grid  +offset,grid))                        # noqa
+        r = np.clip( offset+self.grid,0,np.minimum(self.grid*2+offset,grid))
+        L = np.clip(-offset,          0,np.minimum(grid       -offset,self.grid))
+        R = np.clip(-offset+grid,     0,np.minimum(grid*2     -offset,self.grid))
+
+        canvas[l[0]:r[0],l[1]:r[1],l[2]:r[2]] = self.microstructure[L[0]:R[0],L[1]:R[1],L[2]:R[2]]
+
+        #self.add_comments('geom.py:canvas v{}'.format(version)
+        return self.update(canvas,origin=self.origin+offset*self.size/self.grid,rescale=True)
+
+
+    def substitute(self,from_microstructure,to_microstructure):
+        """Substitude microstructure indices."""
+        substituted = self.get_microstructure()
+        for from_ms,to_ms in zip(from_microstructure,to_microstructure):
+            substituted[self.microstructure==from_ms] = to_ms
+
+        #self.add_comments('geom.py:substitute v{}'.format(version)
+        return self.update(substituted)
diff --git a/python/damask/_lattice.py b/python/damask/_lattice.py
index 5c8166508..0f6cffd04 100644
--- a/python/damask/_lattice.py
+++ b/python/damask/_lattice.py
@@ -229,19 +229,20 @@ class Symmetry:
         Return inverse pole figure color if requested.
         Bases are computed from
 
-        basis = {'cubic' :        np.linalg.inv(np.array([[0.,0.,1.],                               # direction of red
-                                                          [1.,0.,1.]/np.sqrt(2.),                   # direction of green
-                                                          [1.,1.,1.]/np.sqrt(3.)]).T),              # direction of blue
-                 'hexagonal' :    np.linalg.inv(np.array([[0.,0.,1.],                               # direction of red
-                                                          [1.,0.,0.],                               # direction of green
-                                                          [np.sqrt(3.),1.,0.]/np.sqrt(4.)]).T),     # direction of blue
-                 'tetragonal' :   np.linalg.inv(np.array([[0.,0.,1.],                               # direction of red
-                                                          [1.,0.,0.],                               # direction of green
-                                                          [1.,1.,0.]/np.sqrt(2.)]).T),              # direction of blue
-                 'orthorhombic' : np.linalg.inv(np.array([[0.,0.,1.],                               # direction of red
-                                                          [1.,0.,0.],                               # direction of green
-                                                          [0.,1.,0.]]).T),                          # direction of blue
-                }
+        >>> basis = {'cubic' :       np.linalg.inv(np.array([[0.,0.,1.],                            # direction of red
+        ...                                                  [1.,0.,1.]/np.sqrt(2.),                # direction of green
+        ...                                                  [1.,1.,1.]/np.sqrt(3.)]).T),           # direction of blue
+        ...          'hexagonal' :   np.linalg.inv(np.array([[0.,0.,1.],                            # direction of red
+        ...                                                  [1.,0.,0.],                            # direction of green
+        ...                                                  [np.sqrt(3.),1.,0.]/np.sqrt(4.)]).T),  # direction of blue
+        ...          'tetragonal' :  np.linalg.inv(np.array([[0.,0.,1.],                            # direction of red
+        ...                                                  [1.,0.,0.],                            # direction of green
+        ...                                                  [1.,1.,0.]/np.sqrt(2.)]).T),           # direction of blue
+        ...          'orthorhombic': np.linalg.inv(np.array([[0.,0.,1.],                            # direction of red
+        ...                                                  [1.,0.,0.],                            # direction of green
+        ...                                                  [0.,1.,0.]]).T),                       # direction of blue
+        ...         }
+
         """
         if self.lattice == 'cubic':
             basis = {'improper':np.array([ [-1.            ,  0.            ,  1. ],
diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py
index 7c3c614fb..e4dbc7d7c 100644
--- a/python/damask/_orientation.py
+++ b/python/damask/_orientation.py
@@ -8,6 +8,7 @@ class Orientation:
     Crystallographic orientation.
 
     A crystallographic orientation contains a rotation and a lattice.
+
     """
 
     __slots__ = ['rotation','lattice']
@@ -49,8 +50,10 @@ class Orientation:
         Disorientation between myself and given other orientation.
 
         Rotation axis falls into SST if SST == True.
-        (Currently requires same symmetry for both orientations.
-         Look into A. Heinz and P. Neumann 1991 for cases with differing sym.)
+
+        Currently requires same symmetry for both orientations.
+        Look into A. Heinz and P. Neumann 1991 for cases with differing sym.
+
         """
         if self.lattice.symmetry != other.lattice.symmetry:
             raise NotImplementedError('disorientation between different symmetry classes not supported yet.')
diff --git a/python/damask/_result.py b/python/damask/_result.py
index 5a22c7977..3a43850ed 100644
--- a/python/damask/_result.py
+++ b/python/damask/_result.py
@@ -322,9 +322,10 @@ class Result:
         Return groups that contain all requested datasets.
 
         Only groups within
-          - inc?????/constituent/*_*/*
-          - inc?????/materialpoint/*_*/*
-          - inc?????/geometry/*
+          - inc*/constituent/*/*
+          - inc*/materialpoint/*/*
+          - inc*/geometry/*
+
         are considered as they contain user-relevant data.
         Single strings will be treated as list with one entry.
 
@@ -599,9 +600,6 @@ class Result:
 
     @staticmethod
     def _add_deviator(T):
-        if not T['data'].shape[1:] == (3,3):
-            raise ValueError
-
         return {
                 'data':  mechanics.deviatoric_part(T['data']),
                 'label': 's_{}'.format(T['label']),
@@ -832,14 +830,12 @@ class Result:
         pole      = np.array(p)
         unit_pole = pole/np.linalg.norm(pole)
         m         = util.scale_to_coprime(pole)
-        coords    = np.empty((len(q['data']),2))
-
-        for i,qu in enumerate(q['data']):
-            o = Rotation(np.array([qu['w'],qu['x'],qu['y'],qu['z']]))
-            rotatedPole = o*unit_pole                                                               # rotate pole according to crystal orientation
-            (x,y) = rotatedPole[0:2]/(1.+abs(unit_pole[2]))                                         # stereographic projection
-            coords[i] = [np.sqrt(x*x+y*y),np.arctan2(y,x)] if polar else [x,y]
+        rot       = Rotation(q['data'].view(np.double).reshape(-1,4))
 
+        rotatedPole = rot @ np.broadcast_to(unit_pole,rot.shape+(3,))                               # rotate pole according to crystal orientation
+        xy = rotatedPole[:,0:2]/(1.+abs(unit_pole[2]))                                              # stereographic projection
+        coords = xy if not polar else \
+                 np.block([np.sqrt(xy[:,0:1]*xy[:,0:1]+xy[:,1:2]*xy[:,1:2]),np.arctan2(xy[:,1:2],xy[:,0:1])])
         return {
                 'data': coords,
                 'label': 'p^{}_[{} {} {})'.format(u'rφ' if polar else 'xy',*m),
@@ -869,8 +865,6 @@ class Result:
 
     @staticmethod
     def _add_rotational_part(F):
-        if not F['data'].shape[1:] == (3,3):
-            raise ValueError
         return {
                 'data':  mechanics.rotational_part(F['data']),
                 'label': 'R({})'.format(F['label']),
@@ -895,9 +889,6 @@ class Result:
 
     @staticmethod
     def _add_spherical(T):
-        if not T['data'].shape[1:] == (3,3):
-            raise ValueError
-
         return {
                 'data':  mechanics.spherical_part(T['data']),
                 'label': 'p_{}'.format(T['label']),
@@ -922,9 +913,6 @@ class Result:
 
     @staticmethod
     def _add_strain_tensor(F,t,m):
-        if not F['data'].shape[1:] == (3,3):
-            raise ValueError
-
         return {
                 'data':  mechanics.strain_tensor(F['data'],t,m),
                 'label': 'epsilon_{}^{}({})'.format(t,m,F['label']),
@@ -956,9 +944,6 @@ class Result:
 
     @staticmethod
     def _add_stretch_tensor(F,t):
-        if not F['data'].shape[1:] == (3,3):
-            raise ValueError
-
         return {
                 'data':  mechanics.left_stretch(F['data']) if t == 'V' else mechanics.right_stretch(F['data']),
                 'label': '{}({})'.format(t,F['label']),
@@ -1041,7 +1026,7 @@ class Result:
         pool.join()
 
 
-    def write_XMDF(self):
+    def write_XDMF(self):
         """
         Write XDMF file to directly visualize data in DADF5 file.
 
@@ -1049,7 +1034,7 @@ class Result:
         Selection is not taken into account.
         """
         if len(self.constituents) != 1 or not self.structured:
-            raise NotImplementedError
+            raise NotImplementedError('XDMF only available for grid results with 1 constituent.')
 
         xdmf=ET.Element('Xdmf')
         xdmf.attrib={'Version':  '2.0',
@@ -1217,14 +1202,6 @@ class Result:
 ###################################################################################################
 # BEGIN DEPRECATED
 
-    def _time_to_inc(self,start,end):
-        selected = []
-        for i,time in enumerate(self.times):
-            if start <= time <= end:
-                selected.append(self.increments[i])
-        return selected
-
-
     def set_by_time(self,start,end):
         """
         Set active increments based on start and end time.
@@ -1237,4 +1214,4 @@ class Result:
           end time (included)
 
         """
-        self._manage_selection('set','increments',self._time_to_inc(start,end))
+        self._manage_selection('set','times',self.times_in_range(start,end))
diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py
index 387ee84cb..f90c5fcc0 100644
--- a/python/damask/_rotation.py
+++ b/python/damask/_rotation.py
@@ -13,30 +13,31 @@ class Rotation:
     u"""
     Orientation stored with functionality for conversion to different representations.
 
+    The following conventions apply:
+
+    - coordinate frames are right-handed.
+    - a rotation angle ω is taken to be positive for a counterclockwise rotation
+      when viewing from the end point of the rotation axis towards the origin.
+    - rotations will be interpreted in the passive sense.
+    - Euler angle triplets are implemented using the Bunge convention,
+      with the angular ranges as [0, 2π],[0, π],[0, 2π].
+    - the rotation angle ω is limited to the interval [0, π].
+    - the real part of a quaternion is positive, Re(q) > 0
+    - P = -1 (as default).
+
+    Examples
+    --------
+    Rotate vector "a" (defined in coordinate system "A") to
+    coordinates "b" expressed in system "B":
+
+    - b = Q @ a
+    - b = np.dot(Q.asMatrix(),a)
+
     References
     ----------
     D. Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015
     https://doi.org/10.1088/0965-0393/23/8/083501
 
-    Conventions
-    -----------
-    Convention 1: Coordinate frames are right-handed.
-    Convention 2: A rotation angle ω is taken to be positive for a counterclockwise rotation
-                  when viewing from the end point of the rotation axis towards the origin.
-    Convention 3: Rotations will be interpreted in the passive sense.
-    Convention 4: Euler angle triplets are implemented using the Bunge convention,
-                  with the angular ranges as [0, 2π],[0, π],[0, 2π].
-    Convention 5: The rotation angle ω is limited to the interval [0, π].
-    Convention 6: the real part of a quaternion is positive, Re(q) > 0
-    Convention 7: P = -1 (as default).
-
-    Usage
-    -----
-    Vector "a" (defined in coordinate system "A") is passively rotated
-               resulting in new coordinates "b" when expressed in system "B".
-    b = Q @ a
-    b = np.dot(Q.as_matrix(),a)
-
     """
 
     __slots__ = ['quaternion']
@@ -160,10 +161,10 @@ class Rotation:
         if self.shape == ():
             q = np.broadcast_to(self.quaternion,shape+(4,))
         else:
-            q = np.block([np.broadcast_to(self.quaternion[...,0:1],shape+(1,)),
-                          np.broadcast_to(self.quaternion[...,1:2],shape+(1,)),
-                          np.broadcast_to(self.quaternion[...,2:3],shape+(1,)),
-                          np.broadcast_to(self.quaternion[...,3:4],shape+(1,))])
+            q = np.block([np.broadcast_to(self.quaternion[...,0:1],shape).reshape(shape+(1,)),
+                          np.broadcast_to(self.quaternion[...,1:2],shape).reshape(shape+(1,)),
+                          np.broadcast_to(self.quaternion[...,2:3],shape).reshape(shape+(1,)),
+                          np.broadcast_to(self.quaternion[...,3:4],shape).reshape(shape+(1,))])
         return self.__class__(q)
 
 
@@ -537,7 +538,7 @@ class Rotation:
                      )
         # reduce Euler angles to definition range
         eu[np.abs(eu)<1.e-6] = 0.0
-        eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu)
+        eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu)                # needed?
         return eu
 
     @staticmethod
diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py
index 4505a5ebc..58b43dc1f 100644
--- a/python/damask/_vtk.py
+++ b/python/damask/_vtk.py
@@ -120,9 +120,9 @@ class VTK:
         Parameters
         ----------
         fname : str
-            Filename for reading. Valid extensions are *.vtr, *.vtu, *.vtp, and *.vtk.
+            Filename for reading. Valid extensions are .vtr, .vtu, .vtp, and .vtk.
         dataset_type : str, optional
-            Name of the vtk.vtkDataSet subclass when opening an *.vtk file. Valid types are vtkRectilinearGrid,
+            Name of the vtk.vtkDataSet subclass when opening an .vtk file. Valid types are vtkRectilinearGrid,
             vtkUnstructuredGrid, and vtkPolyData.
 
         """
diff --git a/python/tests/reference/Geom/rotate_Eulers=0.0-32.0-240.0.geom b/python/tests/reference/Geom/rotate_Eulers=0.0-32.0-240.0.geom
new file mode 100644
index 000000000..a3d1de0b4
--- /dev/null
+++ b/python/tests/reference/Geom/rotate_Eulers=0.0-32.0-240.0.geom
@@ -0,0 +1,85 @@
+4 header
+grid   a 8 b 10 c 8
+size   x 8e-06 y 1e-05 z 8e-06
+origin x 0.0 y -2.5e-06 z -2e-06
+homogenization 1
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 30 14 18 42 42 42
+42 42 29 13 17 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 31 15 19 42 42 42
+42  6 10  2  2 42 42 42
+42  1  2  2  2  2 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 32 16 20 42 42 42
+42  7 11  2  2 42 42 42
+42  7 11  2  2 42 42 42
+42  2  2  2  2  2 42 42
+42 42  2  2 31 35 42 42
+42 42 22 26 10 14  1 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42  4  2  2  2  2 42 42
+42 42  2  2 32 36 42 42
+42 42 24 28 12 16  1 42
+42 42 42  7  7  1  1 42
+42 42 42  7  7  1  1 42
+42 42 42  1  1  1 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 25 29 13 17  1 42
+42 42 42  8  8  1  1 42
+42 42 42  1  1  1 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42  1  1  1 42 42
+42 42 42  1  1  1 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42
diff --git a/python/tests/reference/Geom/rotate_Eulers=32.0-68.0-21.0.geom b/python/tests/reference/Geom/rotate_Eulers=32.0-68.0-21.0.geom
new file mode 100644
index 000000000..e29cd4266
--- /dev/null
+++ b/python/tests/reference/Geom/rotate_Eulers=32.0-68.0-21.0.geom
@@ -0,0 +1,104 @@
+4 header
+grid   a 11 b 11 c 9
+size   x 1.1e-05 y 1.1000000000000001e-05 z 9e-06
+origin x -1.5e-06 y -3.0000000000000005e-06 z -2.4999999999999998e-06
+homogenization 1
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42  1  2 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42  1 42 42 42 42 42 42
+42 42 42 42  1  5 42 42 42 42 42
+42 42 42  1  7  4 42 42 42 42 42
+42 42 42 42 42 27 42 42 42 42 42
+42 42 42 42 42 42  2 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42  1  1 42 42 42 42 42 42
+42 42 42  1  1  9 29 42 42 42 42
+42 42  1  1 11  8 28  2 42 42 42
+42 42 42  1 10 31  2 42 42 42 42
+42 42 42 42 30  2  2  2 42 42 42
+42 42 42 42 42 42  2  1 42 42 42
+42 42 42 42 42 42 42  1 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42  1 42 42 42 42 42 42 42
+42 42 42  1  1 42 42 42 42 42 42
+42 42  1 16 36 12 32 42 42 42 42
+42 42 42 15 35  2  2  2 42 42 42
+42 42 42 42  2  2  2 11  3 42 42
+42 42 42 42 42 42 10  6 42 42 42
+42 42 42 42 42 42 42  6 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42  1 42 42 42 42 42 42 42
+42 42 42  1 17 42 42 42 42 42 42
+42 42 42 16 36 37  2 42 42 42 42
+42 42 42 42 39  2  2 12 42 42 42
+42 42 42 38  2  2  2 11  8 42 42
+42 42 42 42  2  2 14 30 42 42 42
+42 42 42 42 42 42 13 30 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 40 42 42 42 42 42 42
+42 42 42 42 42  2 42 42 42 42 42
+42 42 42 42 42  2  2 15 42 42 42
+42 42 42 42 42  2 18 42 42 42 42
+42 42 42 42 42 42 17 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42  2 20 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
+42 42 42 42 42 42 42 42 42 42 42
diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py
index 0f8c10066..7a3823ab3 100644
--- a/python/tests/test_Geom.py
+++ b/python/tests/test_Geom.py
@@ -5,12 +5,13 @@ import pytest
 import numpy as np
 
 from damask import Geom
+from damask import Rotation
 
 
 def geom_equal(a,b):
     return np.all(a.get_microstructure() == b.get_microstructure()) and \
-           np.all(a.get_size()           == b.get_size())           and \
-           np.all(a.get_grid()           == b.get_grid())
+           np.all(a.get_grid()           == b.get_grid()) and \
+           np.allclose(a.get_size(), b.get_size())
 
 @pytest.fixture
 def default():
@@ -36,6 +37,7 @@ class TestGeom:
                         default.get_size(),
                         default.get_origin()
                        )
+        print(modified)
         assert geom_equal(modified,default)
 
 
@@ -57,6 +59,22 @@ class TestGeom:
         new = Geom.from_file(tmpdir.join('default.geom'))
         assert geom_equal(new,default)
 
+    def test_invalid_combination(self,default):
+        with pytest.raises(ValueError):
+            default.update(default.microstructure[1:,1:,1:],size=np.ones(3), rescale=True)
+
+    def test_invalid_size(self,default):
+        with pytest.raises(ValueError):
+            default.update(default.microstructure[1:,1:,1:],size=np.ones(2))
+
+    def test_invalid_microstructure(self,default):
+        with pytest.raises(ValueError):
+            default.update(default.microstructure[1])
+
+    def test_invalid_homogenization(self,default):
+        with pytest.raises(TypeError):
+            default.set_homogenization(homogenization=0)
+
     @pytest.mark.parametrize('directions,reflect',[
                                                    (['x'],        False),
                                                    (['x','y','z'],True),
@@ -98,6 +116,52 @@ class TestGeom:
         if update: modified.to_file(reference)
         assert geom_equal(modified,Geom.from_file(reference))
 
+    def test_renumber(self,default):
+        modified = copy.deepcopy(default)
+        microstructure = modified.get_microstructure()
+        for m in np.unique(microstructure):
+            microstructure[microstructure==m] = microstructure.max() + np.random.randint(1,30)
+        modified.update(microstructure)
+        assert not geom_equal(modified,default)
+        modified.renumber()
+        assert geom_equal(modified,default)
+
+    def test_substitute(self,default):
+        modified = copy.deepcopy(default)
+        microstructure = modified.get_microstructure()
+        offset = np.random.randint(1,500)
+        microstructure+=offset
+        modified.update(microstructure)
+        assert not geom_equal(modified,default)
+        modified.substitute(np.arange(default.microstructure.max())+1+offset,
+                            np.arange(default.microstructure.max())+1)
+        assert geom_equal(modified,default)
+
+    @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])])
+    def test_rotate360(self,default,axis_angle):
+        modified = copy.deepcopy(default)
+        for i in range(np.rint(360/axis_angle[3]).astype(int)):
+            modified.rotate(Rotation.from_axis_angle(axis_angle,degrees=True))
+        assert geom_equal(modified,default)
+
+    @pytest.mark.parametrize('Eulers',[[32.0,68.0,21.0],
+                                       [0.0,32.0,240.0]])
+    def test_rotate(self,default,update,reference_dir,Eulers):
+        modified = copy.deepcopy(default)
+        modified.rotate(Rotation.from_Eulers(Eulers,degrees=True))
+        tag = 'Eulers={}-{}-{}'.format(*Eulers)
+        reference = os.path.join(reference_dir,'rotate_{}.geom'.format(tag))
+        if update: modified.to_file(reference)
+        assert geom_equal(modified,Geom.from_file(reference))
+
+    def test_canvas(self,default):
+        grid_add = np.random.randint(0,30,(3))
+        modified = copy.deepcopy(default)
+        modified.canvas(modified.grid + grid_add)
+        e = default.grid
+        assert np.all(modified.microstructure[:e[0],:e[1],:e[2]] == default.microstructure)
+
     @pytest.mark.parametrize('periodic',[True,False])
     def test_tessellation_approaches(self,periodic):
         grid   = np.random.randint(10,20,3)
diff --git a/python/tests/test_Lattice.py b/python/tests/test_Lattice.py
new file mode 100644
index 000000000..fb8ce8ea5
--- /dev/null
+++ b/python/tests/test_Lattice.py
@@ -0,0 +1,41 @@
+import random
+
+import pytest
+import numpy as np
+
+from damask import Symmetry
+
+class TestSymmetry:
+
+    @pytest.mark.parametrize('invalid_symmetry',['fcc','bcc','hello'])
+    def test_invalid_symmetry(self,invalid_symmetry):
+        with pytest.raises(KeyError):
+            s = Symmetry(invalid_symmetry)                                                          # noqa
+
+    def test_equal(self):
+        symmetry = random.choice(Symmetry.lattices)
+        print(symmetry)
+        assert Symmetry(symmetry) == Symmetry(symmetry)
+
+    def test_not_equal(self):
+        symmetries = random.sample(Symmetry.lattices,k=2)
+        assert Symmetry(symmetries[0]) != Symmetry(symmetries[1])
+
+    @pytest.mark.parametrize('lattice',Symmetry.lattices)
+    def test_inFZ(self,lattice):
+          assert Symmetry(lattice).inFZ(np.zeros(3))
+
+    @pytest.mark.parametrize('lattice',Symmetry.lattices)
+    def test_inDisorientationSST(self,lattice):
+          assert Symmetry(lattice).inDisorientationSST(np.zeros(3))
+
+    @pytest.mark.parametrize('lattice',Symmetry.lattices)
+    @pytest.mark.parametrize('proper',[True,False])
+    def test_inSST(self,lattice,proper):
+          assert Symmetry(lattice).inSST(np.zeros(3),proper)
+
+    @pytest.mark.parametrize('function',['inFZ','inDisorientationSST'])
+    def test_invalid_argument(self,function):
+        s = Symmetry()                                                                              # noqa
+        with pytest.raises(ValueError):
+            eval('s.{}(np.ones(4))'.format(function))
diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py
index c117f33f7..4eacdee4c 100644
--- a/python/tests/test_Result.py
+++ b/python/tests/test_Result.py
@@ -4,6 +4,7 @@ import os
 import pytest
 import numpy as np
 
+import damask
 from damask import Result
 from damask import mechanics
 
@@ -13,7 +14,7 @@ def default(tmp_path,reference_dir):
     fname = '12grains6x7x8_tensionY.hdf5'
     shutil.copy(os.path.join(reference_dir,fname),tmp_path)
     f = Result(os.path.join(tmp_path,fname))
-    f.set_by_time(20.0,20.0)
+    f.pick('times',20.0)
     return f
 
 @pytest.fixture
@@ -28,12 +29,60 @@ class TestResult:
         print(default)
 
 
-    def test_time_increments(self,default):
-        shape = default.read_dataset(default.get_dataset_location('F'),0).shape
-        default.set_by_time(0.0,20.0)
-        for i in default.iterate('increments'):
-           assert shape == default.read_dataset(default.get_dataset_location('F'),0).shape
+    def test_pick_all(self,default):
+        default.pick('increments',True)
+        a = default.get_dataset_location('F')
+        default.pick('increments','*')
+        b = default.get_dataset_location('F')
+        default.pick('increments',default.incs_in_range(0,np.iinfo(int).max))
+        c = default.get_dataset_location('F')
 
+        default.pick('times',True)
+        d = default.get_dataset_location('F')
+        default.pick('times','*')
+        e = default.get_dataset_location('F')
+        default.pick('times',default.times_in_range(0.0,np.inf))
+        f = default.get_dataset_location('F')
+        assert a == b == c == d == e ==f
+
+    @pytest.mark.parametrize('what',['increments','times','constituents'])                          # ToDo: discuss materialpoints
+    def test_pick_none(self,default,what):
+        default.pick(what,False)
+        a = default.get_dataset_location('F')
+        default.pick(what,[])
+        b = default.get_dataset_location('F')
+
+        assert a == b == []
+
+    @pytest.mark.parametrize('what',['increments','times','constituents'])                          # ToDo: discuss materialpoints
+    def test_pick_more(self,default,what):
+        default.pick(what,False)
+        default.pick_more(what,'*')
+        a = default.get_dataset_location('F')
+
+        default.pick(what,True)
+        b = default.get_dataset_location('F')
+
+        assert a == b
+
+    @pytest.mark.parametrize('what',['increments','times','constituents'])                          # ToDo: discuss materialpoints
+    def test_pick_less(self,default,what):
+        default.pick(what,True)
+        default.pick_less(what,'*')
+        a = default.get_dataset_location('F')
+
+        default.pick(what,False)
+        b = default.get_dataset_location('F')
+
+        assert a == b == []
+
+    def test_pick_invalid(self,default):
+        with pytest.raises(AttributeError):
+            default.pick('invalid',True)
+
+    def test_add_invalid(self,default):
+        with pytest.raises(Exception):
+            default.add_calculation('#invalid#*2')
 
     def test_add_absolute(self,default):
         default.add_absolute('Fe')
@@ -95,6 +144,20 @@ class TestResult:
         in_file   = default.read_dataset(loc['v(sigma)'],0)
         assert np.allclose(in_memory,in_file)
 
+    @pytest.mark.parametrize('d',[[1,0,0],[0,1,0],[0,0,1]])
+    def test_add_IPFcolor(self,default,d):
+        default.add_IPFcolor('orientation',d)
+        loc = {'orientation': default.get_dataset_location('orientation'),
+               'color':       default.get_dataset_location('IPFcolor_[{} {} {}]'.format(*d))}
+        qu = default.read_dataset(loc['orientation']).view(np.double).reshape(-1,4)
+        crystal_structure = default.get_crystal_structure()
+        in_memory = np.empty((qu.shape[0],3),np.uint8)
+        for i,q in enumerate(qu):
+            o = damask.Orientation(q,crystal_structure).reduced()
+            in_memory[i] = np.uint8(o.IPFcolor(np.array(d))*255)
+        in_file = default.read_dataset(loc['color'])
+        assert np.allclose(in_memory,in_file)
+
     def test_add_maximum_shear(self,default):
         default.add_Cauchy('P','F')
         default.add_maximum_shear('sigma')
@@ -143,6 +206,20 @@ class TestResult:
         in_file   = default.read_dataset(loc['S'],0)
         assert np.allclose(in_memory,in_file)
 
+    @pytest.mark.parametrize('polar',[True,False])
+    def test_add_pole(self,default,polar):
+        pole = np.array([1.,0.,0.])
+        default.add_pole('orientation',pole,polar)
+        loc = {'orientation': default.get_dataset_location('orientation'),
+               'pole':        default.get_dataset_location('p^{}_[1 0 0)'.format(u'rφ' if polar else 'xy'))}
+        rot = damask.Rotation(default.read_dataset(loc['orientation']).view(np.double))
+        rotated_pole = rot * np.broadcast_to(pole,rot.shape+(3,))
+        xy = rotated_pole[:,0:2]/(1.+abs(pole[2]))
+        in_memory = xy if not polar else \
+                    np.block([np.sqrt(xy[:,0:1]*xy[:,0:1]+xy[:,1:2]*xy[:,1:2]),np.arctan2(xy[:,1:2],xy[:,0:1])])
+        in_file = default.read_dataset(loc['pole'])
+        assert np.allclose(in_memory,in_file)
+
     def test_add_rotational_part(self,default):
         default.add_rotational_part('F')
         loc = {'F':    default.get_dataset_location('F'),
@@ -185,3 +262,7 @@ class TestResult:
         in_memory = mechanics.left_stretch(default.read_dataset(loc['F'],0))
         in_file   = default.read_dataset(loc['V(F)'],0)
         assert np.allclose(in_memory,in_file)
+
+    @pytest.mark.parametrize('output',['F',[],['F','P']])
+    def test_vtk(self,default,output):
+        default.to_vtk(output)
diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py
index a9768afb8..5a1cd145d 100644
--- a/python/tests/test_Rotation.py
+++ b/python/tests/test_Rotation.py
@@ -866,12 +866,23 @@ class TestRotation:
         with pytest.raises(ValueError):
             function(invalid_shape)
 
+    @pytest.mark.parametrize('shape',[None,(3,),(4,2)])
+    def test_broadcast(self,shape):
+        rot = Rotation.from_random(shape)
+        new_shape = tuple(np.random.randint(8,32,(3))) if shape is None else \
+                    rot.shape + (np.random.randint(8,32),)
+        rot_broadcast = rot.broadcast_to(tuple(new_shape))
+        for i in range(rot_broadcast.shape[-1]):
+            assert np.allclose(rot_broadcast.quaternion[...,i,:], rot.quaternion)
+
+
     @pytest.mark.parametrize('function,invalid',[(Rotation.from_quaternion, np.array([-1,0,0,0])),
                                                  (Rotation.from_quaternion, np.array([1,1,1,0])),
                                                  (Rotation.from_Eulers,     np.array([1,4,0])),
                                                  (Rotation.from_axis_angle, np.array([1,0,0,4])),
                                                  (Rotation.from_axis_angle, np.array([1,1,0,1])),
                                                  (Rotation.from_matrix,     np.random.rand(3,3)),
+                                                 (Rotation.from_matrix,     np.array([[1,1,0],[1,2,0],[0,0,1]])),
                                                  (Rotation.from_Rodrigues,  np.array([1,0,0,-1])),
                                                  (Rotation.from_Rodrigues,  np.array([1,1,0,1])),
                                                  (Rotation.from_homochoric, np.array([2,2,2])),