diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7b4203ae2..0f488ec37 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -6,7 +6,7 @@ cmake_minimum_required (VERSION 3.10.0 FATAL_ERROR)
 # Find PETSc from system environment
 set(PETSC_DIR $ENV{PETSC_DIR})
 if (PETSC_DIR STREQUAL "")
-  message (FATAL_ERROR "PETSc location (PETSC_DIR) is not defined")
+    message (FATAL_ERROR "PETSc location (PETSC_DIR) is not defined")
 endif ()
 
 set (petsc_conf_variables "${PETSC_DIR}/lib/petsc/conf/variables")
@@ -118,7 +118,7 @@ endif ()
 list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
 
 if (CMAKE_BUILD_TYPE STREQUAL "")
-  set (CMAKE_BUILD_TYPE "RELEASE")
+    set (CMAKE_BUILD_TYPE "RELEASE")
 endif ()
 
 # Predefined sets for OPTIMIZATION/OPENMP based on BUILD_TYPE
@@ -165,13 +165,13 @@ add_definitions (-DDAMASKVERSION="${DAMASK_V}")
 add_definitions (-DPETSc)
 
 if    (CMAKE_Fortran_COMPILER_ID STREQUAL "Intel")
-  include(Compiler-Intel)
+    include (Compiler-Intel)
 elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
-  include(Compiler-GNU)
+    include (Compiler-GNU)
 elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "PGI")
-  include(Compiler-PGI)
+    include (Compiler-PGI)
 else ()
-  message (FATAL_ERROR "Compiler type (CMAKE_Fortran_COMPILER_ID) not recognized")
+    message (FATAL_ERROR "Compiler type (CMAKE_Fortran_COMPILER_ID) not recognized")
 endif ()
 
 
@@ -179,8 +179,8 @@ set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${BUILDCMD_PRE} ${OPENMP_FLAGS} ${
 set (CMAKE_Fortran_LINK_EXECUTABLE           "${BUILDCMD_PRE} ${PETSC_LINKER} ${OPENMP_FLAGS}   ${OPTIMIZATION_FLAGS} ${LINKER_FLAGS}")
 
 if (CMAKE_BUILD_TYPE STREQUAL "DEBUG")
-  set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${DEBUG_FLAGS}")
-  set (CMAKE_Fortran_LINK_EXECUTABLE           "${CMAKE_Fortran_LINK_EXECUTABLE} ${DEBUG_FLAGS}")
+    set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${DEBUG_FLAGS}")
+    set (CMAKE_Fortran_LINK_EXECUTABLE           "${CMAKE_Fortran_LINK_EXECUTABLE} ${DEBUG_FLAGS}")
 endif ()
 
 set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}   "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${PETSC_INCLUDES} ${BUILDCMD_POST}")
diff --git a/examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Ferrite.config b/examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Ferrite.config
index 87b4e2312..5af1eee11 100644
--- a/examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Ferrite.config
+++ b/examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Ferrite.config
@@ -18,4 +18,4 @@ tau0_slip                95.e6  97.e6           # per family, optimization long
 tausat_slip             222.e6 412.7e6          # per family, optimization long simplex 109
 h0_slipslip             1000.0e6
 interaction_slipslip    1 1 1.4 1.4 1.4 1.4
-w0_slip                 2.0
+a_slip                 2.0
diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py
index 36ce3b8e2..db9f243c9 100644
--- a/python/damask/grid_filters.py
+++ b/python/damask/grid_filters.py
@@ -338,7 +338,7 @@ def node_coord0_gridSizeOrigin(coord0,ordered=False):
     Parameters
     ----------
     coord0 : numpy.ndarray
-        array of undeformed nodal coordinates
+        array of undeformed nodal coordinates.
     ordered : bool, optional
         expect coord0 data to be ordered (x fast, z slow).
 
@@ -365,7 +365,19 @@ def node_coord0_gridSizeOrigin(coord0,ordered=False):
 
 
 def regrid(size,F,new_grid):
-    """tbd."""
+    """
+    Return mapping from coordinates in deformed configuration to a regular grid.
+
+    Parameters
+    ----------
+    size : numpy.ndarray
+        physical size
+    F : numpy.ndarray
+        deformation gradient field
+    new_grid : numpy.ndarray
+        new grid for undeformed coordinates
+
+    """
     c = cell_coord0(F.shape[:3][::-1],size) \
       + cell_displacement_avg(size,F) \
       + cell_displacement_fluct(size,F)
@@ -376,4 +388,4 @@ def regrid(size,F,new_grid):
         c[np.where(c[:,:,:,d]>outer[d])] -= outer[d]
 
     tree = spatial.cKDTree(c.reshape((-1,3)),boxsize=outer)
-    return tree.query(cell_coord0(new_grid,outer))[1]
+    return tree.query(cell_coord0(new_grid,outer))[1].flatten()
diff --git a/python/tests/test_grid_filters.py b/python/tests/test_grid_filters.py
index a5455e1ae..ed1336842 100644
--- a/python/tests/test_grid_filters.py
+++ b/python/tests/test_grid_filters.py
@@ -77,3 +77,9 @@ class TestGridFilters:
          grid = np.random.randint(8,32,(3))
          F    = np.broadcast_to(np.random.random((3,3)), tuple(grid)+(3,3))                         # noqa
          assert np.allclose(eval('grid_filters.{}_displacement_fluct(size,F)'.format(mode)),0.0)
+
+    def test_regrid(self):
+         size = np.random.random(3)
+         grid = np.random.randint(8,32,(3))
+         F    = np.broadcast_to(np.eye(3), tuple(grid[::-1])+(3,3))
+         assert all(grid_filters.regrid(size,F,grid) == np.arange(grid.prod()))
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index b2eee5561..2e4336698 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -4,37 +4,37 @@ if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
     SET_SOURCE_FILES_PROPERTIES("lattice.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-240")
 endif()
 
-file(GLOB damask-sources *.f90 *.c)
+file (GLOB damask-sources *.f90 *.c)
 
 # probably we should have subfolders for abaqus and MSC.Marc
-list(FILTER damask-sources EXCLUDE REGEX ".*CPFEM.f90")
-list(FILTER damask-sources EXCLUDE REGEX ".*DAMASK_marc.*.f90")
-list(FILTER damask-sources EXCLUDE REGEX ".*mesh_marc.*.f90")
-list(FILTER damask-sources EXCLUDE REGEX ".*mesh_abaqus.*.f90")
-list(FILTER damask-sources EXCLUDE REGEX ".*commercialFEM_fileList.*.f90")
+list (FILTER damask-sources EXCLUDE REGEX ".*CPFEM.f90")
+list (FILTER damask-sources EXCLUDE REGEX ".*DAMASK_marc.*.f90")
+list (FILTER damask-sources EXCLUDE REGEX ".*mesh_marc.*.f90")
+list (FILTER damask-sources EXCLUDE REGEX ".*mesh_abaqus.*.f90")
+list (FILTER damask-sources EXCLUDE REGEX ".*commercialFEM_fileList.*.f90")
 
 
 if (PROJECT_NAME STREQUAL "damask-grid")
 
