diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 055a40815..001c326a4 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/PRIVATE b/PRIVATE
index 96c32ba42..b898a8b55 160000
--- a/PRIVATE
+++ b/PRIVATE
@@ -1 +1 @@
-Subproject commit 96c32ba4237a51eaad92cd139e1a716ee5b32493
+Subproject commit b898a8b5552bd9d1c555edc3d8134564dd32fe53
diff --git a/examples/config/Phase_Phenopowerlaw_BCC-Martensite.yaml b/examples/config/Phase_Phenopowerlaw_BCC-Martensite.yaml
index 3ad6779eb..f2066b373 100644
--- a/examples/config/Phase_Phenopowerlaw_BCC-Martensite.yaml
+++ b/examples/config/Phase_Phenopowerlaw_BCC-Martensite.yaml
@@ -1,17 +1,12 @@
 # Tasan et.al. 2015 Acta Materalia
 # Tasan et.al. 2015 International Journal of Plasticity
 # Diehl et.al. 2015 Meccanica
-Martensite:
-  lattice: cI
-  mechanical:
-    elastic: {C_11: 417.4e+9, C_12: 242.4e+9, C_44: 211.1e+9, type: Hooke}
-    plastic:
-      N_sl: [12, 12]
-      a_sl: 2.0
-      dot_gamma_0_sl: 0.001
-      h_0_sl-sl: 563.0e+9
-      h_sl-sl: [1, 1.4, 1, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4]
-      n_sl: 20
-      type: phenopowerlaw
-      xi_0_sl: [405.8e+6, 456.7e+6]
-      xi_inf_sl: [872.9e+6, 971.2e+6]
+N_sl: [12, 12]
+a_sl: 2.0
+dot_gamma_0_sl: 0.001
+h_0_sl-sl: 563.0e+9
+h_sl-sl: [1, 1.4, 1, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4]
+n_sl: 20
+type: phenopowerlaw
+xi_0_sl: [405.8e+6, 456.7e+6]
+xi_inf_sl: [872.9e+6, 971.2e+6]
diff --git a/examples/config/phase/AISI304.yaml b/examples/config/phase/AISI304.yaml
new file mode 100644
index 000000000..0ca55635b
--- /dev/null
+++ b/examples/config/phase/AISI304.yaml
@@ -0,0 +1,6 @@
+references:
+  - H.M. Ledbetter
+    physica status solidi (a) 85(1):89-96, 1984
+    https://doi.org/10.1002/pssa.2210850111
+lattice: cF
+rho: 7937.0
diff --git a/examples/config/phase/mechanical/eigen/thermalexpansion_AISI304.yaml b/examples/config/phase/mechanical/eigen/thermalexpansion_AISI304.yaml
new file mode 100644
index 000000000..6a3b339bb
--- /dev/null
+++ b/examples/config/phase/mechanical/eigen/thermalexpansion_AISI304.yaml
@@ -0,0 +1,7 @@
+type: thermalexpansion
+references:
+  - R.H. Bogaard et al.
+    Thermochimica Acta 218:373-393, 1993
+    https://doi.org/10.1016/0040-6031(93)80437-F
+A_11: 15.0e-6
+T_ref: 300.0
diff --git a/examples/config/phase/mechanical/elastic/Hooke_AISI304.yaml b/examples/config/phase/mechanical/elastic/Hooke_AISI304.yaml
new file mode 100644
index 000000000..70fbd25c7
--- /dev/null
+++ b/examples/config/phase/mechanical/elastic/Hooke_AISI304.yaml
@@ -0,0 +1,8 @@
+type: Hooke
+references:
+  - H.M. Ledbetter
+    physica status solidi (a) 85(1):89-96, 1984
+    https://doi.org/10.1002/pssa.2210850111
+C_11: 204.6e+9
+C_12: 137.7e+9
+C_44: 126.2e+9
diff --git a/examples/config/phase/mechanical/elastic/Hooke_SAE1050-martensite.yaml b/examples/config/phase/mechanical/elastic/Hooke_SAE1050-martensite.yaml
new file mode 100644
index 000000000..7ba941066
--- /dev/null
+++ b/examples/config/phase/mechanical/elastic/Hooke_SAE1050-martensite.yaml
@@ -0,0 +1,8 @@
+type: Hooke
+references:
+  - S.A. Kim and W.L. Johnson,
+    Materials Science & Engineering A 452-453:633-639, 2007,
+    https://doi.org/10.1016/j.msea.2006.11.147
+C_11: 268.1e+9
+C_12: 111.2e+9
+C_44: 79.06e+9
diff --git a/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml b/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml
index 362797d0b..4e4ff1a9f 100644
--- a/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml
+++ b/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml
@@ -4,7 +4,8 @@ references:
     International Journal of Plasticity 134:102779, 2020,
     https://doi.org/10.1016/j.ijplas.2020.102779
   - K. Sedighiani et al.,
