diff --git a/CMakeLists.txt b/CMakeLists.txt
index a6686cc..ab59814 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,102 +1,102 @@
 #
 # (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 #
 
 cmake_minimum_required(VERSION 3.0)
 
 # Basic settings
 # ==============
 
 project(GooseFEM)
 
 option(BUILD_TESTS "Build tests" OFF)
 option(BUILD_EXAMPLES "Build examples" OFF)
 
 # Version
 # =======
 
 IF(DEFINED ENV{PKG_VERSION})
     set(GOOSEFEM_VERSION $ENV{PKG_VERSION})
 else()
     execute_process(
         COMMAND python -c "from setuptools_scm import get_version; print(get_version())"
         WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}
         OUTPUT_VARIABLE GOOSEFEM_VERSION
         OUTPUT_STRIP_TRAILING_WHITESPACE)
 endif()
 
-message(STATUS "Building GooseFEM ${GOOSEFEM_VERSION}")
+message(STATUS "Building ${PROJECT_NAME} ${GOOSEFEM_VERSION}")
 
 # Set target
 # ==========
 
 find_package(xtensor REQUIRED)
 
 add_library(GooseFEM INTERFACE)
 
 target_include_directories(GooseFEM INTERFACE
     $<INSTALL_INTERFACE:include>
     $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>)
 
 target_link_libraries(GooseFEM INTERFACE xtensor)
 
 # Installation
 # ============
 
 include(CMakePackageConfigHelpers)
 include(GNUInstallDirs)
 
 install(DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/include/" DESTINATION include)
 
 configure_file("include/${PROJECT_NAME}/version.h"
                "${CMAKE_CURRENT_BINARY_DIR}/version.h"
                 @ONLY)
 
 install(FILES "${CMAKE_CURRENT_BINARY_DIR}/version.h"
         DESTINATION "include/${PROJECT_NAME}/")
 
 install(TARGETS GooseFEM EXPORT GooseFEM-targets)
 
 install(
     EXPORT GooseFEM-targets
     FILE GooseFEMTargets.cmake
     DESTINATION "${CMAKE_INSTALL_LIBDIR}/cmake/GooseFEM")
 
 set(_GOOSEFEM ${CMAKE_SIZEOF_VOID_P})
 unset(CMAKE_SIZEOF_VOID_P)
 
 write_basic_package_version_file("${CMAKE_CURRENT_BINARY_DIR}/GooseFEMConfigVersion.cmake"
     VERSION ${GOOSEFEM_VERSION}
     COMPATIBILITY AnyNewerVersion)
 
 set(CMAKE_SIZEOF_VOID_P ${_GOOSEFEM})
 
 install(FILES "${CMAKE_CURRENT_SOURCE_DIR}/GooseFEMConfig.cmake"
               "${CMAKE_CURRENT_BINARY_DIR}/GooseFEMConfigVersion.cmake"
         DESTINATION "${CMAKE_INSTALL_LIBDIR}/cmake/GooseFEM")
 
 configure_file("${CMAKE_CURRENT_SOURCE_DIR}/GooseFEM.pc.in"
                "${CMAKE_CURRENT_BINARY_DIR}/GooseFEM.pc" @ONLY)
 
 install(FILES "${CMAKE_CURRENT_BINARY_DIR}/GooseFEM.pc"
         DESTINATION "${CMAKE_INSTALL_LIBDIR}/pkgconfig/")
 
 # Add builds
 # ==========
 
 include("GooseFEMConfig.cmake")
 
 include(CTest)
 
 if(BUILD_TESTS)
     enable_testing()
     add_subdirectory(test/basic)
 
     enable_testing()
     add_subdirectory(test/gmat)
 endif()
 
 if(BUILD_EXAMPLES)
     enable_testing()
     add_subdirectory(docs/examples)
 endif()
diff --git a/docs/changelog.rst b/docs/changelog.rst
index 30745d3..3460662 100644
--- a/docs/changelog.rst
+++ b/docs/changelog.rst
@@ -1,186 +1,192 @@
 
 *********
 Changelog
 *********
 
 v0.9.0
 ======
 
 API Changes
 -----------
 
 *   VectorPartitioned::asDofs(dofval_u, dofval_p, dofval) ->
     VectorPartitioned::dofsFromParitioned(dofval_u, dofval_p, dofval)
 
 *   VectorPartitioned::asNode(dofval_u, dofval_p, nodevec) ->
     VectorPartitioned::nodeFromPartitioned(dofval_u, dofval_p, nodevec)
 
 *   VectorPartitioned::asElement(dofval_u, dofval_p, elemvec) ->
     VectorPartitioned::elementFromPartitioned(dofval_u, dofval_p, elemvec)
 
 *   Version defines as replaced by ``#define GOOSEFEM_VERSION``,
     added convenience methods ``GooseFEM::version()`` and ``GooseFEM::version_dependencies()``.
 
 Deprecating in next version
 ----------------------------
 
 *   VectorPartitioned::assembleDofs_u
 *   VectorPartitioned::assembleDofs_p
 *   Mesh::Renumber::get
 *   Mesh::Reordered::get
 
+New functions
+-------------
+
+*   VectorPartitioned::dofs_is_u
+*   VectorPartitioned::dofs_is_p
+
 Changes under the hood
 ----------------------
 
 *   Overloading from Vector (also in Python API)
 *   Overloading from QuadratureBase (also in Python API)
 *   Added doxygen docs (published to GitHub pages)
 
 v0.8.6
 ======
 
 *   String-define safety: stringification + unquoting.
 
 v0.8.2
 ======
 
 *   Using setuptools_scm to manage version (#169)
 
 v0.8.1
 ======
 
 *   Various documentation updates: using doxygen (e.g. #168, #167, #157, #150)
     *  Adding autodocs using doxygen/breathe.
     *  Adding autodocs Python API with references to the C++ docs.
 *   Using GitHub pages for doxygen docs (#156, #155)
 *   Adding version information (incl. git commit hash) (#166)
 *   Adding GooseFEM::Element::Quad4::Quadrature::interp_N_vector
 *   Generalizing GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular
 *   Generalising implementation:
     *   Internally deriving from Vector
     *   Python API: unifying Element
     *   Python API: fixing overloaded methods
     *   Removing internal use of deprecated method
     *   Using "initQuadratureBase" in derived Quadrature classes
     *   Introducing QuadratureBase class -> avoids copies of convenience functions
 *   [CI] Using ctest command to improve output in case of test failure
 *   Restructuring environment (#154)
 *   Fixing readthedocs setup (#153)
 
 v0.8.0
 ======
 
 *   [CI] Using gcc-8
 *   Adding Mesh::Quad4::FineLayer::elementsLayer
 *   Stitch: Adding nodesets to example
 *   Stitch: Adding hybrid example. Adding assertions.
 *   Making API more functional
 *   Adding Mesh::ManualStich
 *   Adding Mesh::Stitch
 *   Minor style update
 *   [CMake] Minor updates in testing
 *   [CI] improve comments (#142)
 *   Combining tests MeshQuad4 (#141)
 *   Using clang on Windows (#139)
 
 v0.7.0
 ======
 
 *   Adding ``Mesh::Quad4::FineLayer::elementgrid_leftright``
 
 v0.6.1
 ======
 
 *   Minor bugfix ``Mesh::Quad4::FineLayer::elementgrid_around_ravel``: allowing huge sizes.
 
 v0.6.0
 ======
 
 *   Adding ``Mesh::Quad4::FineLayer::elementgrid_around_ravel``
 *   ``FineLayer::elementgrid_ravel``: Adding test
 *   Renaming ``elementMatrix`` -> ``elementgrid`` everywhere
 *   Adding ``Mesh::Quad4::FineLayer::elementgrid_ravel``
 *   Adding ``GOOFEM_WIP_ASSERT`` to assert if code needs to be generalized
 *   API change: renaming ``Mesh::Quad4::Regular::elementMatrix``
     -> M``esh::Quad4::Regular::elementgrid``.
 
 v0.5.1
 ======
 
 *   FineLayer - replica: bug-fix in size detection.
 *   Updated examples to new GMat API.
 
 v0.5.0
 ======
 
 *   Renaming ``MatrixDiagonal::AsDiagonal`` -> ``MatrixDiagonal::Todiagonal``
     to maintain API consistency.
 *   Adding ``Mesh::elemmap2nodemap``. Updating Python API.
 *   Adding ``roll`` to FineLayer.
 *   Adding ``Mesh::centers`` and ``Mesh::defaultElementType``.
 *   Mapping connectivity on generating FineLayer-object.
 *   Switching to new GMat API.
 *   Solver: force factorization on the first call.
 *   Sorting output of ``GooseFEM::Mesh::elem2node``. Adding checks.
 *   Switched to GitHub CI.
 *   Adding ``todense`` to sparse matrix classes.
 *   Adding ``dot`` to ``MatrixPartitioned``.
 
 v0.4.2
 ======
 
 *   CMake: using Eigen's CMake target.
 
 v0.4.1
 ======
 
 API additions
 -------------
 
 *   Added  "AllocateElemmat".
 
 v0.4.0
 ======
 
 API additions
 -------------
 
 *   Added "AllocateQtensor", "AllocateQscalar", "AllocateDofval", "AllocateNodevec", "AllocateElemvec".
 
 API changes
 -----------
 
 *   Removing Paraview interface: replaced by external libraries "XDMFWrite_HighFive" and "XDMFWrite_h5py".
 
 *   Element*: "dV" now only returns raw data, the "asTensor" member function (and free function) can be used to convert the 'qscalar' to a 'qtensor'.
 
 *   Separating sparse solver in separate class to offer more flexibility in the future.
 
 *   Adding "dot" to "Matrix".
 
 Other updates
 -------------
 
 *   Applying clang-format to source, python API, tests, and examples..
 
 *   Adding test GMatElastoPlasticQPot.
 
 *   Adding test based on hybrid material definitions.
 
 *   Formatting update: renaming all return variables "out" to "ret".
 
 *   Correction zero allocation to allows for dofval.size() > nodevec.size()
 
 *   Formatting update xt::amax and xt::sum.
 
 *   Renaming private function to begin with caps when the function allocates its return data.
 
 *   Reducing copies when using Eigen.
 
 *   Reducing default size examples.
 
 *   Supporting Windows (#87).
 
 *   Removing xtensor_fixed.
 
 *   Using xt::has_shape.
diff --git a/docs/examples/CMakeLists.txt b/docs/examples/CMakeLists.txt
index fa3fe05..142254a 100644
--- a/docs/examples/CMakeLists.txt
+++ b/docs/examples/CMakeLists.txt
@@ -1,107 +1,107 @@
 
 cmake_minimum_required(VERSION 3.16)
 
 if(CMAKE_CURRENT_SOURCE_DIR STREQUAL CMAKE_SOURCE_DIR)
     project(GooseFEM-examples)
     find_package(GooseFEM REQUIRED CONFIG)
 endif()
 
 set(ASSERT OFF)
 set(DEBUG OFF)
 option(SIMD "Enable xsimd" ON)
 option(WARNINGS "Enable warnings" ON)
 
 set(HIGHFIVE_USE_BOOST 0)
 set(HIGHFIVE_USE_XTENSOR 1)
 find_package(HighFive REQUIRED)
 find_package(XDMFWrite_HighFive REQUIRED)
 find_package(GMatElastic REQUIRED)
 find_package(GMatElastoPlastic REQUIRED)
-# find_package(GMatElastoPlasticFiniteStrainSimo REQUIRED)
+find_package(GMatElastoPlasticFiniteStrainSimo REQUIRED)
 find_package(GMatNonLinearElastic REQUIRED)
 
 add_library(libraries INTERFACE IMPORTED)
 
 target_link_libraries(libraries INTERFACE
     GooseFEM
     HighFive
     XDMFWrite_HighFive
     GMatElastic
     GMatElastoPlastic
-    # GMatElastoPlasticFiniteStrainSimo
+    GMatElastoPlasticFiniteStrainSimo
     GMatNonLinearElastic)
 
 if(SIMD)
     target_link_libraries(libraries INTERFACE xtensor::optimize xtensor::use_xsimd)
 endif()
 
 if(WARNINGS)
     target_link_libraries(libraries INTERFACE GooseFEM::compiler_warnings)
 endif()
 
 if(ASSERT)
     target_link_libraries(libraries INTERFACE GooseFEM::assert)
 endif()
 
 if(DEBUG)
     target_link_libraries(libraries INTERFACE GooseFEM::debug)
 endif()
 
 
 # create executable
 
 set(exec "statics_FixedDisplacements_LinearElastic_example")
 set(source "statics/FixedDisplacements_LinearElastic/example.cpp")
 add_executable(${exec} ${source})
 target_link_libraries(${exec} PRIVATE libraries)
 add_test(NAME ${exec} COMMAND ${exec})
 
 # create executable
 
 set(exec "statics_FixedDisplacements_LinearElastic_manual_partition")
 set(source "statics/FixedDisplacements_LinearElastic/manual_partition.cpp")
 add_executable(${exec} ${source})
 target_link_libraries(${exec} PRIVATE libraries)
 add_test(NAME ${exec} COMMAND ${exec})
 
 # create executable
 
 set(exec "statics_MixedPeriodic_LinearElastic_example")
 set(source "statics/MixedPeriodic_LinearElastic/example.cpp")
 add_executable(${exec} ${source})
 target_link_libraries(${exec} PRIVATE libraries)
 add_test(NAME ${exec} COMMAND ${exec})
 
 # create executable
 
-# set(exec "statics_Periodic_ElastoPlastic_main")
-# set(source "statics/Periodic_ElastoPlastic/main.cpp")
-# add_executable(${exec} ${source})
-# target_link_libraries(${exec} PRIVATE libraries)
-# add_test(NAME ${exec} COMMAND ${exec})
+set(exec "statics_Periodic_ElastoPlastic_main")
+set(source "statics/Periodic_ElastoPlastic/main.cpp")
+add_executable(${exec} ${source})
+target_link_libraries(${exec} PRIVATE libraries)
+add_test(NAME ${exec} COMMAND ${exec})
 
 # create executable
 
-# set(exec "statics_Periodic_ElastoPlasticFiniteStrainSimo_main")
-# set(source "statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp")
-# add_executable(${exec} ${source})
-# target_link_libraries(${exec} PRIVATE libraries)
-# add_test(NAME ${exec} COMMAND ${exec})
+set(exec "statics_Periodic_ElastoPlasticFiniteStrainSimo_main")
+set(source "statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp")
+add_executable(${exec} ${source})
+target_link_libraries(${exec} PRIVATE libraries)
+add_test(NAME ${exec} COMMAND ${exec})
 
 # create executable
 
-# set(exec "statics_Periodic_LinearElastic_main")
-# set(source "statics/Periodic_LinearElastic/main.cpp")
-# add_executable(${exec} ${source})
-# target_link_libraries(${exec} PRIVATE libraries)
-# add_test(NAME ${exec} COMMAND ${exec})
+set(exec "statics_Periodic_LinearElastic_main")
+set(source "statics/Periodic_LinearElastic/main.cpp")
+add_executable(${exec} ${source})
+target_link_libraries(${exec} PRIVATE libraries)
+add_test(NAME ${exec} COMMAND ${exec})
 
 # create executable
 
-# set(exec "statics_Periodic_NonLinearElastic_main")
-# set(source "statics/Periodic_NonLinearElastic/main.cpp")
-# add_executable(${exec} ${source})
-# target_link_libraries(${exec} PRIVATE libraries)
-# add_test(NAME ${exec} COMMAND ${exec})
+set(exec "statics_Periodic_NonLinearElastic_main")
+set(source "statics/Periodic_NonLinearElastic/main.cpp")
+add_executable(${exec} ${source})
+target_link_libraries(${exec} PRIVATE libraries)
+add_test(NAME ${exec} COMMAND ${exec})
 
 
diff --git a/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp b/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp
index 3b98d98..73415f9 100644
--- a/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp
+++ b/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp
@@ -1,237 +1,240 @@
 #include <Eigen/Eigen>
 #include <GMatElastoPlasticFiniteStrainSimo/Cartesian3d.h>
 #include <GooseFEM/GooseFEM.h>
 #include <XDMFWrite_HighFive.hpp>
 #include <highfive/H5Easy.hpp>
 
 namespace GM = GMatElastoPlasticFiniteStrainSimo::Cartesian3d;
 namespace GF = GooseFEM;
 namespace PV = XDMFWrite_HighFive;
 namespace H5 = H5Easy;
 
 int main()
 {
     // mesh
     // ----
 
     // define mesh
     GF::Mesh::Quad4::Regular mesh(5, 5);
 
     // mesh dimensions
     size_t nelem = mesh.nelem();
     size_t nne = mesh.nne();
     size_t ndim = mesh.ndim();
 
     // mesh definitions
     xt::xtensor<double, 2> coor = mesh.coor();
     xt::xtensor<size_t, 2> conn = mesh.conn();
     xt::xtensor<size_t, 2> dofs = mesh.dofs();
     xt::xtensor<size_t, 2> elmat = mesh.elementgrid();
 
     // periodicity and fixed displacements DOFs
     // ----------------------------------------
 
     // add control nodes
     GF::Tyings::Control control(coor, dofs);
     coor = control.coor();
     dofs = control.dofs();
     xt::xtensor<size_t, 2> control_dofs = control.controlDofs();
     xt::xtensor<size_t, 1> control_nodes = control.controlNodes();
 
     // extract fixed DOFs:
     // - all control nodes: to prescribe the deformation gradient
     // - one node of the mesh: to remove rigid body modes
     xt::xtensor<size_t, 1> iip = xt::concatenate(xt::xtuple(
         xt::reshape_view(control_dofs, {ndim * ndim}),
         xt::reshape_view(xt::view(dofs, xt::keep(mesh.nodesOrigin()), xt::all()), {ndim})));
 
     // get DOF-tyings, reorganise system
     GF::Tyings::Periodic tyings(coor, dofs, control_dofs, mesh.nodesPeriodic(), iip);
     dofs = tyings.dofs();
 
     // simulation variables
     // --------------------
 
     // vector definition:
     // provides methods to switch between dofval/nodeval/elemvec, or to manipulate a part of them
     GF::VectorPartitionedTyings vector(conn, dofs, tyings.Cdu(), tyings.Cdp(), tyings.Cdi());
 
     // nodal quantities
     xt::xtensor<double, 2> disp = xt::zeros<double>(coor.shape()); // nodal displacement
     xt::xtensor<double, 2> du = xt::zeros<double>(coor.shape());   // iterative displacement update
     xt::xtensor<double, 2> fint = xt::zeros<double>(coor.shape()); // internal force
     xt::xtensor<double, 2> fext = xt::zeros<double>(coor.shape()); // external force
     xt::xtensor<double, 2> fres = xt::zeros<double>(coor.shape()); // residual force
 
     // element vectors / matrix
     xt::xtensor<double, 3> ue = xt::empty<double>({nelem, nne, ndim});
     xt::xtensor<double, 3> fe = xt::empty<double>({nelem, nne, ndim});
     xt::xtensor<double, 3> Ke = xt::empty<double>({nelem, nne * ndim, nne * ndim});
 
     // DOF values
     xt::xtensor<double, 1> Fext = xt::zeros<double>({tyings.nni()});
     xt::xtensor<double, 1> Fint = xt::zeros<double>({tyings.nni()});
 
     // element/material definition
     // ---------------------------
 
     // FEM quadrature
     GF::Element::Quad4::QuadraturePlanar elem0(vector.AsElement(coor));
     GF::Element::Quad4::QuadraturePlanar elem(vector.AsElement(coor));
     size_t nip = elem0.nip();
 
     // material model
     // even though the problem is 2-d, the material model is 3-d, plane strain is implicitly assumed
-    GM::Matrix mat(nelem, nip);
-    size_t tdim = mat.ndim();
+    GM::Array<2> mat({nelem, nip});
+    size_t tdim = 3;
 
     // some artificial material definition
     xt::xtensor<size_t, 1> ehard = xt::ravel(xt::view(elmat, xt::range(0, 2), xt::range(0, 2)));
     xt::xtensor<size_t, 2> Ihard = xt::zeros<size_t>({nelem, nip});
     xt::view(Ihard, xt::keep(ehard), xt::all()) = 1ul;
     xt::xtensor<size_t, 2> Isoft = xt::ones<size_t>({nelem, nip}) - Ihard;
 
     mat.setLinearHardening(Isoft, 1.0, 1.0, 0.05, 0.05);
     mat.setElastic(Ihard, 1.0, 1.0);
 
     // solve
     // -----
 
     // allocate tensors
     xt::xtensor<double, 4> I2 = mat.I2();
     xt::xtensor<double, 4> F = xt::empty<double>({nelem, nip, tdim, tdim});
     xt::xtensor<double, 4> Eps = xt::empty<double>({nelem, nip, tdim, tdim});
     xt::xtensor<double, 4> Sig = xt::empty<double>({nelem, nip, tdim, tdim});
     xt::xtensor<double, 6> C = xt::empty<double>({nelem, nip, tdim, tdim, tdim, tdim});
 
     // allocate system matrix
     GF::MatrixPartitionedTyings K(conn, dofs, tyings.Cdu(), tyings.Cdp());
     GF::MatrixPartitionedTyingsSolver<> Solver;
 
     // allocate internal variables
     double res;
 
     // some shear strain history
     xt::xtensor<double, 1> dgamma = 0.001 * xt::ones<double>({101});
     dgamma(0) = 0.0;
 
     // initialise output
     // - output-file containing data
     HighFive::File file("main.h5", HighFive::File::Overwrite);
     // - ParaView meta-data
     PV::TimeSeries xdmf;
     // - write mesh
     H5::dump(file, "/conn", conn);
     H5::dump(file, "/coor", coor);
 
     // loop over increments
     for (size_t inc = 0; inc < dgamma.size(); ++inc) {
         // update history
         mat.increment();
 
         // iterate
         for (size_t iter = 0;; ++iter) {
             // deformation gradient tensor
             vector.asElement(disp, ue);
             elem0.gradN_vector_T(ue, F);
             F += I2;
 
             // stress & tangent
-            mat.tangent(F, Sig, C);
+            mat.setDefGrad(F);
+            mat.stress(Sig);
+            mat.tangent(C);
 
             // internal force
             elem.int_gradN_dot_tensor2_dV(Sig, fe);
             vector.assembleNode(fe, fint);
 
             // stiffness matrix
             elem.int_gradN_dot_tensor4_dot_gradNT_dV(C, Ke);
             K.assemble(Ke);
 
             // residual
             xt::noalias(fres) = fext - fint;
 
             // check for convergence (skip the zeroth iteration, as the residual still vanishes)
             if (iter > 0) {
                 // - internal/external force as DOFs (account for periodicity)
                 vector.asDofs_i(fext, Fext);
                 vector.asDofs_i(fint, Fint);
                 // - extract reaction force
                 vector.copy_p(Fint, Fext);
                 // - norm of the residual and the reaction force
                 double nfres = xt::sum(xt::abs(Fext - Fint))[0];
                 double nfext = xt::sum(xt::abs(Fext))[0];
                 // - relative residual, for convergence check
                 if (nfext)
                     res = nfres / nfext;
                 else
                     res = nfres;
                 // - print progress to screen
                 std::cout << "inc = " << inc << ", "
                           << "iter = " << iter << ", "
                           << "res = " << res << std::endl;
                 // - check for convergence
                 if (res < 1.0e-5) {
                     break;
                 }
                 // - safe-guard from infinite loop
                 if (iter > 20) {
                     throw std::runtime_error("Maximal number of iterations exceeded");
                 }
             }
 
             // initialise displacement update
             du.fill(0.0);
 
             // set fixed displacements
             if (iter == 0) {
                 du(control_nodes(0), 1) = dgamma(inc);
             }
 
             // solve
             Solver.solve(K, fres, du);
 
             // add displacement update
             disp += du;
 
             // update shape functions
-            elem.update_x(vector.AsElement(coor + disp));
+            elem.update_x(vector.AsElement(xt::eval(coor + disp)));
         }
 
         // post-process
         // - compute strain and stress
         vector.asElement(disp, ue);
         elem0.gradN_vector_T(ue, F);
         F += I2;
         GM::strain(F, Eps);
-        mat.stress(F, Sig);
+        mat.setDefGrad(F);
+        mat.stress(Sig);
         // - integration point volume
         xt::xtensor<double, 4> dV = elem.AsTensor<2>(elem.dV());
         // - element average stress
         xt::xtensor<double, 3> Sigelem = xt::average(Sig, dV, {1});
         xt::xtensor<double, 3> Epselem = xt::average(Eps, dV, {1});
         // - macroscopic strain and stress
         xt::xtensor_fixed<double, xt::xshape<3, 3>> Sigbar = xt::average(Sig, dV, {0, 1});
         xt::xtensor_fixed<double, xt::xshape<3, 3>> Epsbar = xt::average(Eps, dV, {0, 1});
         // - write to output-file: increment numbers
         H5::dump(file, "/stored", inc, {inc});
         // - write to output-file: macroscopic response
-        H5::dump(file, "/macroscopic/sigeq", GM::Sigeq(Sigbar), {inc});
-        H5::dump(file, "/macroscopic/epseq", GM::Epseq(Epsbar), {inc});
+        H5::dump(file, "/macroscopic/sigeq", GM::Sigeq(Sigbar)(), {inc});
+        H5::dump(file, "/macroscopic/epseq", GM::Epseq(Epsbar)(), {inc});
         // - write to output-file: element quantities
         H5::dump(file, "/sigeq/" + std::to_string(inc), GM::Sigeq(Sigelem));
         H5::dump(file, "/epseq/" + std::to_string(inc), GM::Epseq(Epselem));
         H5::dump(file, "/disp/" + std::to_string(inc), GF::as3d(disp));
         // - update ParaView meta-data
         xdmf.push_back({
             PV::Topology(file, "/conn", mesh.getElementType()),
             PV::Geometry(file, "/coor"),
             PV::Attribute(file, "/disp/" + std::to_string(inc), PV::AttributeCenter::Node, "Displacement"),
             PV::Attribute(file, "/sigeq/" + std::to_string(inc), PV::AttributeCenter::Cell, "Eq. stress"),
             PV::Attribute(file, "/epseq/" + std::to_string(inc), PV::AttributeCenter::Cell, "Eq. strain")});
     }
 
     // write ParaView meta-data
     PV::write("main.xdmf", xdmf.get());
 
     return 0;
 }