diff --git a/CMakeLists.txt b/CMakeLists.txt
index 4ad6309..ec28b31 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,288 +1,288 @@
 # =============================================================================
 # file   CMakeLists.txt
 #
 # @author Till Junge <till.junge@epfl.ch>
 #
 # @date   08 Jan 2018
 #
 # @brief  Main configuration file
 #
 # @section LICENSE
 #
 # Copyright © 2018 Till Junge
 #
 # µSpectre is free software; you can redistribute it and/or
 # modify it under the terms of the GNU Lesser General Public License as
 # published by the Free Software Foundation, either version 3, or (at
 # your option) any later version.
 #
 # µSpectre is distributed in the hope that it will be useful, but
 # WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 # General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with µSpectre; see the file COPYING. If not, write to the
 # Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 # Boston, MA 02111-1307, USA.
 #
 # Additional permission under GNU GPL version 3 section 7
 #
 # If you modify this Program, or any covered work, by linking or combining it
 # with proprietary FFT implementations or numerical libraries, containing parts
 # covered by the terms of those libraries' licenses, the licensors of this
 # Program grant you additional permission to convey the resulting work.
 # =============================================================================
 
 cmake_minimum_required(VERSION 3.0.0)
 
 project(µSpectre)
 
 set(CMAKE_CXX_STANDARD 14)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 set(BUILD_SHARED_LIBS ON)
 set(MUSPECTRE_PYTHON_MAJOR_VERSION 3)
 
 add_compile_options(-Wall -Wextra -Weffc++)
 
 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}/cmake)
 set(MAKE_DOC_TARGET "OFF" CACHE BOOL "If on, a target dev_doc (which builds the documentation) is added")
 set(MAKE_TESTS "ON" CACHE BOOL "If on, several ctest targets will be built automatically")
 set(MAKE_EXAMPLES "ON" CACHE BOOL "If on, the executables in the bin folder will be compiled")
 set(MAKE_BENCHMARKS "ON" CACHE BOOL "If on, the benchmarks will be compiled")
 set(MPI_PARALLEL "OFF" CACHE BOOL "If on, MPI-parallel solvers become available")
 set(RUNNING_IN_CI "OFF" CACHE INTERNAL "changes output format for tests")
 
 if(${MAKE_TESTS})
   enable_testing()
   find_package(Boost COMPONENTS unit_test_framework REQUIRED)
 endif(${MAKE_TESTS})
 
 if(${MPI_PARALLEL})
   add_definitions(-DWITH_MPI)
   find_package(MPI)
   if (NOT ${MPI_FOUND})
     message(SEND_ERROR "You chose MPI but CMake cannot find the MPI package")
   endif(NOT ${MPI_FOUND})
 endif(${MPI_PARALLEL})
 
 include(muspectreTools)
 include(cpplint)
 
 string( TOLOWER "${CMAKE_BUILD_TYPE}" build_type )
 if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang" OR "${CMAKE_CXX_COMPILER_ID}" STREQUAL "AppleClang")
   # using Clang
   add_compile_options(-Wno-missing-braces)
   if ("debug" STREQUAL "${build_type}")
     add_compile_options(-O0)
   endif()
 elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
   # using GCC
   add_compile_options(-Wno-non-virtual-dtor)
   add_compile_options(-march=native)
   if (("relwithdebinfo" STREQUAL "${build_type}") OR ("release" STREQUAL "${build_type}" ))
     add_compile_options(-march=native)
   endif()
   if ("debug" STREQUAL "${build_type}" )
     add_compile_options(-O0)
   endif()
 elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
   # using Intel C++
 elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "MSVC")
   # using Visual Studio C++
 endif()
 
 # Do not trust old gcc. the std::optional has memory bugs
 if(${CMAKE_COMPILER_IS_GNUCC})
   if(${CMAKE_CXX_COMPILER_VERSION} VERSION_LESS 6.0.0)
     add_definitions(-DNO_EXPERIMENTAL)
   endif()
 endif()
 
 add_external_package(Eigen3 VERSION 3.3.0 CONFIG)
+
 add_external_package(pybind11 VERSION 2.2 CONFIG)
 find_package(PythonLibsNew ${MUSPECTRE_PYTHON_MAJOR_VERSION} MODULE REQUIRED)
 
 include_directories(
   ${CMAKE_SOURCE_DIR}/src
   ${CMAKE_SOURCE_DIR}
   )
 
 if(APPLE)
   include_directories(${CMAKE_INSTALL_PREFIX}/include ${Boost_INCLUDE_DIRS})
 endif()
 