-    Mechanics of Materials, submitted
+    Mechanics of Materials, 164:104117, 2022,
+    https://doi.org/10.1016/j.mechmat.2021.104117
 output: [rho_dip, rho_mob]
 N_sl: [12, 12]
 b_sl: [2.49e-10, 2.49e-10]
diff --git a/examples/config/phase/thermal/AISI304.yaml b/examples/config/phase/thermal/AISI304.yaml
new file mode 100644
index 000000000..ef09fd850
--- /dev/null
+++ b/examples/config/phase/thermal/AISI304.yaml
@@ -0,0 +1,9 @@
+references:
+  - B.F. Blackwell et al.
+    Proceedings of 34th National Heat Transfer Conference 2000
+    https://www.osti.gov/servlets/purl/760791
+  - R.H. Bogaard et al.
+    Thermochimica Acta 218:373-393, 1993
+    https://doi.org/10.1016/0040-6031(93)80437-F
+C_p: 470.0
+K_11: 14.34
diff --git a/examples/config/phase/thermal/Steel-0.5C.yaml b/examples/config/phase/thermal/steel-0.5C.yaml
similarity index 100%
rename from examples/config/phase/thermal/Steel-0.5C.yaml
rename to examples/config/phase/thermal/steel-0.5C.yaml
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 e657e6a62..d9749c5a9 100644
--- a/python/damask/VERSION
+++ b/python/damask/VERSION
@@ -1 +1 @@
-v3.0.0-alpha5-333-g01cd92755
+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 ba7617fbb..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,14 +226,14 @@ 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.
 
         Parameters
         ----------
-        field : numpy.array, shape (:,:)
+        field : numpy.ndarray, shape (:,:)
             Data to be shaded.
         bounds : sequence of float, len (2), optional
             Value range (left,right) spanned by 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.
 
@@ -616,7 +613,7 @@ class Colormap(mpl.colors.ListedColormap):
 
 
     @staticmethod
-    def _lab2xyz(lab: np.ndarray, ref_white: np.ndarray = None) -> np.ndarray:
+    def _lab2xyz(lab: np.ndarray, ref_white: np.ndarray = _REF_WHITE) -> np.ndarray:
         """
         CIE Lab to CIE Xyz.
 
@@ -624,6 +621,8 @@ class Colormap(mpl.colors.ListedColormap):
         ----------
         lab : numpy.ndarray, shape (3)
             CIE lab values.
+        ref_white : numpy.ndarray, shape (3)
+            Reference white, default value is the standard 2° observer for D65.
 
         Returns
         -------
@@ -642,10 +641,10 @@ class Colormap(mpl.colors.ListedColormap):
                          f_x**3.                if f_x**3. > _EPS     else (116.*f_x-16.)/_KAPPA,
                          ((lab[0]+16.)/116.)**3 if lab[0]>_KAPPA*_EPS else lab[0]/_KAPPA,
                          f_z**3.                if f_z**3. > _EPS     else (116.*f_z-16.)/_KAPPA
-                        ])*(ref_white if ref_white is not None else _REF_WHITE)
+                        ])*ref_white
 
     @staticmethod
-    def _xyz2lab(xyz: np.ndarray, ref_white: np.ndarray = None) -> np.ndarray:
+    def _xyz2lab(xyz: np.ndarray, ref_white: np.ndarray = _REF_WHITE) -> np.ndarray:
         """
         CIE Xyz to CIE Lab.
 
