diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 6f18e3f1e..1bf5cf8b1 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -45,7 +45,7 @@ variables:
   MPI_INTEL:      "MPI/Intel/2022.0.1/IntelMPI/2021.5.0"
   # ++++++++++++ PETSc ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   PETSC_GNU:       "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1"
-  PETSC_INTELLLVM: "Libraries/PETSc/3.16.2/oneAPI-2022.0.1-IntelMPI-2021.5.0"
+  PETSC_INTELLLVM: "Libraries/PETSc/3.16.3/oneAPI-2022.0.1-IntelMPI-2021.5.0"
   PETSC_INTEL:    "Libraries/PETSc/3.16.2/Intel-2022.0.1-IntelMPI-2021.5.0"
   # ++++++++++++ MSC Marc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   MSC:            "FEM/MSC/2021.3.1"
diff --git a/CMakeLists.txt b/CMakeLists.txt
index e8b774cfc..8c7b129f0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -88,16 +88,12 @@ else()
     message(FATAL_ERROR "Compiler type(CMAKE_Fortran_COMPILER_ID) not recognized")
 endif()
 
-file(STRINGS "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/petsc/conf/petscvariables" PETSC_EXTERNAL_LIB REGEX "PETSC_WITH_EXTERNAL_LIB = .*$?")
-string(REGEX MATCHALL "-[lLW]([^\" ]+)" PETSC_EXTERNAL_LIB "${PETSC_EXTERNAL_LIB}")
-list(REMOVE_DUPLICATES PETSC_EXTERNAL_LIB)
-string(REPLACE ";" " " PETSC_EXTERNAL_LIB "${PETSC_EXTERNAL_LIB}")
+file(STRINGS "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/petsc/conf/petscvariables" PETSC_EXTERNAL_LIB REGEX "PETSC_EXTERNAL_LIB_BASIC = .*$?")
+string(REPLACE "PETSC_EXTERNAL_LIB_BASIC = " "" PETSC_EXTERNAL_LIB "${PETSC_EXTERNAL_LIB}")
 message("PETSC_EXTERNAL_LIB:\n${PETSC_EXTERNAL_LIB}\n")
 
 file(STRINGS "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/petsc/conf/petscvariables" PETSC_INCLUDES REGEX "PETSC_FC_INCLUDES = .*$?")
-string(REGEX MATCHALL "-I([^\" ]+)" PETSC_INCLUDES "${PETSC_INCLUDES}")
-list(REMOVE_DUPLICATES PETSC_INCLUDES)
-string(REPLACE ";" " " PETSC_INCLUDES "${PETSC_INCLUDES}")
+string(REPLACE "PETSC_FC_INCLUDES = " "" PETSC_INCLUDES "${PETSC_INCLUDES}")
 message("PETSC_INCLUDES:\n${PETSC_INCLUDES}\n")
 
 set(CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${BUILDCMD_PRE} ${OPENMP_FLAGS} ${STANDARD_CHECK} ${OPTIMIZATION_FLAGS} ${COMPILE_FLAGS} ${PRECISION_FLAGS}")
@@ -109,7 +105,7 @@ if(CMAKE_BUILD_TYPE STREQUAL "DEBUG")
 endif()
 
 set(CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}   "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${PETSC_INCLUDES} ${BUILDCMD_POST}")
-set(CMAKE_Fortran_LINK_EXECUTABLE             "${CMAKE_Fortran_LINK_EXECUTABLE} <OBJECTS> -o <TARGET> <LINK_LIBRARIES> ${PETSC_EXTERNAL_LIB} -lz ${BUILDCMD_POST}")
+set(CMAKE_Fortran_LINK_EXECUTABLE             "${CMAKE_Fortran_LINK_EXECUTABLE} <OBJECTS> -o <TARGET> <LINK_LIBRARIES> -L${PETSC_LIBRARY_DIRS} -lpetsc ${PETSC_EXTERNAL_LIB} -lz ${BUILDCMD_POST}")
 
 message("Fortran Compiler Flags:\n${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}}\n")
 message("C Compiler Flags:\n${CMAKE_C_FLAGS_${CMAKE_BUILD_TYPE}}\n")
diff --git a/Makefile b/Makefile
index e73205702..d5ba1bdad 100644
--- a/Makefile
+++ b/Makefile
@@ -10,14 +10,12 @@ all: grid mesh
 .PHONY: grid
 grid:
 	@cmake -B build/grid -DDAMASK_SOLVER=grid -DCMAKE_INSTALL_PREFIX=${PWD} -DCMAKE_BUILD_TYPE=${BUILD_TYPE} -DBUILDCMD_POST=${BUILDCMD_POST} -DBUILDCMD_PRE=${BUILDCMD_PRE} -DOPTIMIZATION=${OPTIMIZATION} -DOPENMP=${OPENMP}
-	@cmake --build build/grid --parallel
-	@cmake --install build/grid
+	@cmake --build build/grid --parallel --target install
 
 .PHONY: mesh
 mesh:
 	@cmake -B build/mesh -DDAMASK_SOLVER=mesh -DCMAKE_INSTALL_PREFIX=${PWD} -DCMAKE_BUILD_TYPE=${BUILD_TYPE} -DBUILDCMD_POST=${BUILDCMD_POST} -DBUILDCMD_PRE=${BUILDCMD_PRE} -DOPTIMIZATION=${OPTIMIZATION} -DOPENMP=${OPENMP}
-	@cmake --build build/mesh --parallel
-	@cmake --install build/mesh
+	@cmake --build build/mesh --parallel --target install
 
 .PHONY: clean
 clean:
diff --git a/install/MarcMentat/apply_DAMASK_modifications.py b/install/MarcMentat/apply_DAMASK_modifications.py
index 1aadf5313..dd5a3a413 100755
--- a/install/MarcMentat/apply_DAMASK_modifications.py
+++ b/install/MarcMentat/apply_DAMASK_modifications.py
@@ -67,9 +67,7 @@ os.system(f'xvfb-run -a {executable} -compile {menu_file}')
 
 print('setting file access rights...')
 
-files = (glob.glob(str(marc_root/f'marc{marc_version}/tools/*_damask*')) +
-         glob.glob(str(marc_root/f'mentat{marc_version}/bin/kill[4-6]')) +
-         glob.glob(str(marc_root/f'mentat{marc_version}/bin/submit[4-6]')))
-
-for file in files:
+for file in (glob.glob(str(marc_root/f'marc{marc_version}/tools/*_damask*')) +
+             glob.glob(str(marc_root/f'mentat{marc_version}/bin/kill[4-6]')) +
+             glob.glob(str(marc_root/f'mentat{marc_version}/bin/submit[4-6]'))):
     os.chmod(file , 0o755)