-      list(FILTER damask-sources EXCLUDE REGEX ".*mesh_FEM.*.f90")
-      file(GLOB grid-sources grid/*.f90)
+    list (FILTER damask-sources EXCLUDE REGEX ".*mesh_FEM.*.f90")
+    file (GLOB grid-sources grid/*.f90)
 
-      if(NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
-        add_executable(DAMASK_spectral ${damask-sources} ${grid-sources})
+    if (NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
+        add_executable (DAMASK_spectral ${damask-sources} ${grid-sources})
         install (TARGETS DAMASK_spectral RUNTIME DESTINATION bin)
-      else()
-        add_library(DAMASK_spectral OBJECT ${damask-sources} ${grid-sources})
+    else()
+        add_library (DAMASK_spectral OBJECT ${damask-sources} ${grid-sources})
         exec_program (mktemp OUTPUT_VARIABLE nothing)
         exec_program (mktemp ARGS -d OUTPUT_VARIABLE black_hole)
         install (PROGRAMS ${nothing} DESTINATION ${black_hole})
-      endif()
+    endif()
 
 elseif (PROJECT_NAME STREQUAL "damask-mesh")
       
-      list(FILTER damask-sources EXCLUDE REGEX ".*mesh_grid.*.f90")
-      file(GLOB mesh-sources mesh/*.f90)
+    list (FILTER damask-sources EXCLUDE REGEX ".*mesh_grid.*.f90")
+    file (GLOB mesh-sources mesh/*.f90)
 
-      add_executable(DAMASK_FEM ${damask-sources} ${mesh-sources})
-      install (TARGETS DAMASK_FEM RUNTIME DESTINATION bin)
+    add_executable (DAMASK_FEM ${damask-sources} ${mesh-sources})
+    install (TARGETS DAMASK_FEM RUNTIME DESTINATION bin)
       
 endif()
diff --git a/src/CPFEM.f90 b/src/CPFEM.f90
index 1ef5f1389..fd406219f 100644
--- a/src/CPFEM.f90
+++ b/src/CPFEM.f90
@@ -152,7 +152,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
                                                      i, j, k, l, m, n, ph, homog, mySource
  logical                                             updateJaco                                     ! flag indicating if JAcobian has to be updated
 
- elCP = mesh_FEasCP('elem',elFE)
+ elCP = mesh_FEM2DAMASK_elem(elFE)
 
  if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt &
      .and. elCP == debug_e .and. ip == debug_i) then
diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90
index ba7751372..f7efb5279 100644
--- a/src/DAMASK_marc.f90
+++ b/src/DAMASK_marc.f90
@@ -269,7 +269,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
    if (timinc < theDelta .and. theInc == inc .and. lastLovl /= lovl) &                              ! first after cutback
      computationMode = CPFEM_RESTOREJACOBIAN
  elseif (lovl == 6) then                                                                            ! stress requested by marc
-   cp_en = mesh_FEasCP('elem',m(1))
+   cp_en = mesh_FEM2DAMASK_elem(m(1))
    if (cptim > theTime .or. inc /= theInc) then                                                     ! reached "convergence"
      terminallyIll = .false.
      cycleCounter = -1                                                                              ! first calc step increments this to cycle = 0
@@ -391,7 +391,7 @@ subroutine flux(f,ts,n,time)
   real(pReal), dimension(2),           intent(out) :: &
     f
 
-  call thermal_conduction_getSourceAndItsTangent(f(1), f(2), ts(3), n(3),mesh_FEasCP('elem',n(1)))
+  call thermal_conduction_getSourceAndItsTangent(f(1), f(2), ts(3), n(3),mesh_FEM2DAMASK_elem(n(1)))
 
  end subroutine flux
 
diff --git a/src/FEsolving.f90 b/src/FEsolving.f90
index f81110414..fc04499fa 100644
--- a/src/FEsolving.f90
+++ b/src/FEsolving.f90
@@ -4,10 +4,10 @@
 !> @brief global variables for flow control
 !--------------------------------------------------------------------------------------------------
 module FEsolving
-  use prec
    
   implicit none
- 
+  public
+
   logical :: &
     terminallyIll = .false.                                                                         !< at least one material point is terminally ill
 
diff --git a/src/IO.f90 b/src/IO.f90
index 9406c349b..da3dcbf3c 100644
--- a/src/IO.f90
+++ b/src/IO.f90
@@ -400,15 +400,18 @@ pure function IO_lc(string)
   character(len=*), intent(in) :: string                                                            !< string to convert
   character(len=len(string))   :: IO_lc
 
-  character(26), parameter :: LOWER = 'abcdefghijklmnopqrstuvwxyz'
-  character(26), parameter :: UPPER = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
+  character(len=*),          parameter :: LOWER = 'abcdefghijklmnopqrstuvwxyz'
+  character(len=len(LOWER)), parameter :: UPPER = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
 
-  integer                  :: i,n
+  integer :: i,n
 
   do i=1,len(string)
-    IO_lc(i:i) = string(i:i)
-    n = index(UPPER,IO_lc(i:i))
-    if (n/=0) IO_lc(i:i) = LOWER(n:n)
+    n = index(UPPER,string(i:i))
+    if(n/=0) then
+      IO_lc(i:i) = LOWER(n:n)
+    else
+      IO_lc(i:i) = string(i:i)
+    endif
   enddo
 
 end function IO_lc
@@ -549,6 +552,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
 ! math errors
     case (400)
       msg = 'matrix inversion error'
+    case (401)
+      msg = 'error in Eigenvalue calculation'
     case (402)
       msg = 'invalid orientation specified'
 
@@ -876,11 +881,10 @@ integer function verifyIntValue(string)
  
   character(len=*), intent(in) :: string                                                            !< string for conversion to int value
 
-  integer                      :: readStatus, invalidWhere
+  integer                      :: readStatus
   character(len=*), parameter  :: VALIDCHARS = '0123456789+- '
  
-  invalidWhere = verify(string,VALIDCHARS)
-  valid: if (invalidWhere == 0) then
+  valid: if (verify(string,VALIDCHARS) == 0) then
     read(string,*,iostat=readStatus) verifyIntValue
     if (readStatus /= 0) call IO_error(111,ext_msg=string)
   else valid
@@ -898,11 +902,10 @@ real(pReal) function verifyFloatValue(string)
  
   character(len=*), intent(in) :: string                                                            !< string for conversion to float value
   
-  integer                      :: readStatus, invalidWhere
+  integer                      :: readStatus
   character(len=*), parameter  :: VALIDCHARS = '0123456789eE.+- '
  
-  invalidWhere = verify(string,VALIDCHARS)
-  valid: if (invalidWhere == 0) then
+  valid: if (verify(string,VALIDCHARS) == 0) then
     read(string,*,iostat=readStatus) verifyFloatValue
     if (readStatus /= 0) call IO_error(112,ext_msg=string)
   else valid
diff --git a/src/config.f90 b/src/config.f90
index bcf9b4b3d..58b0d69a1 100644
--- a/src/config.f90
+++ b/src/config.f90
@@ -27,7 +27,7 @@ module config
     config_numerics, &
     config_debug
  
-  character(len=pStringLen), dimension(:), allocatable, public, protected :: &
+  character(len=pStringLen),    public, protected, allocatable, dimension(:) :: &
     config_name_phase, &                                                                            !< name of each phase
     config_name_homogenization, &                                                                   !< name of each homogenization
     config_name_crystallite, &                                                                      !< name of each crystallite setting
diff --git a/src/element.f90 b/src/element.f90
index a38036f93..de4000700 100644
--- a/src/element.f90
+++ b/src/element.f90
@@ -3,7 +3,6 @@
 !> @author Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH
 !--------------------------------------------------------------------------------------------------
 module element
-  use prec
   use IO
  
   implicit none
diff --git a/src/math.f90 b/src/math.f90
index 4be2d4fc9..c63f30881 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -8,14 +8,16 @@
 module math
   use prec
   use IO
-  use debug
   use numerics
+
   implicit none
   public
 #if __INTEL_COMPILER >= 1900
   ! do not make use associated entities available to other modules
   private :: &
-    prec
+    prec, &
+    IO, &
+    numerics
 #endif
 
   real(pReal),    parameter :: PI = acos(-1.0_pReal)                                                !< ratio of a circle's circumference to its diameter
@@ -24,24 +26,24 @@ module math
   complex(pReal), parameter :: TWOPIIMG = cmplx(0.0_pReal,2.0_pReal*PI)                             !< Re(0.0), Im(2xPi)
 
   real(pReal), dimension(3,3), parameter :: &
-    MATH_I3 = reshape([&
+    math_I3 = reshape([&
       1.0_pReal,0.0_pReal,0.0_pReal, &
       0.0_pReal,1.0_pReal,0.0_pReal, &
       0.0_pReal,0.0_pReal,1.0_pReal  &
       ],[3,3])                                                                                      !< 3x3 Identity
 
   real(pReal), dimension(6), parameter, private :: &
-    nrmMandel = [&
+    NRMMANDEL = [&
       1.0_pReal,       1.0_pReal,       1.0_pReal, &
       sqrt(2.0_pReal), sqrt(2.0_pReal), sqrt(2.0_pReal) ]                                           !< weighting for Mandel notation (forward)
 
   real(pReal), dimension(6), parameter , private :: &
-    invnrmMandel = [&
+    INVNRMMANDEL = [&
       1.0_pReal,                 1.0_pReal,                 1.0_pReal, &
       1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal) ]             !< weighting for Mandel notation (backward)
 
   integer, dimension (2,6), parameter, private :: &
-    mapNye = reshape([&
+    MAPNYE = reshape([&
       1,1, &
       2,2, &
       3,3, &
@@ -51,7 +53,7 @@ module math
       ],[2,6])                                                                                      !< arrangement in Nye notation.
 
   integer, dimension (2,6), parameter, private :: &
-    mapVoigt = reshape([&
+    MAPVOIGT = reshape([&
       1,1, &
       2,2, &
       3,3, &
@@ -61,7 +63,7 @@ module math
       ],[2,6])                                                                                      !< arrangement in Voigt notation
 
   integer, dimension (2,9), parameter, private :: &
-    mapPlain = reshape([&
+    MAPPLAIN = reshape([&
       1,1, &
       1,2, &
       1,3, &
@@ -129,13 +131,13 @@ recursive subroutine math_sort(a, istart, iend, sortDim)
   integer, dimension(:,:), intent(inout) :: a
   integer, intent(in),optional :: istart,iend, sortDim
   integer :: ipivot,s,e,d
-  
+
   if(present(istart)) then
     s = istart
   else
     s = lbound(a,2)
   endif
-  
+
   if(present(iend)) then
     e = iend
   else
@@ -147,7 +149,7 @@ recursive subroutine math_sort(a, istart, iend, sortDim)
   else
     d = 1
   endif
-   
+
   if (s < e) then
     ipivot = qsort_partition(a,s, e, d)
     call math_sort(a, s, ipivot-1, d)
@@ -232,14 +234,14 @@ end function math_range
 !--------------------------------------------------------------------------------------------------
 !> @brief second rank identity tensor of specified dimension
 !--------------------------------------------------------------------------------------------------
-pure function math_identity2nd(dimen)
+pure function math_identity2nd(d)
 
-  integer, intent(in) :: dimen                                                                      !< tensor dimension
+  integer, intent(in) :: d                                                                          !< tensor dimension
   integer :: i
-  real(pReal), dimension(dimen,dimen) :: math_identity2nd
+  real(pReal), dimension(d,d) :: math_identity2nd
 
   math_identity2nd = 0.0_pReal
-  do i=1, dimen
+  do i=1,d
     math_identity2nd(i,i) = 1.0_pReal
   enddo
 
@@ -250,16 +252,17 @@ end function math_identity2nd
 !> @brief symmetric fourth rank identity tensor of specified dimension
 !  from http://en.wikipedia.org/wiki/Tensor_derivative_(continuum_mechanics)#Derivative_of_a_second-order_tensor_with_respect_to_itself
 !--------------------------------------------------------------------------------------------------
-pure function math_identity4th(dimen)
+pure function math_identity4th(d)
 
-  integer, intent(in) :: dimen                                                                      !< tensor dimension
+  integer, intent(in) :: d                                                                          !< tensor dimension
   integer :: i,j,k,l
-  real(pReal), dimension(dimen,dimen,dimen,dimen) ::  math_identity4th
-  real(pReal), dimension(dimen,dimen)             ::  identity2nd
+  real(pReal), dimension(d,d,d,d) :: math_identity4th
+  real(pReal), dimension(d,d)     :: identity2nd
 
-  identity2nd = math_identity2nd(dimen)
-  forall(i=1:dimen,j=1:dimen,k=1:dimen,l=1:dimen) &
+  identity2nd = math_identity2nd(d)
+  do i=1,d; do j=1,d; do k=1,d; do l=1,d
     math_identity4th(i,j,k,l) = 0.5_pReal*(identity2nd(i,k)*identity2nd(j,l)+identity2nd(i,l)*identity2nd(j,k))
+  enddo; enddo; enddo; enddo
 
 end function math_identity4th
 
@@ -324,7 +327,9 @@ pure function math_outer(A,B)
   real(pReal), dimension(size(A,1),size(B,1)) ::  math_outer
   integer :: i,j
 
-  forall(i=1:size(A,1),j=1:size(B,1)) math_outer(i,j) = A(i)*B(j)
+  do i=1,size(A,1); do j=1,size(B,1)
+    math_outer(i,j) = A(i)*B(j)
+  enddo; enddo
 
 end function math_outer
 
@@ -347,11 +352,13 @@ end function math_inner
 !--------------------------------------------------------------------------------------------------
 real(pReal) pure function math_mul33xx33(A,B)
 
-  real(pReal), dimension(3,3), intent(in) ::  A,B
+  real(pReal), dimension(3,3), intent(in) :: A,B
   integer :: i,j
-  real(pReal), dimension(3,3) ::              C
+  real(pReal), dimension(3,3) ::             C
 
-  forall(i=1:3,j=1:3) C(i,j) = A(i,j) * B(i,j)
+  do i=1,3; do j=1,3
+    C(i,j) = A(i,j) * B(i,j)
+  enddo; enddo
   math_mul33xx33 = sum(C)
 
 end function math_mul33xx33
@@ -362,12 +369,14 @@ end function math_mul33xx33
 !--------------------------------------------------------------------------------------------------
 pure function math_mul3333xx33(A,B)
 
+  real(pReal), dimension(3,3,3,3), intent(in) :: A
+  real(pReal), dimension(3,3),     intent(in) :: B
   real(pReal), dimension(3,3) :: math_mul3333xx33
-  real(pReal), dimension(3,3,3,3), intent(in) ::  A
-  real(pReal), dimension(3,3), intent(in) ::  B
   integer :: i,j 
 
-  forall(i = 1:3,j = 1:3) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
+  do i=1,3; do j=1,3
+    math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
+  enddo; enddo
 
 end function math_mul3333xx33
 
@@ -378,12 +387,13 @@ end function math_mul3333xx33
 pure function math_mul3333xx3333(A,B)
 
   integer :: i,j,k,l
-  real(pReal), dimension(3,3,3,3), intent(in) ::  A
-  real(pReal), dimension(3,3,3,3), intent(in) ::  B
+  real(pReal), dimension(3,3,3,3), intent(in) :: A
+  real(pReal), dimension(3,3,3,3), intent(in) :: B
   real(pReal), dimension(3,3,3,3) :: math_mul3333xx3333
 
-  forall(i = 1:3,j = 1:3, k = 1:3, l= 1:3) &
+  do i=1,3; do j=1,3; do k=1,3; do l=1,3
     math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
+  enddo; enddo; enddo; enddo
 
 end function math_mul3333xx3333
 
@@ -393,12 +403,11 @@ end function math_mul3333xx3333
 !--------------------------------------------------------------------------------------------------
 pure function math_exp33(A,n)
 
-  integer :: i
-  integer, intent(in), optional :: n
   real(pReal), dimension(3,3), intent(in) :: A
+  integer, intent(in), optional :: n
   real(pReal), dimension(3,3) :: B, math_exp33
   real(pReal) :: invFac
-  integer     :: n_
+  integer     :: n_,i
 
   if (present(n)) then
     n_ = n
@@ -646,7 +655,7 @@ pure function math_33to9(m33)
   integer :: i
   
   do i = 1, 9
-    math_33to9(i) = m33(mapPlain(1,i),mapPlain(2,i))
+    math_33to9(i) = m33(MAPPLAIN(1,i),MAPPLAIN(2,i))
   enddo
 
 end function math_33to9
@@ -663,7 +672,7 @@ pure function math_9to33(v9)
   integer :: i
 
   do i = 1, 9
-    math_9to33(mapPlain(1,i),mapPlain(2,i)) = v9(i)
+    math_9to33(MAPPLAIN(1,i),MAPPLAIN(2,i)) = v9(i)
   enddo
 
 end function math_9to33
@@ -685,13 +694,13 @@ pure function math_sym33to6(m33,weighted)
   integer :: i
   
   if(present(weighted)) then
-    w = merge(nrmMandel,1.0_pReal,weighted)
+    w = merge(NRMMANDEL,1.0_pReal,weighted)
   else
-    w = nrmMandel
+    w = NRMMANDEL
   endif
 
   do i = 1, 6
-    math_sym33to6(i) = w(i)*m33(mapNye(1,i),mapNye(2,i))
+    math_sym33to6(i) = w(i)*m33(MAPNYE(1,i),MAPNYE(2,i))
   enddo 
 
 end function math_sym33to6
@@ -713,14 +722,14 @@ pure function math_6toSym33(v6,weighted)
   integer :: i
   
   if(present(weighted)) then
-    w = merge(invnrmMandel,1.0_pReal,weighted)
+    w = merge(INVNRMMANDEL,1.0_pReal,weighted)
   else
-    w = invnrmMandel
+    w = INVNRMMANDEL
   endif
 
   do i=1,6
-    math_6toSym33(mapNye(1,i),mapNye(2,i)) = w(i)*v6(i)
-    math_6toSym33(mapNye(2,i),mapNye(1,i)) = w(i)*v6(i)
+    math_6toSym33(MAPNYE(1,i),MAPNYE(2,i)) = w(i)*v6(i)
+    math_6toSym33(MAPNYE(2,i),MAPNYE(1,i)) = w(i)*v6(i)
   enddo
 
 end function math_6toSym33
@@ -737,7 +746,7 @@ pure function math_3333to99(m3333)
   integer :: i,j
 
   do i=1,9; do j=1,9
-    math_3333to99(i,j) = m3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j))
+    math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
   enddo; enddo
 
 end function math_3333to99
@@ -754,7 +763,7 @@ pure function math_99to3333(m99)
   integer :: i,j
 
   do i=1,9; do j=1,9
-    math_99to3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) = m99(i,j)
+    math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
   enddo; enddo
 
 end function math_99to3333
@@ -776,13 +785,13 @@ pure function math_sym3333to66(m3333,weighted)
   integer :: i,j
   
   if(present(weighted)) then
-    w = merge(nrmMandel,1.0_pReal,weighted)
+    w = merge(NRMMANDEL,1.0_pReal,weighted)
   else
-    w = nrmMandel
+    w = NRMMANDEL
   endif
 
   do i=1,6; do j=1,6
-    math_sym3333to66(i,j) = w(i)*w(j)*m3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j))
+    math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j))
   enddo; enddo
 
 end function math_sym3333to66
@@ -804,16 +813,16 @@ pure function math_66toSym3333(m66,weighted)
   integer :: i,j
   
   if(present(weighted)) then
-    w = merge(invnrmMandel,1.0_pReal,weighted)
+    w = merge(INVNRMMANDEL,1.0_pReal,weighted)
   else
-    w = invnrmMandel
+    w = INVNRMMANDEL
   endif
 
   do i=1,6; do j=1,6
-    math_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j)
-    math_66toSym3333(mapNye(2,i),mapNye(1,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j)
-    math_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(2,j),mapNye(1,j)) = w(i)*w(j)*m66(i,j)
-    math_66toSym3333(mapNye(2,i),mapNye(1,i),mapNye(2,j),mapNye(1,j)) = w(i)*w(j)*m66(i,j)
+    math_66toSym3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j)) = w(i)*w(j)*m66(i,j)
+    math_66toSym3333(MAPNYE(2,i),MAPNYE(1,i),MAPNYE(1,j),MAPNYE(2,j)) = w(i)*w(j)*m66(i,j)
+    math_66toSym3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(2,j),MAPNYE(1,j)) = w(i)*w(j)*m66(i,j)
+    math_66toSym3333(MAPNYE(2,i),MAPNYE(1,i),MAPNYE(2,j),MAPNYE(1,j)) = w(i)*w(j)*m66(i,j)
   enddo; enddo
 
 end function math_66toSym3333
@@ -829,10 +838,10 @@ pure function math_Voigt66to3333(m66)
   integer :: i,j
 
   do i=1,6; do j=1, 6
-    math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j)
-    math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j)
-    math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(2,j),mapVoigt(1,j)) = m66(i,j)
-    math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(2,j),mapVoigt(1,j)) = m66(i,j)
+    math_Voigt66to3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j)
+    math_Voigt66to3333(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j)
+    math_Voigt66to3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = m66(i,j)
+    math_Voigt66to3333(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = m66(i,j)
   enddo; enddo
 
 end function math_Voigt66to3333
@@ -1201,7 +1210,7 @@ end function math_invariantsSym33
 integer pure function math_factorial(n)
 
   integer, intent(in) :: n
-  
+
   math_factorial = product(math_range(n))
 
 end function math_factorial
@@ -1233,7 +1242,7 @@ integer pure function math_multinomial(alpha)
 
   integer, intent(in), dimension(:) :: alpha
   integer :: i
- 
+
   math_multinomial = 1
   do i = 1, size(alpha)
     math_multinomial = math_multinomial*math_binomial(sum(alpha(1:i)),alpha(i))
@@ -1293,12 +1302,12 @@ end function math_clip
 !> @brief check correctness of (some) math functions
 !--------------------------------------------------------------------------------------------------
 subroutine unitTest
- 
+
   integer, dimension(2,4) :: &
     sort_in_   = reshape([+1,+5,  +5,+6,  -1,-1,  +3,-2],[2,4])
   integer, dimension(2,4), parameter :: &
     sort_out_  = reshape([-1,-1,  +1,+5,  +5,+6,  +3,-2],[2,4])
-    
+
   integer,     dimension(5) :: range_out_ = [1,2,3,4,5]
 
   real(pReal)                 :: det
@@ -1308,8 +1317,12 @@ subroutine unitTest
   real(pReal), dimension(3,3) :: t33,t33_2
   real(pReal), dimension(6,6) :: t66
   real(pReal), dimension(9,9) :: t99,t99_2
+  real(pReal), dimension(:,:), &
+               allocatable    :: txx,txx_2
+  real(pReal)                 :: r
+  integer                     :: d
   logical                     :: e
-  
+
   if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,3.0_pReal,3.0_pReal,3.0_pReal] - &
               math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2,3,0])) > tol_math_check)) &
     call IO_error(0,ext_msg='math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]')
@@ -1322,15 +1335,19 @@ subroutine unitTest
               math_expand([1.0_pReal,2.0_pReal],[1,2,3])) > tol_math_check)) &
     call IO_error(0,ext_msg='math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]')
 
-
   call math_sort(sort_in_,1,3,2)
   if(any(sort_in_ /= sort_out_)) &
     call IO_error(0,ext_msg='math_sort')
-  
+
   if(any(math_range(5) /= range_out_)) &
     call IO_error(0,ext_msg='math_range')
- 
- 
+
+   if(any(dNeq(math_exp33(math_I3,0),math_I3))) &
+    call IO_error(0,ext_msg='math_exp33(math_I3,1)')
+
+  if(any(dNeq(math_exp33(math_I3,256),exp(1.0_pReal)*math_I3))) &
+    call IO_error(0,ext_msg='math_exp33(math_I3,256)')
+
   call random_number(v9)
   if(any(dNeq(math_33to9(math_9to33(v9)),v9))) &
     call IO_error(0,ext_msg='math_33to9/math_9to33')
@@ -1338,56 +1355,66 @@ subroutine unitTest
   call random_number(t99)
   if(any(dNeq(math_3333to99(math_99to3333(t99)),t99))) &
     call IO_error(0,ext_msg='math_3333to99/math_99to3333')
-  
+
   call random_number(v6)
   if(any(dNeq(math_sym33to6(math_6toSym33(v6)),v6))) &
     call IO_error(0,ext_msg='math_sym33to6/math_6toSym33')
-  
+
   call random_number(t66)
   if(any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66))) &
     call IO_error(0,ext_msg='math_sym3333to66/math_66toSym3333')
-    
+
   call random_number(v6)
   if(any(dNeq0(math_6toSym33(v6) - math_symmetric33(math_6toSym33(v6))))) &
     call IO_error(0,ext_msg='math_symmetric33')
- 
+
   call random_number(v3_1)
   call random_number(v3_2)
   call random_number(v3_3)
   call random_number(v3_4)
-  
+
   if(dNeq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0, &
           math_volTetrahedron(v3_1,v3_2,v3_3,v3_4),tol=1.0e-12_pReal)) &
   call IO_error(0,ext_msg='math_volTetrahedron')
-    
+
   call random_number(t33)
   if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) &
     call IO_error(0,ext_msg='math_det33/math_detSym33')
   
+  if(any(dNeq0(math_identity2nd(3),math_inv33(math_I3)))) &
+    call IO_error(0,ext_msg='math_inv33(math_I3)')
+
   do while(abs(math_det33(t33))<1.0e-9_pReal)
     call random_number(t33)
   enddo
   if(any(dNeq0(matmul(t33,math_inv33(t33)) - math_identity2nd(3),tol=1.0e-9_pReal))) &
     call IO_error(0,ext_msg='math_inv33')
-    
+
   call math_invert33(t33_2,det,e,t33)
   if(any(dNeq0(matmul(t33,t33_2) - math_identity2nd(3),tol=1.0e-9_pReal)) .or. e) &
     call IO_error(0,ext_msg='math_invert33: T:T^-1 != I')
   if(dNeq(det,math_det33(t33),tol=1.0e-12_pReal)) &
     call IO_error(0,ext_msg='math_invert33 (determinant)')
-    
+
   call math_invert(t33_2,e,t33)
   if(any(dNeq0(matmul(t33,t33_2) - math_identity2nd(3),tol=1.0e-9_pReal)) .or. e) &
     call IO_error(0,ext_msg='math_invert t33')
-    
+
   t33_2 = transpose(math_rotationalPart33(t33))
   if(any(dNeq0(math_rotationalPart33(matmul(t33_2,t33)) - MATH_I3,tol=5.0e-4_pReal))) &
     call IO_error(0,ext_msg='math_rotationalPart33')
 
+  call random_number(r)
+  d = int(r*5.0_pReal) + 1
+  txx = math_identity2nd(d)
+  allocate(txx_2(d,d))
+  call math_invert(txx_2,e,txx)
+  if(any(dNeq0(txx_2,txx) .or. e)) &
+    call IO_error(0,ext_msg='math_invert(txx)/math_identity2nd')
 
   call math_invert(t99_2,e,t99) ! not sure how likely it is that we get a singular matrix 
   if(any(dNeq0(matmul(t99_2,t99)-math_identity2nd(9),tol=1.0e-9_pReal)) .or. e) &
-    call IO_error(0,ext_msg='math_invert t99')
+    call IO_error(0,ext_msg='math_invert(t99)')
 
   if(any(dNeq(math_clip([4.0_pReal,9.0_pReal],5.0_pReal,6.5_pReal),[5.0_pReal,6.5_pReal]))) &
      call IO_error(0,ext_msg='math_clip')
diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90
index b1ba6b590..bce392d78 100644
--- a/src/mesh_grid.f90
+++ b/src/mesh_grid.f90
@@ -18,10 +18,10 @@ module mesh_grid
   use discretization
   use geometry_plastic_nonlocal
   use FEsolving
- 
+
   implicit none
   private
-  
+
   integer,     dimension(3), public, protected :: &
     grid                                                                                            !< (global) grid
   integer,                   public, protected :: &
@@ -32,10 +32,10 @@ module mesh_grid
   real(pReal),               public, protected :: &
     size3, &                                                                                        !< (local) size in 3rd direction
     size3offset                                                                                     !< (local) size offset in 3rd direction
- 
+
   public :: &
     mesh_init
- 
+
 contains
 
 
@@ -45,7 +45,7 @@ contains
 subroutine mesh_init(ip,el)
 
   integer, intent(in), optional :: el, ip                                                           ! for compatibility reasons
-  
+
   include 'fftw3-mpi.f03'
   real(pReal), dimension(3) :: &
     mySize, &                                                                                       !< domain size of this process
@@ -93,7 +93,7 @@ subroutine mesh_init(ip,el)
   call discretization_init(homogenizationAt,microstructureAt, &
                            IPcoordinates0(myGrid,mySize,grid3Offset), &
                            Nodes0(myGrid,mySize,grid3Offset),&
-                           merge((grid(1)+1) * (grid(2)+1) * (grid3+1),&                            ! write bottom layer 
+                           merge((grid(1)+1) * (grid(2)+1) * (grid3+1),&                            ! write bottom layer
                                  (grid(1)+1) * (grid(2)+1) *  grid3,&                               ! do not write bottom layer (is top of rank-1)
                                  worldrank<1))
 
@@ -113,8 +113,8 @@ subroutine mesh_init(ip,el)
 ! geometry information required by the nonlocal CP model
   call geometry_plastic_nonlocal_setIPvolume(reshape([(product(mySize/real(myGrid,pReal)),j=1,product(myGrid))], &
                                                      [1,product(myGrid)]))
-  call geometry_plastic_nonlocal_setIParea        (cellEdgeArea(mySize,myGrid))
-  call geometry_plastic_nonlocal_setIPareaNormal  (cellEdgeNormal(product(myGrid)))
+  call geometry_plastic_nonlocal_setIParea        (cellSurfaceArea(mySize,myGrid))
+  call geometry_plastic_nonlocal_setIPareaNormal  (cellSurfaceNormal(product(myGrid)))
   call geometry_plastic_nonlocal_setIPneighborhood(IPneighborhood(myGrid))
 
 !--------------------------------------------------------------------------------------------------
@@ -127,7 +127,7 @@ end subroutine mesh_init
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Parses geometry file
-!> @details important variables have an implicit "save" attribute. Therefore, this function is 
+!> @details important variables have an implicit "save" attribute. Therefore, this function is
 ! supposed to be called only once!
 !--------------------------------------------------------------------------------------------------
 subroutine readGeom(grid,geomSize,origin,microstructure,homogenization)
@@ -140,7 +140,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure,homogenization)
   integer,     dimension(:), intent(out), allocatable :: &
     microstructure, &
     homogenization
-   
+
   character(len=:),      allocatable :: rawData
   character(len=65536)               :: line
   integer, allocatable, dimension(:) :: chunkPos
@@ -154,9 +154,9 @@ subroutine readGeom(grid,geomSize,origin,microstructure,homogenization)
     l, &                                                                                            !< line counter
     c, &                                                                                            !< counter for # microstructures in line
     o, &                                                                                            !< order of "to" packing
-    e, &                                                                                            !< "element", i.e. spectral collocation point 
+    e, &                                                                                            !< "element", i.e. spectral collocation point
     i, j
-    
+
   grid = -1
   geomSize = -1.0_pReal
 
@@ -169,7 +169,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure,homogenization)
   allocate(character(len=fileLength)::rawData)
   read(fileUnit) rawData
   close(fileUnit)
-  
+
 !--------------------------------------------------------------------------------------------------
 ! get header length
   endPos = index(rawData,IO_EOL)
@@ -196,7 +196,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure,homogenization)
 
     chunkPos = IO_stringPos(trim(line))
     if (chunkPos(1) < 2) cycle                                                                      ! need at least one keyword value pair
-    
+
     select case (IO_lc(IO_StringValue(trim(line),chunkPos,1)) )
       case ('grid')
         if (chunkPos(1) > 6) then
@@ -211,7 +211,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure,homogenization)
             end select
           enddo
         endif
-        
+
       case ('size')
         if (chunkPos(1) > 6) then
           do j = 2,6,2
@@ -225,7 +225,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure,homogenization)
             end select
           enddo
         endif
-        
+
       case ('origin')
         if (chunkPos(1) > 6) then
           do j = 2,6,2
@@ -258,7 +258,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure,homogenization)
 
   allocate(microstructure(product(grid)), source = -1)                                              ! too large in case of MPI (shrink later, not very elegant)
   allocate(homogenization(product(grid)), source = h)                                               ! too large in case of MPI (shrink later, not very elegant)
-     
+
 !--------------------------------------------------------------------------------------------------
 ! read and interpret content
   e = 1
@@ -269,7 +269,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure,homogenization)
     startPos = endPos + 1
     l = l + 1
     chunkPos = IO_stringPos(trim(line))
-    
+
     noCompression: if (chunkPos(1) /= 3) then
       c = chunkPos(1)
       microstructure(e:e+c-1) =  [(IO_intValue(line,chunkPos,i+1), i=0, c-1)]
@@ -309,7 +309,7 @@ function IPcoordinates0(grid,geomSize,grid3Offset)
   integer :: &
     a,b,c, &
     i
- 
+
   i = 0
   do c = 1, grid(3); do b = 1, grid(2); do a = 1, grid(1)
     i = i + 1
@@ -346,37 +346,37 @@ end function nodes0
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate IP interface areas
 !--------------------------------------------------------------------------------------------------
-pure function cellEdgeArea(geomSize,grid)
-  
+pure function cellSurfaceArea(geomSize,grid)
+
   real(pReal), dimension(3), intent(in) :: geomSize                                                 ! size (for this process!)
   integer,     dimension(3), intent(in) :: grid                                                     ! grid (for this process!)
-  
-  real(pReal), dimension(6,1,product(grid)) :: cellEdgeArea
 
-  cellEdgeArea(1:2,1,:) = geomSize(2)/real(grid(2)) * geomSize(3)/real(grid(3))
-  cellEdgeArea(3:4,1,:) = geomSize(3)/real(grid(3)) * geomSize(1)/real(grid(1))
-  cellEdgeArea(5:6,1,:) = geomSize(1)/real(grid(1)) * geomSize(2)/real(grid(2))
-  
-end function cellEdgeArea
+  real(pReal), dimension(6,1,product(grid)) :: cellSurfaceArea
+
+  cellSurfaceArea(1:2,1,:) = geomSize(2)/real(grid(2)) * geomSize(3)/real(grid(3))
+  cellSurfaceArea(3:4,1,:) = geomSize(3)/real(grid(3)) * geomSize(1)/real(grid(1))
+  cellSurfaceArea(5:6,1,:) = geomSize(1)/real(grid(1)) * geomSize(2)/real(grid(2))
+
+end function cellSurfaceArea
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate IP interface areas normals
 !--------------------------------------------------------------------------------------------------
-pure function cellEdgeNormal(nElems)
+pure function cellSurfaceNormal(nElems)
 
   integer, intent(in) :: nElems
-  
-  real(pReal), dimension(3,6,1,nElems) :: cellEdgeNormal
 
-  cellEdgeNormal(1:3,1,1,:) = spread([+1.0_pReal, 0.0_pReal, 0.0_pReal],2,nElems)
-  cellEdgeNormal(1:3,2,1,:) = spread([-1.0_pReal, 0.0_pReal, 0.0_pReal],2,nElems)
-  cellEdgeNormal(1:3,3,1,:) = spread([ 0.0_pReal,+1.0_pReal, 0.0_pReal],2,nElems)
-  cellEdgeNormal(1:3,4,1,:) = spread([ 0.0_pReal,-1.0_pReal, 0.0_pReal],2,nElems)
-  cellEdgeNormal(1:3,5,1,:) = spread([ 0.0_pReal, 0.0_pReal,+1.0_pReal],2,nElems)
-  cellEdgeNormal(1:3,6,1,:) = spread([ 0.0_pReal, 0.0_pReal,-1.0_pReal],2,nElems)
-  
-end function cellEdgeNormal
+  real(pReal), dimension(3,6,1,nElems) :: cellSurfaceNormal
+
+  cellSurfaceNormal(1:3,1,1,:) = spread([+1.0_pReal, 0.0_pReal, 0.0_pReal],2,nElems)
+  cellSurfaceNormal(1:3,2,1,:) = spread([-1.0_pReal, 0.0_pReal, 0.0_pReal],2,nElems)
+  cellSurfaceNormal(1:3,3,1,:) = spread([ 0.0_pReal,+1.0_pReal, 0.0_pReal],2,nElems)
+  cellSurfaceNormal(1:3,4,1,:) = spread([ 0.0_pReal,-1.0_pReal, 0.0_pReal],2,nElems)
+  cellSurfaceNormal(1:3,5,1,:) = spread([ 0.0_pReal, 0.0_pReal,+1.0_pReal],2,nElems)
+  cellSurfaceNormal(1:3,6,1,:) = spread([ 0.0_pReal, 0.0_pReal,-1.0_pReal],2,nElems)
+
+end function cellSurfaceNormal
 
 
 !--------------------------------------------------------------------------------------------------
@@ -385,9 +385,9 @@ end function cellEdgeNormal
 pure function IPneighborhood(grid)
 
   integer, dimension(3), intent(in) :: grid                                                         ! grid (for this process!)
-  
-  integer, dimension(3,6,1,product(grid)) :: IPneighborhood                                         !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me]
-  
+
+  integer, dimension(3,6,1,product(grid)) :: IPneighborhood                                         !< 6 neighboring IPs as [element ID, IP ID, face ID]
+
   integer :: &
    x,y,z, &
    e
@@ -432,7 +432,7 @@ pure function IPneighborhood(grid)
     IPneighborhood(3,6,1,e) = 5
 
   enddo; enddo; enddo
- 
+
 end function IPneighborhood
 
 
diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90
index 05bb77b8a..14e001eaf 100644
--- a/src/mesh_marc.f90
+++ b/src/mesh_marc.f90
@@ -32,17 +32,12 @@ module mesh
   real(pReal), public, protected :: &
     mesh_unitlength                                                                                 !< physical length of one unit in mesh
  
-  integer, dimension(:,:), allocatable, target :: &
-    mesh_mapFEtoCPelem, &                                                                           !< [sorted FEid,   corresponding CPid]
-    mesh_mapFEtoCPnode                                                                              !< [sorted FEnode, corresponding CPnode]
-
-  integer, dimension(:),   allocatable :: &
-    mapMarc2DAMASK_elem, &                                                                          !< DAMASK element ID for Marc element ID
-    mapMarc2DAMASK_node                                                                             !< DAMASK node ID for Marc node ID
+  integer, dimension(:),   allocatable, public :: &
+    mesh_FEM2DAMASK_elem, &                                                                          !< DAMASK element ID for Marc element ID
+    mesh_FEM2DAMASK_node                                                                             !< DAMASK node ID for Marc node ID
 
   public :: &
-    mesh_init, &
-    mesh_FEasCP
+    mesh_init
  
 contains
 
@@ -89,7 +84,7 @@ subroutine mesh_init(ip,el)
   FEsolving_execIP   = [1,elem%nIPs]
  
   allocate(calcMode(elem%nIPs,nElems),source=.false.)                                               ! pretend to have collected what first call is asking (F = I)
-  calcMode(ip,mesh_FEasCP('elem',el)) = .true.                                                      ! first ip,el needs to be already pingponged to "calc"
+  calcMode(ip,mesh_FEM2DAMASK_elem(el)) = .true.                                                    ! first ip,el needs to be already pingponged to "calc"
  
 
   allocate(cellNodeDefinition(elem%nNodes-1))
@@ -215,13 +210,11 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni
   call inputRead_elemType(elem, &
                           nElems,inputFile)
 
-  allocate (mesh_mapFEtoCPelem(2,nElems), source=0)
-  call inputRead_mapElems(mapMarc2DAMASK_elem,&
-                          elem%nNodes,inputFile)
+  call inputRead_mapElems(mesh_FEM2DAMASK_elem,&
+                          nElems,elem%nNodes,inputFile)
 
-  allocate (mesh_mapFEtoCPnode(2,Nnodes), source=0)
-  call inputRead_mapNodes(mapMarc2DAMASK_node,&
-                          inputFile)
+  call inputRead_mapNodes(mesh_FEM2DAMASK_node,&
+                          nNodes,inputFile)
 
   call inputRead_elemNodes(node0_elem, &
                            Nnodes,inputFile)
@@ -432,14 +425,16 @@ end subroutine inputRead_mapElemSets
 !--------------------------------------------------------------------------------------------------
 !> @brief Maps elements from FE ID to internal (consecutive) representation.
 !--------------------------------------------------------------------------------------------------
-subroutine inputRead_mapElems(map, &
-                              nNodes,fileContent)
+subroutine inputRead_mapElems(FEM2DAMASK, &
+                              nElems,nNodesPerElem,fileContent)
  
-  integer, allocatable, dimension(:), intent(out) :: map
+  integer, allocatable, dimension(:), intent(out) :: FEM2DAMASK
 
-  integer,                            intent(in)  :: nNodes                                         !< number of nodes per element
+  integer,                            intent(in)  :: nElems, &                                      !< number of elements
+                                                     nNodesPerElem                                  !< number of nodes per element
   character(len=*),     dimension(:), intent(in)  :: fileContent                                    !< file content, separated per lines
 
+  integer, dimension(2,nElems)       :: map_unsorted
   integer, allocatable, dimension(:) :: chunkPos
   integer :: i,j,l,nNodesAlreadyRead
 
@@ -448,11 +443,11 @@ subroutine inputRead_mapElems(map, &
     if(chunkPos(1) < 1) cycle
     if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'connectivity') then
       j = 0
-      do i = 1,size(mesh_mapFEtoCPelem,2)
+      do i = 1,nElems
         chunkPos = IO_stringPos(fileContent(l+1+i+j))
-        mesh_mapFEtoCPelem(:,i) = [IO_intValue(fileContent(l+1+i+j),chunkPos,1),i]
+        map_unsorted(:,i) = [IO_intValue(fileContent(l+1+i+j),chunkPos,1),i]
         nNodesAlreadyRead = chunkPos(1) - 2
-        do while(nNodesAlreadyRead < nNodes)                                                        ! read on if not all nodes in one line
+        do while(nNodesAlreadyRead < nNodesPerElem)                                                 ! read on if not all nodes in one line
           j = j + 1
           chunkPos = IO_stringPos(fileContent(l+1+i+j))
           nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1)
@@ -462,10 +457,10 @@ subroutine inputRead_mapElems(map, &
     endif
   enddo
 
-  call math_sort(mesh_mapFEtoCPelem)
-  allocate(map(minval(mesh_mapFEtoCPelem(1,:)):maxval(mesh_mapFEtoCPelem(1,:))),source=-1)
-  do i = 1,size(mesh_mapFEtoCPelem,2)
-     map(mesh_mapFEtoCPelem(1,i)) = mesh_mapFEtoCPelem(2,i)
+  call math_sort(map_unsorted)
+  allocate(FEM2DAMASK(minval(map_unsorted(1,:)):maxval(map_unsorted(1,:))),source=-1)
+  do i = 1, nElems
+    FEM2DAMASK(map_unsorted(1,i)) = map_unsorted(2,i)
   enddo
 
 end subroutine inputRead_mapElems
@@ -474,13 +469,15 @@ end subroutine inputRead_mapElems
 !--------------------------------------------------------------------------------------------------
 !> @brief Maps node from FE ID to internal (consecutive) representation.
 !--------------------------------------------------------------------------------------------------
-subroutine inputRead_mapNodes(map, &
-                              fileContent)
+subroutine inputRead_mapNodes(FEM2DAMASK, &
+                              nNodes,fileContent)
 
-  integer, allocatable, dimension(:), intent(out) :: map
+  integer, allocatable, dimension(:), intent(out) :: FEM2DAMASK
 
+  integer,                            intent(in)  :: nNodes                                         !< number of nodes
   character(len=*),     dimension(:), intent(in)  :: fileContent                                    !< file content, separated per lines
 
+  integer, dimension(2,nNodes)       :: map_unsorted
   integer, allocatable, dimension(:) :: chunkPos
   integer :: i, l
 
@@ -488,18 +485,18 @@ subroutine inputRead_mapNodes(map, &
     chunkPos = IO_stringPos(fileContent(l))
     if(chunkPos(1) < 1) cycle
     if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates') then
-      do i = 1,size(mesh_mapFEtoCPnode,2)
+      do i = 1,nNodes
         chunkPos = IO_stringPos(fileContent(l+1+i))
-        mesh_mapFEtoCPnode(:,i) = [IO_intValue(fileContent(l+1+i),chunkPos,1),i]
+        map_unsorted(:,i) = [IO_intValue(fileContent(l+1+i),chunkPos,1),i]
       enddo
       exit
     endif
   enddo
 
-  call math_sort(mesh_mapFEtoCPnode)
-  allocate(map(minval(mesh_mapFEtoCPnode(1,:)):maxval(mesh_mapFEtoCPnode(1,:))),source=-1)
-  do i = 1,size(mesh_mapFEtoCPnode,2)
-     map(mesh_mapFEtoCPnode(1,i)) = mesh_mapFEtoCPnode(2,i)
+  call math_sort(map_unsorted)
+  allocate(FEM2DAMASK(minval(map_unsorted(1,:)):maxval(map_unsorted(1,:))),source=-1)
+  do i = 1, nNodes
+    FEM2DAMASK(map_unsorted(1,i)) = map_unsorted(2,i)
   enddo
 
 end subroutine inputRead_mapNodes
@@ -526,7 +523,7 @@ subroutine inputRead_elemNodes(nodes, &
     if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates') then
       do i=1,nNode
         chunkPos = IO_stringPos(fileContent(l+1+i))
-        m = mesh_FEasCP('node',IO_intValue(fileContent(l+1+i),chunkPos,1))
+        m = mesh_FEM2DAMASK_node(IO_intValue(fileContent(l+1+i),chunkPos,1))
         do j = 1,3
           nodes(j,m) = mesh_unitlength * IO_floatValue(fileContent(l+1+i),chunkPos,j+1)
         enddo
@@ -650,11 +647,11 @@ function inputRead_connectivityElem(nElem,nNodes,fileContent)
       j = 0
       do i = 1,nElem
         chunkPos = IO_stringPos(fileContent(l+1+i+j))
-        e = mesh_FEasCP('elem',IO_intValue(fileContent(l+1+i+j),chunkPos,1))
+        e = mesh_FEM2DAMASK_elem(IO_intValue(fileContent(l+1+i+j),chunkPos,1))
         if (e /= 0) then                                                                            ! disregard non CP elems
           do k = 1,chunkPos(1)-2
             inputRead_connectivityElem(k,e) = &
-              mesh_FEasCP('node',IO_IntValue(fileContent(l+1+i+j),chunkPos,k+2))
+              mesh_FEM2DAMASK_node(IO_IntValue(fileContent(l+1+i+j),chunkPos,k+2))
           enddo
           nNodesAlreadyRead = chunkPos(1) - 2
           do while(nNodesAlreadyRead < nNodes)                                                      ! read on if not all nodes in one line
@@ -662,7 +659,7 @@ function inputRead_connectivityElem(nElem,nNodes,fileContent)
             chunkPos = IO_stringPos(fileContent(l+1+i+j))
             do k = 1,chunkPos(1)
               inputRead_connectivityElem(nNodesAlreadyRead+k,e) = &
-                mesh_FEasCP('node',IO_IntValue(fileContent(l+1+i+j),chunkPos,k))
+                mesh_FEM2DAMASK_node(IO_IntValue(fileContent(l+1+i+j),chunkPos,k))
             enddo
             nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1)
           enddo
@@ -717,7 +714,7 @@ subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogeniza
           if (initialcondTableStyle == 2) m = m + 2
           contInts = continuousIntValues(fileContent(l+k+m+1:),nElem,nameElemSet,mapElemSet,size(nameElemSet)) ! get affected elements
           do i = 1,contInts(1)
-            e = mesh_FEasCP('elem',contInts(1+i))
+            e = mesh_FEM2DAMASK_elem(contInts(1+i))
             if (sv == 2) microstructureAt(e) = myVal
             if (sv == 3) homogenizationAt(e) = myVal
           enddo
@@ -1050,53 +1047,6 @@ function IPareaNormal(elem,nElem,connectivity,node)
 end function IPareaNormal
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief Gives the FE to CP ID mapping by binary search through lookup array
-!! valid questions (what) are 'elem', 'node'
-!--------------------------------------------------------------------------------------------------
-integer function mesh_FEasCP(what,myID)
-
-  character(len=*), intent(in) :: what
-  integer,          intent(in) :: myID
- 
-  integer, dimension(:,:), pointer :: lookupMap
-  integer :: lower,upper,center
- 
-  mesh_FEasCP = 0
-  select case(IO_lc(what(1:4)))
-    case('elem')
-      lookupMap => mesh_mapFEtoCPelem
-    case('node')
-      lookupMap => mesh_mapFEtoCPnode
-    case default
-      return
-  endselect
- 
-  lower = 1
-  upper = int(size(lookupMap,2),pInt)
- 
-  if (lookupMap(1,lower) == myID) then                                                              ! check at bounds QUESTION is it valid to extend bounds by 1 and just do binary search w/o init check at bounds?
-    mesh_FEasCP = lookupMap(2,lower)
-    return
-  elseif (lookupMap(1,upper) == myID) then
-    mesh_FEasCP = lookupMap(2,upper)
-    return
-  endif
-  binarySearch: do while (upper-lower > 1)
-    center = (lower+upper)/2
-    if (lookupMap(1,center) < myID) then
-      lower = center
-    elseif (lookupMap(1,center) > myID) then
-      upper = center
-    else
-      mesh_FEasCP = lookupMap(2,center)
-      exit
-    endif
-  enddo binarySearch
-
-end function mesh_FEasCP
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief return integer list corresponding to items in consecutive lines.
 !! First integer in array is counter
diff --git a/src/rotations.f90 b/src/rotations.f90
index 16e63c180..c1fb1f2cf 100644
--- a/src/rotations.f90
+++ b/src/rotations.f90
@@ -608,13 +608,13 @@ function om2ax(om) result(ax)
     ax(1:3) = [ 0.0_pReal, 0.0_pReal, 1.0_pReal ]
   else
     call dgeev('N','V',3,om_,3,Wr,Wi,devNull,3,VR,3,work,size(work,1),ierr)
-    if (ierr /= 0) call IO_error(0,ext_msg='Error in om2ax: DGEEV return not zero')
+    if (ierr /= 0) call IO_error(401,ext_msg='Error in om2ax: DGEEV return not zero')
 #if defined(__GFORTRAN__) &&  __GNUC__<9 || defined(__INTEL_COMPILER) && INTEL_COMPILER<1800 || defined(__PGI)
     i = maxloc(merge(1,0,cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal)),dim=1)
 #else
     i = findloc(cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal),.true.,dim=1) !find eigenvalue (1,0)
 #endif
-    if (i == 0) call IO_error(0,ext_msg='Error in om2ax Real: eigenvalue not found')
+    if (i == 0) call IO_error(401,ext_msg='Error in om2ax Real: eigenvalue not found')
     ax(1:3) = VR(1:3,i)
     where (                dNeq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) &
       ax(1:3) = sign(ax(1:3),-P *[om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])