@@ -653,6 +652,8 @@ class Colormap(mpl.colors.ListedColormap):
         ----------
         xyz : numpy.ndarray, shape (3)
             CIE Xyz values.
+        ref_white : numpy.ndarray, shape (3)
+            Reference white, default value is the standard 2° observer for D65.
 
         Returns
         -------
@@ -664,7 +665,6 @@ class Colormap(mpl.colors.ListedColormap):
         http://www.brucelindbloom.com/index.html?Eqn_Lab_to_XYZ.html
 
         """
-        ref_white = ref_white if ref_white is not None else _REF_WHITE
         f = np.where(xyz/ref_white > _EPS,(xyz/ref_white)**(1./3.),(_KAPPA*xyz/ref_white+16.)/116.)
 
         return np.array([
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 42b5a16c4..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: _np.ndarray,
-                       origin: _np.ndarray = _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 = origin        + size/_np.array(cells)*.5
-    end   = origin + 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: _np.ndarray,
-                      origin: _np.ndarray = _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/damask/util.py b/python/damask/util.py
index a39a7b8c7..0581302db 100644
--- a/python/damask/util.py
+++ b/python/damask/util.py
@@ -18,11 +18,11 @@ from . import version
 __all__=[
          'srepr',
          'emph', 'deemph', 'warn', 'strikeout',
-         'run', 
+         'run',
          'natural_sort',
          'show_progress',
          'scale_to_coprime',
-         'project_stereographic',
+         'project_equal_angle', 'project_equal_area',
          'hybrid_IA',
          'execution_stamp',
          'shapeshifter', 'shapeblender',
@@ -267,13 +267,13 @@ def scale_to_coprime(v):
     return m
 
 
-def project_stereographic(vector,direction='z',normalize=True,keepdims=False):
+def project_equal_angle(vector,direction='z',normalize=True,keepdims=False):
     """
-    Apply stereographic projection to vector.
+    Apply equal-angle projection to vector.
 
     Parameters
     ----------
-    vector : numpy.ndarray of shape (...,3)
+    vector : numpy.ndarray, shape (...,3)
         Vector coordinates to be projected.
     direction : str
         Projection direction 'x', 'y', or 'z'.
@@ -281,32 +281,74 @@ def project_stereographic(vector,direction='z',normalize=True,keepdims=False):
     normalize : bool
         Ensure unit length of input vector. Defaults to True.
     keepdims : bool
-        Maintain three-dimensional output coordinates.
-        Default two-dimensional output uses right-handed frame spanned by
+        Maintain three-dimensional output coordinates. Defaults to False.
+        Two-dimensional output uses right-handed frame spanned by
         the next and next-next axis relative to the projection direction,
         e.g. x-y when projecting along z and z-x when projecting along y.
 
     Returns
     -------
-    coordinates : numpy.ndarray of shape (...,2 | 3)
+    coordinates : numpy.ndarray, shape (...,2 | 3)
         Projected coordinates.
 
     Examples
     --------
     >>> import damask
     >>> import numpy as np
-    >>> project_stereographic(np.ones(3))
+    >>> project_equal_angle(np.ones(3))
         [0.3660254, 0.3660254]
-    >>> project_stereographic(np.ones(3),direction='x',normalize=False,keepdims=True)
+    >>> project_equal_angle(np.ones(3),direction='x',normalize=False,keepdims=True)
         [0, 0.5, 0.5]