diff --git a/processing/legacy/addDisplacement.py b/processing/legacy/addDisplacement.py
deleted file mode 100755
index 6fe577ec7..000000000
--- a/processing/legacy/addDisplacement.py
+++ /dev/null
@@ -1,71 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-parser = OptionParser(usage='%prog options [ASCIItable(s)]', description = """
-Add displacments resulting from deformation gradient field.
-Operates on periodic three-dimensional x,y,z-ordered data sets.
-Outputs at cell centers or cell nodes (into separate file).
-
-""", version = scriptID)
-
-parser.add_option('-f',
-                  '--defgrad',
-                  dest    = 'f',
-                  metavar = 'string',
-                  help    = 'label of deformation gradient [%default]')
-parser.add_option('-p',
-                  '--pos', '--position',
-                  dest    = 'pos',
-                  metavar = 'string',
-                  help    = 'label of coordinates [%default]')
-parser.add_option('--nodal',
-                  dest    = 'nodal',
-                  action  = 'store_true',
-                  help    = 'output nodal (instead of cell-centered) displacements')
-
-parser.set_defaults(f   = 'f',
-                    pos = 'pos',
-                   )
-
-(options,filenames) = parser.parse_args()
-
-for name in filenames:
-    damask.util.report(scriptName,name)
-
-    table = damask.Table.load(StringIO(''.join(sys.stdin.read())) if name is None else name)
-    grid,size,origin = damask.grid_filters.cellsSizeOrigin_coordinates0_point(table.get(options.pos))
-
-    F = table.get(options.f).reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+(3,3))
-    if options.nodal:
-        damask.Table(damask.grid_filters.coordinates0_node(grid,size).reshape(-1,3,order='F'),
-                     {'pos':(3,)})\
-              .add('avg({}).{}'.format(options.f,options.pos),
-                   damask.grid_filters.displacement_avg_node(size,F).reshape(-1,3,order='F'),
-                   scriptID+' '+' '.join(sys.argv[1:]))\
-              .add('fluct({}).{}'.format(options.f,options.pos),
-                   damask.grid_filters.displacement_fluct_node(size,F).reshape(-1,3,order='F'),
-                    scriptID+' '+' '.join(sys.argv[1:]))\
-               .save((sys.stdout if name is None else os.path.splitext(name)[0]+'_nodal.txt'))
-    else:
-        table.add('avg({}).{}'.format(options.f,options.pos),
-                  damask.grid_filters.displacement_avg_point(size,F).reshape(-1,3,order='F'),
-                  scriptID+' '+' '.join(sys.argv[1:]))\
-             .add('fluct({}).{}'.format(options.f,options.pos),
-                  damask.grid_filters.displacement_fluct_point(size,F).reshape(-1,3,order='F'),
-                  scriptID+' '+' '.join(sys.argv[1:]))\
-             .save((sys.stdout if name is None else name))
diff --git a/python/damask/VERSION b/python/damask/VERSION
index a467f728c..d9749c5a9 100644
--- a/python/damask/VERSION
+++ b/python/damask/VERSION
@@ -1 +1 @@
-v3.0.0-alpha5-355-gc29428a60
+v3.0.0-alpha5-375-g76fe2d2b3
diff --git a/python/damask/__init__.py b/python/damask/__init__.py
index 231fa8b30..584a97e87 100644
--- a/python/damask/__init__.py
+++ b/python/damask/__init__.py
@@ -8,6 +8,7 @@ with open(_Path(__file__).parent/_Path('VERSION')) as _f:
     version = _re.sub(r'^v','',_f.readline().strip())
     __version__ = version
 
+from .                 import _typehints       # noqa
 from .                 import util             # noqa
 from .                 import seeds            # noqa
 from .                 import tensor           # noqa
diff --git a/python/damask/_colormap.py b/python/damask/_colormap.py
index c2721e0fa..2da92ae3f 100644
--- a/python/damask/_colormap.py
+++ b/python/damask/_colormap.py
@@ -3,13 +3,9 @@ import json
 import functools
 import colorsys
 from pathlib import Path
-from typing import Sequence, Union, TextIO
+from typing import Union, TextIO
 
 import numpy as np
-try:
-    from numpy.typing import ArrayLike
-except ImportError:
-    ArrayLike = Union[np.ndarray,Sequence[float]] # type: ignore
 import scipy.interpolate as interp
 import matplotlib as mpl
 if os.name == 'posix' and 'DISPLAY' not in os.environ:
@@ -18,6 +14,7 @@ import matplotlib.pyplot as plt
 from matplotlib import cm
 from PIL import Image
 
+from ._typehints import FloatSequence, FileHandle
 from . import util
 from . import Table
 
@@ -82,8 +79,8 @@ class Colormap(mpl.colors.ListedColormap):
 
 
     @staticmethod
-    def from_range(low: ArrayLike,
-                   high: ArrayLike,
+    def from_range(low: FloatSequence,
+                   high: FloatSequence,
                    name: str = 'DAMASK colormap',
                    N: int = 256,
                    model: str = 'rgb') -> 'Colormap':
@@ -197,7 +194,7 @@ class Colormap(mpl.colors.ListedColormap):
 
 
     def at(self,
-           fraction : Union[float,Sequence[float]]) -> np.ndarray:
+           fraction : Union[float,FloatSequence]) -> np.ndarray:
         """
         Interpolate color at fraction.
 
@@ -229,7 +226,7 @@ class Colormap(mpl.colors.ListedColormap):
 
     def shade(self,
               field: np.ndarray,
-              bounds: ArrayLike = None,
+              bounds: FloatSequence = None,
               gap: float = None) -> Image:
         """
         Generate PIL image of 2D field using colormap.
@@ -296,7 +293,7 @@ class Colormap(mpl.colors.ListedColormap):
 
 
     def _get_file_handle(self,
-                         fname: Union[TextIO, str, Path, None],
+                         fname: Union[FileHandle, None],
                          suffix: str = '') -> TextIO:
         """
         Provide file handle.
@@ -323,7 +320,7 @@ class Colormap(mpl.colors.ListedColormap):
             return fname
 
 
-    def save_paraview(self, fname: Union[TextIO, str, Path] = None):
+    def save_paraview(self, fname: FileHandle = None):
         """
         Save as JSON file for use in Paraview.
 
@@ -350,7 +347,7 @@ class Colormap(mpl.colors.ListedColormap):
         fhandle.write('\n')
 
 