-
 #build tests (these are before we add -Werror to the compile options)
 if (${MAKE_TESTS})
   ##############################################################################
 
   # build library tests
   file( GLOB TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/test_*.cc")
 
   add_executable(main_test_suite tests/main_test_suite.cc ${TEST_SRCS})
   target_link_libraries(main_test_suite ${Boost_LIBRARIES} muSpectre)
 
 
   muSpectre_add_test(main_test_suite TYPE BOOST main_test_suite --report_level=detailed)
 
   add_executable(mattb_test tests/main_test_suite tests/test_materials_toolbox.cc)
   target_link_libraries(mattb_test ${Boost_LIBRARIES} muSpectre)
   muSpectre_add_test(mattb_test TYPE BOOST mattb_test --report_level=detailed)
 
   # build header tests
   file( GLOB HEADER_TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/header_test_*.cc")
   foreach(header_test ${HEADER_TEST_SRCS})
     get_filename_component(header_test_name ${header_test} NAME_WE)
     string(SUBSTRING ${header_test_name} 12 -1 test_name)
     list(APPEND header_tests ${test_name})
     add_executable(${test_name} tests/main_test_suite.cc ${header_test})
     target_link_libraries(${test_name} ${Boost_LIBRARIES} Eigen3::Eigen)
     target_include_directories(${test_name} INTERFACE ${muSpectre_INCLUDES})
     muSpectre_add_test(${test_name} TYPE BOOST ${test_name}
       --report_level=detailed)
   endforeach(header_test ${HEADER_TEST_SRCS})
 
   add_custom_target(header_tests)
   add_dependencies(header_tests ${header_tests})
 
   ##############################################################################
   # build py_comparison tests
   file (GLOB PY_COMP_TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/py_comparison_*.cc")
 
   find_package(PythonInterp ${MUSPECTRE_PYTHON_MAJOR_VERSION} REQUIRED)
   foreach(py_comp_test ${PY_COMP_TEST_SRCS})
     get_filename_component(py_comp_test_fname ${py_comp_test} NAME_WE)
 
     string (SUBSTRING ${py_comp_test_fname} 19 -1 py_comp_test_name)
 
     pybind11_add_module(${py_comp_test_name} ${py_comp_test})
     target_include_directories(${py_comp_test_name} PUBLIC
       ${PYTHON_INCLUDE_DIRS})
     target_link_libraries(${py_comp_test_name} PRIVATE muSpectre)
     configure_file(
       tests/${py_comp_test_fname}.py
       "${CMAKE_BINARY_DIR}/${py_comp_test_fname}.py"
       COPYONLY)
     muSpectre_add_test(${py_comp_test_fname}
       TYPE PYTHON ${py_comp_test_fname}.py)
 
   endforeach(py_comp_test ${PY_COMP_TEST_SRCS})
 
   ##############################################################################
   # copy python test
   file( GLOB PY_TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/python_*.py")
 
   foreach(pytest ${PY_TEST_SRCS})
     get_filename_component(pytest_name ${pytest} NAME)
     configure_file(
       ${pytest}
       "${CMAKE_BINARY_DIR}/${pytest_name}"
       COPYONLY)
   endforeach(pytest ${PY_TEST_SRCS})
   muSpectre_add_test(python_binding_test  TYPE PYTHON python_binding_tests.py)
 
   if(${MPI_PARALLEL})
     ############################################################################
     # add MPI tests
     file( GLOB TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/mpi_test_*.cc")
 
     add_executable(mpi_main_test_suite tests/mpi_main_test_suite.cc ${TEST_SRCS})
     target_link_libraries(mpi_main_test_suite ${Boost_LIBRARIES} muSpectre)
 
     muSpectre_add_test(mpi_main_test_suite1 TYPE BOOST MPI_NB_PROCS 1
       mpi_main_test_suite --report_level=detailed)
     muSpectre_add_test(mpi_main_test_suite2 TYPE BOOST MPI_NB_PROCS 2
       mpi_main_test_suite --report_level=detailed)
 
     muSpectre_add_test(python_mpi_binding_test1 TYPE PYTHON MPI_NB_PROCS 1
       python_mpi_binding_tests.py)
     muSpectre_add_test(python_mpi_binding_test2 TYPE PYTHON MPI_NB_PROCS 2
       python_mpi_binding_tests.py)
   endif(${MPI_PARALLEL})
 endif(${MAKE_TESTS})
 
 ################################################################################
 # compile the library
 add_compile_options( -Werror)
 add_subdirectory(
   ${CMAKE_SOURCE_DIR}/src/
   )
 add_subdirectory(
   ${CMAKE_SOURCE_DIR}/language_bindings/
   )
 
 if (${MAKE_DOC_TARGET})
   add_subdirectory(
     ${CMAKE_SOURCE_DIR}/doc/
     )
 endif()
 
 ################################################################################
 if (${MAKE_EXAMPLES})
   #compile executables
 
   set(binaries
     ${CMAKE_SOURCE_DIR}/bin/demonstrator1.cc
     ${CMAKE_SOURCE_DIR}/bin/demonstrator_dynamic_solve.cc
     ${CMAKE_SOURCE_DIR}/bin/demonstrator2.cc
     ${CMAKE_SOURCE_DIR}/bin/hyper-elasticity.cc
     ${CMAKE_SOURCE_DIR}/bin/small_case.cc)
 
   if (${MPI_PARALLEL})
     set (binaries
       ${binaries}
       ${CMAKE_SOURCE_DIR}/bin/demonstrator_mpi.cc
       )
   endif (${MPI_PARALLEL})
 
   foreach(binaryfile ${binaries})
     get_filename_component(binaryname  ${binaryfile} NAME_WE)
     add_executable(${binaryname} ${binaryfile})
     target_link_libraries(${binaryname} ${Boost_LIBRARIES} muSpectre)
   endforeach(binaryfile ${binaries})
 
   # or copy them
   file (GLOB pybins "${CMAKE_SOURCE_DIR}/bin/*.py")
   foreach(pybin ${pybins})
 
     get_filename_component(binaryname  ${pybin} NAME_WE)
     configure_file(
       ${pybin}
       "${CMAKE_BINARY_DIR}/${binaryname}.py"
       COPYONLY)
   endforeach(pybin ${pybins})
 
   # additional files to copy
   set (FILES_FOR_COPY
     ${CMAKE_SOURCE_DIR}/bin/odd_image.npz)
 
   foreach(FILE_FOR_COPY ${FILES_FOR_COPY})
 
     get_filename_component(binaryname  ${FILE_FOR_COPY} NAME)
     configure_file(
       ${FILE_FOR_COPY}
       "${CMAKE_BINARY_DIR}/${binaryname}"
       COPYONLY)
   endforeach(FILE_FOR_COPY ${FILES_FOR_COPY})
 
 endif (${MAKE_EXAMPLES})
 
 
 ################################################################################
 # compile benchmarks
 if(${MAKE_BENCHMARKS})
   file(GLOB benchmarks "${CMAKE_SOURCE_DIR}/benchmarks/benchmark*cc")
   foreach(benchmark ${benchmarks})
     get_filename_component(benchmark_name ${benchmark} NAME_WE)
     add_executable(${benchmark_name} ${benchmark})
     target_link_libraries(${benchmark_name} ${BOOST_LIBRARIES} muSpectre)
   endforeach(benchmark ${benchmark})
 endif(${MAKE_BENCHMARKS})
 
 
 cpplint_add_subdirectory("${CMAKE_SOURCE_DIR}/src" "")
 cpplint_add_subdirectory("${CMAKE_SOURCE_DIR}/tests" "")
 cpplint_add_subdirectory("${CMAKE_SOURCE_DIR}/language_bindings" "")
 cpplint_add_subdirectory("${CMAKE_SOURCE_DIR}/bin" "--filter=-build/namespaces")
diff --git a/Jenkinsfile b/Jenkinsfile
index 0cafdd5..475f59f 100644
--- a/Jenkinsfile
+++ b/Jenkinsfile
@@ -1,218 +1,218 @@
 pipeline {
     parameters {string(defaultValue: '', description: 'api-token', name: 'API_TOKEN')
                 string(defaultValue: '', description: 'Token for readthedocs', name: 'RTD_TOKEN')
 	            string(defaultValue: '', description: 'buildable phid', name: 'TARGET_PHID')
 	            string(defaultValue: 'docker_debian_testing', description: 'docker file to use', name: 'DOCKERFILE')
                 string(defaultValue: '', description: 'Commit id', name: 'COMMIT_ID')
     }
 
     agent any
 
     environment {
         OMPI_MCA_plm = 'isolated'
         OMPI_MCA_btl = 'tcp,self'
     }
     options {
         disableConcurrentBuilds()
     }
     stages {
         stage ('wipe build') {
             when {
                 anyOf{
                     changeset glob: "**/*.cmake"
                     changeset glob: "**/CMakeLists.txt"
 
                 }
             }
             steps {
                 sh ' rm -rf build_*'
             }
         }
         stage ('configure') {
             parallel {
                 stage ('docker_debian_testing') {
                     agent {
                         dockerfile {
                             filename 'docker_debian_testing'
                             dir 'dockerfiles'
                         }
                     }
                     steps {
                         configure('docker_debian_testing')
                     }
                 }
                 stage ('docker_debian_stable') {
                     agent {
                         dockerfile {
                             filename 'docker_debian_stable'
                             dir 'dockerfiles'
                         }
                     }
                     steps {
                         configure('docker_debian_stable')
                     }
                 }
             }
         }
         stage ('build') {
             parallel {
                 stage ('docker_debian_testing') {
                     agent {
                         dockerfile {
                             filename 'docker_debian_testing'
                             dir 'dockerfiles'
                         }
                     }
                     steps {
                         build('docker_debian_testing')
                     }
                 }
                 stage ('docker_debian_stable') {
                     agent {
                         dockerfile {
                             filename 'docker_debian_stable'
                             dir 'dockerfiles'
                         }
                     }
                     steps {
                         build('docker_debian_stable')
                     }
                 }
             }
         }
 
 
         stage ('Warnings gcc') {
             steps {
                 warnings(consoleParsers: [[parserName: 'GNU Make + GNU C Compiler (gcc)']])
             }
         }
 
         stage ('test') {
             parallel {
                 ////////////////////////////////////////////////////
                 stage ('docker_debian_testing') {
                     agent {
                         dockerfile {
                             filename 'docker_debian_testing'
                             dir 'dockerfiles'
                         }
                     }
                     steps {
                         run_test('docker_debian_testing', 'g++')
                         run_test('docker_debian_testing', 'clang++')
                     }
                     post {
                         always {
                             collect_test_results('docker_debian_testing', 'g++')
                             collect_test_results('docker_debian_testing', 'clang++')
                         }
                     }
                 }
                 ////////////////////////////////////////////////////
                 stage ('docker_debian_stable') {
                     agent {
                         dockerfile {
                             filename 'docker_debian_stable'
                             dir 'dockerfiles'
                         }
                     }
                     steps {
                         run_test('docker_debian_stable', 'g++')
                         run_test('docker_debian_stable', 'clang++')
                     }
                     post {
                         always {
                             collect_test_results('docker_debian_stable', 'g++')
                             collect_test_results('docker_debian_stable', 'clang++')
                         }
                     }
                 }
             }
         }
     }
 
 
     post {
         always {
             createartifact()
         }
 	    success {
             send_fail_pass('pass')
             trigger_readthedocs()
         }
 
-	    failure {
+	    unsuccessful {
             send_fail_pass('fail')
         }
     }
 }
 
 
 def configure(container_name) {
     def BUILD_DIR = "build_${container_name}"
     for (CXX_COMPILER in ["g++", "clang++"]) {
         sh """
 mkdir -p ${BUILD_DIR}_${CXX_COMPILER}
 cd ${BUILD_DIR}_${CXX_COMPILER}
 CXX=${CXX_COMPILER} cmake -DCMAKE_BUILD_TYPE:STRING=Release -DRUNNING_IN_CI=ON ..
 """
     }
 }
 
 def build(container_name) {
     def BUILD_DIR = "build_${container_name}"
     for (CXX_COMPILER in ["g++", "clang++"]) {
         sh "make -C ${BUILD_DIR}_${CXX_COMPILER}"
     }
 
 }
 
 def run_test(container_name, cxx_compiler) {
     def BUILD_DIR = "build_${container_name}"
     sh "cd ${BUILD_DIR}_${cxx_compiler} && ctest || true"
 }
 
 def send_fail_pass(state) {
     sh """
 set +x
 curl https://c4science.ch/api/harbormaster.sendmessage \
 -d api.token=${API_TOKEN} \
 -d buildTargetPHID=${TARGET_PHID} \
 -d type=${state}
 """
 }
 
 def trigger_readthedocs() {
     sh """
 set +x
 curl -X POST \
 -d "token=${RTD_TOKEN}" \
 -d "branches=master" \
 https://readthedocs.org/api/v2/webhook/muspectre/26537/
 
 curl -X POST \
 -d "token=${RTD_TOKEN}" \
 https://readthedocs.org/api/v2/webhook/muspectre/26537/
 """
 }
 
 
 def collect_test_results(container_name, cxx_compiler) {
     def BUILD_DIR = "build_${container_name}"
     junit "${BUILD_DIR}_${cxx_compiler}/test_results*.xml"
 }
 
 def createartifact() {
     sh """ set +x
 curl https://c4science.ch/api/harbormaster.createartifact \
 -d api.token=${API_TOKEN} \
 -d buildTargetPHID=${TARGET_PHID} \
 -d artifactKey="Jenkins URI" \
 -d artifactType=uri \
 -d artifactData[uri]=${BUILD_URL} \
 -d artifactData[name]="View Jenkins result" \
 -d artifactData[ui.external]=1
 """
 
 }
diff --git a/bin/visualisation_example.py b/bin/visualisation_example.py
new file mode 100644
index 0000000..7af088d
--- /dev/null
+++ b/bin/visualisation_example.py
@@ -0,0 +1,161 @@
+#!/usr/bin/env python3
+# -*- coding:utf-8 -*-
+"""
+@file   visualisation_example.py
+
+@author Richard Leute <richard.leute@imtek.uni-freiburg.de>
+
+@date   14 Jan 2019
+
+@brief  small example of how one can use the visualisation tools:
+        gradient_integration() and vtk_export()
+
+@section LICENSE
+
+Copyright © 2019 Till Junge, Richard Leute
+
+µSpectre is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public License as
+published by the Free Software Foundation, either version 3, or (at
+your option) any later version.
+
+µSpectre is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with µSpectre; see the file COPYING. If not, write to the
+Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.
+
+Additional permission under GNU GPL version 3 section 7
+
+If you modify this Program, or any covered work, by linking or combining it
+with proprietary FFT implementations or numerical libraries, containing parts
+covered by the terms of those libraries' licenses, the licensors of this
+Program grant you additional permission to convey the resulting work.
+"""
+
+
+import numpy as np
+import sys
+import os
+sys.path.append("language_bindings/python/")
+import muSpectre as µ
+import muSpectre.gradient_integration as gi
+import muSpectre.vtk_export as vt_ex
+
+
+### Input parameters ###
+#----------------------#
+#general
+resolutions = [71, 71, 3]
+lengths     = [1.0, 1.0, 0.2]
+Nx, Ny, Nz  = resolutions
+formulation = µ.Formulation.finite_strain
+Young   = [10, 20]   #Youngs modulus for each phase
+Poisson = [0.3, 0.4] #Poissons ratio for each phase
+
+#solver
+newton_tol = 1e-6 #tolerance for newton algo
+cg_tol     = 1e-6 #tolerance for cg algo
+equil_tol  = 1e-6 #tolerance for equilibrium
+maxiter = 100
+verbose = 0
+
+# sinusoidal bump
+d  = 10  #thickness of the phase with higher elasticity in pixles
+l = Nx//3   #length of bump
+h = Ny//4   #height of bump
+
+low_y  = (Ny-d)//2  #lower y-boundary of phase
+high_y = low_y+d    #upper y-boundary of phase
+
+left_x  = (Nx-l)//2  #boundaries in x direction left
+right_x = left_x + l #boundaries in x direction right
+x = np.arange(l)
+p_y = h*np.sin(np.pi/(l-1)*x) #bump function
+
+xy_bump = np.ones((l,h,Nz)) #grid representing bumpx
+for i, threshold in enumerate(np.round(p_y)):
+    xy_bump[i,int(threshold):,:] = 0
+
+phase = np.zeros(resolutions, dtype=int) #0 for surrounding matrix
+phase[:, low_y:high_y, :] = 1
+phase[left_x:right_x, high_y:high_y+h, :] = xy_bump
+
+
+### Run muSpectre ###
+#-------------------#
+cell = µ.Cell(resolutions,
+              lengths,
+              formulation)
+mat = µ.material.MaterialLinearElastic4_3d.make(cell, "material")
+
+for i, pixel in enumerate(cell):
+    #add Young and Poisson depending on the material index
+    m_i = phase.flatten()[i] #m_i = material index / phase index
+    mat.add_pixel(pixel, Young[m_i], Poisson[m_i])
+
+cell.initialise()  #initialization of fft to make faster fft
+DelF = np.array([[0  , 0.7, 0],
+                 [0  , 0  , 0],
+                 [0  , 0  , 0]])
+
+solver_newton = µ.solvers.SolverCG(cell, cg_tol, maxiter, verbose)
+result = µ.solvers.newton_cg(cell, DelF, solver_newton,
+                               newton_tol, equil_tol, verbose)
+
+# print solver results
+print('\n\nstatus messages of OptimizeResults:')
+print('-----------------------------------')
+print('success:            ', result.success)
+print('status:             ', result.status)
+print('message:            ', result.message)
+print('# iterations:       ', result.nb_it)
+print('# cell evaluations: ', result.nb_fev)
+print('formulation:        ', result.formulation)
+
+
+### Visualisation ###
+#-------------------#
+#integration of the deformation gradient field
+placement_n, x = gi.compute_placement(result, lengths, resolutions, order = 0)
+
+### some fields which can be added to the visualisation
+#2-tensor field containing the first Piola Kirchhoff stress
+PK1 = gi.reshape_gradient(result.stress, resolutions)
+#scalar field containing the distance to the origin O
+distance_O = np.linalg.norm(placement_n, axis=-1)
+
+#random fields
+center_shape = tuple(np.array(x.shape[:-1])-1)#shape of the center point grid
+dim = len(resolutions)
+scalar_c  = np.random.random(center_shape)
+vector_c  = np.random.random(center_shape + (dim,))
+tensor2_c = np.random.random(center_shape + (dim,dim))
+scalar_n  = np.random.random(x.shape[:-1])
+vector_n  = np.random.random(x.shape[:-1] + (dim,))
+tensor2_n = np.random.random(x.shape[:-1] + (dim,dim))
+
+#write dictionaries with cell data and point data
+c_data  = {"scalar field"   : scalar_c,
+           "vector field"   : vector_c,
+           "2-tensor field" : tensor2_c,
+           "PK1 stress"     : PK1,
+           "phase"          : phase}
+p_data = {"scalar field"   : scalar_n,
+          "vector field"   : vector_n,
+          "2-tensor field" : tensor2_n,
+          "distance_O"     : distance_O}
+
+vt_ex.vtk_export(fpath = "visualisation_example",
+                 x_n   = x,
+                 placement  = placement_n,
+                 point_data = p_data,
+                 cell_data  = c_data)
+
+print("The file 'visualisation_example.vtr' was successfully written!\n"
+      "You can open it for example with paraview or some other software:\n"
+      "paraview visualisation_example.vtr")
diff --git a/language_bindings/python/CMakeLists.txt b/language_bindings/python/CMakeLists.txt
index fb55a13..b074634 100644
--- a/language_bindings/python/CMakeLists.txt
+++ b/language_bindings/python/CMakeLists.txt
@@ -1,93 +1,93 @@
 #==============================================================================
 # file   CMakeLists.txt
 #
 # @author Till Junge <till.junge@epfl.ch>
 #
 # @date   08 Jan 2018
 #
 # @brief  configuration for python binding using pybind11
 #
 # @section LICENSE
 #
 # Copyright © 2018 Till Junge
 #
 # µSpectre is free software; you can redistribute it and/or
 # modify it under the terms of the GNU Lesser General Public License as
 # published by the Free Software Foundation, either version 3, or (at
 # your option) any later version.
 #
 # µSpectre is distributed in the hope that it will be useful, but
 # WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 # General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with µSpectre; see the file COPYING. If not, write to the
 # Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 # Boston, MA 02111-1307, USA.
 #
 # Additional permission under GNU GPL version 3 section 7
 #
 # If you modify this Program, or any covered work, by linking or combining it
 # with proprietary FFT implementations or numerical libraries, containing parts
 # covered by the terms of those libraries' licenses, the licensors of this
 # Program grant you additional permission to convey the resulting work.
 # =============================================================================
 
 # FIXME! The user should have a choice to configure this path.
 execute_process(COMMAND "${PYTHON_EXECUTABLE}" "-m" "site" "--user-site"
   RESULT_VARIABLE _PYTHON_SUCCESS
   OUTPUT_VARIABLE PYTHON_USER_SITE
   ERROR_VARIABLE _PYTHON_ERROR_VALUE)
 if(NOT _PYTHON_SUCCESS MATCHES 0)
   message(FATAL_ERROR
           "Python config failure:\n${_PYTHON_ERROR_VALUE}")
 endif()
 string(REGEX REPLACE "\n" "" PYTHON_USER_SITE ${PYTHON_USER_SITE})
 
 set (PY_BINDING_SRCS
   ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_module.cc
   ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_common.cc
   ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_cell.cc
   ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_material.cc
   ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_material_linear_elastic1.cc
   ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_material_linear_elastic2.cc
   ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_material_linear_elastic3.cc
   ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_material_linear_elastic4.cc
   ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_material_hyper_elasto_plastic1.cc
   ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_material_linear_elastic_generic.cc
   ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_solvers.cc
   ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_fftengine.cc
   ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_projections.cc
   ${CMAKE_CURRENT_SOURCE_DIR}/bind_py_field_collection.cc
   )
 
 if (${USE_FFTWMPI})
   add_definitions(-DWITH_FFTWMPI)
 endif(${USE_FFTWMPI})
 
 if (${USE_PFFT})
   add_definitions(-DWITH_PFFT)
 endif(${USE_PFFT})
 
 find_package(PythonLibsNew ${MUSPECTRE_PYTHON_MAJOR_VERSION} MODULE REQUIRED)
 
 # On OS X, add -Wno-deprecated-declarations
 IF(${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
   add_compile_options(-Wno-deprecated-declarations)
 ENDIF(${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
 
 pybind11_add_module(pyMuSpectreLib ${PY_BINDING_SRCS})
 target_link_libraries(pyMuSpectreLib PRIVATE muSpectre)
 # Want to rename the output, so that the python module is called muSpectre
 set_target_properties(pyMuSpectreLib PROPERTIES OUTPUT_NAME _muSpectre)
 
 target_include_directories(pyMuSpectreLib PUBLIC ${PYTHON_INCLUDE_DIRS})
 
-add_custom_target(pyMuSpectre ALL SOURCES muSpectre/__init__.py muSpectre/fft.py)
+add_custom_target(pyMuSpectre ALL SOURCES muSpectre/__init__.py muSpectre/fft.py muSpectre/tools.py)
 add_custom_command(TARGET pyMuSpectre POST_BUILD
                    COMMAND ${CMAKE_COMMAND} -E copy_directory
                        ${CMAKE_SOURCE_DIR}/language_bindings/python/muSpectre $<TARGET_FILE_DIR:pyMuSpectreLib>/muSpectre)
 
 install(TARGETS pyMuSpectreLib LIBRARY DESTINATION ${PYTHON_USER_SITE})
 install(FILES muSpectre/__init__.py muSpectre/fft.py DESTINATION ${PYTHON_USER_SITE}/muSpectre)
diff --git a/language_bindings/python/bind_py_solvers.cc b/language_bindings/python/bind_py_solvers.cc
index 9a30c54..40be085 100644
--- a/language_bindings/python/bind_py_solvers.cc
+++ b/language_bindings/python/bind_py_solvers.cc
@@ -1,146 +1,147 @@
 /**
  * @file   bind_py_solver.cc
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   09 Jan 2018
  *
  * @brief  python bindings for the muSpectre solvers
  *
  * Copyright © 2018 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU Lesser General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with µSpectre; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * * Boston, MA 02111-1307, USA.
  *
  * Additional permission under GNU GPL version 3 section 7
  *
  * If you modify this Program, or any covered work, by linking or combining it
  * with proprietary FFT implementations or numerical libraries, containing parts
  * covered by the terms of those libraries' licenses, the licensors of this
  * Program grant you additional permission to convey the resulting work.
  */
 
 #include "common/common.hh"
 #include "solver/solvers.hh"
 #include "solver/solver_cg.hh"
 #include "solver/solver_eigen.hh"
 
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 #include <pybind11/eigen.h>
 
 using muSpectre::Dim_t;
 using muSpectre::OptimizeResult;
 using muSpectre::Real;
 using muSpectre::Uint;
 using pybind11::literals::operator""_a;
 namespace py = pybind11;
 
 /**
  * Solvers instanciated for cells with equal spatial and material dimension
  */
 
 template <class Solver>
 void add_iterative_solver_helper(py::module & mod, std::string name) {
   py::class_<Solver, typename Solver::Parent>(mod, name.c_str())
       .def(py::init<muSpectre::Cell &, Real, Uint, bool>(), "cell"_a, "tol"_a,
            "maxiter"_a, "verbose"_a = false)
       .def("name", &Solver::get_name);
 }
 
 void add_iterative_solver(py::module & mod) {
   std::stringstream name{};
   name << "SolverBase";
   py::class_<muSpectre::SolverBase>(mod, name.str().c_str());
   add_iterative_solver_helper<muSpectre::SolverCG>(mod, "SolverCG");
   add_iterative_solver_helper<muSpectre::SolverCGEigen>(mod, "SolverCGEigen");
   add_iterative_solver_helper<muSpectre::SolverGMRESEigen>(mod,
                                                            "SolverGMRESEigen");
   add_iterative_solver_helper<muSpectre::SolverBiCGSTABEigen>(
       mod, "SolverBiCGSTABEigen");
   add_iterative_solver_helper<muSpectre::SolverDGMRESEigen>(
       mod, "SolverDGMRESEigen");
   add_iterative_solver_helper<muSpectre::SolverMINRESEigen>(
       mod, "SolverMINRESEigen");
 }
 
 void add_newton_cg_helper(py::module & mod) {
   const char name[]{"newton_cg"};
   using solver = muSpectre::SolverBase;
   using grad = py::EigenDRef<Eigen::MatrixXd>;
   using grad_vec = muSpectre::LoadSteps_t;
 
   mod.def(name,
           [](muSpectre::Cell & s, const grad & g, solver & so, Real nt,
              Real eqt, Dim_t verb) -> OptimizeResult {
             Eigen::MatrixXd tmp{g};
             return newton_cg(s, tmp, so, nt, eqt, verb);
           },
           "cell"_a, "ΔF₀"_a, "solver"_a, "newton_tol"_a, "equil_tol"_a,
           "verbose"_a = 0);
   mod.def(name,
           [](muSpectre::Cell & s, const grad_vec & g, solver & so, Real nt,
              Real eqt, Dim_t verb) -> std::vector<OptimizeResult> {
             return newton_cg(s, g, so, nt, eqt, verb);
           },
           "cell"_a, "ΔF₀"_a, "solver"_a, "newton_tol"_a, "equilibrium_tol"_a,
           "verbose"_a = 0);
 }
 
 void add_de_geus_helper(py::module & mod) {
   const char name[]{"de_geus"};
   using solver = muSpectre::SolverBase;
   using grad = py::EigenDRef<Eigen::MatrixXd>;
   using grad_vec = muSpectre::LoadSteps_t;
 
   mod.def(name,
           [](muSpectre::Cell & s, const grad & g, solver & so, Real nt,
              Real eqt, Dim_t verb) -> OptimizeResult {
             Eigen::MatrixXd tmp{g};
             return de_geus(s, tmp, so, nt, eqt, verb);
           },
           "cell"_a, "ΔF₀"_a, "solver"_a, "newton_tol"_a, "equilibrium_tol"_a,
           "verbose"_a = 0);
   mod.def(name,
           [](muSpectre::Cell & s, const grad_vec & g, solver & so, Real nt,
              Real eqt, Dim_t verb) -> std::vector<OptimizeResult> {
             return de_geus(s, g, so, nt, eqt, verb);
           },
           "cell"_a, "ΔF₀"_a, "solver"_a, "newton_tol"_a, "equilibrium_tol"_a,
           "verbose"_a = 0);
 }
 
 void add_solver_helper(py::module & mod) {
   add_newton_cg_helper(mod);
   add_de_geus_helper(mod);
 }
 
 void add_solvers(py::module & mod) {
   auto solvers{mod.def_submodule("solvers")};
   solvers.doc() = "bindings for solvers";
 
   py::class_<OptimizeResult>(mod, "OptimizeResult")
       .def_readwrite("grad", &OptimizeResult::grad)
       .def_readwrite("stress", &OptimizeResult::stress)
       .def_readwrite("success", &OptimizeResult::success)
       .def_readwrite("status", &OptimizeResult::status)
       .def_readwrite("message", &OptimizeResult::message)
       .def_readwrite("nb_it", &OptimizeResult::nb_it)
-      .def_readwrite("nb_fev", &OptimizeResult::nb_fev);
+      .def_readwrite("nb_fev", &OptimizeResult::nb_fev)
+      .def_readwrite("formulation", &OptimizeResult::formulation);
 
   add_iterative_solver(solvers);
 
   add_solver_helper(solvers);
 }
diff --git a/language_bindings/python/muSpectre/__init__.py b/language_bindings/python/muSpectre/__init__.py
index e6d39bc..b778d54 100644
--- a/language_bindings/python/muSpectre/__init__.py
+++ b/language_bindings/python/muSpectre/__init__.py
@@ -1,104 +1,105 @@
 #!/usr/bin/env python3
 # -*- coding:utf-8 -*-
 """
 @file   __init__.py
 
 @author Lars Pastewka <lars.pastewka@imtek.uni-freiburg.de>
 
 @date   21 Mar 2018
 
 @brief  Main entry point for muSpectre Python module
 
 Copyright © 2018 Till Junge
 
 µSpectre is free software; you can redistribute it and/or
 modify it under the terms of the GNU General Lesser Public License as
 published by the Free Software Foundation, either version 3, or (at
 your option) any later version.
 
 µSpectre is distributed in the hope that it will be useful, but
 WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 General Public License for more details.
 
 You should have received a copy of the GNU Lesser General Public License
 along with µSpectre; see the file COPYING. If not, write to the
 Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 Boston, MA 02111-1307, USA.
 
 Additional permission under GNU GPL version 3 section 7
 
 If you modify this Program, or any covered work, by linking or combining it
 with proprietary FFT implementations or numerical libraries, containing parts
 covered by the terms of those libraries' licenses, the licensors of this
 Program grant you additional permission to convey the resulting work.
 """
 
 
 try:
     from mpi4py import MPI
 except ImportError:
     MPI = None
 
 import _muSpectre
 from _muSpectre import (Formulation, get_domain_ccoord, get_domain_index,
                         get_hermitian_sizes, material, solvers,
                         FFT_PlanFlags, FiniteDiff)
 import muSpectre.fft
+import muSpectre.gradient_integration
 
 _factories = {'fftw': ('CellFactory', False),
               'fftwmpi': ('FFTWMPICellFactory', True),
               'pfft': ('PFFTCellFactory', True),
               'p3dfft': ('P3DFFTCellFactory', True)}
 
 
 def Cell(resolutions, lengths, formulation=Formulation.finite_strain,
          fft='fftw', communicator=None):
     """
     Instantiate a muSpectre Cell class.
 
     Parameters
     ----------
     resolutions: list
         Grid resolutions in the Cartesian directions.
     lengths: list
         Physical size of the cell in the Cartesian directions.
     formulation: Formulation
         Formulation for strains and stresses used by the solver. Options are
         `Formulation.finite_strain` and `Formulation.small_strain`. Finite
         strain formulation is the default.
     fft: string
         FFT engine to use. Options are 'fftw', 'fftwmpi', 'pfft' and 'p3dfft'.
         Default is 'fftw'.
     communicator: mpi4py communicator
         mpi4py communicator object passed to parallel FFT engines. Note that
         the default 'fftw' engine does not support parallel execution.
 
 
     Returns
     -------
     cell: object
         Return a muSpectre Cell object.
     """
     try:
         factory_name, is_parallel = _factories[fft]
     except KeyError:
         raise KeyError("Unknown FFT engine '{}'.".format(fft))
     try:
         factory = _muSpectre.__dict__[factory_name]
     except KeyError:
         raise KeyError("FFT engine '{}' has not been compiled into the "
                        "muSpectre library.".format(fft))
     if is_parallel:
         if MPI is None:
             raise RuntimeError('Parallel solver requested but mpi4py could'
                                ' not be imported.')
         if communicator is None:
             communicator = MPI.COMM_SELF
         return factory(resolutions, lengths, formulation,
                        MPI._handleof(communicator))
     else:
         if communicator is not None:
             raise ValueError("FFT engine '{}' does not support parallel "
                              "execution.".format(fft))
         return factory(resolutions, lengths, formulation)
diff --git a/language_bindings/python/muSpectre/gradient_integration.py b/language_bindings/python/muSpectre/gradient_integration.py
new file mode 100644
index 0000000..528a5d6
--- /dev/null
+++ b/language_bindings/python/muSpectre/gradient_integration.py
@@ -0,0 +1,365 @@
+#!/usr/bin/env python3
+# -*- coding:utf-8 -*-
+"""
+@file   gradient_integration.py
+
+@author Till Junge <till.junge@epfl.ch>
+
+@date   22 Nov 2018
+
+@brief  Functions for the integration of periodic first- and second-rank
+        tensor fields on an n-dimensional rectangular grid
+
+Copyright © 2018 Till Junge
+
+µSpectre is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public License as
+published by the Free Software Foundation, either version 3, or (at
+your option) any later version.
+
+µSpectre is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with µSpectre; see the file COPYING. If not, write to the
+Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.
+
+Additional permission under GNU GPL version 3 section 7
+
+If you modify this Program, or any covered work, by linking or combining it
+with proprietary FFT implementations or numerical libraries, containing parts
+covered by the terms of those libraries' licenses, the licensors of this
+Program grant you additional permission to convey the resulting work.
+"""
+
+import numpy as np
+import sys
+
+def compute_wave_vectors(lengths, resolutions):
+    """Computes the wave vectors for a dim-dimensional rectangular or
+    cubic (or hypercubic) grid as a function of the edge lengths and
+    resolutions.
+
+    Note: the norm of the wave vectors corresponds to the angular
+    velocity, not the frequency.
+
+    Keyword Arguments:
+    lengths     -- np.ndarray of length dim with the edge lengths in each
+                   spatial dimension (dtype = float)
+    resolutions -- np.ndarray of length dim with the resolutions in each
+                   spatial dimension (dtype = int)
+
+    Returns:
+    np.ndarary of shape resolutions + [dim]. The wave vector for a
+    given pixel/voxel is given in the last dimension
+
+    """
+    return np.moveaxis(np.meshgrid(
+        *[2*np.pi*np.fft.fftfreq(r, l/r) for l,r in zip(lengths, resolutions)],
+        indexing="ij"), 0, -1)
+
+def compute_grid(lengths, resolutions):
+    """For a dim-dimensional pixel/voxel grid, computes the pixel/voxel
+    centre and corner positions as a function of the grid's edge
+    lengths and resolutions
+
+    Keyword Arguments:
+    lengths     -- np.ndarray of length dim with the edge lengths in each
+                   spatial dimension (dtype = float)
+    resolutions -- np.ndarray of length dim with the resolutions in each
+                   spatial dimension (dtype = int)
+    Returns:
+    tuple((x_n, x_c)) two ndarrays with nodal/corner positions and
+    centre positions respectively. x_n has one more entry in every
+    direction than the resolution of the grid (added points correspond
+    to the periodic repetitions)
+
+    """
+    x_n = np.moveaxis(np.meshgrid(
+        *[np.linspace(0, l, r+1) for l,r in zip(lengths, resolutions)],
+        indexing="ij"), 0 ,-1)
+    dx = lengths/resolutions
+
+    x_c = np.moveaxis(np.meshgrid(
+        *[np.linspace(0, l, r, endpoint=False) for l,r in
+          zip(lengths, resolutions)],
+        indexing="ij"), 0, -1) + .5*dx
+
+    return x_n, x_c
+
+def reshape_gradient(F, resolutions):
+    """reshapes a flattened second-rank tensor into a multidimensional array of
+    shape resolutions + [dim, dim].
+
+    Note: this reshape entails a copy, because of the column-major to
+    row-major transposition between Eigen and numpy
+
+    Keyword Arguments:
+    F           -- flattenen array of gradients as in OptimizeResult
+    resolutions -- np.ndarray of length dim with the resolutions in each
+                   spatial dimension (dtype = int)
+
+    Returns:
+    np.ndarray
+    """
+
+    dim = len(resolutions)
+    if not isinstance(resolutions, list):
+        raise Exception("resolutions needs to be in list form, "+
+                        "for concatenation")
+    expected_input_shape = [np.prod(resolutions) * dim**2]
+    output_shape = resolutions + [dim, dim]
+    if not ((F.shape[0] == expected_input_shape[0]) and
+            (F.size == expected_input_shape[0])):
+        raise Exception("expected gradient of shape {}, got {}".format(
+            expected_input_shape, F.shape))
+
+    order = list(range(dim+2))
+    order[-2:] = reversed(order[-2:])
+    return F.reshape(output_shape).transpose(*order)
+
+def complement_periodically(array, dim):
+    """Takes an arbitrary multidimensional array of at least dimension dim
+    and returns an augmented copy with periodic copies of the
+    left/lower entries in the added right/upper boundaries
+
+    Keyword Arguments:
+    array -- arbitrary np.ndarray of at least dim dimensions
+    dim   -- nb of dimension to complement periodically
+
+    """
+    shape = list(array.shape)
+    shape[:dim] = [d+1 for d in shape[:dim]]
+    out_arr = np.empty(shape, dtype = array.dtype)
+    sl = tuple([slice(0, s) for s in array.shape])
+    out_arr[sl] = array
+
+    for i in range(dim):
+        lower_slice = tuple([slice(0,s) if (d != i) else  0 for (d,s) in
+                             enumerate(shape)])
+        upper_slice = tuple([slice(0,s) if (d != i) else -1 for (d,s) in
+                             enumerate(shape)])
+        out_arr[upper_slice] = out_arr[lower_slice]
+
+    return out_arr
+
+
+def get_integrator(x, freqs, order=0):
+    """returns the discrete Fourier-space integration operator as a
+    function of the position grid (used to determine the spatial
+    dimension and resolution), the wave vectors, and the integration
+    order
+
+    Keyword Arguments:
+    x     -- np.ndarray of pixel/voxel centre positons in shape
+             resolution + [dim]
+    freqs -- wave vectors as computed by compute_wave_vectors
+    order -- (default 0) integration order. 0 stands for exact integration
+
+    Returns:
+    (dim, shape, integrator)
+    """
+    dim = x.shape[-1]
+    shape = x.shape[:-1]
+    delta_x = x[tuple([1]*dim)] - x[tuple([0]*dim)]
+    def correct_denom(denom):
+        """Corrects the denominator to avoid division by zero
+        for freqs = 0 and on even grids for freqs*delta_x=pi."""
+        denom[tuple([0]*dim)] = 1.
+        for i, n in enumerate(shape):
+            if n%2 == 0:
+                denom[tuple([np.s_[:]]*i + [n//2] + [np.s_[:]]*(dim-1-i))] = 1.
+        return denom[..., np.newaxis]
+    if order == 0:
+        freqs_norm_square = np.einsum("...i, ...i -> ...", freqs, freqs)
+        freqs_norm_square.reshape(-1)[0] = 1.
+        integrator = 1j * freqs/freqs_norm_square[...,np.newaxis]
+    # Higher order corrections after:
+    # A. Vidyasagar et al., Computer Methods in Applied Mechanics and
+    # Engineering 106 (2017) 133-151, sec. 3.4 and Appendix C.
+    # and
+    # P. Eisenlohr et al., International Journal of Plasticity 46 (2013)
+    # 37-53, Appendix B, eq. 23.
+    elif order == 1:
+        sin_1 = np.sin(freqs*delta_x)
+        denom = correct_denom((sin_1**2).sum(axis=-1))
+        integrator = 1j*delta_x*sin_1 / denom
+    elif order == 2:
+        sin_1, sin_2 = 8*np.sin(freqs*delta_x), np.sin(2*freqs*delta_x)
+        denom = correct_denom((sin_1**2).sum(axis=-1) - (sin_2**2).sum(axis=-1))
+        integrator = 1j*6*delta_x*(sin_1+sin_2) / denom
+    else:
+        raise Exception("order '{}' is not implemented".format(order))
+    return dim, shape, integrator
+
+
+def integrate_tensor_2(grad, x, freqs, staggered_grid=False, order=0):
+    """Integrates a second-rank tensor gradient field to a chosen order as
+    a function of the given field, the grid positions, and wave
+    vectors. Optionally, the integration can be performed on the
+    pixel/voxel corners (staggered grid).
+
+    Keyword Arguments:
+    grad           -- np.ndarray of shape resolution + [dim, dim] containing the
+                      second-rank gradient to be integrated
+    x              -- np.ndarray of shape resolution + [dim] (or augmented
+                      resolution + [dim]) containing the pixel/voxel centre
+                      positions (for un-staggered grid integration) or the pixel
+                      /voxel corner positions (for staggered grid integration)
+    freqs          -- wave vectors as computed by compute_wave_vectors
+    staggered_grid -- (default False) if set to True, the integration is
+                      performed on the pixel/voxel corners, rather than the
+                      centres. This leads to a different integration scheme
+    order          -- (default 0) integration order.
+                      0 stands for exact integration
+
+    Returns:
+    np.ndarray contaning the integrated field
+    """
+
+    dim, shape, integrator = get_integrator(x, freqs, order)
+
+    axes = range(dim)
+    grad_k = np.fft.fftn(grad, axes=axes)
+    f_k = np.einsum("...j, ...ij -> ...i", integrator, grad_k)
+    normalisation = np.prod(grad.shape[:dim])
+    grad_k_0 = grad_k[tuple((0 for _ in range(dim)))].real/normalisation
+    homogeneous = np.einsum("ij, ...j -> ...i",
+                            grad_k_0, x)
+    if not staggered_grid:
+        fluctuation = -np.fft.ifftn(f_k, axes=axes).real
+    else:
+        del_x = (x[tuple((1 for _ in range(dim)))] -
+                 x[tuple((0 for _ in range(dim)))])
+        k_del_x = np.einsum("...i, i ->...", freqs, del_x)[...,np.newaxis]
+        if order == 0:
+            shift = np.exp(-1j * k_del_x/2)
+        elif order == 1:
+            shift = (np.exp(-1j * k_del_x) + 1) / 2
+        elif order == 2:
+            shift = np.exp(-1j*k_del_x/2) * np.cos(k_del_x/2) *\
+                    (np.cos(k_del_x) - 4) / (np.cos(k_del_x/2) - 4)
+        fluctuation = complement_periodically(
+            -np.fft.ifftn(shift*f_k, axes=axes).real, dim)
+
+    return fluctuation + homogeneous
+
+
+def integrate_vector(df, x, freqs, staggered_grid=False, order=0):
+    """Integrates a first-rank tensor gradient field to a chosen order as
+    a function of the given field, the grid positions, and wave
+    vectors. Optionally, the integration can be performed on the
+    pixel/voxel corners (staggered_grid)
+
+    Keyword Arguments:
+    df             -- np.ndarray of shape resolution + [dim] containing the
+                      first-rank tensor gradient to be integrated
+    x              -- np.ndarray of shape resolution + [dim] (or augmented
+                      resolution + [dim]) containing the pixel/voxel centre
+                      positions (for un-staggered grid integration) or the pixel
+                      /voxel corner positions (for staggered grid integration)
+    freqs          -- wave vectors as computed by compute_wave_vectors
+    staggered_grid -- (default False) if set to True, the integration is
+                      performed on the pixel/voxel corners, rather than the
+                      centres. This leads to a different integration scheme...
+    order          -- (default 0) integration order.
+                      0 stands for exact integration
+
+    Returns:
+    np.ndarray contaning the integrated field
+    """
+    dim, shape, integrator = get_integrator(x, freqs, order)
+
+    axes = range(dim)
+    df_k = np.fft.fftn(df, axes=axes)
+    f_k = np.einsum("...i, ...i -> ...", df_k, integrator)
+    df_k_0 = df_k[tuple((0 for _ in range(dim)))].real
+    homogeneous = np.einsum("i, ...i -> ...", df_k_0, x/np.prod(shape))
+
+    if not staggered_grid:
+        fluctuation = -np.fft.ifftn(f_k, axes=axes).real
+    else:
+        del_x = x[tuple((1 for _ in range(dim)))] \
+                - x[tuple((0 for _ in range(dim)))]
+        k_del_x = np.einsum("...i, i ->...", freqs, del_x)
+        if order == 0:
+            shift = np.exp(-1j * k_del_x/2)
+        elif order == 1:
+            shift = (np.exp(-1j * k_del_x) + 1) / 2
+        elif order == 2:
+            shift = np.exp(-1j*k_del_x/2) * np.cos(k_del_x/2) *\
+                    (np.cos(k_del_x) - 4) / (np.cos(k_del_x/2) - 4)
+        fluctuation = complement_periodically(
+            -np.fft.ifftn(shift*f_k, axes=axes).real, dim)
+
+    return fluctuation + homogeneous
+
+
+def compute_placement(result, lengths, resolutions, order=0, formulation=None):
+    """computes the placement (the sum of original position and
+    displacement) as a function of a OptimizeResult, domain edge
+    lengths, domain discretisation resolutions, the chosen
+    integration order and the continuum mechanics description(small or finite
+    strain description)
+
+    Keyword Arguments:
+    result      -- OptimiseResult, or just the grad field of an OptimizeResult
+    lengths     -- np.ndarray of length dim with the edge lengths in each
+                   spatial dimension (dtype = float)
+    resolutions -- np.ndarray of length dim with the resolutions in each
+                   spatial dimension (dtype = int)
+    order       -- (default 0) integration order. 0 stands for exact integration
+    formulation -- (default None) the formulation is derived from the
+                   OptimiseResult argument. If this is not possible you have to
+                   fix the formulation to either "small_strain" or
+                   "finite_strain". (dtype = str)
+    Returns:
+    (placement, x_n) tuple of ndarrays containing the placement and the
+                     corresponding original positions
+
+    """
+
+    lengths = np.array(lengths)
+    resolutions = np.array(resolutions)
+
+    #Check whether result is a np.array or an OptimiseResult object
+    if isinstance(result, np.ndarray):
+        if formulation == None:
+            #exit the program, if the formulation is unknown!
+            raise ValueError('\n'
+                'You have to specify your continuum mechanics description.\n'
+                'Either you use a formulation="small_strain" or '
+                '"finite_strain" description.\n'
+                'Otherwise you can give a result=OptimiseResult object, which '
+                'tells me the formulation.')
+        form = formulation
+        grad = reshape_gradient(result, resolutions.tolist())
+    else:
+        form = str(result.formulation)[12:]
+        if form != formulation and formulation != None:
+            #exit the program, if the formulation is ambiguous!
+            raise ValueError('\nThe given formulation "{}" differs from the '
+                             'one saved in your result "{}"!'
+                             .format(formulation, form))
+        grad = reshape_gradient(result.grad, resolutions.tolist())
+
+    #reshape the gradient depending on the formulation
+    if form == "small_strain" or form == "small_strain_sym":
+        raise NotImplementedError('\nIntegration of small strains'
+                                  'is not implemented yet!')
+    elif form == "finite_strain":
+        grad = grad
+    else:
+        raise ValueError('\nThe formulation: "{}" is unknown!'
+                         .format(formulation))
+
+    #compute the placement by integrating
+    x_n, x_c = compute_grid(lengths, resolutions)
+    freqs = compute_wave_vectors(lengths, resolutions)
+    placement = integrate_tensor_2(grad, x_n, freqs,
+                                   staggered_grid=True, order=order)
+
+    return placement, x_n
diff --git a/language_bindings/python/muSpectre/vtk_export.py b/language_bindings/python/muSpectre/vtk_export.py
new file mode 100644
index 0000000..c771fdf
--- /dev/null
+++ b/language_bindings/python/muSpectre/vtk_export.py
@@ -0,0 +1,178 @@
+#!/usr/bin/env python3
+# -*- coding:utf-8 -*-
+"""
+@file   vtk_export.py
+
+@author Till Junge <till.junge@epfl.ch>
+
+@date   22 Nov 2018
+
+@brief  function for export of vtk files
+
+Copyright © 2018 Till Junge
+
+µSpectre is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Lesser Public License as
+published by the Free Software Foundation, either version 3, or (at
+your option) any later version.
+
+µSpectre is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with µSpectre; see the file COPYING. If not, write to the
+Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.
+
+Additional permission under GNU GPL version 3 section 7
+
+If you modify this Program, or any covered work, by linking or combining it
+with proprietary FFT implementations or numerical libraries, containing parts
+covered by the terms of those libraries' licenses, the licensors of this
+Program grant you additional permission to convey the resulting work.
+"""
+
+import numpy as np
+from tvtk.api import tvtk, write_data
+
+def vtk_export(fpath, x_n, placement, point_data=None, cell_data=None,
+               legacy_ascii=False):
+    """write a vtr (vtk for rectilinear grids) file for visualisation of
+    µSpectre results. Optionally, a legacy ascii file can be written
+    (useful for debugging)
+
+    Keyword Arguments:
+    fpath        -- file name for output path *WITHOUT EXTENSION*
+    x_n          -- nodal positions, as computed by
+                    gradient_integration.compute_grid
+    placement    -- nodal deformed placement as computed by
+                    gradient_integration.compute_placement
+    point_data   -- (default None) dictionary of point data. These must have
+                    either of the following shapes:
+                    a) tuple(x_n.shape[:dim])
+                       interpreted as scalar field,
+                    b) tuple(x_n.shape[:dim] + [dim])
+                       interpreted as vector field, or
+                    c) tuple(x_n.shape[:dim] + [dim, dim])
+                       interpreted as second-rank tensor field
+    cell_data    -- (default None) dictionary of cell data. These must have
+                    either of the following shapes:
+                    a) tuple(x_c.shape[:dim])
+                       interpreted as scalar field,
+                    b) tuple(x_c.shape[:dim] + [dim])
+                       interpreted as vector field, or
+                    c) tuple(x_c.shape[:dim] + [dim, dim])
+                       interpreted as second-rank tensor field
+                    where x_c is the array of pixel/voxel positions as computed
+                    by gradient_integration.compute_grid
+    legacy_ascii -- (default False) If set to True, a human-readable, but larger
+                     ascii file is written
+
+    """
+    dim = len(x_n.shape[:-1])
+    if dim not in (2, 3):
+        raise Exception(
+            ("should be two- or three-dimensional, got positions for a {}-"
+             "dimensional problem").format(dim))
+    res_n = list(x_n.shape[:-1])
+    vtk_res = res_n if dim == 3 else res_n + [1]
+    res_c = [max(1, r-1) for r in res_n]
+
+    # setting up the geometric grid
+    vtk_obj = tvtk.RectilinearGrid()
+    vtk_obj.dimensions = vtk_res
+    vtk_obj.x_coordinates = x_n[:,0,0] if dim == 2 else x_n[:,0,0,0]
+    vtk_obj.y_coordinates = x_n[0,:,1] if dim == 2 else x_n[0,:,0,1]
+    vtk_obj.z_coordinates = np.zeros_like(
+        vtk_obj.x_coordinates) if dim == 2 else x_n[0,0,:,2]
+
+    # displacements are mandatory, so they get added independently of
+    # any additional point data
+    disp = np.zeros([np.prod(vtk_res), 3])
+    disp[:,:dim] = (placement - x_n).reshape(-1, dim, order="F")
+    vtk_obj.point_data.vectors = disp
+    vtk_obj.point_data.vectors.name = "displacement"
+
+    # check name clashes:
+    if point_data is None:
+        point_data = dict()
+    if cell_data is None:
+        cell_data = dict()
+    if "displacement" in point_data.keys():
+        raise Exception("Name 'displacement' is reserved")
+    if "displacement" in cell_data.keys():
+        raise Exception("Name 'displacement' is reserved")
+    #clash_set = set(point_data.keys()) & set(cell_data.keys())
+    #if clash_set:
+    #    clash_names = ["'{}'".format(name) for name in clash_set]
+    #    clash_names_fmt = ", ".join(clash_names)
+    #    raise Exception(
+    #        ("Only unique names are allowed, but the names {} appear in both "
+    #         "point- and cell data").format(clash_names_fmt))
+    #Note: is this necessary -> only same name in same dictionary doesn't make sense.
+    #      I think it's ok to have e.g. the material phase as cell and point data?
+    #Richard: I think the same, so I commented the check and would throw it away
+
+    # helper functions to add data
+    def add_data(value, name, point=True):
+        data = vtk_obj.point_data if point else vtk_obj.cell_data
+        data.add_array(value)
+        data.get_array(data._get_number_of_arrays()-1).name = name
+        return
+
+    def add_scalar(value, name, point=True):
+        add_data(value.reshape(-1, order="F"), name, point)
+        return
+
+    def add_vector(value, name, point=True):
+        res = vtk_res if point else res_c
+        vec = np.zeros([np.prod(res), 3])
+        vec[:,:dim] = value.reshape(-1, dim, order="F")
+        add_data(vec, name, point)
+        return
+
+    def add_tensor(value, name, point=True):
+        res = vtk_res if point else res_c
+        tens = np.zeros([np.prod(res), 3, 3])
+        tens[:,:dim, :dim] = value.reshape(-1, dim, dim, order="F")
+        add_data(tens.reshape(-1, 3*3), name, point)
+        return
+
+    adders = {()        : add_scalar,
+              (dim,)    : add_vector,
+              (dim, dim): add_tensor}
+
+    def shape_checker(value, reference):
+        """
+        checks whether values have the right shape and determines the
+        appropriate function to add them to the output file
+        """
+        res = value.shape[:dim]
+        shape = value.shape[dim:]
+        if not res == tuple(reference[:dim]):
+            raise Exception(
+                ("the first dim dimensions of dataset '{}' have the wrong size,"
+                 " should be {}, but got {}").format(
+                     key, reference[:dim], res))
+        if not shape in (adders.keys()):
+            raise Exception(
+                ("Can only handle scalar, vectorial, and tensorial fields, but "
+                 "got a field of shape {}").format(shape))
+        return res, shape, adders[shape]
+
+    # add point data
+    for key, value in point_data.items():
+        res, shape, adder = shape_checker(value, res_n)
+        adder(value, key, point=True)
+
+    # add cell data
+    for key, value in cell_data.items():
+        res, shape, adder = shape_checker(value, res_c)
+        adder(value, key, point=False)
+
+    path = fpath + (".vtk" if legacy_ascii else "")
+
+    write_data(vtk_obj, path)
+    return vtk_obj
diff --git a/src/common/field.hh b/src/common/field.hh
index f8d35b4..2d3ca1f 100644
--- a/src/common/field.hh
+++ b/src/common/field.hh
@@ -1,668 +1,666 @@
 /**
  * @file   field.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   07 Sep 2017
  *
  * @brief  header-only implementation of a field for field collections
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU Lesser General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with µSpectre; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * * Boston, MA 02111-1307, USA.
  *
  * Additional permission under GNU GPL version 3 section 7
  *
  * If you modify this Program, or any covered work, by linking or combining it
  * with proprietary FFT implementations or numerical libraries, containing parts
  * covered by the terms of those libraries' licenses, the licensors of this
  * Program grant you additional permission to convey the resulting work.
  */
 
 #ifndef SRC_COMMON_FIELD_HH_
 #define SRC_COMMON_FIELD_HH_
 
 #include "common/T4_map_proxy.hh"
 #include "common/field_typed.hh"
 
 #include <Eigen/Dense>
 
 #include <algorithm>
 #include <cmath>
 #include <sstream>
 #include <string>
 #include <type_traits>
 #include <typeinfo>
 #include <utility>
 #include <vector>
 
 namespace muSpectre {
 
   namespace internal {
 
     /* ---------------------------------------------------------------------- */
     //! declaraton for friending
     template <class FieldCollection, typename T, Dim_t NbComponents,
               bool isConst>
     class FieldMap;
 
     /* ---------------------------------------------------------------------- */
     /**
      * A `TypedSizedFieldBase` is the base class for fields that contain a
      * statically known number of scalars of a statically known type per pixel
      * in a `FieldCollection`. The actual data for all pixels is
      * stored in `TypedSizeFieldBase::values`.
      * `TypedSizedFieldBase` is the base class for `MatrixField` and
      * `TensorField`.
      */
     template <class FieldCollection, typename T, Dim_t NbComponents>
     class TypedSizedFieldBase : public TypedField<FieldCollection, T> {
       friend class FieldMap<FieldCollection, T, NbComponents, true>;
       friend class FieldMap<FieldCollection, T, NbComponents, false>;
 
      public:
       //! for compatibility checks
       constexpr static auto nb_components{NbComponents};
       using Parent = TypedField<FieldCollection, T>;  //!< base class
       using Scalar = T;                               //!< for type checking
       using Base = typename Parent::Base;             //!< root base class
 
       //! storage container
       using Storage_t = typename Parent::Storage_t;
 
       //! Plain type that is being mapped (Eigen lingo)
       using EigenRep_t = Eigen::Array<T, NbComponents, Eigen::Dynamic>;
       //! maps returned when iterating over field
       using EigenMap_t = Eigen::Map<EigenRep_t>;
 
       //! maps returned when iterating over field
       using ConstEigenMap_t = Eigen::Map<const EigenRep_t>;
 
       //! constructor
       TypedSizedFieldBase(std::string unique_name,
                           FieldCollection & collection);
       virtual ~TypedSizedFieldBase() = default;
 
       //! add a new value at the end of the field
       template <class Derived>
       inline void push_back(const Eigen::DenseBase<Derived> & value);
 
       //! add a new scalar value at the end of the field
       template <bool scalar_store = NbComponents == 1>
       inline std::enable_if_t<scalar_store> push_back(const T & value);
 
       /**
        * returns an upcasted reference to a field, or throws an
        * exception if the field is incompatible
        */
       static TypedSizedFieldBase & check_ref(Base & other);
       /**
        * returns an upcasted reference to a field, or throws an
        * exception if the field is incompatible
        */
       static const TypedSizedFieldBase & check_ref(const Base & other);
 
       //! return a map representing the entire field as a single `Eigen::Array`
       inline EigenMap_t eigen();
       //! return a map representing the entire field as a single `Eigen::Array`
       inline ConstEigenMap_t eigen() const;
       /**
        * return a map representing the entire field as a single
        * dynamically sized `Eigen::Array` (for python bindings)
        */
       inline typename Parent::EigenMap_t dyn_eigen() { return Parent::eigen(); }
 
       //! inner product between compatible fields
       template <typename T2>
       inline Real inner_product(
           const TypedSizedFieldBase<FieldCollection, T2, NbComponents> & other)
           const;
 
      protected:
       //! returns a raw pointer to the entry, for `Eigen::Map`
       inline T * get_ptr_to_entry(const size_t && index);
 
       //! returns a raw pointer to the entry, for `Eigen::Map`
       inline const T * get_ptr_to_entry(const size_t && index) const;
     };
 
   }  // namespace internal
 
   /* ---------------------------------------------------------------------- */
   /**
    * The `TensorField` is a subclass of
    * `muSpectre::internal::TypedSizedFieldBase` that represents tensorial
    * fields, i.e. arbitrary-dimensional arrays with identical number of
    * rows/columns (that typically correspond to the spatial cartesian
    * dimensions). It is defined by the stored scalar type @a T, the tensorial
    * order @a order (often also called degree or rank) and the number of spatial
    * dimensions @a dim.
    */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
   class TensorField : public internal::TypedSizedFieldBase<FieldCollection, T,
                                                            ipow(dim, order)> {
    public:
     //! base class
     using Parent =
         internal::TypedSizedFieldBase<FieldCollection, T, ipow(dim, order)>;
     using Base = typename Parent::Base;  //!< root base class
     //! polymorphic base class
     using Field_p = typename FieldCollection::Field_p;
     using Scalar = typename Parent::Scalar;  //!< for type checking
     //! Copy constructor
     TensorField(const TensorField & other) = delete;
 
     //! Move constructor
     TensorField(TensorField && other) = delete;
 
     //! Destructor
     virtual ~TensorField() = default;
 
     //! Copy assignment operator
     TensorField & operator=(const TensorField & other) = delete;
 
     //! Move assignment operator
     TensorField & operator=(TensorField && other) = delete;
 
     //! return the order of the stored tensor
     inline Dim_t get_order() const;
     //! return the dimension of the stored tensor
     inline Dim_t get_dim() const;
 
     //! factory function
     template <class FieldType, class CollectionType, typename... Args>
     friend FieldType & make_field(std::string unique_name,
                                   CollectionType & collection, Args &&... args);
 
     //! return a reference or throw an exception if `other` is incompatible
     static TensorField & check_ref(Base & other) {
       return static_cast<TensorField &>(Parent::check_ref(other));
     }
     //! return a reference or throw an exception if `other` is incompatible
     static const TensorField & check_ref(const Base & other) {
       return static_cast<const TensorField &>(Parent::check_ref(other));
     }
 
     /**
      * Convenience functions to return a map onto this field. A map allows
      * iteration over all pixels. The map's iterator returns an object that
      * represents the underlying mathematical structure of the field and
      * implements common linear algebra operations on it.
      * Specifically, this function returns
      * - A `MatrixFieldMap` with @a dim rows and one column if the tensorial
      * order @a order is unity.
      * - A `MatrixFieldMap` with @a dim rows and @a dim columns if the tensorial
      * order @a order is 2.
      * - A `T4MatrixFieldMap` if the tensorial order is 4.
      */
     inline decltype(auto) get_map();
     /**
      * Convenience functions to return a map onto this field. A map allows
      * iteration over all pixels. The map's iterator returns an object that
      * represents the underlying mathematical structure of the field and
      * implements common linear algebra operations on it.
      * Specifically, this function returns
      * - A `MatrixFieldMap` with @a dim rows and one column if the tensorial
      * order @a order is unity.
      * - A `MatrixFieldMap` with @a dim rows and @a dim columns if the tensorial
      * order @a order is 2.
      * - A `T4MatrixFieldMap` if the tensorial order is 4.
      */
     inline decltype(auto) get_const_map();
     /**
      * Convenience functions to return a map onto this field. A map allows
      * iteration over all pixels. The map's iterator returns an object that
      * represents the underlying mathematical structure of the field and
      * implements common linear algebra operations on it.
      * Specifically, this function returns
      * - A `MatrixFieldMap` with @a dim rows and one column if the tensorial
      * order @a order is unity.
      * - A `MatrixFieldMap` with @a dim rows and @a dim columns if the tensorial
      * order @a order is 2.
      * - A `T4MatrixFieldMap` if the tensorial order is 4.
      */
     inline decltype(auto) get_map() const;
 
     /**
      * creates a `TensorField` same size and type as this, but all
      * entries are zero. Convenience function
      */
     inline TensorField & get_zeros_like(std::string unique_name) const;
 
    protected:
     //! constructor protected!
     TensorField(std::string unique_name, FieldCollection & collection);
 
    private:
   };
 
   /* ---------------------------------------------------------------------- */
   /**
    * The `MatrixField` is subclass of `muSpectre::internal::TypedSizedFieldBase`
    * that represents matrix fields, i.e. a two dimensional arrays, defined by
    * the stored scalar type @a T and the number of rows @a NbRow and columns
    * @a NbCol of the matrix.
    */
   template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol = NbRow>
   class MatrixField : public internal::TypedSizedFieldBase<FieldCollection, T,
                                                            NbRow * NbCol> {
    public:
     //! base class
     using Parent =
         internal::TypedSizedFieldBase<FieldCollection, T, NbRow * NbCol>;
     using Base = typename Parent::Base;  //!< root base class
     //! polymorphic base field ptr to store
     using Field_p = std::unique_ptr<internal::FieldBase<FieldCollection>>;
     //! Copy constructor
     MatrixField(const MatrixField & other) = delete;
 
     //! Move constructor
     MatrixField(MatrixField && other) = delete;
 
     //! Destructor
     virtual ~MatrixField() = default;
 
     //! Copy assignment operator
     MatrixField & operator=(const MatrixField & other) = delete;
 
     //! Move assignment operator
     MatrixField & operator=(MatrixField && other) = delete;
 
     //! returns the number of rows
     inline Dim_t get_nb_row() const;
     //! returns the number of columns
     inline Dim_t get_nb_col() const;
 
     //! factory function
     template <class FieldType, class CollectionType, typename... Args>
     friend FieldType & make_field(std::string unique_name,
                                   CollectionType & collection, Args &&... args);
 
     //! returns a `MatrixField` reference if `other` is a compatible field
     static MatrixField & check_ref(Base & other) {
       return static_cast<MatrixField &>(Parent::check_ref(other));
     }
     //! returns a `MatrixField` reference if `other` is a compatible field
     static const MatrixField & check_ref(const Base & other) {
       return static_cast<const MatrixField &>(Parent::check_ref(other));
     }
 
     /**
      * Convenience functions to return a map onto this field. A map allows
      * iteration over all pixels. The map's iterator returns an object that
      * represents the underlying mathematical structure of the field and
      * implements common linear algebra operations on it.
      * Specifically, this function returns
      * - A `ScalarFieldMap` if @a NbRows and @a NbCols are unity.
      * - A `MatrixFieldMap` with @a NbRows rows and @a NbCols columns
      * otherwise.
      */
     inline decltype(auto) get_map();
     /**
      * Convenience functions to return a map onto this field. A map allows
      * iteration over all pixels. The map's iterator returns an object that
      * represents the underlying mathematical structure of the field and
      * implements common linear algebra operations on it.
      * Specifically, this function returns
      * - A `ScalarFieldMap` if @a NbRows and @a NbCols are unity.
      * - A `MatrixFieldMap` with @a NbRows rows and @a NbCols columns
      * otherwise.
      */
     inline decltype(auto) get_const_map();
     /**
      * Convenience functions to return a map onto this field. A map allows
      * iteration over all pixels. The map's iterator returns an object that
      * represents the underlying mathematical structure of the field and
      * implements common linear algebra operations on it.
      * Specifically, this function returns
      * - A `ScalarFieldMap` if @a NbRows and @a NbCols are unity.
      * - A `MatrixFieldMap` with @a NbRows rows and @a NbCols columns
      * otherwise.
      */
     inline decltype(auto) get_map() const;
 
     /**
      * creates a `MatrixField` same size and type as this, but all
      * entries are zero. Convenience function
      */
     inline MatrixField & get_zeros_like(std::string unique_name) const;
 
    protected:
     //! constructor protected!
     MatrixField(std::string unique_name, FieldCollection & collection);
-
-   private:
   };
 
   /* ---------------------------------------------------------------------- */
   //! convenience alias (
   template <class FieldCollection, typename T>
   using ScalarField = MatrixField<FieldCollection, T, 1, 1>;
   /* ---------------------------------------------------------------------- */
   // Implementations
   /* ---------------------------------------------------------------------- */
   namespace internal {
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents>
     TypedSizedFieldBase<FieldCollection, T, NbComponents>::TypedSizedFieldBase(
         std::string unique_name, FieldCollection & collection)
         : Parent(unique_name, collection, NbComponents) {
       static_assert(
           (std::is_arithmetic<T>::value || std::is_same<T, Complex>::value),
           "Use TypedSizedFieldBase for integer, real or complex scalars for T");
       static_assert(NbComponents > 0,
                     "Only fields with more than 0 components");
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents>
     TypedSizedFieldBase<FieldCollection, T, NbComponents> &
     TypedSizedFieldBase<FieldCollection, T, NbComponents>::check_ref(
         Base & other) {
       if (typeid(T).hash_code() != other.get_stored_typeid().hash_code()) {
         std::stringstream err_str{};
         err_str << "Cannot create a reference of type '" << typeid(T).name()
                 << "' for field '" << other.get_name() << "' of type '"
                 << other.get_stored_typeid().name() << "'";
         throw std::runtime_error(err_str.str());
       }
       // check size compatibility
       if (NbComponents != other.get_nb_components()) {
         std::stringstream err_str{};
         err_str << "Cannot create a reference to a field with " << NbComponents
                 << " components "
                 << "for field '" << other.get_name() << "' with "
                 << other.get_nb_components() << " components";
         throw std::runtime_error{err_str.str()};
       }
       return static_cast<TypedSizedFieldBase &>(other);
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents>
     const TypedSizedFieldBase<FieldCollection, T, NbComponents> &
     TypedSizedFieldBase<FieldCollection, T, NbComponents>::check_ref(
         const Base & other) {
       if (typeid(T).hash_code() != other.get_stored_typeid().hash_code()) {
         std::stringstream err_str{};
         err_str << "Cannot create a reference of type '" << typeid(T).name()
                 << "' for field '" << other.get_name() << "' of type '"
                 << other.get_stored_typeid().name() << "'";
         throw std::runtime_error(err_str.str());
       }
       // check size compatibility
       if (NbComponents != other.get_nb_components()) {
         std::stringstream err_str{};
         err_str << "Cannot create a reference toy a field with " << NbComponents
                 << " components "
                 << "for field '" << other.get_name() << "' with "
                 << other.get_nb_components() << " components";
         throw std::runtime_error{err_str.str()};
       }
       return static_cast<const TypedSizedFieldBase &>(other);
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents>
     auto TypedSizedFieldBase<FieldCollection, T, NbComponents>::eigen()
         -> EigenMap_t {
       return EigenMap_t(this->data(), NbComponents, this->size());
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents>
     auto TypedSizedFieldBase<FieldCollection, T, NbComponents>::eigen() const
         -> ConstEigenMap_t {
       return ConstEigenMap_t(this->data(), NbComponents, this->size());
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents>
     template <typename otherT>
     Real TypedSizedFieldBase<FieldCollection, T, NbComponents>::inner_product(
         const TypedSizedFieldBase<FieldCollection, otherT, NbComponents> &
             other) const {
       return (this->eigen() * other.eigen()).sum();
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents>
     T * TypedSizedFieldBase<FieldCollection, T, NbComponents>::get_ptr_to_entry(
         const size_t && index) {
       return this->data_ptr + NbComponents * std::move(index);
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents>
     const T *
     TypedSizedFieldBase<FieldCollection, T, NbComponents>::get_ptr_to_entry(
         const size_t && index) const {
       return this->data_ptr + NbComponents * std::move(index);
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents>
     template <class Derived>
     void TypedSizedFieldBase<FieldCollection, T, NbComponents>::push_back(
         const Eigen::DenseBase<Derived> & value) {
       static_assert(Derived::SizeAtCompileTime == NbComponents,
                     "You provided an array with the wrong number of entries.");
       static_assert((Derived::RowsAtCompileTime == 1) or
                         (Derived::ColsAtCompileTime == 1),
                     "You have not provided a column or row vector.");
       static_assert(not FieldCollection::Global,
                     "You can only push_back data into local field "
                     "collections.");
       for (Dim_t i = 0; i < NbComponents; ++i) {
         this->values.push_back(value(i));
       }
       ++this->current_size;
       this->data_ptr = &this->values.front();
     }
 
     /* ---------------------------------------------------------------------- */
     template <class FieldCollection, typename T, Dim_t NbComponents>
     template <bool scalar_store>
     std::enable_if_t<scalar_store>
     TypedSizedFieldBase<FieldCollection, T, NbComponents>::push_back(
         const T & value) {
       static_assert(scalar_store, "SFINAE");
       this->values.push_back(value);
       ++this->current_size;
       this->data_ptr = &this->values.front();
     }
 
   }  // namespace internal
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
   TensorField<FieldCollection, T, order, dim>::TensorField(
       std::string unique_name, FieldCollection & collection)
       : Parent(unique_name, collection) {}
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
   Dim_t TensorField<FieldCollection, T, order, dim>::get_order() const {
     return order;
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
   Dim_t TensorField<FieldCollection, T, order, dim>::get_dim() const {
     return dim;
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
   MatrixField<FieldCollection, T, NbRow, NbCol>::MatrixField(
       std::string unique_name, FieldCollection & collection)
       : Parent(unique_name, collection) {}
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
   Dim_t MatrixField<FieldCollection, T, NbRow, NbCol>::get_nb_col() const {
     return NbCol;
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
   Dim_t MatrixField<FieldCollection, T, NbRow, NbCol>::get_nb_row() const {
     return NbRow;
   }
 
 }  // namespace muSpectre
 
 #include "common/field_map.hh"
 
 namespace muSpectre {
 
   namespace internal {
 
     /* ---------------------------------------------------------------------- */
     /**
      * defines the default mapped type obtained when calling
      * `muSpectre::TensorField::get_map()`
      */
     template <class FieldCollection, typename T, size_t order, Dim_t dim,
               bool ConstMap>
     struct tensor_map_type {};
 
     /// specialisation for vectors
     template <class FieldCollection, typename T, Dim_t dim, bool ConstMap>
     struct tensor_map_type<FieldCollection, T, firstOrder, dim, ConstMap> {
       //! use this type
       using type = MatrixFieldMap<FieldCollection, T, dim, 1, ConstMap>;
     };
 
     /// specialisation to second-order tensors (matrices)
     template <class FieldCollection, typename T, Dim_t dim, bool ConstMap>
     struct tensor_map_type<FieldCollection, T, secondOrder, dim, ConstMap> {
       //! use this type
       using type = MatrixFieldMap<FieldCollection, T, dim, dim, ConstMap>;
     };
 
     /// specialisation to fourth-order tensors
     template <class FieldCollection, typename T, Dim_t dim, bool ConstMap>
     struct tensor_map_type<FieldCollection, T, fourthOrder, dim, ConstMap> {
       //! use this type
       using type = T4MatrixFieldMap<FieldCollection, T, dim, ConstMap>;
     };
 
     /* ---------------------------------------------------------------------- */
     /**
      * defines the default mapped type obtained when calling
      * `muSpectre::MatrixField::get_map()`
      */
     template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol,
               bool ConstMap>
     struct matrix_map_type {
       //! mapping type
       using type = MatrixFieldMap<FieldCollection, T, NbRow, NbCol, ConstMap>;
     };
 
     //! specialisation to scalar fields
     template <class FieldCollection, typename T, bool ConstMap>
     struct matrix_map_type<FieldCollection, T, oneD, oneD, ConstMap> {
       //! mapping type
       using type = ScalarFieldMap<FieldCollection, T, ConstMap>;
     };
 
   }  // namespace internal
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
   auto TensorField<FieldCollection, T, order, dim>::get_map()
       -> decltype(auto) {
     constexpr bool map_constness{false};
     using RawMap_t =
         typename internal::tensor_map_type<FieldCollection, T, order, dim,
                                            map_constness>::type;
     return RawMap_t(*this);
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
   auto TensorField<FieldCollection, T, order, dim>::get_const_map()
       -> decltype(auto) {
     constexpr bool map_constness{true};
     using RawMap_t =
         typename internal::tensor_map_type<FieldCollection, T, order, dim,
                                            map_constness>::type;
     return RawMap_t(*this);
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
   auto TensorField<FieldCollection, T, order, dim>::get_map() const
       -> decltype(auto) {
     constexpr bool map_constness{true};
     using RawMap_t =
         typename internal::tensor_map_type<FieldCollection, T, order, dim,
                                            map_constness>::type;
     return RawMap_t(*this);
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
   auto TensorField<FieldCollection, T, order, dim>::get_zeros_like(
       std::string unique_name) const -> TensorField & {
     return make_field<TensorField>(unique_name, this->collection);
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
   auto MatrixField<FieldCollection, T, NbRow, NbCol>::get_map()
       -> decltype(auto) {
     constexpr bool map_constness{false};
     using RawMap_t =
         typename internal::matrix_map_type<FieldCollection, T, NbRow, NbCol,
                                            map_constness>::type;
     return RawMap_t(*this);
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
   auto MatrixField<FieldCollection, T, NbRow, NbCol>::get_const_map()
       -> decltype(auto) {
     constexpr bool map_constness{true};
     using RawMap_t =
         typename internal::matrix_map_type<FieldCollection, T, NbRow, NbCol,
                                            map_constness>::type;
     return RawMap_t(*this);
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
   auto MatrixField<FieldCollection, T, NbRow, NbCol>::get_map() const
       -> decltype(auto) {
     constexpr bool map_constness{true};
     using RawMap_t =
         typename internal::matrix_map_type<FieldCollection, T, NbRow, NbCol,
                                            map_constness>::type;
     return RawMap_t(*this);
   }
 
   /* ---------------------------------------------------------------------- */
   template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
   auto MatrixField<FieldCollection, T, order, dim>::get_zeros_like(
       std::string unique_name) const -> MatrixField & {
     return make_field<MatrixField>(unique_name, this->collection);
   }
 
 }  // namespace muSpectre
 
 #endif  // SRC_COMMON_FIELD_HH_
diff --git a/src/materials/material_linear_elastic4.cc b/src/materials/material_linear_elastic4.cc
index 73be1db..4637363 100644
--- a/src/materials/material_linear_elastic4.cc
+++ b/src/materials/material_linear_elastic4.cc
@@ -1,79 +1,79 @@
 /**
  * @file   material_linear_elastic4.cc
  *
  * @author Richard Leute <richard.leute@imtek.uni-freiburg.de
  *
  * @date   15 March 2018
  *
  * @brief linear elastic material with distribution of stiffness properties.
  *        In difference to material_linear_elastic3 two Lame constants are
  *        stored per pixel instead of the whole elastic matrix C.
  *        Uses the MaterialMuSpectre facilities to keep it simple.
  *
  * Copyright © 2018 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU Lesser General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with µSpectre; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * * Boston, MA 02111-1307, USA.
  *
  * Additional permission under GNU GPL version 3 section 7
  *
  * If you modify this Program, or any covered work, by linking or combining it
  * with proprietary FFT implementations or numerical libraries, containing parts
  * covered by the terms of those libraries' licenses, the licensors of this
  * Program grant you additional permission to convey the resulting work.
  */
 
 #include "materials/material_linear_elastic4.hh"
 
 namespace muSpectre {
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   MaterialLinearElastic4<DimS, DimM>::MaterialLinearElastic4(std::string name)
       : Parent{name}, lambda_field{make_field<Field_t>(
                           "local first Lame constant", this->internal_fields)},
         mu_field{
             make_field<Field_t>("local second Lame constant(shear modulus)",
                                 this->internal_fields)},
         internal_variables{lambda_field.get_const_map(),
                            mu_field.get_const_map()} {}
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
-  void MaterialLinearElastic4<DimS, DimM>::add_pixel(
-      const Ccoord_t<DimS> & /*pixel*/) {
-    throw std::runtime_error(
-        "this material needs pixels with Youngs modulus and Poisson ratio.");
+  void MaterialLinearElastic4<DimS, DimM>::
+  add_pixel(const Ccoord_t<DimS> & /*pixel*/) {
+    throw std::runtime_error
+      ("This material needs pixels with Youngs modulus and Poisson ratio.");
   }
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   void
   MaterialLinearElastic4<DimS, DimM>::add_pixel(const Ccoord_t<DimS> & pixel,
                                                 const Real & Young_modulus,
                                                 const Real & Poisson_ratio) {
     this->internal_fields.add_pixel(pixel);
     // store the first(lambda) and second(mu) Lame constant in the field
     Real lambda = Hooke::compute_lambda(Young_modulus, Poisson_ratio);
     Real mu = Hooke::compute_mu(Young_modulus, Poisson_ratio);
     this->lambda_field.push_back(lambda);
     this->mu_field.push_back(mu);
   }
 
   template class MaterialLinearElastic4<twoD, twoD>;
   template class MaterialLinearElastic4<twoD, threeD>;
   template class MaterialLinearElastic4<threeD, threeD>;
 
 }  // namespace muSpectre
diff --git a/src/materials/material_linear_elastic4.hh b/src/materials/material_linear_elastic4.hh
index d9ae9e8..ba27004 100644
--- a/src/materials/material_linear_elastic4.hh
+++ b/src/materials/material_linear_elastic4.hh
@@ -1,207 +1,208 @@
 /**
  * @file   material_linear_elastic4.hh
  *
  * @author Richard Leute <richard.leute@imtek.uni-freiburg.de>
  *
  * @date   15 March 2018
  *
  * @brief linear elastic material with distribution of stiffness properties.
  *        In difference to material_linear_elastic3 two Lame constants are
  *        stored per pixel instead of the whole elastic matrix C.
  *        Uses the MaterialMuSpectre facilities to keep it simple.
  *
  * Copyright © 2018 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU Lesser General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with µSpectre; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * * Boston, MA 02111-1307, USA.
  *
  * Additional permission under GNU GPL version 3 section 7
  *
  * If you modify this Program, or any covered work, by linking or combining it
  * with proprietary FFT implementations or numerical libraries, containing parts
  * covered by the terms of those libraries' licenses, the licensors of this
  * Program grant you additional permission to convey the resulting work.
  */
 
 #ifndef SRC_MATERIALS_MATERIAL_LINEAR_ELASTIC4_HH_
 #define SRC_MATERIALS_MATERIAL_LINEAR_ELASTIC4_HH_
 
 #include "materials/material_linear_elastic1.hh"
 #include "common/field.hh"
 #include "common/tensor_algebra.hh"
 
 #include <Eigen/Dense>
 
 namespace muSpectre {
 
   template <Dim_t DimS, Dim_t DimM>
   class MaterialLinearElastic4;
 
   /**
    * traits for objective linear elasticity with eigenstrain
    */
   template <Dim_t DimS, Dim_t DimM>
   struct MaterialMuSpectre_traits<MaterialLinearElastic4<DimS, DimM>> {
     //! global field collection
     using GFieldCollection_t =
         typename MaterialBase<DimS, DimM>::GFieldCollection_t;
 
     //! expected map type for strain fields
     using StrainMap_t =
         MatrixFieldMap<GFieldCollection_t, Real, DimM, DimM, true>;
     //! expected map type for stress fields
     using StressMap_t = MatrixFieldMap<GFieldCollection_t, Real, DimM, DimM>;
     //! expected map type for tangent stiffness fields
     using TangentMap_t = T4MatrixFieldMap<GFieldCollection_t, Real, DimM>;
 
     //! declare what type of strain measure your law takes as input
     constexpr static auto strain_measure{StrainMeasure::GreenLagrange};
     //! declare what type of stress measure your law yields as output
     constexpr static auto stress_measure{StressMeasure::PK2};
 
     //! local field_collections used for internals
     using LFieldColl_t = LocalFieldCollection<DimS>;
     //! local Lame constant type
     using LLameConstantMap_t = ScalarFieldMap<LFieldColl_t, Real, true>;
     //! elasticity without internal variables
     using InternalVariables =
         std::tuple<LLameConstantMap_t, LLameConstantMap_t>;
   };
 
   /**
    * implements objective linear elasticity with an eigenstrain per pixel
    */
   template <Dim_t DimS, Dim_t DimM>
   class MaterialLinearElastic4
       : public MaterialMuSpectre<MaterialLinearElastic4<DimS, DimM>, DimS,
                                  DimM> {
    public:
     //! base class
     using Parent = MaterialMuSpectre<MaterialLinearElastic4, DimS, DimM>;
     /**
      * type used to determine whether the
      * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only
      * stresses or also tangent stiffnesses
      */
     using NeedTangent = typename Parent::NeedTangent;
     //! global field collection
 
     using Stiffness_t =
         Eigen::TensorFixedSize<Real, Eigen::Sizes<DimM, DimM, DimM, DimM>>;
 
     //! traits of this material
     using traits = MaterialMuSpectre_traits<MaterialLinearElastic4>;
 
     //! Type of container used for storing eigenstrain
     using InternalVariables = typename traits::InternalVariables;
 
     //! Hooke's law implementation
     using Hooke =
         typename MatTB::Hooke<DimM, typename traits::StrainMap_t::reference,
                               typename traits::TangentMap_t::reference>;
 
     //! Default constructor
     MaterialLinearElastic4() = delete;
 
     //! Construct by name
     explicit MaterialLinearElastic4(std::string name);
 
     //! Copy constructor
     MaterialLinearElastic4(const MaterialLinearElastic4 & other) = delete;
 
     //! Move constructor
     MaterialLinearElastic4(MaterialLinearElastic4 && other) = delete;
 
     //! Destructor
     virtual ~MaterialLinearElastic4() = default;
 
     //! Copy assignment operator
     MaterialLinearElastic4 &
     operator=(const MaterialLinearElastic4 & other) = delete;
 
     //! Move assignment operator
     MaterialLinearElastic4 &
     operator=(MaterialLinearElastic4 && other) = delete;
 
     /**
      * evaluates second Piola-Kirchhoff stress given the Green-Lagrange
      * strain (or Cauchy stress if called with a small strain tensor), the first
      * Lame constant (lambda) and the second Lame constant (shear modulus/mu).
      */
     template <class s_t>
     inline decltype(auto) evaluate_stress(s_t && E, const Real & lambda,
                                           const Real & mu);
 
     /**
      * evaluates both second Piola-Kirchhoff stress and stiffness given
      * the Green-Lagrange strain (or Cauchy stress and stiffness if
      * called with a small strain tensor), the first Lame constant (lambda) and
      * the second Lame constant (shear modulus/mu).
      */
     template <class s_t>
     inline decltype(auto) evaluate_stress_tangent(s_t && E, const Real & lambda,
                                                   const Real & mu);
 
     /**
      * return the empty internals tuple
      */
     InternalVariables & get_internals() { return this->internal_variables; }
 
     /**
      * overload add_pixel to write into loacal stiffness tensor
      */
     void add_pixel(const Ccoord_t<DimS> & pixel) final;
 
     /**
      * overload add_pixel to write into local stiffness tensor
      */
-    void add_pixel(const Ccoord_t<DimS> & pixel, const Real & Poisson_ratio,
-                   const Real & Youngs_modulus);
+    void add_pixel(const Ccoord_t<DimS> & pixel,
+                   const Real & Youngs_modulus, const Real & Poisson_ratio);
 
    protected:
     //! storage for first Lame constant 'lambda'
     //! and second Lame constant(shear modulus) 'mu'
     using Field_t = MatrixField<LocalFieldCollection<DimS>, Real, oneD, oneD>;
     Field_t & lambda_field;
     Field_t & mu_field;
     //! tuple for iterable eigen_field
     InternalVariables internal_variables;
 
    private:
   };
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   template <class s_t>
   auto MaterialLinearElastic4<DimS, DimM>::evaluate_stress(s_t && E,
                                                            const Real & lambda,
                                                            const Real & mu)
       -> decltype(auto) {
     auto C = Hooke::compute_C_T4(lambda, mu);
     return Matrices::tensmult(C, E);
   }
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   template <class s_t>
   auto MaterialLinearElastic4<DimS, DimM>::evaluate_stress_tangent(
       s_t && E, const Real & lambda, const Real & mu) -> decltype(auto) {
-    auto C = Hooke::compute_C_T4(lambda, mu);
-    return std::make_tuple(Matrices::tensmult(C, E), C);
+    T4Mat<Real, DimM> C = Hooke::compute_C_T4(lambda, mu);
+    return std::make_tuple(
+        this->evaluate_stress(std::forward<s_t>(E), lambda, mu), C);
   }
 
 }  // namespace muSpectre
 
 #endif  // SRC_MATERIALS_MATERIAL_LINEAR_ELASTIC4_HH_
diff --git a/src/materials/material_muSpectre_base.hh b/src/materials/material_muSpectre_base.hh
index 491c3bb..0cc3384 100644
--- a/src/materials/material_muSpectre_base.hh
+++ b/src/materials/material_muSpectre_base.hh
@@ -1,721 +1,721 @@
 /**
  * @file   material_muSpectre_base.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   25 Oct 2017
  *
  * @brief  Base class for materials written for µSpectre specifically. These
  *         can take full advantage of the configuration-change utilities of
  *         µSpectre. The user can inherit from them to define new constitutive
  *         laws and is merely required to provide the methods for computing the
  *         second Piola-Kirchhoff stress and Tangent. This class uses the
  *         "curiously recurring template parameter" to avoid virtual calls.
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU Lesser General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with µSpectre; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * * Boston, MA 02111-1307, USA.
+ * Boston, MA 02111-1307, USA.
  *
  * Additional permission under GNU GPL version 3 section 7
  *
  * If you modify this Program, or any covered work, by linking or combining it
  * with proprietary FFT implementations or numerical libraries, containing parts
  * covered by the terms of those libraries' licenses, the licensors of this
  * Program grant you additional permission to convey the resulting work.
  */
 
 #ifndef SRC_MATERIALS_MATERIAL_MUSPECTRE_BASE_HH_
 #define SRC_MATERIALS_MATERIAL_MUSPECTRE_BASE_HH_
 
 #include "common/common.hh"
 #include "materials/material_base.hh"
 #include "materials/materials_toolbox.hh"
 #include "materials/material_evaluator.hh"
 #include "common/field_collection.hh"
 #include "common/field.hh"
 #include "common//utilities.hh"
 
 #include <tuple>
 #include <type_traits>
 #include <iterator>
 #include <stdexcept>
 
 namespace muSpectre {
 
   // Forward declaration for factory function
   template <Dim_t DimS, Dim_t DimM>
   class CellBase;
 
   /**
    * material traits are used by `muSpectre::MaterialMuSpectre` to
    * break the circular dependence created by the curiously recurring
    * template parameter. These traits must define
    * - these `muSpectre::FieldMap`s:
    *   - `StrainMap_t`: typically a `muSpectre::MatrixFieldMap` for a
    *                    constant second-order `muSpectre::TensorField`
    *   - `StressMap_t`: typically a `muSpectre::MatrixFieldMap` for a
    *                    writable secord-order `muSpectre::TensorField`
    *   - `TangentMap_t`: typically a `muSpectre::T4MatrixFieldMap` for a
    *                     writable fourth-order `muSpectre::TensorField`
    * - `strain_measure`: the expected strain type (will be replaced by the
    *                     small-strain tensor ε
    *                     `muspectre::StrainMeasure::Infinitesimal` in small
    *                     strain computations)
    * - `stress_measure`: the measure of the returned stress. Is used by
    *                     `muspectre::MaterialMuSpectre` to transform it into
    *                     Cauchy stress (`muspectre::StressMeasure::Cauchy`) in
    *                     small-strain computations and into first
    *                     Piola-Kirchhoff stress `muspectre::StressMeasure::PK1`
    *                     in finite-strain computations
    * - `InternalVariables`: a tuple of `muSpectre::FieldMap`s containing
    *                        internal variables
    */
   template <class Material>
   struct MaterialMuSpectre_traits {};
 
   template <class Material, Dim_t DimS, Dim_t DimM>
   class MaterialMuSpectre;
 
   /**
    * Base class for most convenient implementation of materials
    */
   template <class Material, Dim_t DimS, Dim_t DimM>
   class MaterialMuSpectre : public MaterialBase<DimS, DimM> {
    public:
     /**
      * type used to determine whether the
      * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only
      * stresses or also tangent stiffnesses
      */
     using NeedTangent = MatTB::NeedTangent;
     using Parent = MaterialBase<DimS, DimM>;  //!< base class
     //! global field collection
     using GFieldCollection_t = typename Parent::GFieldCollection_t;
     //! expected type for stress fields
     using StressField_t = typename Parent::StressField_t;
     //! expected type for strain fields
     using StrainField_t = typename Parent::StrainField_t;
     //! expected type for tangent stiffness fields
     using TangentField_t = typename Parent::TangentField_t;
 
     //! traits for the CRTP subclass
     using traits = MaterialMuSpectre_traits<Material>;
 
     //! Default constructor
     MaterialMuSpectre() = delete;
 
     //! Construct by name
     explicit MaterialMuSpectre(std::string name);
 
     //! Copy constructor
     MaterialMuSpectre(const MaterialMuSpectre & other) = delete;
 
     //! Move constructor
     MaterialMuSpectre(MaterialMuSpectre && other) = delete;
 
     //! Destructor
     virtual ~MaterialMuSpectre() = default;
 
     //! Factory
     template <class... ConstructorArgs>
     static Material & make(CellBase<DimS, DimM> & cell,
                            ConstructorArgs &&... args);
     /** Factory
      * takes all arguments after the name of the underlying
      * Material's constructor. E.g., if the underlying material is a
      * `muSpectre::MaterialLinearElastic1<threeD>`, these would be Young's
      * modulus and Poisson's ratio.
      */
     template <class... ConstructorArgs>
     static std::tuple<std::shared_ptr<Material>, MaterialEvaluator<DimM>>
     make_evaluator(ConstructorArgs &&... args);
 
     //! Copy assignment operator
     MaterialMuSpectre & operator=(const MaterialMuSpectre & other) = delete;
 
     //! Move assignment operator
     MaterialMuSpectre & operator=(MaterialMuSpectre && other) = delete;
 
     //! allocate memory, etc
     void initialise() override;
 
     using Parent::compute_stresses;
     using Parent::compute_stresses_tangent;
     //! computes stress
     void compute_stresses(const StrainField_t & F, StressField_t & P,
                           Formulation form) final;
     //! computes stress and tangent modulus
     void compute_stresses_tangent(const StrainField_t & F, StressField_t & P,
                                   TangentField_t & K, Formulation form) final;
 
    protected:
     //! computes stress with the formulation available at compile time
     //! __attribute__ required by g++-6 and g++-7 because of this bug:
     //! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=80947
     template <Formulation Form>
     inline void compute_stresses_worker(const StrainField_t & F,
                                         StressField_t & P)
         __attribute__((visibility("default")));
 
     //! computes stress with the formulation available at compile time
     //! __attribute__ required by g++-6 and g++-7 because of this bug:
     //! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=80947
     template <Formulation Form>
     inline void compute_stresses_worker(const StrainField_t & F,
                                         StressField_t & P, TangentField_t & K)
         __attribute__((visibility("default")));
     //! this iterable class is a default for simple laws that just take a strain
     //! the iterable is just a templated wrapper to provide a range to iterate
     //! over that does or does not include tangent moduli
     template <NeedTangent need_tgt = NeedTangent::no>
     class iterable_proxy;
 
     /**
      * inheriting classes with internal variables need to overload this function
      */
     typename traits::InternalVariables & get_internals() {
       return static_cast<Material &>(*this).get_internals();
     }
 
     bool is_initialised{false};  //!< to handle double initialisation right
 
    private:
   };
 
   /* ---------------------------------------------------------------------- */
   template <class Material, Dim_t DimS, Dim_t DimM>
   MaterialMuSpectre<Material, DimS, DimM>::MaterialMuSpectre(std::string name)
       : Parent(name) {
     using stress_compatible =
         typename traits::StressMap_t::template is_compatible<StressField_t>;
     using strain_compatible =
         typename traits::StrainMap_t::template is_compatible<StrainField_t>;
     using tangent_compatible =
         typename traits::TangentMap_t::template is_compatible<TangentField_t>;
 
     static_assert((stress_compatible::value && stress_compatible::explain()),
                   "The material's declared stress map is not compatible "
                   "with the stress field. More info in previously shown "
                   "assert.");
 
     static_assert((strain_compatible::value && strain_compatible::explain()),
                   "The material's declared strain map is not compatible "
                   "with the strain field. More info in previously shown "
                   "assert.");
 
     static_assert((tangent_compatible::value && tangent_compatible::explain()),
                   "The material's declared tangent map is not compatible "
                   "with the tangent field. More info in previously shown "
                   "assert.");
   }
 
   /* ---------------------------------------------------------------------- */
   template <class Material, Dim_t DimS, Dim_t DimM>
   template <class... ConstructorArgs>
   Material &
   MaterialMuSpectre<Material, DimS, DimM>::make(CellBase<DimS, DimM> & cell,
                                                 ConstructorArgs &&... args) {
     auto mat = std::make_unique<Material>(args...);
     auto & mat_ref = *mat;
     cell.add_material(std::move(mat));
     return mat_ref;
   }
 
   /* ---------------------------------------------------------------------- */
   template <class Material, Dim_t DimS, Dim_t DimM>
   template <class... ConstructorArgs>
   std::tuple<std::shared_ptr<Material>, MaterialEvaluator<DimM>>
   MaterialMuSpectre<Material, DimS, DimM>::make_evaluator(
       ConstructorArgs &&... args) {
     auto mat = std::make_shared<Material>("name", args...);
     using Ret_t =
         std::tuple<std::shared_ptr<Material>, MaterialEvaluator<DimM>>;
     return Ret_t(mat, MaterialEvaluator<DimM>{mat});
   }
 
   /* ---------------------------------------------------------------------- */
   template <class Material, Dim_t DimS, Dim_t DimM>
   void MaterialMuSpectre<Material, DimS, DimM>::initialise() {
     if (!this->is_initialised) {
       this->internal_fields.initialise();
       this->is_initialised = true;
     }
   }
 
   /* ---------------------------------------------------------------------- */
   template <class Material, Dim_t DimS, Dim_t DimM>
   void MaterialMuSpectre<Material, DimS, DimM>::compute_stresses(
       const StrainField_t & F, StressField_t & P, Formulation form) {
     switch (form) {
     case Formulation::finite_strain: {
       this->template compute_stresses_worker<Formulation::finite_strain>(F, P);
       break;
     }
     case Formulation::small_strain: {
       this->template compute_stresses_worker<Formulation::small_strain>(F, P);
       break;
     }
     default:
       throw std::runtime_error("Unknown formulation");
       break;
     }
   }
 
   /* ---------------------------------------------------------------------- */
   template <class Material, Dim_t DimS, Dim_t DimM>
   void MaterialMuSpectre<Material, DimS, DimM>::compute_stresses_tangent(
       const StrainField_t & F, StressField_t & P, TangentField_t & K,
       Formulation form) {
     switch (form) {
     case Formulation::finite_strain: {
       this->template compute_stresses_worker<Formulation::finite_strain>(F, P,
                                                                          K);
       break;
     }
     case Formulation::small_strain: {
       this->template compute_stresses_worker<Formulation::small_strain>(F, P,
                                                                         K);
       break;
     }
     default:
       throw std::runtime_error("Unknown formulation");
       break;
     }
   }
 
   /* ---------------------------------------------------------------------- */
   template <class Material, Dim_t DimS, Dim_t DimM>
   template <Formulation Form>
   void MaterialMuSpectre<Material, DimS, DimM>::compute_stresses_worker(
       const StrainField_t & F, StressField_t & P, TangentField_t & K) {
     /* These lambdas are executed for every integration point.
 
        F contains the transformation gradient for finite strain calculations and
        the infinitesimal strain tensor in small strain problems
 
        The internal_variables tuple contains whatever internal variables
        Material declared (e.g., eigenstrain, strain rate, etc.)
     */
     using Strains_t = std::tuple<typename traits::StrainMap_t::reference>;
     using Stresses_t = std::tuple<typename traits::StressMap_t::reference,
                                   typename traits::TangentMap_t::reference>;
     auto constitutive_law_small_strain = [this](Strains_t Strains,
                                                 Stresses_t Stresses,
                                                 auto && internal_variables) {
       constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)};
       constexpr StrainMeasure expected_strain_m{
           get_formulation_strain_type(Form, traits::strain_measure)};
 
       auto & this_mat = static_cast<Material &>(*this);
 
       // Transformation gradient is first in the strains tuple
       auto & F = std::get<0>(Strains);
       auto && strain =
           MatTB::convert_strain<stored_strain_m, expected_strain_m>(F);
       // return value contains a tuple of rvalue_refs to both stress and tangent
       // moduli
       Stresses = apply(
           [&strain, &this_mat](auto &&... internals) {
             return this_mat.evaluate_stress_tangent(std::move(strain),
                                                     internals...);
           },
           internal_variables);
     };
 
     auto constitutive_law_finite_strain = [this](Strains_t Strains,
                                                  Stresses_t Stresses,
                                                  auto && internal_variables) {
       constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)};
       constexpr StrainMeasure expected_strain_m{
           get_formulation_strain_type(Form, traits::strain_measure)};
       auto & this_mat = static_cast<Material &>(*this);
 
       // Transformation gradient is first in the strains tuple
       auto & grad = std::get<0>(Strains);
       auto && strain =
           MatTB::convert_strain<stored_strain_m, expected_strain_m>(grad);
       // TODO(junge): Figure this out: I can't std::move(internals...),
       // because if there are no internals, compilation fails with "no
       // matching function for call to ‘move()’'. These are tuples of
       // lvalue references, so it shouldn't be too bad, but still
       // irksome.
 
       // return value contains a tuple of rvalue_refs to both stress
       // and tangent moduli
       auto stress_tgt = apply(
           [&strain, &this_mat](auto &&... internals) {
             return this_mat.evaluate_stress_tangent(std::move(strain),
                                                     internals...);
           },
           internal_variables);
       auto & stress = std::get<0>(stress_tgt);
       auto & tangent = std::get<1>(stress_tgt);
       Stresses =
           MatTB::PK1_stress<traits::stress_measure, traits::strain_measure>(
               std::move(grad), std::move(stress), std::move(tangent));
     };
 
     iterable_proxy<NeedTangent::yes> fields{*this, F, P, K};
     for (auto && arglist : fields) {
       /**
        * arglist is a tuple of three tuples containing only Lvalue
        * references (see value_tye in the class definition of
        * iterable_proxy::iterator). Tuples contain strains, stresses
        * and internal variables, respectively,
        */
 
       // auto && stress_tgt = std::get<0>(tuples);
       // auto && inputs = std::get<1>(tuples);TODO:clean this
       static_assert(std::is_same<typename traits::StrainMap_t::reference,
                                  std::remove_reference_t<decltype(std::get<0>(
                                      std::get<0>(arglist)))>>::value,
                     "Type mismatch for strain reference, check iterator "
                     "value_type");
       static_assert(std::is_same<typename traits::StressMap_t::reference,
                                  std::remove_reference_t<decltype(std::get<0>(
                                      std::get<1>(arglist)))>>::value,
                     "Type mismatch for stress reference, check iterator"
                     "value_type");
       static_assert(std::is_same<typename traits::TangentMap_t::reference,
                                  std::remove_reference_t<decltype(std::get<1>(
                                      std::get<1>(arglist)))>>::value,
                     "Type mismatch for tangent reference, check iterator"
                     "value_type");
 
       switch (Form) {
       case Formulation::small_strain: {
         apply(constitutive_law_small_strain, std::move(arglist));
         break;
       }
       case Formulation::finite_strain: {
         apply(constitutive_law_finite_strain, std::move(arglist));
         break;
       }
       }
     }
   }
 
   /* ---------------------------------------------------------------------- */
   template <class Material, Dim_t DimS, Dim_t DimM>
   template <Formulation Form>
   void MaterialMuSpectre<Material, DimS, DimM>::compute_stresses_worker(
       const StrainField_t & F, StressField_t & P) {
     /* These lambdas are executed for every integration point.
 
        F contains the transformation gradient for finite strain calculations and
        the infinitesimal strain tensor in small strain problems
 
        The internal_variables tuple contains whatever internal variables
        Material declared (e.g., eigenstrain, strain rate, etc.)
     */
 
     using Strains_t = std::tuple<typename traits::StrainMap_t::reference>;
     using Stresses_t = std::tuple<typename traits::StressMap_t::reference>;
     auto constitutive_law_small_strain = [this](Strains_t Strains,
                                                 Stresses_t Stresses,
                                                 auto && internal_variables) {
       constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)};
       constexpr StrainMeasure expected_strain_m{
           get_formulation_strain_type(Form, traits::strain_measure)};
 
       auto & this_mat = static_cast<Material &>(*this);
 
       // Transformation gradient is first in the strains tuple
       auto & F = std::get<0>(Strains);
       auto && strain =
           MatTB::convert_strain<stored_strain_m, expected_strain_m>(F);
       // return value contains a tuple of rvalue_refs to both stress and tangent
       // moduli
       auto & sigma = std::get<0>(Stresses);
       sigma = apply(
           [&strain, &this_mat](auto &&... internals) {
             return this_mat.evaluate_stress(std::move(strain), internals...);
           },
           internal_variables);
     };
 
     auto constitutive_law_finite_strain = [this](Strains_t Strains,
                                                  Stresses_t && Stresses,
                                                  auto && internal_variables) {
       constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)};
       constexpr StrainMeasure expected_strain_m{
           get_formulation_strain_type(Form, traits::strain_measure)};
       auto & this_mat = static_cast<Material &>(*this);
 
       // Transformation gradient is first in the strains tuple
       auto & F = std::get<0>(Strains);
       auto && strain =
           MatTB::convert_strain<stored_strain_m, expected_strain_m>(F);
 
       // TODO(junge): Figure this out: I can't std::move(internals...),
       // because if there are no internals, compilation fails with "no
       // matching function for call to ‘move()’'. These are tuples of
       // lvalue references, so it shouldn't be too bad, but still
       // irksome.
 
       // return value contains a tuple of rvalue_refs to both stress
       // and tangent moduli
       auto && stress = apply(
           [&strain, &this_mat](auto &&... internals) {
             return this_mat.evaluate_stress(std::move(strain), internals...);
           },
           internal_variables);
       auto & P = get<0>(Stresses);
       P = MatTB::PK1_stress<traits::stress_measure, traits::strain_measure>(
           F, stress);
     };
 
     iterable_proxy<NeedTangent::no> fields{*this, F, P};
 
     for (auto && arglist : fields) {
       /**
        * arglist is a tuple of three tuples containing only Lvalue
        * references (see value_tye in the class definition of
        * iterable_proxy::iterator). Tuples contain strains, stresses
        * and internal variables, respectively,
        */
 
       // auto && stress_tgt = std::get<0>(tuples);
       // auto && inputs = std::get<1>(tuples);TODO:clean this
       static_assert(std::is_same<typename traits::StrainMap_t::reference,
                                  std::remove_reference_t<decltype(std::get<0>(
                                      std::get<0>(arglist)))>>::value,
                     "Type mismatch for strain reference, check iterator "
                     "value_type");
       static_assert(std::is_same<typename traits::StressMap_t::reference,
                                  std::remove_reference_t<decltype(std::get<0>(
                                      std::get<1>(arglist)))>>::value,
                     "Type mismatch for stress reference, check iterator"
                     "value_type");
 
       switch (Form) {
       case Formulation::small_strain: {
         apply(constitutive_law_small_strain, std::move(arglist));
         break;
       }
       case Formulation::finite_strain: {
         apply(constitutive_law_finite_strain, std::move(arglist));
         break;
       }
       }
     }
   }
 
   /* ---------------------------------------------------------------------- */
   //! this iterator class is a default for simple laws that just take a strain
   template <class Material, Dim_t DimS, Dim_t DimM>
   template <MatTB::NeedTangent NeedTgt>
   class MaterialMuSpectre<Material, DimS, DimM>::iterable_proxy {
    public:
     //! Default constructor
     iterable_proxy() = delete;
     /**
      * type used to determine whether the
      * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only
      * stresses or also tangent stiffnesses
      */
     using NeedTangent =
         typename MaterialMuSpectre<Material, DimS, DimM>::NeedTangent;
     /** Iterator uses the material's internal variables field
         collection to iterate selectively over the global fields
         (such as the transformation gradient F and first
         Piola-Kirchhoff stress P.
     **/
     template <bool DoNeedTgt = (NeedTgt == NeedTangent::yes)>
     iterable_proxy(MaterialMuSpectre & mat, const StrainField_t & F,
                    StressField_t & P,
                    std::enable_if_t<DoNeedTgt, TangentField_t> & K)
         : material{mat}, strain_field{F},
           stress_tup{P, K}, internals{material.get_internals()} {};
 
     /** Iterator uses the material's internal variables field
         collection to iterate selectively over the global fields
         (such as the transformation gradient F and first
         Piola-Kirchhoff stress P.
     **/
     template <bool DontNeedTgt = (NeedTgt == NeedTangent::no)>
     iterable_proxy(MaterialMuSpectre & mat, const StrainField_t & F,
                    std::enable_if_t<DontNeedTgt, StressField_t> & P)
         : material{mat}, strain_field{F},
           stress_tup{P}, internals{material.get_internals()} {};
 
     //! Expected type for strain fields
     using StrainMap_t = typename traits::StrainMap_t;
     //! Expected type for stress fields
     using StressMap_t = typename traits::StressMap_t;
     //! Expected type for tangent stiffness fields
     using TangentMap_t = typename traits::TangentMap_t;
     //! expected type for strain values
     using Strain_t = typename traits::StrainMap_t::reference;
     //! expected type for stress values
     using Stress_t = typename traits::StressMap_t::reference;
     //! expected type for tangent stiffness values
     using Tangent_t = typename traits::TangentMap_t::reference;
 
     //! tuple of intervnal variables, depends on the material
     using InternalVariables = typename traits::InternalVariables;
     //! tuple containing a stress and possibly a tangent stiffness field
     using StressFieldTup =
         std::conditional_t<(NeedTgt == NeedTangent::yes),
                            std::tuple<StressField_t &, TangentField_t &>,
                            std::tuple<StressField_t &>>;
     //! tuple containing a stress and possibly a tangent stiffness field map
     using StressMapTup =
         std::conditional_t<(NeedTgt == NeedTangent::yes),
                            std::tuple<StressMap_t, TangentMap_t>,
                            std::tuple<StressMap_t>>;
     //! tuple containing a stress and possibly a tangent stiffness value ref
     using Stress_tTup = std::conditional_t<(NeedTgt == NeedTangent::yes),
                                            std::tuple<Stress_t, Tangent_t>,
                                            std::tuple<Stress_t>>;
 
     //! Copy constructor
     iterable_proxy(const iterable_proxy & other) = default;
 
     //! Move constructor
     iterable_proxy(iterable_proxy && other) = default;
 
     //! Destructor
     virtual ~iterable_proxy() = default;
 
     //! Copy assignment operator
     iterable_proxy & operator=(const iterable_proxy & other) = default;
 
     //! Move assignment operator
     iterable_proxy & operator=(iterable_proxy && other) = default;
 
     /**
      * dereferences into a tuple containing strains, and internal
      * variables, as well as maps to the stress and potentially
      * stiffness maps where to write the response of a pixel
      */
     class iterator {
      public:
       //! type to refer to internal variables owned by a CRTP material
       using InternalReferences = MatTB::ReferenceTuple_t<InternalVariables>;
       //! return type to be unpacked per pixel my the constitutive law
       using value_type =
           std::tuple<std::tuple<Strain_t>, Stress_tTup, InternalReferences>;
       using iterator_category = std::forward_iterator_tag;  //!< stl conformance
 
       //! Default constructor
       iterator() = delete;
 
       /** Iterator uses the material's internal variables field
           collection to iterate selectively over the global fields
           (such as the transformation gradient F and first
           Piola-Kirchhoff stress P.
       **/
       explicit iterator(const iterable_proxy & it, bool begin = true)
           : it{it}, strain_map{it.strain_field}, stress_map(it.stress_tup),
             index{begin ? 0 : it.material.internal_fields.size()} {}
 
       //! Copy constructor
       iterator(const iterator & other) = default;
 
       //! Move constructor
       iterator(iterator && other) = default;
 
       //! Destructor
       virtual ~iterator() = default;
 
       //! Copy assignment operator
       iterator & operator=(const iterator & other) = default;
 
       //! Move assignment operator
       iterator & operator=(iterator && other) = default;
 
       //! pre-increment
       inline iterator & operator++();
       //! dereference
       inline value_type operator*();
       //! inequality
       inline bool operator!=(const iterator & other) const;
 
      protected:
       const iterable_proxy & it;  //!< ref to the proxy
       StrainMap_t strain_map;     //!< map onto the global strain field
       //! map onto the global stress field and possibly tangent stiffness
       StressMapTup stress_map;
       size_t index;  //!< index or pixel currently referred to
     };
 
     //! returns iterator to first pixel if this material
     iterator begin() { return std::move(iterator(*this)); }
     //! returns iterator past the last pixel in this material
     iterator end() { return std::move(iterator(*this, false)); }
 
    protected:
     MaterialMuSpectre & material;        //!< reference to the proxied material
     const StrainField_t & strain_field;  //!< cell's global strain field
     //! references to the global stress field and perhaps tangent stiffness
     StressFieldTup stress_tup;
     //! references to the internal variables
     InternalVariables & internals;
 
    private:
   };
 
   /* ---------------------------------------------------------------------- */
   template <class Material, Dim_t DimS, Dim_t DimM>
   template <MatTB::NeedTangent NeedTgt>
   bool
   MaterialMuSpectre<Material, DimS, DimM>::iterable_proxy<NeedTgt>::iterator::
   operator!=(const iterator & other) const {
     return (this->index != other.index);
   }
 
   /* ---------------------------------------------------------------------- */
   template <class Material, Dim_t DimS, Dim_t DimM>
   template <MatTB::NeedTangent NeedTgt>
   typename MaterialMuSpectre<Material, DimS,
                              DimM>::template iterable_proxy<NeedTgt>::iterator &
   MaterialMuSpectre<Material, DimS, DimM>::iterable_proxy<NeedTgt>::iterator::
   operator++() {
     this->index++;
     return *this;
   }
 
   /* ---------------------------------------------------------------------- */
   template <class Material, Dim_t DimS, Dim_t DimM>
   template <MatTB::NeedTangent NeedTgT>
   typename MaterialMuSpectre<Material, DimS, DimM>::template iterable_proxy<
       NeedTgT>::iterator::value_type
       MaterialMuSpectre<Material, DimS,
                         DimM>::iterable_proxy<NeedTgT>::iterator::operator*() {
     const Ccoord_t<DimS> pixel{
         this->it.material.internal_fields.get_ccoord(this->index)};
     auto && strain = std::make_tuple(this->strain_map[pixel]);
 
     auto && stresses = apply(
         [&pixel](auto &&... stress_tgt) {
           return std::make_tuple(stress_tgt[pixel]...);
         },
         this->stress_map);
     auto && internal = this->it.material.get_internals();
     const auto id{this->index};
     auto && internals = apply(
         [id](auto &&... internals_) {
           return InternalReferences{internals_[id]...};
         },
         internal);
     return std::make_tuple(std::move(strain), std::move(stresses),
                            std::move(internals));
   }
 
 }  // namespace muSpectre
 
 #endif  // SRC_MATERIALS_MATERIAL_MUSPECTRE_BASE_HH_
diff --git a/src/materials/materials_toolbox.hh b/src/materials/materials_toolbox.hh
index 474c477..383dcc6 100644
--- a/src/materials/materials_toolbox.hh
+++ b/src/materials/materials_toolbox.hh
@@ -1,577 +1,591 @@
 /**
  * @file   materials_toolbox.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   02 Nov 2017
  *
  * @brief  collection of common continuum mechanics tools
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU Lesser General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with µSpectre; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * * Boston, MA 02111-1307, USA.
  *
  * Additional permission under GNU GPL version 3 section 7
  *
  * If you modify this Program, or any covered work, by linking or combining it
  * with proprietary FFT implementations or numerical libraries, containing parts
  * covered by the terms of those libraries' licenses, the licensors of this
  * Program grant you additional permission to convey the resulting work.
  */
 
 #ifndef SRC_MATERIALS_MATERIALS_TOOLBOX_HH_
 #define SRC_MATERIALS_MATERIALS_TOOLBOX_HH_
 
 #include "common/common.hh"
 #include "common/tensor_algebra.hh"
 #include "common/eigen_tools.hh"
 #include "common/T4_map_proxy.hh"
 
 #include <Eigen/Dense>
 #include <unsupported/Eigen/MatrixFunctions>
 
 #include <exception>
 #include <sstream>
 #include <iostream>
 #include <tuple>
 #include <type_traits>
 
 namespace muSpectre {
 
   namespace MatTB {
     /**
      * thrown when generic materials-related runtime errors occur
      * (mostly continuum mechanics problems)
      */
     class MaterialsToolboxError : public std::runtime_error {
      public:
       //! constructor
       explicit MaterialsToolboxError(const std::string & what)
           : std::runtime_error(what) {}
       //! constructor
       explicit MaterialsToolboxError(const char * what)
           : std::runtime_error(what) {}
     };
 
     /* ---------------------------------------------------------------------- */
     /**
      * Flag used to designate whether the material should compute both stress
      * and tangent moduli or only stress
      */
     enum class NeedTangent {
       yes,  //!< compute both stress and tangent moduli
       no    //!< compute only stress
     };
 
     /**
      * struct used to determine the exact type of a tuple of references obtained
      * when a bunch of iterators over fiel_maps are dereferenced and their
      * results are concatenated into a tuple
      */
     template <class... T>
     struct ReferenceTuple {
       //! use this type
       using type = std::tuple<typename T::reference...>;
     };
 
     /**
      * specialisation for tuples
      */
     // template <>
     template <class... T>
     struct ReferenceTuple<std::tuple<T...>> {
       //! use this type
       using type = typename ReferenceTuple<T...>::type;
     };
 
     /**
      * helper type for ReferenceTuple
      */
     template <class... T>
     using ReferenceTuple_t = typename ReferenceTuple<T...>::type;
 
     /* ---------------------------------------------------------------------- */
     namespace internal {
 
       /** Structure for functions returning one strain measure as a
        *  function of another
        **/
       template <StrainMeasure In, StrainMeasure Out = In>
       struct ConvertStrain {
         static_assert((In == StrainMeasure::Gradient) ||
                           (In == StrainMeasure::Infinitesimal),
                       "This situation makes me suspect that you are not using "
                       "MatTb as intended. Disable this assert only if you are "
                       "sure about what you are doing.");
 
         //! returns the converted strain
         template <class Strain_t>
         inline static decltype(auto) compute(Strain_t && input) {
           // transparent case, in which no conversion is required:
           // just a perfect forwarding
           static_assert((In == Out),
                         "This particular strain conversion is not implemented");
           return std::forward<Strain_t>(input);
         }
       };
 
       /* ----------------------------------------------------------------------
        */
       /** Specialisation for getting Green-Lagrange strain from the
           transformation gradient
           E = ¹/₂ (C - I) = ¹/₂ (Fᵀ·F - I)
       **/
       template <>
       struct ConvertStrain<StrainMeasure::Gradient,
                            StrainMeasure::GreenLagrange> {
         //! returns the converted strain
         template <class Strain_t>
         inline static decltype(auto) compute(Strain_t && F) {
           return .5 * (F.transpose() * F - Strain_t::PlainObject::Identity());
         }
       };
 
       /* ----------------------------------------------------------------------
        */
       /** Specialisation for getting Left Cauchy-Green strain from the
           transformation gradient
           B = F·Fᵀ = V²
       **/
       template <>
       struct ConvertStrain<StrainMeasure::Gradient,
                            StrainMeasure::LCauchyGreen> {
         //! returns the converted strain
         template <class Strain_t>
         inline static decltype(auto) compute(Strain_t && F) {
           return F * F.transpose();
         }
       };
 
       /* ----------------------------------------------------------------------
        */
       /** Specialisation for getting Right Cauchy-Green strain from the
           transformation gradient
           C = Fᵀ·F = U²
       **/
       template <>
       struct ConvertStrain<StrainMeasure::Gradient,
                            StrainMeasure::RCauchyGreen> {
         //! returns the converted strain
         template <class Strain_t>
         inline static decltype(auto) compute(Strain_t && F) {
           return F.transpose() * F;
         }
       };
 
       /* ----------------------------------------------------------------------
        */
       /** Specialisation for getting logarithmic (Hencky) strain from the
           transformation gradient
           E₀ = ¹/₂ ln C = ¹/₂ ln (Fᵀ·F)
       **/
       template <>
       struct ConvertStrain<StrainMeasure::Gradient, StrainMeasure::Log> {
         //! returns the converted strain
         template <class Strain_t>
         inline static decltype(auto) compute(Strain_t && F) {
           constexpr Dim_t dim{EigenCheck::tensor_dim<Strain_t>::value};
           return (.5 * logm(Eigen::Matrix<Real, dim, dim>{F.transpose() * F}))
               .eval();
         }
       };
 
     }  // namespace internal
 
     /* ---------------------------------------------------------------------- */
     //! set of functions returning one strain measure as a function of
     //! another
     template <StrainMeasure In, StrainMeasure Out, class Strain_t>
     decltype(auto) convert_strain(Strain_t && strain) {
       return internal::ConvertStrain<In, Out>::compute(std::move(strain));
     }
 
     namespace internal {
 
       //! Base template for elastic modulus conversion
       template <ElasticModulus Out, ElasticModulus In1, ElasticModulus In2>
       struct Converter {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real & /*in1*/,
                                              const Real & /*in2*/) {
           // if no return has happened until now, the conversion is not
           // implemented yet
           static_assert(
               (In1 == In2),
               "This conversion has not been implemented yet, please add "
               "it here below as a specialisation of this function "
               "template. Check "
               "https://en.wikipedia.org/wiki/Lam%C3%A9_parameters for "
               "the formula.");
           return 0;
         }
       };
 
       /**
        * Spectialisation for when the output is the first input
        */
       template <ElasticModulus Out, ElasticModulus In>
       struct Converter<Out, Out, In> {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real & A,
                                              const Real & /*B*/) {
           return A;
         }
       };
 
       /**
        * Spectialisation for when the output is the second input
        */
       template <ElasticModulus Out, ElasticModulus In>
       struct Converter<Out, In, Out> {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real & /*A*/,
                                              const Real & B) {
           return B;
         }
       };
 
       /**
        * Specialisation μ(E, ν)
        */
       template <>
       struct Converter<ElasticModulus::Shear, ElasticModulus::Young,
                        ElasticModulus::Poisson> {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real & E, const Real & nu) {
           return E / (2 * (1 + nu));
         }
       };
 
       /**
        * Specialisation λ(E, ν)
        */
       template <>
       struct Converter<ElasticModulus::lambda, ElasticModulus::Young,
                        ElasticModulus::Poisson> {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real & E, const Real & nu) {
           return E * nu / ((1 + nu) * (1 - 2 * nu));
         }
       };
 
       /**
        * Specialisation K(E, ν)
        */
       template <>
       struct Converter<ElasticModulus::Bulk, ElasticModulus::Young,
                        ElasticModulus::Poisson> {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real & E, const Real & nu) {
           return E / (3 * (1 - 2 * nu));
         }
       };
 
       /**
        * Specialisation E(K, µ)
        */
       template <>
       struct Converter<ElasticModulus::Young, ElasticModulus::Bulk,
                        ElasticModulus::Shear> {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real & K, const Real & G) {
           return 9 * K * G / (3 * K + G);
         }
       };
 
       /**
        * Specialisation ν(K, µ)
        */
       template <>
       struct Converter<ElasticModulus::Poisson, ElasticModulus::Bulk,
                        ElasticModulus::Shear> {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real & K, const Real & G) {
           return (3 * K - 2 * G) / (2 * (3 * K + G));
         }
       };
 
       /**
        * Specialisation E(λ, µ)
        */
       template <>
       struct Converter<ElasticModulus::Young, ElasticModulus::lambda,
                        ElasticModulus::Shear> {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real & lambda,
                                              const Real & G) {
           return G * (3 * lambda + 2 * G) / (lambda + G);
         }
       };
 
       /**
        * Specialisation λ(K, µ)
        */
       template <>
       struct Converter<ElasticModulus::lambda, ElasticModulus::Bulk,
                        ElasticModulus::Shear> {
         //! wrapped function (raison d'être)
         inline constexpr static Real compute(const Real & K, const Real & mu) {
           return K - 2. * mu / 3.;
         }
       };
 
+      /**
+       * Specialisation K(λ, µ)
+       */
+      template <>
+      struct Converter<ElasticModulus::Bulk,
+                       ElasticModulus::lambda,
+                       ElasticModulus::Shear> {
+        //! wrapped function (raison d'être)
+        inline constexpr static Real compute(const Real & lambda,
+                                             const Real & G) {
+          return lambda + (2 * G) / 3;
+        }
+      };
+
     }  // namespace internal
 
     /**
      * allows the conversion from any two distinct input moduli to a
      * chosen output modulus
      */
     template <ElasticModulus Out, ElasticModulus In1, ElasticModulus In2>
     inline constexpr Real convert_elastic_modulus(const Real & in1,
                                                   const Real & in2) {
       // enforcing sanity
       static_assert((In1 != In2),
                     "The input modulus types cannot be identical");
 
       // enforcing independence from order in which moduli are supplied
       constexpr bool inverted{In1 > In2};
       using Converter =
           std::conditional_t<inverted, internal::Converter<Out, In2, In1>,
                              internal::Converter<Out, In1, In2>>;
       if (inverted) {
         return Converter::compute(std::move(in2), std::move(in1));
       } else {
         return Converter::compute(std::move(in1), std::move(in2));
       }
     }
 
     //! static inline implementation of Hooke's law
     template <Dim_t Dim, class Strain_t, class Tangent_t>
     struct Hooke {
       /**
        * compute Lamé's first constant
        * @param young: Young's modulus
        * @param poisson: Poisson's ratio
        */
       inline static constexpr Real compute_lambda(const Real & young,
                                                   const Real & poisson) {
         return convert_elastic_modulus<ElasticModulus::lambda,
                                        ElasticModulus::Young,
                                        ElasticModulus::Poisson>(young, poisson);
       }
 
       /**
        * compute Lamé's second constant (i.e., shear modulus)
        * @param young: Young's modulus
        * @param poisson: Poisson's ratio
        */
       inline static constexpr Real compute_mu(const Real & young,
                                               const Real & poisson) {
         return convert_elastic_modulus<ElasticModulus::Shear,
                                        ElasticModulus::Young,
                                        ElasticModulus::Poisson>(young, poisson);
       }
 
       /**
        * compute the bulk modulus
        * @param young: Young's modulus
        * @param poisson: Poisson's ratio
        */
       inline static constexpr Real compute_K(const Real & young,
                                              const Real & poisson) {
         return convert_elastic_modulus<ElasticModulus::Bulk,
                                        ElasticModulus::Young,
                                        ElasticModulus::Poisson>(young, poisson);
       }
 
       /**
        * compute the stiffness tensor
        * @param lambda: Lamé's first constant
        * @param mu: Lamé's second constant (i.e., shear modulus)
        */
       inline static Eigen::TensorFixedSize<Real,
                                            Eigen::Sizes<Dim, Dim, Dim, Dim>>
       compute_C(const Real & lambda, const Real & mu) {
         return lambda *
                    Tensors::outer<Dim>(Tensors::I2<Dim>(), Tensors::I2<Dim>()) +
                2 * mu * Tensors::I4S<Dim>();
       }
 
       /**
        * compute the stiffness tensor
        * @param lambda: Lamé's first constant
        * @param mu: Lamé's second constant (i.e., shear modulus)
        */
       inline static T4Mat<Real, Dim> compute_C_T4(const Real & lambda,
                                                   const Real & mu) {
         return lambda * Matrices::Itrac<Dim>() +
                2 * mu * Matrices::Isymm<Dim>();
       }
 
       /**
        * return stress
        * @param lambda: First Lamé's constant
        * @param mu: Second Lamé's constant (i.e. shear modulus)
        * @param E: Green-Lagrange or small strain tensor
        */
       template <class s_t>
       inline static decltype(auto) evaluate_stress(const Real & lambda,
                                                    const Real & mu, s_t && E) {
         return E.trace() * lambda * Strain_t::Identity() + 2 * mu * E;
       }
 
       /**
        * return stress and tangent stiffness
        * @param lambda: First Lamé's constant
        * @param mu: Second Lamé's constant (i.e. shear modulus)
        * @param E: Green-Lagrange or small strain tensor
        * @param C: stiffness tensor (Piola-Kirchhoff 2 (or σ) w.r.t to `E`)
        */
       template <class s_t>
       inline static decltype(auto) evaluate_stress(const Real & lambda,
                                                    const Real & mu,
                                                    Tangent_t && C, s_t && E) {
         return std::make_tuple(
             std::move(evaluate_stress(lambda, mu, std::move(E))), std::move(C));
       }
     };
 
     namespace internal {
 
       /* ----------------------------------------------------------------------
        */
       template <Dim_t Dim, FiniteDiff FinDif>
       struct NumericalTangentHelper {
         using T4_t = T4Mat<Real, Dim>;
         using T2_t = Eigen::Matrix<Real, Dim, Dim>;
         using T2_vec = Eigen::Map<Eigen::Matrix<Real, Dim * Dim, 1>>;
 
         template <class FunType, class Derived>
         static inline T4_t compute(FunType && fun,
                                    const Eigen::MatrixBase<Derived> & strain,
                                    Real delta);
       };
 
       /* ----------------------------------------------------------------------
        */
       template <Dim_t Dim, FiniteDiff FinDif>
       template <class FunType, class Derived>
       auto NumericalTangentHelper<Dim, FinDif>::compute(
           FunType && fun, const Eigen::MatrixBase<Derived> & strain, Real delta)
           -> T4_t {
         static_assert((FinDif == FiniteDiff::forward) or
                           (FinDif == FiniteDiff::backward),
                       "Not implemented");
         T4_t tangent{T4_t::Zero()};
 
         const T2_t fun_val{fun(strain)};
         for (Dim_t i{}; i < Dim * Dim; ++i) {
           T2_t strain2{strain};
           T2_vec strain_vec{strain2.data()};
           switch (FinDif) {
           case FiniteDiff::forward: {
             strain_vec(i) += delta;
 
             T2_t del_f_del{(fun(strain2) - fun_val) / delta};
 
             tangent.col(i) = T2_vec(del_f_del.data());
             break;
           }
           case FiniteDiff::backward: {
             strain_vec(i) -= delta;
 
             T2_t del_f_del{(fun_val - fun(strain2)) / delta};
 
             tangent.col(i) = T2_vec(del_f_del.data());
             break;
           }
           }
           static_assert(Int(decltype(tangent.col(i))::SizeAtCompileTime) ==
                             Int(T2_t::SizeAtCompileTime),
                         "wrong column size");
         }
         return tangent;
       }
 
       /* ----------------------------------------------------------------------
        */
       template <Dim_t Dim>
       struct NumericalTangentHelper<Dim, FiniteDiff::centred> {
         using T4_t = T4Mat<Real, Dim>;
         using T2_t = Eigen::Matrix<Real, Dim, Dim>;
         using T2_vec = Eigen::Map<Eigen::Matrix<Real, Dim * Dim, 1>>;
 
         template <class FunType, class Derived>
         static inline T4_t compute(FunType && fun,
                                    const Eigen::MatrixBase<Derived> & strain,
                                    Real delta) {
           T4_t tangent{T4_t::Zero()};
 
           for (Dim_t i{}; i < Dim * Dim; ++i) {
             T2_t strain1{strain};
             T2_t strain2{strain};
             T2_vec strain1_vec{strain1.data()};
             T2_vec strain2_vec{strain2.data()};
             strain1_vec(i) += delta;
             strain2_vec(i) -= delta;
 
             T2_t del_f_del{(fun(strain1).eval() - fun(strain2).eval()) /
                            (2 * delta)};
 
             tangent.col(i) = T2_vec(del_f_del.data());
             static_assert(Int(decltype(tangent.col(i))::SizeAtCompileTime) ==
                               Int(T2_t::SizeAtCompileTime),
                           "wrong column size");
           }
           return tangent;
         }
       };
 
     }  // namespace internal
     /**
      * Helper function to numerically determine tangent, intended for
      * testing, rather than as a replacement for analytical tangents
      */
     template <Dim_t Dim, FiniteDiff FinDif = FiniteDiff::centred, class FunType,
               class Derived>
     inline T4Mat<Real, Dim> compute_numerical_tangent(
         FunType && fun, const Eigen::MatrixBase<Derived> & strain, Real delta) {
       static_assert(Derived::RowsAtCompileTime == Dim,
                     "can't handle dynamic matrix");
       static_assert(Derived::ColsAtCompileTime == Dim,
                     "can't handle dynamic matrix");
 
       using T2_t = Eigen::Matrix<Real, Dim, Dim>;
       using T2_vec = Eigen::Map<Eigen::Matrix<Real, Dim * Dim, 1>>;
 
       static_assert(
           std::is_convertible<FunType, std::function<T2_t(T2_t)>>::value,
           "Function argument 'fun' needs to be a function taking "
           "one second-rank tensor as input and returning a "
           "second-rank tensor");
 
       static_assert(Dim_t(T2_t::SizeAtCompileTime) ==
                         Dim_t(T2_vec::SizeAtCompileTime),
                     "wrong map size");
       return internal::NumericalTangentHelper<Dim, FinDif>::compute(
           std::forward<FunType>(fun), strain, delta);
     }
 
   }  // namespace MatTB
 
 }  // namespace muSpectre
 
 #endif  // SRC_MATERIALS_MATERIALS_TOOLBOX_HH_
diff --git a/src/solver/deprecated_solvers.cc b/src/solver/deprecated_solvers.cc
index 349e07d..eee0b94 100644
--- a/src/solver/deprecated_solvers.cc
+++ b/src/solver/deprecated_solvers.cc
@@ -1,396 +1,397 @@
 /**
  * @file   solvers.cc
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   20 Dec 2017
  *
  * @brief  implementation of solver functions
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU Lesser General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with µSpectre; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * * Boston, MA 02111-1307, USA.
  *
  * Additional permission under GNU GPL version 3 section 7
  *
  * If you modify this Program, or any covered work, by linking or combining it
  * with proprietary FFT implementations or numerical libraries, containing parts
  * covered by the terms of those libraries' licenses, the licensors of this
  * Program grant you additional permission to convey the resulting work.
  */
 
 #include "solver/deprecated_solvers.hh"
 #include "solver/deprecated_solver_cg.hh"
 #include "common/iterators.hh"
 
 #include <Eigen/IterativeLinearSolvers>
 
 #include <iomanip>
 #include <cmath>
 
 namespace muSpectre {
 
   template <Dim_t DimS, Dim_t DimM>
   std::vector<OptimizeResult>
   deprecated_de_geus(CellBase<DimS, DimM> & cell,
                      const GradIncrements<DimM> & delFs,
                      DeprecatedSolverBase<DimS, DimM> & solver, Real newton_tol,
                      Real equil_tol, Dim_t verbose) {
     using Field_t = typename MaterialBase<DimS, DimM>::StrainField_t;
     const Communicator & comm = cell.get_communicator();
     auto solver_fields{std::make_unique<GlobalFieldCollection<DimS>>()};
     solver_fields->initialise(cell.get_subdomain_resolutions(),
                               cell.get_subdomain_locations());
 
     // Corresponds to symbol δF or δε
     auto & incrF{make_field<Field_t>("δF", *solver_fields)};
 
     // Corresponds to symbol ΔF or Δε
     auto & DeltaF{make_field<Field_t>("ΔF", *solver_fields)};
 
     // field to store the rhs for cg calculations
     auto & rhs{make_field<Field_t>("rhs", *solver_fields)};
 
     solver.initialise();
 
     if (solver.get_maxiter() == 0) {
       solver.set_maxiter(cell.size() * DimM * DimM * 10);
     }
 
     size_t count_width{};
     const auto form{cell.get_formulation()};
     std::string strain_symb{};
     if (verbose > 0 && comm.rank() == 0) {
       // setup of algorithm 5.2 in Nocedal, Numerical Optimization (p. 111)
       std::cout << "de Geus-" << solver.name() << " for ";
       switch (form) {
       case Formulation::small_strain: {
         strain_symb = "ε";
         std::cout << "small";
         break;
       }
       case Formulation::finite_strain: {
         strain_symb = "F";
         std::cout << "finite";
         break;
       }
       default:
         throw SolverError("unknown formulation");
         break;
       }
       std::cout << " strain with" << std::endl
                 << "newton_tol = " << newton_tol
                 << ", cg_tol = " << solver.get_tol()
                 << " maxiter = " << solver.get_maxiter() << " and Δ"
                 << strain_symb << " =" << std::endl;
       for (auto && tup : akantu::enumerate(delFs)) {
         auto && counter{std::get<0>(tup)};
         auto && grad{std::get<1>(tup)};
         std::cout << "Step " << counter + 1 << ":" << std::endl
                   << grad << std::endl;
       }
       count_width = size_t(std::log10(solver.get_maxiter())) + 1;
     }
 
     // initialise F = I or ε = 0
     auto & F{cell.get_strain()};
     switch (form) {
     case Formulation::finite_strain: {
       F.get_map() = Matrices::I2<DimM>();
       break;
     }
     case Formulation::small_strain: {
       F.get_map() = Matrices::I2<DimM>().Zero();
       for (const auto & delF : delFs) {
         if (!check_symmetry(delF)) {
           throw SolverError("all Δε must be symmetric!");
         }
       }
       break;
     }
     default:
       throw SolverError("Unknown formulation");
       break;
     }
 
     // initialise return value
     std::vector<OptimizeResult> ret_val{};
 
     Grad_t<DimM> previous_grad{Grad_t<DimM>::Zero()};
     for (const auto & delF : delFs) {  // incremental loop
       std::string message{"Has not converged"};
       Real incrNorm{2 * newton_tol}, gradNorm{1};
       Real stressNorm{2 * equil_tol};
       bool has_converged{false};
       auto convergence_test = [&incrNorm, &gradNorm, &newton_tol, &stressNorm,
                                &equil_tol, &message, &has_converged]() {
         bool incr_test = incrNorm / gradNorm <= newton_tol;
         bool stress_test = stressNorm < equil_tol;
         if (incr_test) {
           message = "Residual  tolerance reached";
         } else if (stress_test) {
           message = "Reached stress divergence tolerance";
         }
         has_converged = incr_test || stress_test;
         return has_converged;
       };
       Uint newt_iter{0};
       for (; (newt_iter < solver.get_maxiter()) &&
              (!has_converged || (newt_iter == 1));
            ++newt_iter) {
         // obtain material response
         auto res_tup{cell.evaluate_stress_tangent(F)};
         auto & P{std::get<0>(res_tup)};
         auto & K{std::get<1>(res_tup)};
 
         auto tangent_effect = [&cell, &K](const Field_t & dF, Field_t & dP) {
           cell.directional_stiffness(K, dF, dP);
         };
 
         if (newt_iter == 0) {
           DeltaF.get_map() = -(delF - previous_grad);  // neg sign because rhs
           tangent_effect(DeltaF, rhs);
           stressNorm = std::sqrt(comm.sum(rhs.eigen().matrix().squaredNorm()));
           if (convergence_test()) {
             break;
           }
           incrF.eigenvec() = solver.solve(rhs.eigenvec(), incrF.eigenvec());
           F.eigen() -= DeltaF.eigen();
         } else {
           rhs.eigen() = -P.eigen();
           cell.project(rhs);
           stressNorm = std::sqrt(comm.sum(rhs.eigen().matrix().squaredNorm()));
           if (convergence_test()) {
             break;
           }
           incrF.eigen() = 0;
           incrF.eigenvec() = solver.solve(rhs.eigenvec(), incrF.eigenvec());
         }
 
         F.eigen() += incrF.eigen();
 
         incrNorm = std::sqrt(comm.sum(incrF.eigen().matrix().squaredNorm()));
         gradNorm = std::sqrt(comm.sum(F.eigen().matrix().squaredNorm()));
         if (verbose > 0 && comm.rank() == 0) {
           std::cout << "at Newton step " << std::setw(count_width) << newt_iter
                     << ", |δ" << strain_symb << "|/|Δ" << strain_symb
                     << "| = " << std::setw(17) << incrNorm / gradNorm
                     << ", tol = " << newton_tol << std::endl;
           if (verbose - 1 > 1) {
             std::cout << "<" << strain_symb << "> =" << std::endl
                       << F.get_map().mean() << std::endl;
           }
         }
         convergence_test();
       }
       // update previous gradient
       previous_grad = delF;
 
       ret_val.push_back(OptimizeResult{
           F.eigen(), cell.get_stress().eigen(), has_converged,
-          Int(has_converged), message, newt_iter, solver.get_counter()});
+          Int(has_converged), message, newt_iter, solver.get_counter(), form});
 
       // store history variables here
       cell.save_history_variables();
     }
 
     return ret_val;
   }
 
   //! instantiation for two-dimensional cells
   template std::vector<OptimizeResult>
   deprecated_de_geus(CellBase<twoD, twoD> & cell,
                      const GradIncrements<twoD> & delF0,
                      DeprecatedSolverBase<twoD, twoD> & solver, Real newton_tol,
                      Real equil_tol, Dim_t verbose);
 
   //! instantiation for three-dimensional cells
   template std::vector<OptimizeResult>
   deprecated_de_geus(CellBase<threeD, threeD> & cell,
                      const GradIncrements<threeD> & delF0,
                      DeprecatedSolverBase<threeD, threeD> & solver,
                      Real newton_tol, Real equil_tol, Dim_t verbose);
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   std::vector<OptimizeResult>
   deprecated_newton_cg(CellBase<DimS, DimM> & cell,
                        const GradIncrements<DimM> & delFs,
                        DeprecatedSolverBase<DimS, DimM> & solver,
                        Real newton_tol, Real equil_tol, Dim_t verbose) {
     using Field_t = typename MaterialBase<DimS, DimM>::StrainField_t;
     const Communicator & comm = cell.get_communicator();
     auto solver_fields{std::make_unique<GlobalFieldCollection<DimS>>()};
     solver_fields->initialise(cell.get_subdomain_resolutions(),
                               cell.get_subdomain_locations());
 
     // Corresponds to symbol δF or δε
     auto & incrF{make_field<Field_t>("δF", *solver_fields)};
 
     // field to store the rhs for cg calculations
     auto & rhs{make_field<Field_t>("rhs", *solver_fields)};
 
     solver.initialise();
 
     if (solver.get_maxiter() == 0) {
       solver.set_maxiter(cell.size() * DimM * DimM * 10);
     }
 
     size_t count_width{};
     const auto form{cell.get_formulation()};
     std::string strain_symb{};
     if (verbose > 0 && comm.rank() == 0) {
       // setup of algorithm 5.2 in Nocedal, Numerical Optimization (p. 111)
       std::cout << "Newton-" << solver.name() << " for ";
       switch (form) {
       case Formulation::small_strain: {
         strain_symb = "ε";
         std::cout << "small";
         break;
       }
       case Formulation::finite_strain: {
         strain_symb = "F";
         std::cout << "finite";
         break;
       }
       default:
         throw SolverError("unknown formulation");
         break;
       }
       std::cout << " strain with" << std::endl
                 << "newton_tol = " << newton_tol
                 << ", cg_tol = " << solver.get_tol()
                 << " maxiter = " << solver.get_maxiter() << " and Δ"
                 << strain_symb << " =" << std::endl;
       for (auto && tup : akantu::enumerate(delFs)) {
         auto && counter{std::get<0>(tup)};
         auto && grad{std::get<1>(tup)};
         std::cout << "Step " << counter + 1 << ":" << std::endl
                   << grad << std::endl;
       }
       count_width = size_t(std::log10(solver.get_maxiter())) + 1;
     }
 
     // initialise F = I or ε = 0
     auto & F{cell.get_strain()};
     switch (cell.get_formulation()) {
     case Formulation::finite_strain: {
       F.get_map() = Matrices::I2<DimM>();
       break;
     }
     case Formulation::small_strain: {
       F.get_map() = Matrices::I2<DimM>().Zero();
       for (const auto & delF : delFs) {
         if (!check_symmetry(delF)) {
           throw SolverError("all Δε must be symmetric!");
         }
       }
       break;
     }
     default:
       throw SolverError("Unknown formulation");
       break;
     }
 
     // initialise return value
     std::vector<OptimizeResult> ret_val{};
 
     Grad_t<DimM> previous_grad{Grad_t<DimM>::Zero()};
     for (const auto & delF : delFs) {  // incremental loop
       // apply macroscopic strain increment
       for (auto && grad : F.get_map()) {
         grad += delF - previous_grad;
       }
 
       std::string message{"Has not converged"};
       Real incrNorm{2 * newton_tol}, gradNorm{1};
       Real stressNorm{2 * equil_tol};
       bool has_converged{false};
       auto convergence_test = [&incrNorm, &gradNorm, &newton_tol, &stressNorm,
                                &equil_tol, &message, &has_converged]() {
         bool incr_test = incrNorm / gradNorm <= newton_tol;
         bool stress_test = stressNorm < equil_tol;
         if (incr_test) {
           message = "Residual  tolerance reached";
         } else if (stress_test) {
           message = "Reached stress divergence tolerance";
         }
         has_converged = incr_test || stress_test;
         return has_converged;
       };
       Uint newt_iter{0};
 
       for (; newt_iter < solver.get_maxiter() && !has_converged; ++newt_iter) {
         // obtain material response
         auto res_tup{cell.evaluate_stress_tangent(F)};
         auto & P{std::get<0>(res_tup)};
 
         rhs.eigen() = -P.eigen();
         cell.project(rhs);
         stressNorm = std::sqrt(comm.sum(rhs.eigen().matrix().squaredNorm()));
         if (convergence_test()) {
           break;
         }
         incrF.eigen() = 0;
 
         incrF.eigenvec() = solver.solve(rhs.eigenvec(), incrF.eigenvec());
 
         F.eigen() += incrF.eigen();
 
         incrNorm = std::sqrt(comm.sum(incrF.eigen().matrix().squaredNorm()));
         gradNorm = std::sqrt(comm.sum(F.eigen().matrix().squaredNorm()));
         if (verbose > 0 && comm.rank() == 0) {
           std::cout << "at Newton step " << std::setw(count_width) << newt_iter
                     << ", |δ" << strain_symb << "|/|Δ" << strain_symb
                     << "| = " << std::setw(17) << incrNorm / gradNorm
                     << ", tol = " << newton_tol << std::endl;
 
           if (verbose - 1 > 1) {
             std::cout << "<" << strain_symb << "> =" << std::endl
                       << F.get_map().mean() << std::endl;
           }
         }
         convergence_test();
       }
       // update previous gradient
       previous_grad = delF;
 
-      ret_val.push_back(OptimizeResult{
-          F.eigen(), cell.get_stress().eigen(), convergence_test(),
-          Int(convergence_test()), message, newt_iter, solver.get_counter()});
+      ret_val.push_back(OptimizeResult{F.eigen(), cell.get_stress().eigen(),
+                                       convergence_test(),
+                                       Int(convergence_test()), message,
+                                       newt_iter, solver.get_counter(), form});
 
       // store history variables for next load increment
       cell.save_history_variables();
     }
 
     return ret_val;
   }
 
   //! instantiation for two-dimensional cells
   template std::vector<OptimizeResult>
   deprecated_newton_cg(CellBase<twoD, twoD> & cell,
                        const GradIncrements<twoD> & delF0,
                        DeprecatedSolverBase<twoD, twoD> & solver,
                        Real newton_tol, Real equil_tol, Dim_t verbose);
 
   //! instantiation for three-dimensional cells
   template std::vector<OptimizeResult>
   deprecated_newton_cg(CellBase<threeD, threeD> & cell,
                        const GradIncrements<threeD> & delF0,
                        DeprecatedSolverBase<threeD, threeD> & solver,
                        Real newton_tol, Real equil_tol, Dim_t verbose);
 
 }  // namespace muSpectre
diff --git a/src/solver/solver_common.hh b/src/solver/solver_common.hh
index 812278b..6744387 100644
--- a/src/solver/solver_common.hh
+++ b/src/solver/solver_common.hh
@@ -1,101 +1,103 @@
 /**
  * @file   solver_common.hh
  *
  * @author Till Junge <till.junge@altermail.ch>
  *
  * @date   28 Dec 2017
  *
  * @brief  Errors raised by solvers and other common utilities
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU Lesser General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with µSpectre; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * * Boston, MA 02111-1307, USA.
  *
  * Additional permission under GNU GPL version 3 section 7
  *
  * If you modify this Program, or any covered work, by linking or combining it
  * with proprietary FFT implementations or numerical libraries, containing parts
  * covered by the terms of those libraries' licenses, the licensors of this
  * Program grant you additional permission to convey the resulting work.
  */
 
 #ifndef SRC_SOLVER_SOLVER_COMMON_HH_
 #define SRC_SOLVER_SOLVER_COMMON_HH_
 
 #include "common/common.hh"
 #include "common/tensor_algebra.hh"
 
 #include <Eigen/Dense>
 
 #include <stdexcept>
 
 namespace muSpectre {
 
   /**
    * emulates scipy.optimize.OptimizeResult
    */
   struct OptimizeResult {
     //! Strain ε or Gradient F at solution
     Eigen::ArrayXXd grad;
     //! Cauchy stress σ or first Piola-Kirchhoff stress P at solution
     Eigen::ArrayXXd stress;
     //! whether or not the solver exited successfully
     bool success;
     //! Termination status of the optimizer. Its value depends on the
     //! underlying solver. Refer to message for details.
     Int status;
     //! Description of the cause of the termination.
     std::string message;
     //! number of iterations
     Uint nb_it;
     //! number of cell evaluations
     Uint nb_fev;
+    //! continuum mechanic flag
+    Formulation formulation;
   };
 
   /**
    * Field type that solvers expect gradients to be expressed in
    */
   template <Dim_t Dim>
   using Grad_t = Matrices::Tens2_t<Dim>;
   /**
    * multiple increments can be submitted at once (useful for
    * path-dependent materials)
    */
   template <Dim_t Dim>
   using GradIncrements =
       std::vector<Grad_t<Dim>, Eigen::aligned_allocator<Grad_t<Dim>>>;
 
   /* ---------------------------------------------------------------------- */
   class SolverError : public std::runtime_error {
     using runtime_error::runtime_error;
   };
 
   /* ---------------------------------------------------------------------- */
   class ConvergenceError : public SolverError {
     using SolverError::SolverError;
   };
 
   /* ---------------------------------------------------------------------- */
   /**
    * check whether a strain is symmetric, for the purposes of small
    * strain problems
    */
   bool check_symmetry(const Eigen::Ref<const Eigen::ArrayXXd> & eps,
                       Real rel_tol = 1e-8);
 
 }  // namespace muSpectre
 
 #endif  // SRC_SOLVER_SOLVER_COMMON_HH_
diff --git a/src/solver/solvers.cc b/src/solver/solvers.cc
index cd072e8..56f01b0 100644
--- a/src/solver/solvers.cc
+++ b/src/solver/solvers.cc
@@ -1,480 +1,482 @@
 /**
  * file   solvers.cc
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   24 Apr 2018
  *
  * @brief  implementation of dynamic newton-cg solver
  *
  * Copyright © 2018 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU Lesser General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with µSpectre; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * * Boston, MA 02111-1307, USA.
  *
  * Additional permission under GNU GPL version 3 section 7
  *
  * If you modify this Program, or any covered work, by linking or combining it
  * with proprietary FFT implementations or numerical libraries, containing parts
  * covered by the terms of those libraries' licenses, the licensors of this
  * Program grant you additional permission to convey the resulting work.
  */
 
 #include "solver/solvers.hh"
 #include "common/iterators.hh"
 
 #include <Eigen/Dense>
 
 #include <iomanip>
 
 namespace muSpectre {
 
   Eigen::IOFormat format(Eigen::FullPrecision, 0, ", ", ",\n", "[", "]", "[",
                          "]");
 
   //--------------------------------------------------------------------------//
   std::vector<OptimizeResult> newton_cg(Cell & cell,
                                         const LoadSteps_t & load_steps,
                                         SolverBase & solver, Real newton_tol,
                                         Real equil_tol, Dim_t verbose) {
     const Communicator & comm = cell.get_communicator();
 
     using Vector_t = Eigen::Matrix<Real, Eigen::Dynamic, 1>;
     using Matrix_t = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>;
 
     // Corresponds to symbol δF or δε
     Vector_t incrF(cell.get_nb_dof());
 
     // field to store the rhs for cg calculations
     Vector_t rhs(cell.get_nb_dof());
 
     solver.initialise();
 
     size_t count_width{};
     const auto form{cell.get_formulation()};
     std::string strain_symb{};
     if (verbose > 0 && comm.rank() == 0) {
       // setup of algorithm 5.2 in Nocedal, Numerical Optimization (p. 111)
       std::cout << "Newton-" << solver.get_name() << " for ";
       switch (form) {
       case Formulation::small_strain: {
         strain_symb = "ε";
         std::cout << "small";
         break;
       }
       case Formulation::finite_strain: {
         strain_symb = "F";
         std::cout << "finite";
         break;
       }
       default:
         throw SolverError("unknown formulation");
         break;
       }
       std::cout << " strain with" << std::endl
                 << "newton_tol = " << newton_tol
                 << ", cg_tol = " << solver.get_tol()
                 << " maxiter = " << solver.get_maxiter() << " and "
                 << strain_symb << " from " << strain_symb << "₁ =" << std::endl
                 << load_steps.front() << std::endl
                 << " to " << strain_symb << "ₙ =" << std::endl
                 << load_steps.back() << std::endl
                 << "in increments of Δ" << strain_symb << " =" << std::endl
                 << (load_steps.back() - load_steps.front()) / load_steps.size()
                 << std::endl;
       count_width = size_t(std::log10(solver.get_maxiter())) + 1;
     }
 
     auto shape{cell.get_strain_shape()};
 
     switch (form) {
     case Formulation::finite_strain: {
       cell.set_uniform_strain(Matrix_t::Identity(shape[0], shape[1]));
       for (const auto & delF : load_steps) {
         if (not((delF.rows() == shape[0]) and (delF.cols() == shape[1]))) {
           std::stringstream err{};
           err << "Load increments need to be given in " << shape[0] << "×"
               << shape[1] << " matrices, but I got a " << delF.rows() << "×"
               << delF.cols() << " matrix:" << std::endl
               << delF;
           throw SolverError(err.str());
         }
       }
       break;
     }
     case Formulation::small_strain: {
       cell.set_uniform_strain(Matrix_t::Zero(shape[0], shape[1]));
       for (const auto & delF : load_steps) {
         if (not((delF.rows() == shape[0]) and (delF.cols() == shape[1]))) {
           std::stringstream err{};
           err << "Load increments need to be given in " << shape[0] << "×"
               << shape[1] << " matrices, but I got a " << delF.rows() << "×"
               << delF.cols() << " matrix:" << std::endl
               << delF;
           throw SolverError(err.str());
         }
         if (not check_symmetry(delF)) {
           throw SolverError("all Δε must be symmetric!");
         }
       }
       break;
     }
     default:
       throw SolverError("Unknown strain measure");
       break;
     }
 
     // initialise return value
     std::vector<OptimizeResult> ret_val{};
 
     // storage for the previous mean strain (to compute ΔF or Δε)
     Matrix_t previous_macro_strain{load_steps.back().Zero(shape[0], shape[1])};
 
     auto F{cell.get_strain_vector()};
     //! incremental loop
     for (const auto & tup : akantu::enumerate(load_steps)) {
       const auto & strain_step{std::get<0>(tup)};
       const auto & macro_strain{std::get<1>(tup)};
       if ((verbose > 0) and (comm.rank() == 0)) {
         std::cout << "at Load step " << std::setw(count_width)
                   << strain_step + 1 << std::endl;
       }
       using StrainMap_t = RawFieldMap<Eigen::Map<Eigen::MatrixXd>>;
       for (auto && strain : StrainMap_t(F, shape[0], shape[1])) {
         strain += macro_strain - previous_macro_strain;
       }
 
       std::string message{"Has not converged"};
       Real incr_norm{2 * newton_tol}, grad_norm{1};
       Real stress_norm{2 * equil_tol};
       bool has_converged{false};
 
       auto convergence_test = [&incr_norm, &grad_norm, &newton_tol,
                                &stress_norm, &equil_tol, &message,
                                &has_converged]() {
         bool incr_test = incr_norm / grad_norm <= newton_tol;
         bool stress_test = stress_norm < equil_tol;
         if (incr_test) {
           message = "Residual  tolerance reached";
         } else if (stress_test) {
           message = "Reached stress divergence tolerance";
         }
         has_converged = incr_test || stress_test;
         return has_converged;
       };
       Uint newt_iter{0};
 
       for (; newt_iter < solver.get_maxiter() && !has_converged; ++newt_iter) {
         // obtain material response
         auto res_tup{cell.evaluate_stress_tangent()};
         auto & P{std::get<0>(res_tup)};
 
         rhs = -P;
         cell.apply_projection(rhs);
         stress_norm = std::sqrt(comm.sum(rhs.squaredNorm()));
 
         if (convergence_test()) {
           break;
         }
 
         try {
           incrF = solver.solve(rhs);
         } catch (ConvergenceError & error) {
           std::stringstream err{};
           err << "Failure at load step " << strain_step + 1 << " of "
               << load_steps.size() << ". In Newton-Raphson step " << newt_iter
               << ":" << std::endl
               << error.what() << std::endl
               << "The applied boundary condition is Δ" << strain_symb << " ="
               << std::endl
               << macro_strain << std::endl
               << "and the load increment is Δ" << strain_symb << " ="
               << std::endl
               << macro_strain - previous_macro_strain << std::endl;
           throw ConvergenceError(err.str());
         }
 
         F += incrF;
 
         incr_norm = std::sqrt(comm.sum(incrF.squaredNorm()));
         grad_norm = std::sqrt(comm.sum(F.squaredNorm()));
 
         if ((verbose > 1) and (comm.rank() == 0)) {
           std::cout << "at Newton step " << std::setw(count_width) << newt_iter
                     << ", |δ" << strain_symb << "|/|Δ" << strain_symb
                     << "| = " << std::setw(17) << incr_norm / grad_norm
                     << ", tol = " << newton_tol << std::endl;
 
           if (verbose - 1 > 1) {
             std::cout << "<" << strain_symb << "> =" << std::endl
                       << StrainMap_t(F, shape[0], shape[1]).mean() << std::endl;
           }
         }
         convergence_test();
       }
       if (newt_iter == solver.get_maxiter()) {
         std::stringstream err{};
         err << "Failure at load step " << strain_step + 1 << " of "
             << load_steps.size() << ". Newton-Raphson failed to converge. "
             << "The applied boundary condition is Δ" << strain_symb << " ="
             << std::endl
             << macro_strain << std::endl
             << "and the load increment is Δ" << strain_symb << " =" << std::endl
             << macro_strain - previous_macro_strain << std::endl;
         throw ConvergenceError(err.str());
       }
 
       // update previous macroscopic strain
       previous_macro_strain = macro_strain;
 
       // store results
-      ret_val.emplace_back(OptimizeResult{
-          F, cell.get_stress_vector(), convergence_test(),
-          Int(convergence_test()), message, newt_iter, solver.get_counter()});
+      ret_val.emplace_back(
+          OptimizeResult{F, cell.get_stress_vector(), convergence_test(),
+                         Int(convergence_test()), message, newt_iter,
+                         solver.get_counter(), form});
 
       // store history variables for next load increment
       cell.save_history_variables();
     }
 
     return ret_val;
   }
 
   //----------------------------------------------------------------------------//
   std::vector<OptimizeResult> de_geus(Cell & cell,
                                       const LoadSteps_t & load_steps,
                                       SolverBase & solver, Real newton_tol,
                                       Real equil_tol, Dim_t verbose) {
     const Communicator & comm = cell.get_communicator();
 
     using Vector_t = Eigen::Matrix<Real, Eigen::Dynamic, 1>;
     using Matrix_t = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>;
 
     // Corresponds to symbol δF or δε
     Vector_t incrF(cell.get_nb_dof());
 
     // Corresponds to symbol ΔF or Δε
     Vector_t DeltaF(cell.get_nb_dof());
 
     // field to store the rhs for cg calculations
     Vector_t rhs(cell.get_nb_dof());
 
     solver.initialise();
 
     size_t count_width{};
     const auto form{cell.get_formulation()};
     std::string strain_symb{};
     if (verbose > 0 && comm.rank() == 0) {
       // setup of algorithm 5.2 in Nocedal, Numerical Optimization (p. 111)
       std::cout << "de Geus-" << solver.get_name() << " for ";
       switch (form) {
       case Formulation::small_strain: {
         strain_symb = "ε";
         std::cout << "small";
         break;
       }
       case Formulation::finite_strain: {
         strain_symb = "F";
         std::cout << "finite";
         break;
       }
       default:
         throw SolverError("unknown formulation");
         break;
       }
       std::cout << " strain with" << std::endl
                 << "newton_tol = " << newton_tol
                 << ", cg_tol = " << solver.get_tol()
                 << " maxiter = " << solver.get_maxiter() << " and "
                 << strain_symb << " from " << strain_symb << "₁ =" << std::endl
                 << load_steps.front() << std::endl
                 << " to " << strain_symb << "ₙ =" << std::endl
                 << load_steps.back() << std::endl
                 << "in increments of Δ" << strain_symb << " =" << std::endl
                 << (load_steps.back() - load_steps.front()) / load_steps.size()
                 << std::endl;
       count_width = size_t(std::log10(solver.get_maxiter())) + 1;
     }
 
     auto shape{cell.get_strain_shape()};
     Matrix_t default_strain_val{};
 
     switch (form) {
     case Formulation::finite_strain: {
       cell.set_uniform_strain(Matrix_t::Identity(shape[0], shape[1]));
       for (const auto & delF : load_steps) {
         auto rows = delF.rows();
         auto cols = delF.cols();
         if (not((rows == shape[0]) and (cols == shape[1]))) {
           std::stringstream err{};
           err << "Load increments need to be given in " << shape[0] << "×"
               << shape[1] << " matrices, but I got a " << delF.rows() << "×"
               << delF.cols() << " matrix:" << std::endl
               << delF;
           throw SolverError(err.str());
         }
       }
       break;
     }
     case Formulation::small_strain: {
       cell.set_uniform_strain(Matrix_t::Zero(shape[0], shape[1]));
       for (const auto & delF : load_steps) {
         if (not((delF.rows() == shape[0]) and (delF.cols() == shape[1]))) {
           std::stringstream err{};
           err << "Load increments need to be given in " << shape[0] << "×"
               << shape[1] << " matrices, but I got a " << delF.rows() << "×"
               << delF.cols() << " matrix:" << std::endl
               << delF;
           throw SolverError(err.str());
         }
         if (not check_symmetry(delF)) {
           throw SolverError("all Δε must be symmetric!");
         }
       }
       break;
     }
     default:
       throw SolverError("Unknown strain measure");
       break;
     }
 
     // initialise return value
     std::vector<OptimizeResult> ret_val{};
 
     // storage for the previous mean strain (to compute ΔF or Δε)
     Matrix_t previous_macro_strain{load_steps.back().Zero(shape[0], shape[1])};
 
     auto F{cell.get_strain_vector()};
     //! incremental loop
     for (const auto & tup : akantu::enumerate(load_steps)) {
       const auto & strain_step{std::get<0>(tup)};
       const auto & macro_strain{std::get<1>(tup)};
       using StrainMap_t = RawFieldMap<Eigen::Map<Eigen::MatrixXd>>;
       if ((verbose > 0) and (comm.rank() == 0)) {
         std::cout << "at Load step " << std::setw(count_width)
                   << strain_step + 1 << ", " << strain_symb << " =" << std::endl
                   << (macro_strain + default_strain_val).format(format)
                   << std::endl;
       }
 
       std::string message{"Has not converged"};
       Real incr_norm{2 * newton_tol}, grad_norm{1};
       Real stress_norm{2 * equil_tol};
       bool has_converged{false};
 
       auto convergence_test = [&incr_norm, &grad_norm, &newton_tol,
                                &stress_norm, &equil_tol, &message,
                                &has_converged]() {
         bool incr_test = incr_norm / grad_norm <= newton_tol;
         bool stress_test = stress_norm < equil_tol;
         if (incr_test) {
           message = "Residual  tolerance reached";
         } else if (stress_test) {
           message = "Reached stress divergence tolerance";
         }
         has_converged = incr_test || stress_test;
         return has_converged;
       };
       Uint newt_iter{0};
 
       for (; ((newt_iter < solver.get_maxiter()) and (!has_converged)) or
              (newt_iter < 2);
            ++newt_iter) {
         // obtain material response
         auto res_tup{cell.evaluate_stress_tangent()};
         auto & P{std::get<0>(res_tup)};
 
         try {
           if (newt_iter == 0) {
             for (auto && strain : StrainMap_t(DeltaF, shape[0], shape[1])) {
               strain = macro_strain - previous_macro_strain;
             }
             rhs = -cell.evaluate_projected_directional_stiffness(DeltaF);
             F += DeltaF;
             stress_norm = std::sqrt(comm.sum(rhs.matrix().squaredNorm()));
             if (stress_norm < equil_tol) {
               incrF.setZero();
             } else {
               incrF = solver.solve(rhs);
             }
           } else {
             rhs = -P;
             cell.apply_projection(rhs);
             stress_norm = std::sqrt(comm.sum(rhs.matrix().squaredNorm()));
             if (stress_norm < equil_tol) {
               incrF.setZero();
             } else {
               incrF = solver.solve(rhs);
             }
           }
         } catch (ConvergenceError & error) {
           std::stringstream err{};
           err << "Failure at load step " << strain_step + 1 << " of "
               << load_steps.size() << ". In Newton-Raphson step " << newt_iter
               << ":" << std::endl
               << error.what() << std::endl
               << "The applied boundary condition is Δ" << strain_symb << " ="
               << std::endl
               << macro_strain << std::endl
               << "and the load increment is Δ" << strain_symb << " ="
               << std::endl
               << macro_strain - previous_macro_strain << std::endl;
           throw ConvergenceError(err.str());
         }
 
         F += incrF;
 
         incr_norm = std::sqrt(comm.sum(incrF.squaredNorm()));
         grad_norm = std::sqrt(comm.sum(F.squaredNorm()));
 
         if ((verbose > 0) and (comm.rank() == 0)) {
           std::cout << "at Newton step " << std::setw(count_width) << newt_iter
                     << ", |δ" << strain_symb << "|/|Δ" << strain_symb
                     << "| = " << std::setw(17) << incr_norm / grad_norm
                     << ", tol = " << newton_tol << std::endl;
 
           if (verbose - 1 > 1) {
             std::cout << "<" << strain_symb << "> =" << std::endl
                       << StrainMap_t(F, shape[0], shape[1]).mean() << std::endl;
           }
         }
         convergence_test();
       }
       if (newt_iter == solver.get_maxiter()) {
         std::stringstream err{};
         err << "Failure at load step " << strain_step + 1 << " of "
             << load_steps.size() << ". Newton-Raphson failed to converge. "
             << "The applied boundary condition is Δ" << strain_symb << " ="
             << std::endl
             << macro_strain << std::endl
             << "and the load increment is Δ" << strain_symb << " =" << std::endl
             << macro_strain - previous_macro_strain << std::endl;
         throw ConvergenceError(err.str());
       }
 
       // update previous macroscopic strain
       previous_macro_strain = macro_strain;
 
       // store results
-      ret_val.emplace_back(OptimizeResult{
-          F, cell.get_stress_vector(), convergence_test(),
-          Int(convergence_test()), message, newt_iter, solver.get_counter()});
+      ret_val.emplace_back(
+          OptimizeResult{F, cell.get_stress_vector(), convergence_test(),
+                         Int(convergence_test()), message, newt_iter,
+                         solver.get_counter(), form});
 
       // store history variables for next load increment
       cell.save_history_variables();
     }
 
     return ret_val;
   }
 }  // namespace muSpectre
diff --git a/tests/py_comparison_test_material_linear_elastic1.cc b/tests/py_comparison_test_material_linear_elastic1.cc
deleted file mode 100644
index ce3c1ea..0000000
--- a/tests/py_comparison_test_material_linear_elastic1.cc
+++ /dev/null
@@ -1,108 +0,0 @@
-/**
- * @file   py_comparison_test_material_linear_elastic_1.cc
- *
- * @author Till Junge <till.junge@epfl.ch>
- *
- * @date   05 Dec 2018
- *
- * @brief  simple wrapper around the MaterialLinearElastic1 to test it
- *         with arbitrary input
- *
- * Copyright © 2018 Till Junge
- *
- * µSpectre is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License as
- * published by the Free Software Foundation, either version 3, or (at
- * your option) any later version.
- *
- * µSpectre is distributed in the hope that it will be useful, but
- * WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with µSpectre; see the file COPYING. If not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA 02111-1307, USA.
- *
- * Additional permission under GNU GPL version 3 section 7
- *
- * If you modify this Program, or any covered work, by linking or combining it
- * with proprietary FFT implementations or numerical libraries, containing parts
- * covered by the terms of those libraries' licenses, the licensors of this
- * Program grant you additional permission to convey the resulting work.
- */
-
-#include "materials/stress_transformations_PK1.hh"
-#include "materials/material_linear_elastic1.hh"
-#include "materials/materials_toolbox.hh"
-
-#include <pybind11/pybind11.h>
-#include <pybind11/eigen.h>
-
-#include <tuple>
-
-using pybind11::literals::operator""_a;
-namespace py = pybind11;
-
-namespace muSpectre {
-  template <Dim_t Dim>
-  std::tuple<Eigen::Matrix<Real, Dim, Dim>,
-             Eigen::Matrix<Real, Dim * Dim, Dim * Dim>>
-  material_wrapper(Real Young, Real Poisson,
-                   py::EigenDRef<const Eigen::MatrixXd> F) {
-    using Mat_t = MaterialLinearElastic1<Dim, Dim>;
-    Mat_t mat("Name", Young, Poisson);
-
-    auto & coll{mat.get_collection()};
-    coll.add_pixel({0});
-    coll.initialise();
-
-    Eigen::Matrix<Real, Dim, Dim> F_mat(F);
-    Eigen::Matrix<Real, Dim, Dim> E{
-        MatTB::convert_strain<StrainMeasure::Gradient,
-                              StrainMeasure::GreenLagrange>(F_mat)};
-    return mat.evaluate_stress_tangent(std::move(E));
-  }
-
-  template <Dim_t Dim>
-  py::tuple PK2_fun(Real Young, Real Poisson,
-                    py::EigenDRef<const Eigen::MatrixXd> F) {
-    auto && tup{material_wrapper<Dim>(Young, Poisson, F)};
-    auto && S{std::get<0>(tup)};
-    Eigen::MatrixXd C{std::get<1>(tup)};
-
-    return py::make_tuple(Eigen::MatrixXd{S}, C);
-  }
-
-  template <Dim_t Dim>
-  py::tuple PK1_fun(Real Young, Real Poisson,
-                    py::EigenDRef<const Eigen::MatrixXd> F) {
-    auto && tup{material_wrapper<Dim>(Young, Poisson, F)};
-    auto && S{std::get<0>(tup)};
-    auto && C{std::get<1>(tup)};
-
-    using Mat_t = MaterialLinearElastic1<Dim, Dim>;
-    constexpr auto StrainM{Mat_t::traits::strain_measure};
-    constexpr auto StressM{Mat_t::traits::stress_measure};
-
-    auto && PK_tup{MatTB::PK1_stress<StressM, StrainM>(
-        Eigen::Matrix<Real, Dim, Dim>{F}, S, C)};
-    auto && P{std::get<0>(PK_tup)};
-    auto && K{std::get<1>(PK_tup)};
-
-    return py::make_tuple(std::move(P), std::move(K));
-  }
-
-  PYBIND11_MODULE(material_linear_elastic1, mod) {
-    mod.doc() = "Comparison provider for MaterialLinearelastic1";
-    auto adder{[&mod](auto name, auto & fun) {
-      mod.def(name, &fun, "Young"_a, "Poisson"_a, "F"_a);
-    }};
-    adder("PK2_fun_2d", PK2_fun<twoD>);
-    adder("PK2_fun_3d", PK2_fun<threeD>);
-    adder("PK1_fun_2d", PK1_fun<twoD>);
-    adder("PK1_fun_3d", PK1_fun<threeD>);
-  }
-
-}  // namespace muSpectre
diff --git a/tests/py_comparison_test_material_linear_elastic1.py b/tests/py_comparison_test_material_linear_elastic1.py
deleted file mode 100644
index cb59e8c..0000000
--- a/tests/py_comparison_test_material_linear_elastic1.py
+++ /dev/null
@@ -1,149 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding:utf-8 -*-
-"""
-@file   py_comparison_test_material_linear_elastic_1.py
-
-@author Till Junge <till.junge@epfl.ch>
-
-@date   05 Dec 2018
-
-@brief  compares MaterialLinearElastic1 to de Geus's python implementation
-
-Copyright © 2018 Till Junge
-
-µSpectre is free software; you can redistribute it and/or
-modify it under the terms of the GNU General Lesser Public License as
-published by the Free Software Foundation, either version 3, or (at
-your option) any later version.
-
-µSpectre is distributed in the hope that it will be useful, but
-WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-General Public License for more details.
-
-You should have received a copy of the GNU Lesser General Public License
-along with µSpectre; see the file COPYING. If not, write to the
-Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-Boston, MA 02111-1307, USA.
-
-Additional permission under GNU GPL version 3 section 7
-
-If you modify this Program, or any covered work, by linking or combining it
-with proprietary FFT implementations or numerical libraries, containing parts
-covered by the terms of those libraries' licenses, the licensors of this
-Program grant you additional permission to convey the resulting work.
-"""
-
-from  material_linear_elastic1 import *
-import itertools
-import numpy as np
-np.set_printoptions(linewidth=180)
-
-import unittest
-
-#####################
-dyad22 = lambda A2,B2: np.einsum('ij  ,kl  ->ijkl',A2,B2)
-dyad11 = lambda A1,B1: np.einsum('i   ,j   ->ij  ',A1,B1)
-dot22  = lambda A2,B2: np.einsum('ij  ,jk  ->ik  ',A2,B2)
-dot24  = lambda A2,B4: np.einsum('ij  ,jkmn->ikmn',A2,B4)
-dot42  = lambda A4,B2: np.einsum('ijkl,lm  ->ijkm',A4,B2)
-inv2 = np.linalg.inv
-trans2 = np.transpose
-ddot22 = lambda A2,B2: np.einsum('ij  ,ji  ->    ',A2,B2)
-ddot42 = lambda A4,B2: np.einsum('ijkl,lk  ->ij  ',A4,B2)
-ddot44 = lambda A4,B4: np.einsum('ijkl,lkmn->ijmn',A4,B4)
-
-class MatTest(unittest.TestCase):
-    def constitutive(self, F, dim):
-        I  = np.eye(dim)
-        II = dyad22(I,I)
-        I4 = np.einsum('il,jk',I,I)
-        I4rt = np.einsum('ik,jl',I,I)
-        I4s  = (I4+I4rt)/2.
-
-        C4 = self.K*II+2.*self.mu*(I4s-1./3.*II)
-        S  = ddot42(C4,.5*(dot22(trans2(F),F)-I))
-        P  = dot22(F,S)
-        K4 = dot24(S,I4)+ddot44(ddot44(I4rt,dot42(dot24(F,C4),trans2(F))),I4rt)
-        return P, S, K4, C4
-
-    def setUp(self):
-        pass
-    def prep(self, dimension):
-        self.dim=dimension
-        self.Young = 200e9+100*np.random.rand()
-        self.Poisson = .3 + .1*(np.random.rand()-.5)
-        self.K = self.Young/(3*(1-2*self.Poisson))
-        self.mu = self.Young/(2*(1+self.Poisson))
-        self.F = np.eye(self.dim) + (np.random.random((self.dim, self.dim))-.5)/10
-        self.tol = 1e-13
-        self.verbose=True
-
-    def test_equivalence_S_C(self):
-        for dim in (2, 3):
-            self.runner_equivalence_S_C(dim)
-
-    def runner_equivalence_S_C(self, dimension):
-        self.prep(dimension)
-        fun = PK2_fun_2d if self.dim == 2 else PK2_fun_3d
-        S_µ, C_µ_s = fun(self.Young, self.Poisson, self.F)
-        shape = (self.dim, self.dim, self.dim, self.dim)
-        C_µ = C_µ_s.reshape(shape).transpose((0,1,3,2))
-
-        response_p = self.constitutive(self.F, self.dim)
-        S_p, C_p = response_p[1], response_p[3]
-
-        S_error = np.linalg.norm(S_µ- S_p)/np.linalg.norm(S_µ)
-        if not S_error < self.tol:
-            print("Error(S) = {}".format(S_error))
-            print("S_µ:\n{}".format(S_µ))
-            print("S_p:\n{}".format(S_p))
-        self.assertLess(S_error,
-                        self.tol)
-
-        C_error = np.linalg.norm(C_µ- C_p)/np.linalg.norm(C_µ)
-        if not C_error < self.tol:
-            print("Error(C) = {}".format(C_error))
-            flat_shape = (self.dim**2, self.dim**2)
-            print("C_µ:\n{}".format(C_µ.reshape(flat_shape)))
-            print("C_p:\n{}".format(C_p.reshape(flat_shape)))
-        self.assertLess(C_error,
-                        self.tol)
-
-    def test_equivalence_P_K(self):
-        for dim in (2, 3):
-            self.runner_equivalence_P_K(dim)
-
-    def runner_equivalence_P_K(self, dimension):
-        self.prep(dimension)
-        fun = PK1_fun_2d if self.dim == 2 else PK1_fun_3d
-        P_µ, K_µ_s = fun(self.Young, self.Poisson, self.F)
-        shape = (self.dim, self.dim, self.dim, self.dim)
-        K_µ = K_µ_s.reshape(shape).transpose((0,1,3,2))
-
-        response_p = self.constitutive(self.F, self.dim)
-        P_p, K_p = response_p[0], response_p[2]
-
-        P_error = np.linalg.norm(P_µ- P_p)/np.linalg.norm(P_µ)
-        if not P_error < self.tol:
-            print("Error(P) = {}".format(P_error))
-            print("P_µ:\n{}".format(P_µ))
-            print("P_p:\n{}".format(P_p))
-        K_error = np.linalg.norm(K_µ- K_p)/np.linalg.norm(K_µ)
-        if not K_error < self.tol:
-            print("Error(K) = {}".format(K_error))
-            flat_shape = (self.dim**2, self.dim**2)
-            print("K_µ:\n{}".format(K_µ.reshape(flat_shape)))
-            print("K_p:\n{}".format(K_p.reshape(flat_shape)))
-            print("diff:\n{}".format(K_p.reshape(flat_shape)-
-                                     K_µ.reshape(flat_shape)))
-        self.assertLess(P_error,
-                        self.tol)
-
-        self.assertLess(K_error,
-                        self.tol)
-
-
-if __name__ == "__main__":
-    unittest.main()
-
diff --git a/tests/python_binding_tests.py b/tests/python_binding_tests.py
index 14e4294..4606781 100755
--- a/tests/python_binding_tests.py
+++ b/tests/python_binding_tests.py
@@ -1,210 +1,220 @@
 #!/usr/bin/env python3
 """
 file   python_binding_tests.py
 
 @author Till Junge <till.junge@epfl.ch>
 
 @date   09 Jan 2018
 
 @brief  Unit tests for python bindings
 
 @section LICENCE
 
 Copyright © 2018 Till Junge
 
 µSpectre is free software; you can redistribute it and/or
 modify it under the terms of the GNU Lesser General Public License as
 published by the Free Software Foundation, either version 3, or (at
 your option) any later version.
 
 µSpectre is distributed in the hope that it will be useful, but
 WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 General Public License for more details.
 
 You should have received a copy of the GNU Lesser General Public License
 along with µSpectre; see the file COPYING. If not, write to the
 Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 Boston, MA 02111-1307, USA.
 
 Additional permission under GNU GPL version 3 section 7
 
 If you modify this Program, or any covered work, by linking or combining it
 with proprietary FFT implementations or numerical libraries, containing parts
 covered by the terms of those libraries' licenses, the licensors of this
 Program grant you additional permission to convey the resulting work.
 """
 
 import unittest
 import numpy as np
 
 from python_test_imports import µ
 
 from python_fft_tests import FFT_Check
 from python_projection_tests import *
 from python_material_linear_elastic3_test import MaterialLinearElastic3_Check
 from python_material_linear_elastic4_test import MaterialLinearElastic4_Check
 from python_material_linear_elastic_generic1_test import MaterialLinearElasticGeneric1_Check
 from python_material_linear_elastic_generic2_test import MaterialLinearElasticGeneric2_Check
 from python_field_tests import FieldCollection_Check
 
 from python_exact_reference_elastic_test import LinearElastic_Check
 
 from python_field_tests import FieldCollection_Check
+from python_muSpectre_gradient_integration_test import \
+    MuSpectre_gradient_integration_Check
 
 from python_material_evaluator_test import MaterialEvaluator_Check
 
 class CellCheck(unittest.TestCase):
     def test_Construction(self):
         """
         Simple check for cell constructors
         """
         resolution = [5,7]
         lengths = [5.2, 8.3]
         formulation = µ.Formulation.small_strain
         try:
             sys = µ.Cell(resolution,
                          lengths,
                          formulation)
             mat = µ.material.MaterialLinearElastic1_2d.make(sys, "material",
                                                             210e9, .33)
 
         except Exception as err:
             print(err)
             raise err
 
 class MaterialLinearElastic1_2dCheck(unittest.TestCase):
     def setUp(self):
         self.resolution = [5,7]
         self.lengths = [5.2, 8.3]
         self.formulation = µ.Formulation.small_strain
         self.sys = µ.Cell(self.resolution,
                           self.lengths,
                           self.formulation)
         self.mat = µ.material.MaterialLinearElastic1_2d.make(
             self.sys, "material", 210e9, .33)
 
     def test_add_material(self):
         self.mat.add_pixel([2,1])
 
 
 class SolverCheck(unittest.TestCase):
     def setUp(self):
         self.resolution = [3, 3]#[5,7]
         self.lengths = [3., 3.]#[5.2, 8.3]
         self.formulation = µ.Formulation.finite_strain
         self.sys = µ.Cell(self.resolution,
                           self.lengths,
                           self.formulation)
         self.hard = µ.material.MaterialLinearElastic1_2d.make(
             self.sys, "hard", 210e9, .33)
         self.soft = µ.material.MaterialLinearElastic1_2d.make(
             self.sys, "soft",  70e9, .33)
 
     def test_solve(self):
         for i, pixel in enumerate(self.sys):
             if i < 3:
                 self.hard.add_pixel(pixel)
             else:
                 self.soft.add_pixel(pixel)
 
         self.sys.initialise()
         tol = 1e-6
         Del0 = np.array([[0, .1],
                          [0,  0]])
         maxiter = 100
         verbose = 0
 
         solver=µ.solvers.SolverCG(self.sys, tol, maxiter, verbose)
         r = µ.solvers.de_geus(self.sys, Del0,
                               solver,tol, verbose)
         #print(r)
 
 
 class EigenStrainCheck(unittest.TestCase):
     def setUp(self):
         self.resolution = [3, 3]#[5,7]
         self.lengths = [3., 3.]#[5.2, 8.3]
         self.formulation = µ.Formulation.small_strain
         self.cell1 = µ.Cell(self.resolution,
                             self.lengths,
                             self.formulation)
         self.cell2 = µ.Cell(self.resolution,
                             self.lengths,
                             self.formulation)
         self.mat1 = µ.material.MaterialLinearElastic1_2d.make(
             self.cell1, "simple", 210e9, .33)
         self.mat2 = µ.material.MaterialLinearElastic2_2d.make(
             self.cell2, "eigen", 210e9, .33)
         self.mat3 = µ.material.MaterialLinearElastic2_2d.make(
             self.cell2, "eigen2", 120e9, .33)
 
     def test_globalisation(self):
         for pixel in self.cell2:
             self.mat2.add_pixel(pixel, np.random.rand(2,2))
         loc_eigenstrain = self.mat2.collection.get_real_field("Eigenstrain").array
         glo_eigenstrain = self.cell2.get_globalised_internal_real_array("Eigenstrain")
         error = np.linalg.norm(loc_eigenstrain-glo_eigenstrain)
         self.assertEqual(error, 0)
 
     def test_globalisation_constant(self):
         for i, pixel in enumerate(self.cell2):
             if i%2 == 0:
                 self.mat2.add_pixel(pixel, np.ones((2,2)))
             else:
                 self.mat3.add_pixel(pixel, np.ones((2,2)))
         glo_eigenstrain = self.cell2.get_globalised_internal_real_array("Eigenstrain")
         error = np.linalg.norm(glo_eigenstrain-1)
         self.assertEqual(error, 0)
 
+    def test_globalisation(self):
+        for pixel in self.cell2:
+            self.mat2.add_pixel(pixel, np.random.rand(2,2))
+        loc_eigenstrain = self.mat2.collection.get_real_field("Eigenstrain").array
+        glo_eigenstrain = self.cell2.get_globalised_internal_real_array("Eigenstrain")
+        error = np.linalg.norm(loc_eigenstrain-glo_eigenstrain)
+        self.assertEqual(error, 0)
+
     def test_solve(self):
         verbose_test = False
         if verbose_test:
             print("start test_solve")
         grad = np.array([[1.1,  .2],
                          [ .3, 1.5]])
         gl_strain = -0.5*(grad.T.dot(grad) - np.eye(2))
         gl_strain = -0.5*(grad.T + grad - 2*np.eye(2))
         grad = -gl_strain
         if verbose_test:
             print("grad =\n{}\ngl_strain =\n{}".format(grad, gl_strain))
         for i, pixel in enumerate(self.cell1):
             self.mat1.add_pixel(pixel)
             self.mat2.add_pixel(pixel, gl_strain)
         self.cell1.initialise()
         self.cell2.initialise()
         tol = 1e-6
         Del0_1 = grad
         Del0_2 = np.zeros_like(grad)
         maxiter = 2
         verbose = 0
 
         def solve(cell, grad):
             solver=µ.solvers.SolverCG(cell, tol, maxiter, verbose)
             r = µ.solvers.newton_cg(cell, grad,
                                     solver, tol, tol, verbose)
             return r
         results = [solve(cell, del0) for (cell, del0)
                    in zip((self.cell1, self.cell2),
                           (Del0_1, Del0_2))]
         P1 = results[0].stress
         P2 = results[1].stress
         error = np.linalg.norm(P1-P2)/np.linalg.norm(.5*(P1+P2))
 
         if verbose_test:
             print("cell 1, no eigenstrain")
             print("P1:\n{}".format(P1[:,0]))
             print("F1:\n{}".format(results[0].grad[:,0]))
 
             print("cell 2, with eigenstrain")
             print("P2:\n{}".format(P2[:,0]))
             print("F2:\n{}".format(results[1].grad[:,0]))
             print("end test_solve")
         self.assertLess(error, tol)
 
 
 
 
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/python_mpi_material_linear_elastic4_test.py b/tests/python_mpi_material_linear_elastic4_test.py
index f9b88cb..e2677b9 100644
--- a/tests/python_mpi_material_linear_elastic4_test.py
+++ b/tests/python_mpi_material_linear_elastic4_test.py
@@ -1,103 +1,103 @@
 #!/usr/bin/env python3
 # -*- coding:utf-8 -*-
 """
 @file   python_mpi_material_linear_elastic4_test.py
 
 @author Richard Leute <richard.leute@imtek.uni-freiburg.de>
 
 @date   27 Mar 2018
 
 @brief  test MPI-parallel linear elastic material
 
 @section LICENSE
 
 Copyright © 2018 Till Junge
 
 µSpectre is free software; you can redistribute it and/or
 modify it under the terms of the GNU Lesser General Public License as
 published by the Free Software Foundation, either version 3, or (at
 your option) any later version.
 
 µSpectre is distributed in the hope that it will be useful, but
 WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 General Public License for more details.
 
 You should have received a copy of the GNU Lesser General Public License
 along with µSpectre; see the file COPYING. If not, write to the
 Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 Boston, MA 02111-1307, USA.
 
 Additional permission under GNU GPL version 3 section 7
 
 If you modify this Program, or any covered work, by linking or combining it
 with proprietary FFT implementations or numerical libraries, containing parts
 covered by the terms of those libraries' licenses, the licensors of this
 Program grant you additional permission to convey the resulting work.
 """
 
 try:
     from mpi4py import MPI
 except ImportError:
     MPI = None
 
 import unittest
 import numpy as np
 
 from python_test_imports import µ
 
 
 def build_test_classes(fft):
     class MaterialLinearElastic4_Check(unittest.TestCase):
         """
         Check the implementation of storing the first and second Lame constant in
         each cell. Assign the same Youngs modulus and Poisson ratio to each cell,
         from which the two Lame constants are internally computed. Then calculate
         the stress and compare the result with stress=2*mu*Del0 (Hooke law for small
         symmetric strains).
         """
         def setUp(self):
             self.resolution = [7,7]
             self.lengths = [2.3, 3.9]
             self.formulation = µ.Formulation.small_strain
             self.sys = µ.Cell(self.resolution,
                               self.lengths,
                               self.formulation,
                               fft=fft,
                               communicator=MPI.COMM_WORLD)
             self.mat = µ.material.MaterialLinearElastic4_2d.make(
                 self.sys, "material")
 
         def test_solver(self):
             Youngs_modulus = 10.
             Poisson_ratio  = 0.3
 
             for i, pixel in enumerate(self.sys):
                 self.mat.add_pixel(pixel, Youngs_modulus, Poisson_ratio)
-    
+
             self.sys.initialise()
             tol = 1e-6
             Del0 = np.array([[0, 0.025],
                              [0.025,  0]])
             maxiter = 100
             verbose = 1
-    
+
             solver=µ.solvers.SolverCG(self.sys, tol, maxiter, verbose)
             r = µ.solvers.newton_cg(self.sys, Del0,
                                     solver, tol, tol, verbose)
-    
+
             #compare the computed stress with the trivial by hand computed stress
             mu = (Youngs_modulus/(2*(1+Poisson_ratio)))
             stress = 2*mu*Del0
-    
+
             self.assertLess(np.linalg.norm(r.stress-stress.reshape(-1,1)), 1e-8)
 
     return MaterialLinearElastic4_Check
 
 linear_elastic4 = {}
 for fft, is_parallel in µ.fft.fft_engines:
     if is_parallel:
         linear_elastic4[fft] = build_test_classes(fft)
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/tests/python_muSpectre_gradient_integration_test.py b/tests/python_muSpectre_gradient_integration_test.py
new file mode 100644
index 0000000..7fcc399
--- /dev/null
+++ b/tests/python_muSpectre_gradient_integration_test.py
@@ -0,0 +1,661 @@
+#!/usr/bin/env python3
+# -*- coding:utf-8 -*-
+"""
+@file   python_muSpectre_gradient_integration_test.py
+
+@author Richard Leute <richard.leute@imtek.uni-freiburg.de>
+
+@date   23 Nov 2018
+
+@brief  test the functionality of gradient_integration.py
+
+@section LICENSE
+
+Copyright © 2018 Till Junge, Richard Leute
+
+µSpectre is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public License as
+published by the Free Software Foundation, either version 3, or (at
+your option) any later version.
+
+µSpectre is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with µSpectre; see the file COPYING. If not, write to the
+Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.
+
+Additional permission under GNU GPL version 3 section 7
+
+If you modify this Program, or any covered work, by linking or combining it
+with proprietary FFT implementations or numerical libraries, containing parts
+covered by the terms of those libraries' licenses, the licensors of this
+Program grant you additional permission to convey the resulting work.
+"""
+
+import unittest
+import numpy as np
+import scipy.misc as sm
+import itertools
+from python_test_imports import µ
+
+### Helper functions
+def init_X_F_Chi(lens, res, rank=2):
+    """
+    Setup all the needed parameters for initialization of the deformation
+    gradient F and the corresponding deformation map/field Chi_X.
+
+    Keyword Arguments:
+    lens -- list [Lx, Ly, ...] of box lengths in each direction (dtype=float)
+    res  -- list [Nx, Ny, ...] of grid resoultions (dtype = int)
+    rank -- int (default=2), rank of the deformation gradient tensor F.
+            (dtype = int)
+
+    Returns:
+    d    : np.array of grid spacing for each spatial direction (dtype = float)
+    dim  : int dimension of the structure, derived from len(res).
+    x_n  : np.ndarray shape=(res.shape+1, dim) initial nodal/corner positions
+           as created by gradient_integration.compute_grid (dtype = float)
+    x_c  : np.ndarray shape=(res.shape+1, dim) initial cell center positions
+           as created by gradient_integration.compute_grid (dtype = float)
+    F    : np.zeros shape=(res.shape, dim*rank) initialise deformation gradient
+           (dtype = float)
+    Chi_n: np.zeros shape=((res+1).shape, dim) initialise deformation field
+           (dtype = float)
+    Chi_c: np.zeros shape=(res.shape, dim) initialise deformation field
+           (dtype = float)
+    freqs: np.ndarray as returned by compute_wave_vectors(). (dtype = float)
+    """
+    lens = np.array(lens)
+    res  = np.array(res)
+    d    = lens / res
+    dim  = len(res)
+    x_n, x_c = µ.gradient_integration.compute_grid(lens, res)
+    F     = np.zeros(x_c.shape + (dim,)*(rank-1))
+    Chi_n = np.zeros(x_n.shape)
+    Chi_c = np.zeros(x_c.shape)
+    freqs = µ.gradient_integration.compute_wave_vectors(lens, res)
+
+    return d, dim, x_n, x_c, F, Chi_n, Chi_c, freqs
+
+def correct_to_zero_mean(Chi, d, nodal=False):
+    """
+    corrects the displacement field such that it's integral is zero. By this one
+    can get rid of a constant factor in the deformation gradient. This function
+    is specialized for this file and should not be used somewhere else.
+
+    Keywords:
+    Chi   : np.ndarray of the uncorrected analytic placements (dtype = float)
+    d     : np.array of the gridspacing in each spatial direction (dtype = float)
+    nodal : bool (default False) specifies if the input are nodal or cell/center
+            values. Default interpretation are cell/center values.
+
+    Returns:
+    Chi   : np.ndarray of zero mean corrected analytic placements (dtype = float)
+    """
+    Chi_zm = np.copy(Chi)
+    res = np.array(Chi_zm.shape[:-1])
+    dim = res.size
+    Chi_zm -= (Chi_zm.sum(axis=tuple(range(dim)))/np.prod(res))\
+              .reshape((1,)*dim + (dim,))
+    return Chi_zm
+
+def test_integrate(order, F, Chi_n, Chi_c, tol):
+    """
+    make the convergence tests for the integration
+
+    Keywords:
+    order : list of integration orders which are tested (dtype = int)
+    F     : np.ndarray applied deformation gradient (dtype = float)
+    Chi_n : np.ndarray expected node positions (dtype = float)
+    Chi_c : np.ndarray expected cell positions (dtype = float)
+    tol   : list of tolerances for each order. If it is a single value the same
+            tolerance is used for each order. (dtype = float)
+    """
+    print('Maybe implement a function like this...')
+
+def central_diff_derivative(data, d, order, rank=1):
+    """
+    Compute the first derivative of a function with values 'data' with the
+    central difference approximation to the order 'order'. The function values
+    are on a rectangualar grid with constant grid spacing 'd' in each direction.
+
+    CAUTION:
+    The function is assumed to be periodic (pbc)!
+    Thus, if there is a discontinuity at the boundaries you have to expect
+    errors in the derivative at the vicinity close to the discontinuity.
+
+    Keyword Arguments:
+    data  -- np.ndarray of shape=(resolution, dim*rank) function values on an
+             equally spaced grid, with grid spacing 'd' (dtype = float)
+    d     -- scalar or np.array of grid spacing in each direction. Scalar is
+             interpreted as equal spacing in each direction (dtype = float)
+    order -- int >= 1, gives the accuracy order of the central difference
+             approximation (dtype = int)
+    rank  -- int, rank of the data tensor
+
+    Returns:
+    deriv: np.ndarray of shape=(resolution, dim, dim) central difference
+           derivative of given order (dtype = float)
+    """
+    dim = len(data.shape)-rank
+    weights = sm.central_diff_weights(2*order + 1)
+    deriv = np.zeros(data.shape + (dim,))
+    for ax in range(dim):
+        for i in range(2*order + 1):
+            deriv[...,ax] += weights[i]*np.roll(data, order-i, axis=ax) / d[ax]
+
+    return deriv
+
+
+class MuSpectre_gradient_integration_Check(unittest.TestCase):
+    """
+    Check the implementation of all muSpectre.gradient_integration functions.
+    """
+
+    def setUp(self):
+        self.lengths    = np.array([2.4, 3.7, 4.1])
+        self.resolution = np.array([5, 3, 5])
+
+        self.norm_tol = 1e-8
+
+    def test_central_diff_derivative(self):
+        """
+        Test of the central difference approximation by central_diff_derivative
+        of the first derivative of a function on a grid.
+        """
+        res  = self.resolution * 15
+        lens = self.lengths
+        for j in range(1, len(res)):
+            d, dim, x_n, x_c, deriv, f_n, f_c, freqs = init_X_F_Chi(lens[:j],
+                                                                    res[:j])
+            f_c = np.sin(2*np.pi/lens[:j] * x_c)
+            for i in range(j):
+                deriv[...,i,i] =2*np.pi/lens[i]*np.cos(2*np.pi*
+                                                       x_c[...,i]/lens[i])
+            approx_deriv = central_diff_derivative(f_c, d[:j+1], order=5)
+            self.assertLess(np.linalg.norm(deriv-approx_deriv),
+                            self.norm_tol)
+
+    def test_compute_wave_vectors(self):
+        """
+        Test the construction of a wave vector grid by compute_wave_vectors
+        for an arbitrary dimension.
+        """
+        lens  = [4, 6, 7]
+        res   = [3, 4, 5]
+        Nx, Ny, Nz = res
+        q_1d = lambda i: 2*np.pi/lens[i] * \
+               np.append(np.arange(res[i]-res[i]//2),
+                         -np.arange(1, res[i]//2 + 1)[::-1])
+        qx = q_1d(0)
+        qy = q_1d(1)
+        qz = q_1d(2)
+        q  = np.zeros(tuple(res) + (3,))
+        for i,j,k in itertools.product(range(Nx), range(Ny), range(Nz)):
+            q[i,j,k,:] = np.array([qx[i], qy[j], qz[k]])
+        for n in range(1,4):
+            comp_q = µ.gradient_integration.compute_wave_vectors(lens[:n],
+                                                                 res[:n])
+            s = (np.s_[:],)*n + (0,)*(3-n) + (np.s_[:n],)
+            self.assertLess(np.linalg.norm(comp_q - q[s]), self.norm_tol)
+
+    def test_compute_grid(self):
+        """
+        Test the function compute_grid which creates an orthogonal
+        equally spaced grid of the given resolution and lengths.
+        """
+        lens = self.lengths
+        res  = self.resolution
+        d    = np.array(lens)/np.array(res)
+        grid_n = np.zeros(tuple(res+1) + (len(res),))
+        Nx, Ny, Nz = res+1
+        for i,j,k in itertools.product(range(Nx), range(Ny), range(Nz)):
+            grid_n[i,j,k,:] = np.array([i*d[0], j*d[1], k*d[2]])
+        grid_c = (grid_n - d/2)[1:,1:,1:,:]
+        for n in range(1,4):
+            x_n, x_c  = µ.gradient_integration.compute_grid(lens[:n], res[:n])
+            s = (np.s_[:],)*n + (0,)*(3-n) + (np.s_[:n],)
+            self.assertLess(np.linalg.norm(x_c - grid_c[s]), self.norm_tol)
+            self.assertLess(np.linalg.norm(x_n - grid_n[s]), self.norm_tol)
+
+    def test_reshape_gradient(self):
+        """
+        Test if reshape gradient transforms a flattend second order tensor in
+        the right way to a shape resolution + [dim, dim].
+        """
+        lens = list(self.lengths)
+        res  = list(self.resolution)
+        tol  = 1e-5
+        formulation = µ.Formulation.finite_strain
+        DelF = np.array([[0   , 0.01, 0.02],
+                         [0.03, 0   , 0.04],
+                         [0.05, 0.06, 0   ]])
+        one  = np.eye(3,3)
+        for n in range(2,4):
+            sys = µ.Cell(res[:n], lens[:n], formulation)
+            if n == 2:
+                mat = µ.material.MaterialLinearElastic1_2d.make(sys, "material",
+                                                                10, 0.3)
+            if n == 3:
+                mat = µ.material.MaterialLinearElastic1_3d.make(sys, "material",
+                                                                10, 0.3)
+            for pixel in sys:
+                mat.add_pixel(pixel)
+            solver = µ.solvers.SolverCG(sys, tol, maxiter=100, verbose=0)
+            r = µ.solvers.newton_cg(sys, DelF[:n, :n],
+                                    solver, tol, tol , verbose=0)
+            grad = µ.gradient_integration.reshape_gradient(r.grad,list(res[:n]))
+            grad_theo = (DelF[:n, :n] + one[:n, :n]).reshape((1,)*n+(n,n,))
+            self.assertEqual(grad.shape, tuple(res[:n])+(n,n,))
+            self.assertLess(np.linalg.norm(grad - grad_theo), self.norm_tol)
+
+    def test_complement_periodically(self):
+        """
+        Test the periodic reconstruction of an array. Lower left entries are
+        added into the upper right part of the array.
+        """
+        #1D grid scalars
+        x_test   = np.array([0,1,2,3])
+        x_test_p = np.array([0,1,2,3, 0])
+        x_p      = µ.gradient_integration.complement_periodically(x_test, 1)
+        self.assertLess(np.linalg.norm(x_p-x_test_p), self.norm_tol)
+
+        #2D grid scalars
+        x_test   = np.array([[1,2,3,4],
+                             [5,6,7,8]])
+        x_test_p = np.array([[1,2,3,4,1],
+                             [5,6,7,8,5],
+                             [1,2,3,4,1]])
+        x_p      = µ.gradient_integration.complement_periodically(x_test, 2)
+        self.assertLess(np.linalg.norm(x_p-x_test_p), self.norm_tol)
+
+        #2D grid vectors
+        x_test   = np.array([[[1,2,3]   , [3,4,5]]  ,
+                             [[6,7,8]   , [9,10,11]],
+                             [[12,13,14], [15,6,17]] ])
+        x_test_p = np.array([[[1,2,3]   , [3,4,5]  , [1,2,3]]   ,
+                             [[6,7,8]   , [9,10,11], [6,7,8]]   ,
+                             [[12,13,14], [15,6,17], [12,13,14]],
+                             [[1,2,3]   , [3,4,5]  , [1,2,3]]    ])
+        x_p      = µ.gradient_integration.complement_periodically(x_test, 2)
+        self.assertLess(np.linalg.norm(x_p-x_test_p), self.norm_tol)
+
+    def test_get_integrator(self):
+        """
+        Test if the right integrator is computed.
+        """
+        #even grid
+        lens_e = np.array([1,1,1])
+        res_e  = np.array([2,2,2])
+        x_n_e, x_c_e = µ.gradient_integration.compute_grid(lens_e, res_e)
+        freqs_e = µ.gradient_integration.compute_wave_vectors(lens_e, res_e)
+        #odd grid
+        lens_o = np.array([1,1])
+        res_o  = np.array([3,3])
+        x_n_o, x_c_o = µ.gradient_integration.compute_grid(lens_o, res_o)
+        delta_x = 1/3
+        freqs_o = µ.gradient_integration.compute_wave_vectors(lens_o, res_o)
+
+        ### order=0
+        int_ana = 1j/(2*np.pi)*np.array([[[[ 0  , 0  , 0], [ 0  , 0  ,-1  ]] ,
+                                          [[ 0  ,-1  , 0], [ 0  ,-1/2,-1/2]]],
+                                         [[[-1  , 0  , 0], [-1/2, 0  ,-1/2]] ,
+                                          [[-1/2,-1/2, 0], [-1/3,-1/3,-1/3]]] ])
+        dim,shape,integrator = µ.gradient_integration.\
+                               get_integrator(x_c_e, freqs_e, order=0)
+        self.assertEqual(dim, len(res_e))
+        self.assertEqual(shape, tuple(res_e))
+        self.assertLess(np.linalg.norm(integrator-int_ana), self.norm_tol)
+
+        ### order=1
+        #even grid
+        int_ana = np.zeros(res_e.shape)
+        dim,shape,integrator = µ.gradient_integration.\
+                               get_integrator(x_c_e, freqs_e, order=1)
+        self.assertEqual(dim, len(res_e))
+        self.assertEqual(shape, tuple(res_e))
+        self.assertLess(np.linalg.norm(integrator-int_ana), self.norm_tol)
+
+        #odd grid
+        s  = lambda q: np.sin(2*np.pi*q*delta_x)
+        sq = lambda q: (np.sin(2*np.pi*np.array(q)*delta_x)**2).sum()
+        int_ana = 1j * delta_x *\
+                  np.array([[[ 0                 , 0                 ],
+                             [ 0                 , s(1)/sq([0,1])    ],
+                             [ 0                 , s(-1)/sq([0,-1])  ]],
+                            [[ s(1)/sq([1,0])    , 0                 ],
+                             [ s(1)/sq([1,1])    , s(1)/sq([1,1])    ],
+                             [ s(1)/sq([1,-1])   , s(-1)/sq([1,-1])  ]],
+                            [[ s(-1)/sq([-1,0])  , 0                 ],
+                             [ s(-1)/sq([-1,1])  , s(1)/sq([-1,1])   ],
+                             [ s(-1)/sq([-1,-1]) , s(-1)/sq([-1,-1]) ]]])
+
+        dim,shape,integrator = µ.gradient_integration.\
+                               get_integrator(x_c_o, freqs_o, order=1)
+        self.assertEqual(dim, len(res_o))
+        self.assertEqual(shape, tuple(res_o))
+        self.assertLess(np.linalg.norm(integrator-int_ana), self.norm_tol)
+
+        ### order=2
+        #even grid
+        int_ana = np.zeros(res_e.shape)
+        dim,shape,integrator = µ.gradient_integration.\
+                               get_integrator(x_c_e, freqs_e, order=2)
+        self.assertEqual(dim, len(res_e))
+        self.assertEqual(shape, tuple(res_e))
+        self.assertLess(np.linalg.norm(integrator-int_ana), self.norm_tol)
+
+        #odd grid
+        lens_o = np.array([1,1])
+        res_o  = np.array([3,3])
+        x_n, x_c = µ.gradient_integration.compute_grid(lens_o, res_o)
+        delta_x = 1/3
+        freqs    = µ.gradient_integration.compute_wave_vectors(lens_o, res_o)
+        s  = lambda q: 8*np.sin(2*np.pi*q*delta_x) + np.sin(2*2*np.pi*q*delta_x)
+        sq = lambda q: ( (64*np.sin(2*np.pi*np.array(q)*delta_x)**2).sum() -
+                         (np.sin(2*2*np.pi*np.array(q)*delta_x)**2).sum() )
+        int_ana = 6 * 1j * delta_x *\
+                  np.array([[[ 0                 , 0                 ],
+                             [ 0                 , s(1)/sq([0,1])    ],
+                             [ 0                 , s(-1)/sq([0,-1])  ]],
+                            [[ s(1)/sq([1,0])    , 0                 ],
+                             [ s(1)/sq([1,1])    , s(1)/sq([1,1])    ],
+                             [ s(1)/sq([1,-1])   , s(-1)/sq([1,-1])  ]],
+                            [[ s(-1)/sq([-1,0])  , 0                 ],
+                             [ s(-1)/sq([-1,1])  , s(1)/sq([-1,1])   ],
+                             [ s(-1)/sq([-1,-1]) , s(-1)/sq([-1,-1]) ]]])
+
+        dim,shape,integrator = µ.gradient_integration.\
+                               get_integrator(x_c_o, freqs_o, order=2)
+        self.assertEqual(dim, len(res_o))
+        self.assertEqual(shape, tuple(res_o))
+        self.assertLess(np.linalg.norm(integrator-int_ana), self.norm_tol)
+
+    def test_integrate_tensor_2(self):
+        """
+        Test the correct integration of a second-rank tensor gradient field,
+        like the deformation gradient.
+        """
+        order = [1, 2] #list of higher order finite difference integration which
+                       #will be checked.
+
+        ### cosinus, diagonal deformation gradient
+        res  = [15, 15, 14]
+        lens = [7, 1.4, 3]
+        d, dim, x_n, x_c, F, Chi_n, Chi_c, freqs = init_X_F_Chi(lens, res)
+        for i in range(dim):
+            F[:,:,:,i,i] = 0.8*np.pi/lens[i]*np.cos(4*np.pi*
+                                                    x_c[:,:,:,i]/lens[i])
+        Chi_n = 0.2 * np.sin(4*np.pi*x_n/lens)
+        Chi_c = 0.2 * np.sin(4*np.pi*x_c/lens)
+        # zeroth order correction
+        placement_n = µ.gradient_integration.integrate_tensor_2(F, x_n, freqs,
+                        staggered_grid=True, order=0)
+        placement_c = µ.gradient_integration.integrate_tensor_2(F, x_c, freqs,
+                        staggered_grid=False, order=0)
+        self.assertLess(np.linalg.norm(Chi_n - placement_n), self.norm_tol)
+        self.assertLess(np.linalg.norm(Chi_c - placement_c), self.norm_tol)
+        # higher order correction
+        for n in order:
+            tol_n = [1.334, 0.2299] #adjusted tolerances for node points
+            F_c = central_diff_derivative(Chi_c, d, order=n)
+            placement_n = µ.gradient_integration.integrate_tensor_2(F_c, x_n,
+                            freqs, staggered_grid=True, order=n)
+            placement_c = µ.gradient_integration.integrate_tensor_2(F_c, x_c,
+                            freqs, staggered_grid=False,order=n)
+            self.assertLess(np.linalg.norm(Chi_n - placement_n), tol_n[n-1])
+            self.assertLess(np.linalg.norm(Chi_c - placement_c), self.norm_tol)
+
+        ### cosinus, non-diagonal deformation gradient
+        res  = [15, 12, 11]
+        lens = [8, 8, 8]
+        d, dim, x_n, x_c, F, Chi_n, Chi_c, freqs = init_X_F_Chi(lens, res)
+        F[:,:,:,0,0] = 4*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*x_c[:,:,:,0])
+        F[:,:,:,1,1] = 2*np.pi/lens[1]*np.cos(2*np.pi/lens[1]*x_c[:,:,:,1])
+        F[:,:,:,2,2] = 2*np.pi/lens[2]*np.cos(2*np.pi/lens[2]*x_c[:,:,:,2])
+        F[:,:,:,1,0] = 2*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*x_c[:,:,:,0])
+        F[:,:,:,2,0] = 2*np.pi/lens[0]*np.cos(2*np.pi/lens[0]*x_c[:,:,:,0])
+        for i in range(dim):
+            Chi_c[:,:,:,i]= np.sin(2*np.pi*x_c[:,:,:,i]/lens[i])  \
+                            + np.sin(2*np.pi*x_c[:,:,:,0]/lens[0])
+            Chi_n[:,:,:,i]= np.sin(2*np.pi*x_n[:,:,:,i]/lens[i])  \
+                            + np.sin(2*np.pi*x_n[:,:,:,0]/lens[0])
+        # zeroth order correction
+        placement_n = µ.gradient_integration.integrate_tensor_2(F, x_n, freqs,
+                        staggered_grid=True, order=0)
+        placement_c = µ.gradient_integration.integrate_tensor_2(F, x_c, freqs,
+                        staggered_grid=False, order=0)
+        self.assertLess(np.linalg.norm(Chi_n - placement_n), self.norm_tol)
+        self.assertLess(np.linalg.norm(Chi_c - placement_c), self.norm_tol)
+        # higher order correction
+        for n in order:
+            tol_n = [2.563, 0.1544] #adjusted tolerances for node points
+            F_c = central_diff_derivative(Chi_c, d, order=n)
+            placement_n = µ.gradient_integration.integrate_tensor_2(F_c, x_n,
+                            freqs, staggered_grid=True, order=n)
+            placement_c = µ.gradient_integration.integrate_tensor_2(F_c, x_c,
+                            freqs, staggered_grid=False,order=n)
+            self.assertLess(np.linalg.norm(Chi_n - placement_n), tol_n[n-1])
+            self.assertLess(np.linalg.norm(Chi_c - placement_c), self.norm_tol)
+
+        ### polynomial, diagonal deformation gradient
+        # Choose the prefactors of the polynomial such that at least Chi_X and F
+        # have respectively the same values at the boundaries (here X_i=0 and
+        # X_i=4).
+        res  = [13, 14, 11]
+        lens = [4, 4, 4]
+        d, dim, x_n, x_c, F, Chi_n, Chi_c, freqs = init_X_F_Chi(lens, res)
+        for i in range(dim):
+            F[:,:,:,i,i] = 32*x_c[:,:,:,i] -24*x_c[:,:,:,i]**2+4*x_c[:,:,:,i]**3
+        Chi_n = -128/15 + 16*x_n**2 -8*x_n**3 +x_n**4
+        Chi_c = -128/15 + 16*x_c**2 -8*x_c**3 +x_c**4
+        #subtract the mean of Chi_c, because the numeric integration is done to
+        #give a zero mean fluctuating deformation field.
+        mean_c = Chi_c.sum(axis=tuple(range(dim)))/ \
+                 np.array(Chi_c.shape[:-1]).prod()
+        Chi_n -= mean_c.reshape((1,)*dim + (dim,))
+        Chi_c -= mean_c.reshape((1,)*dim + (dim,))
+        # zeroth order correction
+        placement_n = µ.gradient_integration.integrate_tensor_2(F, x_n, freqs,
+                        staggered_grid=True, order=0)
+        placement_c = µ.gradient_integration.integrate_tensor_2(F, x_c, freqs,
+                        staggered_grid=False, order=0)
+        self.assertLess(np.linalg.norm(Chi_n - placement_n), 0.19477)
+        self.assertLess(np.linalg.norm(Chi_c - placement_c), 0.67355)
+        # higher order correction
+        for n in order:
+            tol_n = [18.266, 2.9073] #adjusted tolerances for node points
+            F_c = central_diff_derivative(Chi_c, d, order=n)
+            placement_n = µ.gradient_integration.integrate_tensor_2(F_c, x_n,
+                            freqs, staggered_grid=True, order=n)
+            placement_c = µ.gradient_integration.integrate_tensor_2(F_c, x_c,
+                            freqs, staggered_grid=False,order=n)
+            self.assertLess(np.linalg.norm(Chi_n - placement_n), tol_n[n-1])
+            self.assertLess(np.linalg.norm(Chi_c - placement_c), self.norm_tol)
+
+        ### Realistic test:
+        #   shear of a two dimensional material with two different Young moduli.
+        order_all = [0]+order #orders for which the test is run
+        #initialize material structure
+        res  = [ 9, 21] #resolution
+        lens = [ 9, 21] #lengths
+        d, dim, x_n, x_c, _, _, _, freqs = init_X_F_Chi(lens, res)
+        formulation = µ.Formulation.finite_strain
+        Young   = [10, 20]   #Youngs modulus for each phase (soft, hard)
+        Poisson = [0.3, 0.3] #Poissons ratio for each phase
+
+        #geometry (two slabs stacked in y-direction with,
+        #hight h (soft material) and hight res[1]-h (hard material))
+        h            = res[1]//2
+        phase        = np.zeros(tuple(res), dtype=int)
+        phase[:, h:] = 1
+        phase        = phase.flatten()
+        cell = µ.Cell(res, lens, formulation)
+        mat  = µ.material.MaterialLinearElastic4_2d.make(cell, "material")
+        for i, pixel in enumerate(cell):
+            mat.add_pixel(pixel, Young[phase[i]], Poisson[phase[i]])
+        cell.initialise()
+        DelF = np.array([[0 , 0.01],
+                         [0 , 0   ]])
+
+        # µSpectre solution
+        solver = µ.solvers.SolverCG(cell, tol=1e-6, maxiter=100, verbose=0)
+        result = µ.solvers.newton_cg(cell, DelF, solver, newton_tol=1e-6,
+                                     equil_tol=1e-6, verbose=0)
+        F = µ.gradient_integration.reshape_gradient(result.grad, res)
+        fin_pos = {} #µSpectre computed center and node positions for all orders
+        for n in order_all:
+            placement_n = µ.gradient_integration.integrate_tensor_2(F, x_n,
+                            freqs, staggered_grid=True, order=n)
+            placement_c = µ.gradient_integration.integrate_tensor_2(F, x_c,
+                            freqs, staggered_grid=False, order=n)
+            fin_pos[str(n)+'_n'] = placement_n
+            fin_pos[str(n)+'_c'] = placement_c
+
+        # analytic solution, "placement_ana" (node and center)
+        l_soft = d[1] * h           #height soft material
+        l_hard = d[1] * (res[1]-h)  #height hard material
+        Shear_modulus = np.array(Young) / (2 * (1+np.array(Poisson)))
+        mean_shear_strain = 2*DelF[0,1]
+        shear_strain_soft = (lens[1]*mean_shear_strain) / (l_soft
+                                + l_hard * Shear_modulus[0]/Shear_modulus[1])
+        shear_strain_hard = (lens[1]*mean_shear_strain) / (l_soft
+                                * Shear_modulus[1]/Shear_modulus[0] + l_hard)
+        placement_ana_n = np.zeros(x_n.shape)
+        placement_ana_c = np.zeros(x_c.shape)
+
+        #x-coordinate
+        #soft material
+        placement_ana_n[:,:h+1,0] = shear_strain_soft/2 * x_n[:, :h+1, 1]
+        placement_ana_c[:,:h  ,0] = shear_strain_soft/2 * x_c[:, :h  , 1]
+        #hard material
+        placement_ana_n[:,h+1:,0] =shear_strain_hard/2 * (x_n[:,h+1:,1]-l_soft)\
+                                   + shear_strain_soft/2 * l_soft
+        placement_ana_c[:,h:  ,0] =shear_strain_hard/2 * (x_c[:,h:  ,1]-l_soft)\
+                                    + shear_strain_soft/2 * l_soft
+        #y-coordinate
+        placement_ana_n[:, :, 1] = 0
+        placement_ana_c[:, :, 1] = 0
+
+        #shift the analytic solution such that the average nonaffine deformation
+        #is zero (integral of the nonaffine deformation gradient + N*const != 0)
+        F_homo    = (1./(np.prod(res)) * F.sum(axis=tuple(np.arange(dim))))\
+                    .reshape((1,)*dim+(dim,)*2)
+        #integration constant = integral of the nonaffine deformation gradient/N
+        int_const = - ((placement_ana_c[:,:,0] - F_homo[:,:,0,1] * x_c[:,:,1])
+                       .sum(axis=1))[0] / res[1]
+        ana_sol_n = placement_ana_n + x_n + \
+                    np.array([int_const, 0]).reshape((1,)*dim+(dim,))
+        ana_sol_c = placement_ana_c + x_c + \
+                    np.array([int_const, 0]).reshape((1,)*dim+(dim,))
+
+        # check the numeric vs the analytic solution
+        tol_n = [2.2112e-3, 1.3488e-3, 1.8124e-3]
+        tol_c = [3.1095e-3, 3.2132e-2, 1.8989e-2]
+        for n in order_all:
+            norm_n = np.linalg.norm(fin_pos[str(n)+'_n'] - ana_sol_n)
+            norm_c = np.linalg.norm(fin_pos[str(n)+'_c'] - ana_sol_c)
+            self.assertLess(norm_n, tol_n[n])
+            self.assertLess(norm_c, tol_c[n])
+
+
+    def test_integrate_vector(self):
+        """Test the integration of a first-rank tensor gradient field."""
+        order = [1,2]
+
+        ### cosinus deformation gradient vector field
+        res  = [13, 14, 13]
+        lens = [ 7,  4, 5]
+        d, dim, x_n, x_c, df, Chi_n, Chi_c, freqs = init_X_F_Chi(lens, res, 1)
+        for i in range(dim):
+            df[:,:,:,i] = 0.8*np.pi/lens[i]*np.cos(4*np.pi*x_c[:,:,:,i]/lens[i])
+        Chi_n = 0.2 * np.sin(4*np.pi*x_n/lens).sum(axis=-1)
+        Chi_c = 0.2 * np.sin(4*np.pi*x_c/lens).sum(axis=-1)
+        # zeroth order correction
+        placement_n = µ.gradient_integration.integrate_vector(
+            df, x_n, freqs, staggered_grid=True, order=0)
+        placement_c = µ.gradient_integration.integrate_vector(
+            df, x_c, freqs, staggered_grid=False, order=0)
+        self.assertLess(np.linalg.norm(Chi_n - placement_n), self.norm_tol)
+        self.assertLess(np.linalg.norm(Chi_c - placement_c), self.norm_tol)
+        # higher order correction
+        for n in order:
+            tol_n = [1.404, 0.2882] #adjusted tolerances for node points
+            df_c = central_diff_derivative(Chi_c, d, order=n, rank=0)
+            placement_n = µ.gradient_integration.integrate_vector(
+                df_c, x_n, freqs, staggered_grid=True, order=n)
+            placement_c = µ.gradient_integration.integrate_vector(
+                df_c, x_c, freqs, staggered_grid=False, order=n)
+            self.assertLess(np.linalg.norm(Chi_n - placement_n), tol_n[n-1])
+            self.assertLess(np.linalg.norm(Chi_c - placement_c), self.norm_tol)
+
+        ### polynomial deformation gradient vector field
+        # Choose the prefactors of the polynomial such that at least Chi_X and F
+        # have respectively the same values at the boundaries (here X_i=0 and
+        # X_i=4).
+        res  = [12, 11, 13]
+        lens = [4, 4, 4]
+        d, dim, x_n, x_c, df, Chi_n, Chi_c, freqs = init_X_F_Chi(lens, res, 1)
+        for i in range(dim):
+            df[:,:,:,i] = 32*x_c[:,:,:,i]-24*x_c[:,:,:,i]**2+4*x_c[:,:,:,i]**3
+        Chi_n = (-128/15 + 16*x_n**2 -8*x_n**3 +x_n**4).sum(axis=-1)
+        Chi_c = (-128/15 + 16*x_c**2 -8*x_c**3 +x_c**4).sum(axis=-1)
+        #subtract the mean of Chi_c, because the numeric integration is done to
+        #give a zero mean fluctuating deformation field.
+        mean_c = Chi_c.sum() / np.array(Chi_c.shape).prod()
+        Chi_n -= mean_c
+        Chi_c -= mean_c
+        # zeroth order correction
+        placement_n = µ.gradient_integration.integrate_vector(
+            df, x_n, freqs, staggered_grid=True, order=0)
+        placement_c = µ.gradient_integration.integrate_vector(
+            df, x_c, freqs, staggered_grid=False, order=0)
+        self.assertLess(np.linalg.norm(Chi_n - placement_n), 0.20539)
+        self.assertLess(np.linalg.norm(Chi_c - placement_c), 0.67380)
+        # higher order correction
+        for n in order:
+            tol_n = [18.815, 3.14153] #adjusted tolerances for node points
+            df_c = central_diff_derivative(Chi_c, d, order=n, rank=0)
+            placement_n = µ.gradient_integration.integrate_vector(
+                df_c, x_n, freqs, staggered_grid=True, order=n)
+            placement_c = µ.gradient_integration.integrate_vector(
+                df_c, x_c, freqs, staggered_grid=False, order=n)
+            self.assertLess(np.linalg.norm(Chi_n - placement_n), tol_n[n-1])
+            self.assertLess(np.linalg.norm(Chi_c - placement_c), self.norm_tol)
+
+
+    def test_compute_placement(self):
+        """Test the computation of placements and the original positions."""
+        ### shear of a homogeneous material ###
+        res   = [ 3, 11] #resolution
+        lens  = [10, 10] #lengths
+        dim   = len(res) #dimension
+        x_n=µ.gradient_integration.compute_grid(np.array(lens),np.array(res))[0]
+
+        ### finite strain
+        formulation = µ.Formulation.finite_strain
+        cell = µ.Cell(res, lens, formulation)
+        mat  = µ.material.MaterialLinearElastic1_2d.make(cell, "material",
+                                                         Young=10, Poisson=0.3)
+        for pixel in cell:
+            mat.add_pixel(pixel)
+        cell.initialise()
+        DelF = np.array([[0 , 0.05],
+                         [0 , 0   ]])
+        # analytic
+        placement_ana = np.copy(x_n)
+        placement_ana[:,:,0] += DelF[0,1]*x_n[:,:,1]
+        # µSpectre solution
+        solver = µ.solvers.SolverCG(cell, tol=1e-6, maxiter=100, verbose=0)
+        result = µ.solvers.newton_cg(cell, DelF, solver, newton_tol=1e-6,
+                                     equil_tol=1e-6, verbose=0)
+        for r in [result, result.grad]:
+            #check input of result=OptimiseResult and result=np.ndarray
+            placement, x = µ.gradient_integration.compute_placement(
+                r, lens, res, order=0, formulation='finite_strain')
+            self.assertLess(np.linalg.norm(placement_ana - placement), 1e-12)
+            self.assertTrue((x_n == x).all())
diff --git a/tests/python_muSpectre_vtk_export_test.py b/tests/python_muSpectre_vtk_export_test.py
new file mode 100644
index 0000000..b4e8206
--- /dev/null
+++ b/tests/python_muSpectre_vtk_export_test.py
@@ -0,0 +1,82 @@
+#!/usr/bin/env python3
+# -*- coding:utf-8 -*-
+"""
+@file   python_muSpectre_vtk_export_test.py
+
+@author Richard Leute <richard.leute@imtek.uni-freiburg.de>
+
+@date   10 Jan 2019
+
+@brief  test the functionality of vtk_export.py
+
+@section LICENSE
+
+Copyright © 2019 Till Junge, Richard Leute
+
+µSpectre is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public License as
+published by the Free Software Foundation, either version 3, or (at
+your option) any later version.
+
+µSpectre is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with µSpectre; see the file COPYING. If not, write to the
+Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.
+
+Additional permission under GNU GPL version 3 section 7
+
+If you modify this Program, or any covered work, by linking or combining it
+with proprietary FFT implementations or numerical libraries, containing parts
+covered by the terms of those libraries' licenses, the licensors of this
+Program grant you additional permission to convey the resulting work.
+"""
+
+import unittest
+import numpy as np
+import tempfile
+import os
+from python_test_imports import µ
+
+class MuSpectre_vtk_export_Check(unittest.TestCase):
+    def setUp(self):
+        self.lengths    = np.array([1.1, 2.2, 3])
+        self.resolution = np.array([3, 5, 7])
+    def test_vtk_export(self):
+        """
+        Check the possibility to write scalar-, vector- and second rank tensor-
+        fields on the cell and node grid. The throw of exceptions is not checked
+        """
+        grad = np.array([[0, 0.01, 0],
+                         [0,    0, 0],
+                         [0,    0, 0]])
+        for dim in [2,3]:
+            #test if files are written for 2D and 3D
+            lens = self.lengths[:dim]
+            res  = self.resolution[:dim]
+            F = grad[:dim, :dim].reshape((1,)*dim + (dim,)*2)
+            x_n, x_c = µ.gradient_integration.compute_grid(lens, res)
+            freqs = µ.gradient_integration.compute_wave_vectors(lens, res)
+            placement_n = µ.gradient_integration.integrate_tensor_2(F, x_n,
+                            freqs, staggered_grid=True, order=0)
+            p_d = {'scalar'  : np.random.random(x_n.shape[:-1]),
+                   'vector'  : np.random.random(x_n.shape[:-1] + (dim,)),
+                   '2-tensor': np.random.random(x_n.shape[:-1] + (dim,)*2)}
+            c_d = {'scalar'  : np.random.random(x_c.shape[:-1]),
+                   'vector'  : np.random.random(x_c.shape[:-1] + (dim,)),
+                   '2-tensor': np.random.random(x_c.shape[:-1] + (dim,)*2)}
+            #The temporary directory is atomatically cleand up after one is
+            #exiting the block.
+            with tempfile.TemporaryDirectory(dir=os.getcwd()) as dir_name:
+                os.chdir(dir_name)
+                file_name = 'vtk_export_'+str(dim)+'D_test'
+                µ.vtk_export.vtk_export(file_name, x_n, placement_n,
+                                        point_data = p_d, cell_data = c_d)
+                assert os.path.exists(file_name+'.vtr') == 1,\
+                "vtk_export() was not able to write the {}D output file "\
+                "'{}.vtr'.".format(dim, file_name)
+                os.chdir('../')