-    >>> project_stereographic([0,1,1],direction='y',normalize=True,keepdims=False)
+    >>> project_equal_angle([0,1,1],direction='y',normalize=True,keepdims=False)
         [0.41421356, 0]
 
     """
     shift = 'zyx'.index(direction)
-    v_ = np.roll(vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector,
-                 shift,axis=-1)
-    return np.roll(np.block([v_[...,:2]/(1+np.abs(v_[...,2:3])),np.zeros_like(v_[...,2:3])]),
+    v = np.roll(vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector,
+                shift,axis=-1)
+    return np.roll(np.block([v[...,:2]/(1.0+np.abs(v[...,2:3])),np.zeros_like(v[...,2:3])]),
+                   -shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2]
+
+def project_equal_area(vector,direction='z',normalize=True,keepdims=False):
+    """
+    Apply equal-area projection to vector.
+
+    Parameters
+    ----------
+    vector : numpy.ndarray, shape (...,3)
+        Vector coordinates to be projected.
+    direction : str
+        Projection direction 'x', 'y', or 'z'.
+        Defaults to 'z'.
+    normalize : bool
+        Ensure unit length of input vector. Defaults to True.
+    keepdims : bool
+        Maintain three-dimensional output coordinates. Defaults to False.
+        Two-dimensional output uses right-handed frame spanned by
+        the next and next-next axis relative to the projection direction,
+        e.g. x-y when projecting along z and z-x when projecting along y.
+
+    Returns
+    -------
+    coordinates : numpy.ndarray, shape (...,2 | 3)
+        Projected coordinates.
+
+    Examples
+    --------
+    >>> import damask
+    >>> import numpy as np
+    >>> project_equal_area(np.ones(3))
+        [0.45970084, 0.45970084]
+    >>> project_equal_area(np.ones(3),direction='x',normalize=False,keepdims=True)
+        [0.0, 0.70710678, 0.70710678]
+    >>> project_equal_area([0,1,1],direction='y',normalize=True,keepdims=False)
+        [0.5411961, 0.0]
+
+    """
+    shift = 'zyx'.index(direction)
+    v = np.roll(vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector,
+                shift,axis=-1)
+    return np.roll(np.block([v[...,:2]/np.sqrt(1.0+np.abs(v[...,2:3])),np.zeros_like(v[...,2:3])]),
                    -shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2]
 
 
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/python/tests/test_util.py b/python/tests/test_util.py
index 13731eda2..ee345605f 100644
--- a/python/tests/test_util.py
+++ b/python/tests/test_util.py
@@ -59,9 +59,22 @@ class TestUtil:
                               ([1,1,0],'x',False,False,[0.5,0]),
                               ([1,1,1],'y',True, True, [0.3660254, 0,0.3660254]),
                              ])
-    def test_project_stereographic(self,point,direction,normalize,keepdims,answer):
-        assert np.allclose(util.project_stereographic(np.array(point),direction=direction,
-                                                      normalize=normalize,keepdims=keepdims),answer)
+    def test_project_equal_angle(self,point,direction,normalize,keepdims,answer):
+        assert np.allclose(util.project_equal_angle(np.array(point),direction=direction,
+                                                    normalize=normalize,keepdims=keepdims),answer)
+
+    @pytest.mark.parametrize('point,direction,normalize,keepdims,answer',
+                             [
+                              ([1,0,0],'z',False,True, [1,0,0]),
+                              ([1,0,0],'z',True, False,[1,0]),
+                              ([0,1,1],'z',False,True, [0,0.70710678,0]),
+                              ([0,1,1],'y',True, False,[0.5411961,0]),
+                              ([1,1,0],'x',False,False,[0.70710678,0]),
+                              ([1,1,1],'y',True, True, [0.45970084,0,0.45970084]),
+                             ])
+    def test_project_equal_area(self,point,direction,normalize,keepdims,answer):
+        assert np.allclose(util.project_equal_area(np.array(point),direction=direction,
+                                                   normalize=normalize,keepdims=keepdims),answer)
 
     @pytest.mark.parametrize('fro,to,mode,answer',
                              [
diff --git a/src/LAPACK_interface.f90 b/src/LAPACK_interface.f90
index 7d3043ed0..e11f74875 100644
--- a/src/LAPACK_interface.f90
+++ b/src/LAPACK_interface.f90
@@ -6,7 +6,7 @@
 module LAPACK_interface
   interface
 
-    subroutine dgeev(jobvl,jobvr,n,a,lda,wr,wi,vl,ldvl,vr,ldvr,work,lwork,info)
+    pure subroutine dgeev(jobvl,jobvr,n,a,lda,wr,wi,vl,ldvl,vr,ldvr,work,lwork,info)
       use prec
       character,   intent(in)                             :: jobvl,jobvr
       integer,     intent(in)                             :: n,lda,ldvl,ldvr,lwork
@@ -18,16 +18,16 @@ module LAPACK_interface
       integer,     intent(out)                            :: info
     end subroutine dgeev
 
-    subroutine dgesv(n,nrhs,a,lda,ipiv,b,ldb,info)
+    pure subroutine dgesv(n,nrhs,a,lda,ipiv,b,ldb,info)
       use prec
       integer,     intent(in)                             :: n,nrhs,lda,ldb
       real(pReal), intent(inout), dimension(lda,n)        :: a
       integer,     intent(out),   dimension(n)            :: ipiv
-      real(pReal), intent(out),   dimension(ldb,nrhs)     :: b
+      real(pReal), intent(inout), dimension(ldb,nrhs)     :: b
       integer,     intent(out)                            :: info
     end subroutine dgesv
 
-    subroutine dgetrf(m,n,a,lda,ipiv,info)
+    pure subroutine dgetrf(m,n,a,lda,ipiv,info)
       use prec
       integer,     intent(in)                             :: m,n,lda
       real(pReal), intent(inout), dimension(lda,n)        :: a
@@ -35,16 +35,16 @@ module LAPACK_interface
       integer,     intent(out)                            :: info
     end subroutine dgetrf
 
-    subroutine dgetri(n,a,lda,ipiv,work,lwork,info)
+    pure subroutine dgetri(n,a,lda,ipiv,work,lwork,info)
       use prec
       integer,     intent(in)                             :: n,lda,lwork
       real(pReal), intent(inout), dimension(lda,n)        :: a
-      integer,     intent(out),   dimension(n)            :: ipiv
+      integer,     intent(in),    dimension(n)            :: ipiv
       real(pReal), intent(out),   dimension(max(1,lwork)) :: work
       integer,     intent(out)                            :: info
     end subroutine dgetri
 
-    subroutine dsyev(jobz,uplo,n,a,lda,w,work,lwork,info)
+    pure subroutine dsyev(jobz,uplo,n,a,lda,w,work,lwork,info)
       use prec
       character,   intent(in)                             :: jobz,uplo
       integer,     intent(in)                             :: n,lda,lwork
diff --git a/src/constants.f90 b/src/constants.f90
index dd26ce78c..3b99763de 100644
--- a/src/constants.f90
+++ b/src/constants.f90
@@ -9,7 +9,8 @@ module constants
   public
 
   real(pReal), parameter :: &
-    T_ROOM = 300.0_pReal, &                                                                         !< Room temperature in K
-    K_B = 1.38e-23_pReal                                                                             !< Boltzmann constant in J/Kelvin
+    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.02214076e23_pReal                                                                      !< Avogadro constant in 1/mol
 
 end module constants
diff --git a/src/lattice.f90 b/src/lattice.f90
index 50d4b6c51..3089aa30b 100644
--- a/src/lattice.f90
+++ b/src/lattice.f90
@@ -2070,7 +2070,7 @@ end function getlabels
 !> @brief Equivalent Poisson's ratio (ν)
 !> @details https://doi.org/10.1143/JPSJ.20.635
 !--------------------------------------------------------------------------------------------------
-function lattice_equivalent_nu(C,assumption) result(nu)
+pure function lattice_equivalent_nu(C,assumption) result(nu)
 
   real(pReal), dimension(6,6), intent(in) :: C                                                      !< Stiffness tensor (Voigt notation)
   character(len=5),            intent(in) :: assumption                                             !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
@@ -2103,7 +2103,7 @@ end function lattice_equivalent_nu
 !> @brief Equivalent shear modulus (μ)
 !> @details https://doi.org/10.1143/JPSJ.20.635
 !--------------------------------------------------------------------------------------------------
-function lattice_equivalent_mu(C,assumption) result(mu)
+pure function lattice_equivalent_mu(C,assumption) result(mu)
 
   real(pReal), dimension(6,6), intent(in) :: C                                                      !< Stiffness tensor (Voigt notation)
   character(len=5),            intent(in) :: assumption                                             !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
diff --git a/src/math.f90 b/src/math.f90
index 2d8f8fb3b..3ea4d6a08 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -512,7 +512,7 @@ end subroutine math_invert33
 !--------------------------------------------------------------------------------------------------
 !> @brief Inversion of symmetriced 3x3x3x3 matrix
 !--------------------------------------------------------------------------------------------------
-function math_invSym3333(A)
+pure function math_invSym3333(A)
 
   real(pReal),dimension(3,3,3,3)            :: math_invSym3333
 
@@ -538,7 +538,7 @@ end function math_invSym3333
 !--------------------------------------------------------------------------------------------------
 !> @brief invert quadratic matrix of arbitrary dimension
 !--------------------------------------------------------------------------------------------------
-subroutine math_invert(InvA, error, A)
+pure subroutine math_invert(InvA, error, A)
 
   real(pReal), dimension(:,:),                 intent(in)  :: A
   real(pReal), dimension(size(A,1),size(A,1)), intent(out) :: invA
@@ -996,7 +996,7 @@ end subroutine math_normal
 !--------------------------------------------------------------------------------------------------
 !> @brief eigenvalues and eigenvectors of symmetric matrix
 !--------------------------------------------------------------------------------------------------
-subroutine math_eigh(w,v,error,m)
+pure subroutine math_eigh(w,v,error,m)
 
   real(pReal), dimension(:,:),                  intent(in)  :: m                                    !< quadratic matrix to compute eigenvectors and values of
   real(pReal), dimension(size(m,1)),            intent(out) :: w                                    !< eigenvalues
@@ -1021,7 +1021,7 @@ end subroutine math_eigh
 !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
 !> @details See http://arxiv.org/abs/physics/0610206 (DSYEVH3)
 !--------------------------------------------------------------------------------------------------
-subroutine math_eigh33(w,v,m)
+pure subroutine math_eigh33(w,v,m)
 
   real(pReal), dimension(3,3),intent(in)  :: m                                                      !< 3x3 matrix to compute eigenvectors and values of
   real(pReal), dimension(3),  intent(out) :: w                                                      !< eigenvalues
@@ -1114,7 +1114,7 @@ end function math_rotationalPart
 !> @brief Eigenvalues of symmetric matrix
 ! will return NaN on error
 !--------------------------------------------------------------------------------------------------
-function math_eigvalsh(m)
+pure function math_eigvalsh(m)
 
   real(pReal), dimension(:,:),                  intent(in)  :: m                                    !< symmetric matrix to compute eigenvalues of
   real(pReal), dimension(size(m,1))                         :: math_eigvalsh
@@ -1137,7 +1137,7 @@ end function math_eigvalsh
 !> but apparently more stable solution and has general LAPACK powered version for arbritrary sized
 !> matrices as fallback
 !--------------------------------------------------------------------------------------------------
-function math_eigvalsh33(m)
+pure function math_eigvalsh33(m)
 
   real(pReal), intent(in), dimension(3,3) :: m                                                      !< 3x3 symmetric matrix to compute eigenvalues of
   real(pReal), dimension(3) :: math_eigvalsh33,I
@@ -1432,9 +1432,11 @@ subroutine selfTest
     error stop 'math_LeviCivita'
 
   normal_distribution: block
-    real(pReal), dimension(500000) :: r
+    integer, parameter :: N = 1000000
+    real(pReal), dimension(:), allocatable :: r
     real(pReal) :: mu, sigma
 
+    allocate(r(N))
     call random_number(mu)
     call random_number(sigma)
 
@@ -1443,11 +1445,11 @@ subroutine selfTest
 
     call math_normal(r,mu,sigma)
 
-    if (abs(mu -sum(r)/real(size(r),pReal))>5.0e-2_pReal) &
+    if (abs(mu -sum(r)/real(N,pReal))>5.0e-2_pReal) &
       error stop 'math_normal(mu)'
 
-    mu = sum(r)/real(size(r),pReal)
-    if (abs(sigma**2 -1.0_pReal/real(size(r)-1,pReal) * sum((r-mu)**2))/sigma > 5.0e-2_pReal) &
+    mu = sum(r)/real(N,pReal)
+    if (abs(sigma**2 -1.0_pReal/real(N-1,pReal) * sum((r-mu)**2))/sigma > 5.0e-2_pReal) &
       error stop 'math_normal(sigma)'
   end block normal_distribution
 
diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90
index 37e6fe7bf..1917d81c9 100644
--- a/src/phase_mechanical.f90
+++ b/src/phase_mechanical.f90
@@ -168,17 +168,17 @@ submodule(phase) mechanical
       integer,     intent(in) :: ph,en
     end function plastic_dislotwin_homogenizedC
 
-    module function elastic_C66(ph,en) result(C66)
+    pure module function elastic_C66(ph,en) result(C66)
       real(pReal), dimension(6,6) :: C66
       integer,     intent(in) :: ph, en
     end function elastic_C66
 
-    module function elastic_mu(ph,en) result(mu)
+    pure module function elastic_mu(ph,en) result(mu)
       real(pReal) :: mu
       integer, intent(in) :: ph, en
     end function elastic_mu
 
-    module function elastic_nu(ph,en) result(nu)
+    pure module function elastic_nu(ph,en) result(nu)
       real(pReal) :: nu
       integer, intent(in) :: ph, en
     end function elastic_nu
diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90
index acbbac236..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  -+>>>'
@@ -86,7 +86,7 @@ end subroutine elastic_init
 !--------------------------------------------------------------------------------------------------
 !> @brief return 6x6 elasticity tensor
 !--------------------------------------------------------------------------------------------------
-module function elastic_C66(ph,en) result(C66)
+pure module function elastic_C66(ph,en) result(C66)
 
   integer, intent(in) :: &
     ph, &
@@ -140,7 +140,7 @@ end function elastic_C66
 !--------------------------------------------------------------------------------------------------
 !> @brief return shear modulus
 !--------------------------------------------------------------------------------------------------
-module function elastic_mu(ph,en) result(mu)
+pure module function elastic_mu(ph,en) result(mu)
 
   integer, intent(in) :: &
     ph, &
@@ -157,7 +157,7 @@ end function elastic_mu
 !--------------------------------------------------------------------------------------------------
 !> @brief return Poisson ratio
 !--------------------------------------------------------------------------------------------------
-module function elastic_nu(ph,en) result(nu)
+pure module function elastic_nu(ph,en) result(nu)
 
   integer, intent(in) :: &
     ph, &
diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90
index 6c0e1e0cc..03eb27f31 100644
--- a/src/phase_mechanical_plastic_kinehardening.f90
+++ b/src/phase_mechanical_plastic_kinehardening.f90
@@ -86,6 +86,8 @@ module function plastic_kinehardening_init() result(myPlasticity)
   print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:kinehardening init  -+>>>'
   print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
+  print'(/,1x,a)', 'J.A. Wollmershauser et al., International Journal of Fatigue 36:181–193, 2012'
+  print'(  1x,a)', 'https://doi.org/10.1016/j.ijfatigue.2011.07.008'
 
   phases => config_material%get('phase')
   allocate(param(phases%length))
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), &