-    def save_ASCII(self, fname: Union[TextIO, str, Path] = None):
+    def save_ASCII(self, fname: FileHandle = None):
         """
         Save as ASCII file.
 
@@ -365,7 +362,7 @@ class Colormap(mpl.colors.ListedColormap):
         t.save(self._get_file_handle(fname,'.txt'))
 
 
-    def save_GOM(self, fname: Union[TextIO, str, Path] = None):
+    def save_GOM(self, fname: FileHandle = None):
         """
         Save as ASCII file for use in GOM Aramis.
 
@@ -385,7 +382,7 @@ class Colormap(mpl.colors.ListedColormap):
         self._get_file_handle(fname,'.legend').write(GOM_str)
 
 
-    def save_gmsh(self, fname: Union[TextIO, str, Path] = None):
+    def save_gmsh(self, fname: FileHandle = None):
         """
         Save as ASCII file for use in gmsh.
 
diff --git a/python/damask/_crystal.py b/python/damask/_crystal.py
index 689398cf6..7a21f9245 100644
--- a/python/damask/_crystal.py
+++ b/python/damask/_crystal.py
@@ -114,12 +114,13 @@ class Crystal():
 
     def __repr__(self):
         """Represent."""
-        return '\n'.join([f'Crystal family {self.family}']
-                       + ([] if self.lattice is None else [f'Bravais lattice {self.lattice}']+
-                                                           list(map(lambda x:f'{x[0]}: {x[1]:.5g}',
-                                                                    zip(['a','b','c','α','β','γ',],
-                                                                         self.parameters))))
-                        )
+        family = f'Crystal family: {self.family}'
+        return family if self.lattice is None else \
+               '\n'.join([family,
+                          f'Bravais lattice: {self.lattice}',
+                          'a={:.5g}m, b={:.5g}m, c={:.5g}m'.format(*self.parameters[:3]),
+                          'α={:.5g}°, β={:.5g}°, γ={:.5g}°'.format(*np.degrees(self.parameters[3:]))])
+
 
     def __eq__(self,other):
         """
@@ -378,7 +379,7 @@ class Crystal():
         """
         _kinematics = {
             'cF': {
-                'slip' :[np.array([
+                'slip': [np.array([
                            [+0,+1,-1, +1,+1,+1],
                            [-1,+0,+1, +1,+1,+1],
                            [+1,-1,+0, +1,+1,+1],
@@ -398,7 +399,7 @@ class Crystal():
                            [+1,+0,-1, +1,+0,+1],
                            [+0,+1,+1, +0,+1,-1],
                            [+0,+1,-1, +0,+1,+1]])],
-                'twin' :[np.array([
+                'twin': [np.array([
                            [-2, 1, 1,  1, 1, 1],
                            [ 1,-2, 1,  1, 1, 1],
                            [ 1, 1,-2,  1, 1, 1],
@@ -413,7 +414,7 @@ class Crystal():
                            [-1, 1, 2, -1, 1,-1]])]
             },
             'cI': {
-                'slip' :[np.array([
+                'slip': [np.array([
                            [+1,-1,+1, +0,+1,+1],
                            [-1,-1,+1, +0,+1,+1],
                            [+1,+1,+1, +0,-1,+1],
@@ -464,7 +465,7 @@ class Crystal():
                            [+1,+1,+1, -3,+2,+1],
                            [+1,+1,-1, +3,-2,+1],
                            [+1,-1,+1, +3,+2,-1]])],
-                'twin' :[np.array([
+                'twin': [np.array([
                            [-1, 1, 1,  2, 1, 1],
                            [ 1, 1, 1, -2, 1, 1],
                            [ 1, 1,-1,  2,-1, 1],
@@ -479,7 +480,7 @@ class Crystal():
                            [ 1, 1, 1,  1, 1,-2]])]
             },
             'hP': {
-                'slip' :[np.array([
+                'slip': [np.array([
                            [+2,-1,-1,+0, +0,+0,+0,+1],
                            [-1,+2,-1,+0, +0,+0,+0,+1],
                            [-1,-1,+2,+0, +0,+0,+0,+1]]),
@@ -514,7 +515,7 @@ class Crystal():
                            [+1,+1,-2,+3, -1,-1,+2,+2],
                            [-1,+2,-1,+3, +1,-2,+1,+2],
                            [-2,+1,+1,+3, +2,-1,-1,+2]])],
-                'twin' :[np.array([
+                'twin': [np.array([
                            [-1, 0, 1, 1,  1, 0,-1, 2],   # shear = (3-(c/a)^2)/(sqrt(3) c/a) <-10.1>{10.2}
                            [ 0,-1, 1, 1,  0, 1,-1, 2],
                            [ 1,-1, 0, 1, -1, 1, 0, 2],
@@ -542,7 +543,74 @@ class Crystal():
                            [-1,-1, 2,-3, -1,-1, 2, 2],
                            [ 1,-2, 1,-3,  1,-2, 1, 2],
                            [ 2,-1,-1,-3,  2,-1,-1, 2]])]
-                },
+            },
+            'tI': {
+                'slip': [np.array([
+                           [+0,+0,+1, +1,+0,+0],
+                           [+0,+0,+1, +0,+1,+0]]),
+                         np.array([
+                           [+0,+0,+1, +1,+1,+0],
+                           [+0,+0,+1, -1,+1,+0]]),
+                         np.array([
+                           [+0,+1,+0, +1,+0,+0],
+                           [+1,+0,+0, +0,+1,+0]]),
+                         np.array([
+                           [+1,-1,+1, +1,+1,+0],
+                           [+1,-1,-1, +1,+1,+0],
+                           [-1,-1,-1, -1,+1,+0],
+                           [-1,-1,+1, -1,+1,+0]]),
+                         np.array([
+                           [+1,-1,+0, +1,+1,+0],
+                           [+1,+1,+0, +1,-1,+0]]),
+                         np.array([
+                           [+0,+1,+1, +1,+0,+0],
+                           [+0,-1,+1, +1,+0,+0],
+                           [-1,+0,+1, +0,+1,+0],
+                           [+1,+0,+1, +0,+1,+0]]),
+                         np.array([
+                           [+0,+1,+0, +0,+0,+1],
+                           [+1,+0,+0, +0,+0,+1]]),
+                         np.array([
+                           [+1,+1,+0, +0,+0,+1],
+                           [-1,+1,+0, +0,+0,+1]]),
+                         np.array([
+                           [+0,+1,-1, +0,+1,+1],
+                           [+0,-1,-1, +0,-1,+1],
+                           [-1,+0,-1, -1,+0,+1],
+                           [+1,+0,-1, +1,+0,+1]]),
+                         np.array([
+                           [+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]]),
+                         np.array([
+                           [+1,+0,+0, +0,+1,+1],
+                           [+1,+0,+0, +0,+1,-1],
+                           [+0,+1,+0, +1,+0,+1],
+                           [+0,+1,+0, +1,+0,-1]]),
+                         np.array([
+                           [+0,+1,-1, +2,+1,+1],
+                           [+0,-1,-1, +2,-1,+1],
+                           [+1,+0,-1, +1,+2,+1],
+                           [-1,+0,-1, -1,+2,+1],
+                           [+0,+1,-1, -2,+1,+1],
+                           [+0,-1,-1, -2,-1,+1],
+                           [-1,+0,-1, -1,-2,+1],
+                           [+1,+0,-1, +1,-2,+1]]),
+                         np.array([
+                           [-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, -2,+1,+1],
+                           [+1,-1,+1, -2,-1,+1],
+                           [-1,+1,+1, -1,-2,+1],
+                           [+1,+1,+1, +1,-2,+1]])]
+                }
         }
         master = _kinematics[self.lattice][mode]
         if self.lattice == 'hP':
diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py
index e727c54ae..2fc2a7d4a 100644
--- a/python/damask/_orientation.py
+++ b/python/damask/_orientation.py
@@ -514,6 +514,17 @@ class Orientation(Rotation,Crystal):
          [ 0.07359167 -0.36505797  0.92807163]]
         Bunge Eulers / deg: (11.40, 21.86, 0.60)
 
+        Plot a sample from the Mackenzie distribution.
+
+        >>> import matplotlib.pyplot as plt
+        >>> import damask
+        >>> N = 10000
+        >>> a = damask.Orientation.from_random(shape=N,family='cubic')
+        >>> b = damask.Orientation.from_random(shape=N,family='cubic')
+        >>> d = a.disorientation(b).as_axis_angle(degrees=True,pair=True)[1]
+        >>> plt.hist(d,25)
+        >>> plt.show()
+
         """
         if self.family != other.family:
             raise NotImplementedError('disorientation between different crystal families')
diff --git a/python/damask/_typehints.py b/python/damask/_typehints.py
new file mode 100644
index 000000000..67e920957
--- /dev/null
+++ b/python/damask/_typehints.py
@@ -0,0 +1,11 @@
+"""Functionality for typehints."""
+
+from typing import Sequence, Union, TextIO
+from pathlib import Path
+
+import numpy as np
+
+
+FloatSequence = Union[np.ndarray,Sequence[float]]
+IntSequence = Union[np.ndarray,Sequence[int]]
+FileHandle = Union[TextIO, str, Path]
diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py
index 9d0e055c8..d9d658d0b 100644
--- a/python/damask/grid_filters.py
+++ b/python/damask/grid_filters.py
@@ -12,21 +12,23 @@ the following operations are required for tensorial data:
 
 """
 
-from typing import Sequence, Tuple, Union
+from typing import Tuple as _Tuple
 
 from scipy import spatial as _spatial
 import numpy as _np
 
+from ._typehints import FloatSequence as _FloatSequence, IntSequence as _IntSequence
 
-def _ks(size: _np.ndarray, cells: Union[_np.ndarray,Sequence[int]], first_order: bool = False) -> _np.ndarray:
+
+def _ks(size: _FloatSequence, cells: _IntSequence, first_order: bool = False) -> _np.ndarray:
     """
     Get wave numbers operator.
 
     Parameters
     ----------
-    size : numpy.ndarray of shape (3)
+    size : sequence of float, len (3)
         Physical size of the periodic field.
-    cells : numpy.ndarray of shape (3)
+    cells : sequence of int, len (3)
         Number of cells.
     first_order : bool, optional
         Correction for first order derivatives, defaults to False.
@@ -45,20 +47,20 @@ def _ks(size: _np.ndarray, cells: Union[_np.ndarray,Sequence[int]], first_order:
     return _np.stack(_np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij'), axis=-1)
 
 
-def curl(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray:
+def curl(size: _FloatSequence, f: _np.ndarray) -> _np.ndarray:
     u"""
     Calculate curl of a vector or tensor field in Fourier space.
 
     Parameters
     ----------
-    size : numpy.ndarray of shape (3)
+    size : sequence of float, len (3)
         Physical size of the periodic field.
-    f : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3)
+    f : numpy.ndarray, shape (:,:,:,3) or (:,:,:,3,3)
         Periodic field of which the curl is calculated.
 
     Returns
     -------
-    ∇ × f : numpy.ndarray
+    ∇ × f : numpy.ndarray, shape (:,:,:,3) or (:,:,:,3,3)
         Curl of f.
 
     """
@@ -76,20 +78,20 @@ def curl(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray:
     return _np.fft.irfftn(curl_,axes=(0,1,2),s=f.shape[:3])
 
 
-def divergence(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray:
+def divergence(size: _FloatSequence, f: _np.ndarray) -> _np.ndarray:
     u"""
     Calculate divergence of a vector or tensor field in Fourier space.
 
     Parameters
     ----------
-    size : numpy.ndarray of shape (3)
+    size : sequence of float, len (3)
         Physical size of the periodic field.
-    f : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3)
+    f : numpy.ndarray, shape (:,:,:,3) or (:,:,:,3,3)
         Periodic field of which the divergence is calculated.
 
     Returns
     -------
-    ∇ · f : numpy.ndarray
+    ∇ · f : numpy.ndarray, shape (:,:,:,1) or (:,:,:,3)
         Divergence of f.
 
     """
@@ -103,20 +105,20 @@ def divergence(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray:
     return _np.fft.irfftn(div_,axes=(0,1,2),s=f.shape[:3])
 
 
-def gradient(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray:
+def gradient(size: _FloatSequence, f: _np.ndarray) -> _np.ndarray:
     u"""
-    Calculate gradient of a scalar or vector fieldin Fourier space.
+    Calculate gradient of a scalar or vector field in Fourier space.
 
     Parameters
     ----------
-    size : numpy.ndarray of shape (3)
+    size : sequence of float, len (3)
         Physical size of the periodic field.
-    f : numpy.ndarray of shape (:,:,:,1) or (:,:,:,3)
+    f : numpy.ndarray, shape (:,:,:,1) or (:,:,:,3)
         Periodic field of which the gradient is calculated.
 
     Returns
     -------
-    ∇ f : numpy.ndarray
+    ∇ f : numpy.ndarray, shape (:,:,:,3) or (:,:,:,3,3)
         Divergence of f.
 
     """
@@ -130,29 +132,30 @@ def gradient(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray:
     return _np.fft.irfftn(grad_,axes=(0,1,2),s=f.shape[:3])
 
 
-def coordinates0_point(cells: Union[_np.ndarray, Sequence[int]],
-                       size: Union[_np.ndarray, Sequence[float]],
-                       origin: Union[_np.ndarray, Sequence[float]] = _np.zeros(3)) -> _np.ndarray:
+def coordinates0_point(cells: _IntSequence,
+                       size: _FloatSequence,
+                       origin: _FloatSequence = _np.zeros(3)) -> _np.ndarray:
     """
     Cell center positions (undeformed).
 
     Parameters
     ----------
-    cells : numpy.ndarray of shape (3)
+    cells : sequence of int, len (3)
         Number of cells.
-    size : numpy.ndarray of shape (3)
+    size : sequence of float, len (3)
         Physical size of the periodic field.
-    origin : numpy.ndarray, optional
+    origin : sequence of float, len(3), optional
         Physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
 
     Returns
     -------
-    x_p_0 : numpy.ndarray
+    x_p_0 : numpy.ndarray, shape (:,:,:,3)
         Undeformed cell center coordinates.
 
     """
-    start = _np.array(origin)                   + size/_np.array(cells)*.5
-    end   = _np.array(origin) + _np.array(size) - size/_np.array(cells)*.5
+    size_ = _np.array(size,float)
+    start = origin         + size_/_np.array(cells,int)*.5
+    end   = origin + size_ - size_/_np.array(cells,int)*.5
 
     return _np.stack(_np.meshgrid(_np.linspace(start[0],end[0],cells[0]),
                                   _np.linspace(start[1],end[1],cells[1]),
@@ -160,24 +163,24 @@ def coordinates0_point(cells: Union[_np.ndarray, Sequence[int]],
                      axis = -1)
 
 
-def displacement_fluct_point(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
+def displacement_fluct_point(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
     """
     Cell center displacement field from fluctuation part of the deformation gradient field.
 
     Parameters
     ----------
-    size : numpy.ndarray of shape (3)
+    size : sequence of float, len (3)
         Physical size of the periodic field.
-    F : numpy.ndarray
+    F : numpy.ndarray, shape (:,:,:,3,3)
         Deformation gradient field.
 
     Returns
     -------
-    u_p_fluct : numpy.ndarray
+    u_p_fluct : numpy.ndarray, shape (:,:,:,3)
         Fluctuating part of the cell center displacements.
 
     """
-    integrator = 0.5j*size/_np.pi
+    integrator = 0.5j*_np.array(size,float)/_np.pi
 
     k_s = _ks(size,F.shape[:3],False)
     k_s_squared = _np.einsum('...l,...l',k_s,k_s)
@@ -192,20 +195,20 @@ def displacement_fluct_point(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
     return _np.fft.irfftn(displacement,axes=(0,1,2),s=F.shape[:3])
 
 
-def displacement_avg_point(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
+def displacement_avg_point(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
     """
     Cell center displacement field from average part of the deformation gradient field.
 
     Parameters
     ----------
-    size : numpy.ndarray of shape (3)
+    size : sequence of float, len (3)
         Physical size of the periodic field.
-    F : numpy.ndarray
+    F : numpy.ndarray, shape (:,:,:,3,3)
         Deformation gradient field.
 
     Returns
     -------
-    u_p_avg : numpy.ndarray
+    u_p_avg : numpy.ndarray, shape (:,:,:,3)
         Average part of the cell center displacements.
 
     """
@@ -213,42 +216,42 @@ def displacement_avg_point(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
     return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),coordinates0_point(F.shape[:3],size))
 
 
-def displacement_point(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
+def displacement_point(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
     """
     Cell center displacement field from deformation gradient field.
 
     Parameters
     ----------
-    size : numpy.ndarray of shape (3)
+    size : sequence of float, len (3)
         Physical size of the periodic field.
-    F : numpy.ndarray
+    F : numpy.ndarray, shape (:,:,:,3,3)
         Deformation gradient field.
 
     Returns
     -------
-    u_p : numpy.ndarray
+    u_p : numpy.ndarray, shape (:,:,:,3)
         Cell center displacements.
 
     """
     return displacement_avg_point(size,F) + displacement_fluct_point(size,F)
 
 
-def coordinates_point(size: _np.ndarray, F: _np.ndarray, origin: _np.ndarray = _np.zeros(3)) -> _np.ndarray:
+def coordinates_point(size: _FloatSequence, F: _np.ndarray, origin: _FloatSequence = _np.zeros(3)) -> _np.ndarray:
     """
     Cell center positions.
 
     Parameters
     ----------
-    size : numpy.ndarray of shape (3)
+    size : sequence of float, len (3)
         Physical size of the periodic field.
-    F : numpy.ndarray
+    F : numpy.ndarray, shape (:,:,:,3,3)
         Deformation gradient field.
-    origin : numpy.ndarray of shape (3), optional
+    origin : sequence of float, len(3), optional
         Physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
 
     Returns
     -------
-    x_p : numpy.ndarray
+    x_p : numpy.ndarray, shape (:,:,:,3)
         Cell center coordinates.
 
     """
@@ -256,14 +259,14 @@ def coordinates_point(size: _np.ndarray, F: _np.ndarray, origin: _np.ndarray = _
 
 
 def cellsSizeOrigin_coordinates0_point(coordinates0: _np.ndarray,
-                                       ordered: bool = True) -> Tuple[_np.ndarray,_np.ndarray,_np.ndarray]:
+                                       ordered: bool = True) -> _Tuple[_np.ndarray,_np.ndarray,_np.ndarray]:
     """
     Return grid 'DNA', i.e. cells, size, and origin from 1D array of point positions.
 
     Parameters
     ----------
-    coordinates0 : numpy.ndarray of shape (:,3)
-        Undeformed cell coordinates.
+    coordinates0 : numpy.ndarray, shape (:,3)
+        Undeformed cell center coordinates.
     ordered : bool, optional
         Expect coordinates0 data to be ordered (x fast, z slow).
         Defaults to True.
@@ -277,7 +280,7 @@ def cellsSizeOrigin_coordinates0_point(coordinates0: _np.ndarray,
     coords    = [_np.unique(coordinates0[:,i]) for i in range(3)]
     mincorner = _np.array(list(map(min,coords)))
     maxcorner = _np.array(list(map(max,coords)))
-    cells     = _np.array(list(map(len,coords)),'i')
+    cells     = _np.array(list(map(len,coords)),int)
     size      = cells/_np.maximum(cells-1,1) * (maxcorner-mincorner)
     delta     = size/cells
     origin    = mincorner - delta*.5
@@ -305,24 +308,24 @@ def cellsSizeOrigin_coordinates0_point(coordinates0: _np.ndarray,
     return (cells,size,origin)
 
 
-def coordinates0_node(cells: Union[_np.ndarray, Sequence[int]],
-                      size: Union[_np.ndarray, Sequence[int]],
-                      origin: Union[_np.ndarray, Sequence[int]] = _np.zeros(3)) -> _np.ndarray:
+def coordinates0_node(cells: _IntSequence,
+                      size: _FloatSequence,
+                      origin: _FloatSequence = _np.zeros(3)) -> _np.ndarray:
     """
     Nodal positions (undeformed).
 
     Parameters
     ----------
-    cells : numpy.ndarray of shape (3)
+    cells : sequence of int, len (3)
         Number of cells.
-    size : numpy.ndarray of shape (3)
+    size : sequence of float, len (3)
         Physical size of the periodic field.
-    origin : numpy.ndarray of shape (3), optional
+    origin : sequence of float, len(3), optional
         Physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
 
     Returns
     -------
-    x_n_0 : numpy.ndarray
+    x_n_0 : numpy.ndarray, shape (:,:,:,3)
         Undeformed nodal coordinates.
 
     """
@@ -332,40 +335,40 @@ def coordinates0_node(cells: Union[_np.ndarray, Sequence[int]],
                      axis = -1)
 
 
-def displacement_fluct_node(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
+def displacement_fluct_node(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
     """
     Nodal displacement field from fluctuation part of the deformation gradient field.
 
     Parameters
     ----------
-    size : numpy.ndarray of shape (3)
+    size : sequence of float, len (3)
         Physical size of the periodic field.
-    F : numpy.ndarray
+    F : numpy.ndarray, shape (:,:,:,3,3)
         Deformation gradient field.
 
     Returns
     -------
-    u_n_fluct : numpy.ndarray
+    u_n_fluct : numpy.ndarray, shape (:,:,:,3)
         Fluctuating part of the nodal displacements.
 
     """
     return point_to_node(displacement_fluct_point(size,F))
 
 
-def displacement_avg_node(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
+def displacement_avg_node(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
     """
     Nodal displacement field from average part of the deformation gradient field.
 
     Parameters
     ----------
-    size : numpy.ndarray of shape (3)
+    size : sequence of float, len (3)
         Physical size of the periodic field.
-    F : numpy.ndarray
+    F : numpy.ndarray, shape (:,:,:,3,3)
         Deformation gradient field.
 
     Returns
     -------
-    u_n_avg : numpy.ndarray
+    u_n_avg : numpy.ndarray, shape (:,:,:,3)
         Average part of the nodal displacements.
 
     """
@@ -373,42 +376,42 @@ def displacement_avg_node(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
     return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),coordinates0_node(F.shape[:3],size))
 
 
-def displacement_node(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray:
+def displacement_node(size: _FloatSequence, F: _np.ndarray) -> _np.ndarray:
     """
     Nodal displacement field from deformation gradient field.
 
     Parameters
     ----------
-    size : numpy.ndarray of shape (3)
+    size : sequence of float, len (3)
         Physical size of the periodic field.
-    F : numpy.ndarray
+    F : numpy.ndarray, shape (:,:,:,3,3)
         Deformation gradient field.
 
     Returns
     -------
-    u_p : numpy.ndarray
+    u_p : numpy.ndarray, shape (:,:,:,3)
         Nodal displacements.
 
     """
     return displacement_avg_node(size,F) + displacement_fluct_node(size,F)
 
 
-def coordinates_node(size: _np.ndarray, F: _np.ndarray, origin: _np.ndarray = _np.zeros(3)) -> _np.ndarray:
+def coordinates_node(size: _FloatSequence, F: _np.ndarray, origin: _FloatSequence = _np.zeros(3)) -> _np.ndarray:
     """
     Nodal positions.
 
     Parameters
     ----------
-    size : numpy.ndarray of shape (3)
+    size : sequence of float, len (3)
         Physical size of the periodic field.
-    F : numpy.ndarray
+    F : numpy.ndarray, shape (:,:,:,3,3)
         Deformation gradient field.
-    origin : numpy.ndarray of shape (3), optional
+    origin : sequence of float, len(3), optional
         Physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
 
     Returns
     -------
-    x_n : numpy.ndarray
+    x_n : numpy.ndarray, shape (:,:,:,3)
         Nodal coordinates.
 
     """
@@ -416,13 +419,13 @@ def coordinates_node(size: _np.ndarray, F: _np.ndarray, origin: _np.ndarray = _n
 
 
 def cellsSizeOrigin_coordinates0_node(coordinates0: _np.ndarray,
-                                      ordered: bool = True) -> Tuple[_np.ndarray,_np.ndarray,_np.ndarray]:
+                                      ordered: bool = True) -> _Tuple[_np.ndarray,_np.ndarray,_np.ndarray]:
     """
     Return grid 'DNA', i.e. cells, size, and origin from 1D array of nodal positions.
 
     Parameters
     ----------
-    coordinates0 : numpy.ndarray of shape (:,3)
+    coordinates0 : numpy.ndarray, shape (:,3)
         Undeformed nodal coordinates.
     ordered : bool, optional
         Expect coordinates0 data to be ordered (x fast, z slow).
@@ -437,7 +440,7 @@ def cellsSizeOrigin_coordinates0_node(coordinates0: _np.ndarray,
     coords    = [_np.unique(coordinates0[:,i]) for i in range(3)]
     mincorner = _np.array(list(map(min,coords)))
     maxcorner = _np.array(list(map(max,coords)))
-    cells     = _np.array(list(map(len,coords)),'i') - 1
+    cells     = _np.array(list(map(len,coords)),int) - 1
     size      = maxcorner-mincorner
     origin    = mincorner
 
@@ -463,12 +466,12 @@ def point_to_node(cell_data: _np.ndarray) -> _np.ndarray:
 
     Parameters
     ----------
-    cell_data : numpy.ndarray of shape (:,:,:,...)
+    cell_data : numpy.ndarray, shape (:,:,:,...)
         Data defined on the cell centers of a periodic grid.
 
     Returns
     -------
-    node_data : numpy.ndarray of shape (:,:,:,...)
+    node_data : numpy.ndarray, shape (:,:,:,...)
         Data defined on the nodes of a periodic grid.
 
     """
@@ -485,12 +488,12 @@ def node_to_point(node_data: _np.ndarray) -> _np.ndarray:
 
     Parameters
     ----------
-    node_data : numpy.ndarray of shape (:,:,:,...)
+    node_data : numpy.ndarray, shape (:,:,:,...)
         Data defined on the nodes of a periodic grid.
 
     Returns
     -------
-    cell_data : numpy.ndarray of shape (:,:,:,...)
+    cell_data : numpy.ndarray, shape (:,:,:,...)
         Data defined on the cell centers of a periodic grid.
 
     """
@@ -507,7 +510,7 @@ def coordinates0_valid(coordinates0: _np.ndarray) -> bool:
 
     Parameters
     ----------
-    coordinates0 : numpy.ndarray
+    coordinates0 : numpy.ndarray, shape (:,3)
         Array of undeformed cell coordinates.
 
     Returns
@@ -523,17 +526,17 @@ def coordinates0_valid(coordinates0: _np.ndarray) -> bool:
         return False
 
 
-def regrid(size: _np.ndarray, F: _np.ndarray, cells: Union[_np.ndarray,Sequence[int]]) -> _np.ndarray:
+def regrid(size: _FloatSequence, F: _np.ndarray, cells: _IntSequence) -> _np.ndarray:
     """
     Return mapping from coordinates in deformed configuration to a regular grid.
 
     Parameters
     ----------
-    size : numpy.ndarray of shape (3)
+    size : sequence of float, len (3)
         Physical size.
-    F : numpy.ndarray of shape (:,:,:,3,3)
+    F : numpy.ndarray, shape (:,:,:,3,3), shape (:,:,:,3,3)
         Deformation gradient field.
-    cells : numpy.ndarray of shape (3)
+    cells : sequence of int, len (3)
         Cell count along x,y,z of remapping grid.
 
     """
diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py
index 22e3aeabf..7c1af6c5f 100644
--- a/python/damask/mechanics.py
+++ b/python/damask/mechanics.py
@@ -5,7 +5,7 @@ All routines operate on numpy.ndarrays of shape (...,3,3).
 
 """
 
-from typing import Sequence
+from typing import Sequence as _Sequence
 
 import numpy as _np
 
@@ -243,7 +243,7 @@ def stretch_right(T: _np.ndarray) -> _np.ndarray:
     return _polar_decomposition(T,'U')[0]
 
 
-def _polar_decomposition(T: _np.ndarray, requested: Sequence[str]) -> tuple:
+def _polar_decomposition(T: _np.ndarray, requested: _Sequence[str]) -> tuple:
     """
     Perform singular value decomposition.
 
diff --git a/python/damask/seeds.py b/python/damask/seeds.py
index 4d5a8c624..7e01a40e6 100644
--- a/python/damask/seeds.py
+++ b/python/damask/seeds.py
@@ -1,25 +1,27 @@
 """Functionality for generation of seed points for Voronoi or Laguerre tessellation."""
 
-from typing import Sequence,Tuple
+from typing import Tuple as _Tuple
 
 from scipy import spatial as _spatial
 import numpy as _np
 
+from ._typehints import FloatSequence as _FloatSequence, IntSequence as _IntSequence
 from . import util as _util
 from . import grid_filters as _grid_filters
 
 
-def from_random(size: _np.ndarray, N_seeds: int, cells: _np.ndarray = None, rng_seed=None) -> _np.ndarray:
+def from_random(size: _FloatSequence, N_seeds: int, cells: _IntSequence = None,
+                rng_seed=None) -> _np.ndarray:
     """
     Place seeds randomly in space.
 
     Parameters
     ----------
-    size : numpy.ndarray of shape (3)
+    size : sequence of float, len (3)
         Physical size of the seeding domain.
     N_seeds : int
         Number of seeds.
-    cells : numpy.ndarray of shape (3), optional.
+    cells : sequence of int, len (3), optional.
         If given, ensures that each seed results in a grain when a standard Voronoi
         tessellation is performed using the given grid resolution (i.e. size/cells).
     rng_seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
@@ -28,29 +30,30 @@ def from_random(size: _np.ndarray, N_seeds: int, cells: _np.ndarray = None, rng_
 
     Returns
     -------
-    coords : numpy.ndarray of shape (N_seeds,3)
+    coords : numpy.ndarray, shape (N_seeds,3)
         Seed coordinates in 3D space.
 
     """
+    size_ = _np.array(size,float)
     rng = _np.random.default_rng(rng_seed)
     if cells is None:
-        coords = rng.random((N_seeds,3)) * size
+        coords = rng.random((N_seeds,3)) * size_
     else:
         grid_coords = _grid_filters.coordinates0_point(cells,size).reshape(-1,3,order='F')
         coords = grid_coords[rng.choice(_np.prod(cells),N_seeds, replace=False)] \
-               + _np.broadcast_to(size/cells,(N_seeds,3))*(rng.random((N_seeds,3))*.5-.25)          # wobble without leaving cells
+               + _np.broadcast_to(size_/_np.array(cells,int),(N_seeds,3))*(rng.random((N_seeds,3))*.5-.25) # wobble w/o leaving grid
 
     return coords
 
 
-def from_Poisson_disc(size: _np.ndarray, N_seeds: int, N_candidates: int, distance: float,
+def from_Poisson_disc(size: _FloatSequence, N_seeds: int, N_candidates: int, distance: float,
                       periodic: bool = True, rng_seed=None) -> _np.ndarray:
     """
     Place seeds according to a Poisson disc distribution.
 
     Parameters
     ----------
-    size : numpy.ndarray of shape (3)
+    size : sequence of float, len (3)
         Physical size of the seeding domain.
     N_seeds : int
         Number of seeds.
@@ -66,13 +69,13 @@ def from_Poisson_disc(size: _np.ndarray, N_seeds: int, N_candidates: int, distan
 
     Returns
     -------
-    coords : numpy.ndarray of shape (N_seeds,3)
+    coords : numpy.ndarray, shape (N_seeds,3)
         Seed coordinates in 3D space.
 
     """
     rng = _np.random.default_rng(rng_seed)
     coords = _np.empty((N_seeds,3))
-    coords[0] = rng.random(3) * size
+    coords[0] = rng.random(3) * _np.array(size,float)
 
     s = 1
     i = 0
@@ -96,8 +99,8 @@ def from_Poisson_disc(size: _np.ndarray, N_seeds: int, N_candidates: int, distan
     return coords
 
 
-def from_grid(grid, selection: Sequence[int] = None,
-              invert: bool = False, average: bool = False, periodic: bool = True) -> Tuple[_np.ndarray, _np.ndarray]:
+def from_grid(grid, selection: _IntSequence = None,
+              invert: bool = False, average: bool = False, periodic: bool = True) -> _Tuple[_np.ndarray, _np.ndarray]:
     """
     Create seeds from grid description.
 
@@ -105,7 +108,7 @@ def from_grid(grid, selection: Sequence[int] = None,
     ----------
     grid : damask.Grid
         Grid from which the material IDs are used as seeds.
-    selection : iterable of integers, optional
+    selection : sequence of int, optional
         Material IDs to consider.
     invert : boolean, false
         Consider all material IDs except those in selection. Defaults to False.
@@ -116,7 +119,7 @@ def from_grid(grid, selection: Sequence[int] = None,
 
     Returns
     -------
-    coords, materials : numpy.ndarray of shape (:,3), numpy.ndarray of shape (:)
+    coords, materials : numpy.ndarray, shape (:,3); numpy.ndarray, shape (:)
         Seed coordinates in 3D space, material IDs.
 
     """
diff --git a/python/tests/test_Crystal.py b/python/tests/test_Crystal.py
index 90fbdce74..9c6de965d 100644
--- a/python/tests/test_Crystal.py
+++ b/python/tests/test_Crystal.py
@@ -79,3 +79,23 @@ class TestCrystal:
                     a=a,b=b,c=c,
                     alpha=alpha,beta=beta,gamma=gamma)
         assert np.allclose(points,c.lattice_points)
+
+    @pytest.mark.parametrize('crystal,length',
+                            [(Crystal(lattice='cF'),[12,6]),
+                             (Crystal(lattice='cI'),[12,12,24]),
+                             (Crystal(lattice='hP'),[3,3,6,12,6]),
+                             (Crystal(lattice='tI',c=1.2),[2,2,2,4,2,4,2,2,4,8,4,8,8])
+                            ])
+    def test_N_slip(self,crystal,length):
+        assert [len(s) for s in crystal.kinematics('slip')['direction']] == length
+        assert [len(s) for s in crystal.kinematics('slip')['plane']] == length
+
+    @pytest.mark.parametrize('crystal,length',
+                            [(Crystal(lattice='cF'),[12]),
+                             (Crystal(lattice='cI'),[12]),
+                             (Crystal(lattice='hP'),[6,6,6,6]),
+                            ])
+    def test_N_twin(self,crystal,length):
+        assert [len(s) for s in crystal.kinematics('twin')['direction']] == length
+        assert [len(s) for s in crystal.kinematics('twin')['plane']] == length
+
diff --git a/python/tests/test_grid_filters.py b/python/tests/test_grid_filters.py
index d9c074c11..d5458f0eb 100644
--- a/python/tests/test_grid_filters.py
+++ b/python/tests/test_grid_filters.py
@@ -2,6 +2,8 @@ import pytest
 import numpy as np
 
 from damask import grid_filters
+from damask import Grid
+from damask import seeds
 
 class TestGridFilters:
 
@@ -139,12 +141,19 @@ class TestGridFilters:
         else:
             function(unordered,mode)
 
-    def test_regrid(self):
+    def test_regrid_identity(self):
          size = np.random.random(3)
          cells = np.random.randint(8,32,(3))
-         F    = np.broadcast_to(np.eye(3), tuple(cells)+(3,3))
+         F = np.broadcast_to(np.eye(3), tuple(cells)+(3,3))
          assert all(grid_filters.regrid(size,F,cells) == np.arange(cells.prod()))
 
+    def test_regrid_double_cells(self):
+         size = np.random.random(3)
+         cells = np.random.randint(8,32,(3))
+         g = Grid.from_Voronoi_tessellation(cells,size,seeds.from_random(size,10))
+         F = np.broadcast_to(np.eye(3), tuple(cells)+(3,3))
+         assert all(g.scale(cells*2).material.flatten() ==
+                    g.material.flatten()[grid_filters.regrid(size,F,cells*2)])
 
     @pytest.mark.parametrize('differential_operator',[grid_filters.curl,
                                                       grid_filters.divergence,
diff --git a/src/constants.f90 b/src/constants.f90
index 43bb60953..3b99763de 100644
--- a/src/constants.f90
+++ b/src/constants.f90
@@ -11,6 +11,6 @@ module constants
   real(pReal), parameter :: &
     T_ROOM = 300.0_pReal, &                                                                         !< Room temperature in K. ToDo: IUPAC: 298.15
     K_B = 1.38e-23_pReal, &                                                                         !< Boltzmann constant in J/Kelvin
-    N_A = 6.02214076e-23_pReal                                                                      !< Avogadro constant in 1/mol
+    N_A = 6.02214076e23_pReal                                                                      !< Avogadro constant in 1/mol
 
 end module constants
diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90
index 24797a880..bb530747b 100644
--- a/src/phase_mechanical_elastic.f90
+++ b/src/phase_mechanical_elastic.f90
@@ -30,7 +30,7 @@ module subroutine elastic_init(phases)
     phase, &
     mech, &
     elastic
-  logical :: thermal_active
+
 
   print'(/,1x,a)', '<<<+-  phase:mechanical:elastic init  -+>>>'
   print'(/,1x,a)', '<<<+-  phase:mechanical:elastic:Hooke init  -+>>>'
diff --git a/src/rotations.f90 b/src/rotations.f90
index 97e8104d5..b4728e38b 100644
--- a/src/rotations.f90
+++ b/src/rotations.f90
@@ -372,7 +372,7 @@ end function rotTensor4
 
 
 !---------------------------------------------------------------------------------------------------
-!> @brief Rotate a rank-4 tensor in Voigt 6x6 notation passively (default) or actively.
+!> @brief Rotate a rank-4 stiffness tensor in Voigt 6x6 notation passively (default) or actively.
 !> @details: https://scicomp.stackexchange.com/questions/35600
 !! ToDo: Need to check active/passive !!!
 !---------------------------------------------------------------------------------------------------
@@ -393,11 +393,11 @@ pure function rotStiffness(self,C,active) result(cRot)
     R = self%asMatrix()
   endif
 
-  M = reshape([R(1,1)**2.0_pReal,           R(2,1)**2.0_pReal,           R(3,1)**2.0_pReal, &
+  M = reshape([R(1,1)**2,                   R(2,1)**2,                   R(3,1)**2, &
                R(2,1)*R(3,1),               R(1,1)*R(3,1),               R(1,1)*R(2,1), &
-               R(1,2)**2.0_pReal,           R(2,2)**2.0_pReal,           R(3,2)**2.0_pReal, &
+               R(1,2)**2,                   R(2,2)**2,                   R(3,2)**2, &
                R(2,2)*R(3,2),               R(1,2)*R(3,2),               R(1,2)*R(2,2), &
-               R(1,3)**2.0_pReal,           R(2,3)**2.0_pReal,           R(3,3)**2.0_pReal, &
+               R(1,3)**2,                   R(2,3)**2,                   R(3,3)**2, &
                R(2,3)*R(3,3),               R(1,3)*R(3,3),               R(1,3)*R(2,3), &
                2.0_pReal*R(1,2)*R(1,3),     2.0_pReal*R(2,2)*R(2,3),     2.0_pReal*R(3,2)*R(3,3), &
                R(2,2)*R(3,3)+R(2,3)*R(3,2), R(1,2)*R(3,3)+R(1,3)*R(3,2), R(1,2)*R(2,3)+R(1,3)*R(2,2), &