diff --git a/CMakeLists.txt b/CMakeLists.txt
index b4a3a8e..4bd237d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,145 +1,145 @@
 #==================================================================================================#
 #                                                                                                  #
 # (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM       #
 #                                                                                                  #
 #==================================================================================================#
 
 # initialization
 # --------------
 
 # required to specify the c++ standard
 cmake_minimum_required(VERSION 3.0)
 
 # required for install
 include(CMakePackageConfigHelpers)
 include(GNUInstallDirs)
 
 # options
 option(PKGCONFIG "Build pkg-config" ON)
 option(BUILD_TESTS "${PROJECT_NAME} test suite" OFF)
 
 # project settings
 # ----------------
 
 # name
 project(xGooseFEM)
 
 # file that contains the version information
 set(parse_version include/xGooseFEM/GooseFEM.h)
 
 # header files
 set(headers
   include/xGooseFEM/GooseFEM.h
   include/xGooseFEM/Mesh.hpp
   include/xGooseFEM/Mesh.h
   include/xGooseFEM/MeshTri3.hpp
   include/xGooseFEM/MeshTri3.h
   include/xGooseFEM/MeshQuad4.hpp
   include/xGooseFEM/MeshQuad4.h
   include/xGooseFEM/MeshHex8.hpp
   include/xGooseFEM/MeshHex8.h
   include/xGooseFEM/Element.hpp
   include/xGooseFEM/Element.h
   include/xGooseFEM/ElementQuad4.hpp
   include/xGooseFEM/ElementQuad4.h
   include/xGooseFEM/ElementQuad4Planar.hpp
   include/xGooseFEM/ElementQuad4Planar.h
   include/xGooseFEM/ElementHex8.hpp
   include/xGooseFEM/ElementHex8.h
   include/xGooseFEM/Vector.hpp
   include/xGooseFEM/Vector.h
   include/xGooseFEM/VectorPartitioned.hpp
   include/xGooseFEM/VectorPartitioned.h
   include/xGooseFEM/MatrixPartitioned.hpp
   include/xGooseFEM/MatrixPartitioned.h
   include/xGooseFEM/MatrixDiagonal.hpp
   include/xGooseFEM/MatrixDiagonal.h
   include/xGooseFEM/MatrixDiagonalPartitioned.hpp
   include/xGooseFEM/MatrixDiagonalPartitioned.h
   include/xGooseFEM/Iterate.hpp
   include/xGooseFEM/Iterate.h
   include/xGooseFEM/Dynamics.hpp
   include/xGooseFEM/Dynamics.h
 )
 
 # automatically parse the version number
 file(READ "${parse_version}" version)
 string(REGEX MATCH "define[ \t]+GOOSEFEM_WORLD_VERSION[ \t]+([0-9]+)" _ "${version}")
 set(world_version "${CMAKE_MATCH_1}")
 string(REGEX MATCH "define[ \t]+GOOSEFEM_MAJOR_VERSION[ \t]+([0-9]+)" _ "${version}")
 set(major_version "${CMAKE_MATCH_1}")
 string(REGEX MATCH "define[ \t]+GOOSEFEM_MINOR_VERSION[ \t]+([0-9]+)" _ "${version}")
 set(minor_version "${CMAKE_MATCH_1}")
 set(GOOSEFEM_VERSION ${world_version}.${major_version}.${minor_version})
 
 # print information to screen
 message(STATUS "Building ${PROJECT_NAME} v${GOOSEFEM_VERSION}")
 
 # paths
 # -----
 
 set(GOOSEFEM_ROOT_DIR        "${CMAKE_INSTALL_PREFIX}")
 set(GOOSEFEM_INCLUDE_DIR     "${CMAKE_INSTALL_PREFIX}/${INCLUDE_INSTALL_DIR}")
 set(INCLUDE_INSTALL_DIR      "${CMAKE_INSTALL_INCLUDEDIR}/${PROJECT_NAME}")
 set(CMAKEPACKAGE_INSTALL_DIR "${CMAKE_INSTALL_DATADIR}/${PROJECT_NAME}")
 set(PKGCONFIG_INSTALL_DIR    "${CMAKE_INSTALL_DATADIR}/pkgconfig")
 set(fcmake                   "GooseFEMConfig.cmake")
 set(fpkg                     "GooseFEM.pc")
 
 # configure CMake
 # ---------------
 
 configure_package_config_file(
   ${CMAKE_CURRENT_SOURCE_DIR}/GooseFEMConfig.cmake.in
   ${CMAKE_CURRENT_BINARY_DIR}/GooseFEMConfig.cmake
   PATH_VARS GOOSEFEM_INCLUDE_DIR GOOSEFEM_ROOT_DIR
   INSTALL_DESTINATION ${CMAKEPACKAGE_INSTALL_DIR}
   NO_CHECK_REQUIRED_COMPONENTS_MACRO
 )
 
 # install/build
 # -------------
 
 # build test cases
 if(BUILD_TESTS)
   add_subdirectory(test)
 endif()
 
 # disable pkg-config for native Windows builds
 if(WIN32 OR CMAKE_HOST_SYSTEM_NAME MATCHES Windows)
   option(PKGCONFIG "Build pkg-config ${fpkg} file" OFF)
 endif()
 
 # pkg-config
 if(PKGCONFIG)
   configure_file(${fpkg}.in ${fpkg} @ONLY)
   install(FILES ${CMAKE_CURRENT_BINARY_DIR}/${fpkg} DESTINATION ${PKGCONFIG_INSTALL_DIR})
 endif()
 
 # CMake
 install(FILES ${CMAKE_CURRENT_BINARY_DIR}/${fcmake} DESTINATION ${CMAKEPACKAGE_INSTALL_DIR})
 
 # header files
 install(FILES ${CMAKE_CURRENT_SOURCE_DIR}/${headers} DESTINATION ${INCLUDE_INSTALL_DIR})
 
 # print information to screen
 # ---------------------------
 
 message(STATUS "")
 message(STATUS "+---------------------------------------------------------------------------------")
 message(STATUS "|")
 message(STATUS "| Use 'make install' to install in")
 message(STATUS "|   ${CMAKE_INSTALL_PREFIX}")
 if(BUILD_TESTS)
   message(STATUS "|")
-  message(STATUS "| Use 'make' and './test' to run the tests")
+  message(STATUS "| Use 'make' and './test/test' to run the tests")
 endif()
 message(STATUS "|")
 message(STATUS "| To specify a custom directory call")
 message(STATUS "|   cmake /path/to/${PROJECT_NAME} -DCMAKE_INSTALL_PREFIX=yourprefix")
 message(STATUS "|")
 message(STATUS "| For custom paths, add the following line to your '~/.bashrc'")
 message(STATUS "|   export PKG_CONFIG_PATH=${CMAKE_INSTALL_PREFIX}/share/pkgconfig:$PKG_CONFIG_PATH")
 message(STATUS "|")
 message(STATUS "+---------------------------------------------------------------------------------")
 message(STATUS "")
diff --git a/docs/examples/Dynamics/laminateElastic/noDamping/CMakeLists.txt b/docs/examples/Dynamics/laminateElastic/noDamping/CMakeLists.txt
index bd14491..e9ed08c 100644
--- a/docs/examples/Dynamics/laminateElastic/noDamping/CMakeLists.txt
+++ b/docs/examples/Dynamics/laminateElastic/noDamping/CMakeLists.txt
@@ -1,54 +1,55 @@
 cmake_minimum_required(VERSION 3.1)
 
 # basic info
 # ----------
 
 # define a project name
 project(example)
 
 # define empty list of libraries to link
 set(PROJECT_LIBS "")
 
 
 # basic compiler options
 # ----------------------
 
 # set optimization level
 set(CMAKE_BUILD_TYPE Release)
 
 # set C++ standard
 set(CMAKE_CXX_STANDARD 14)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 
 # set warnings on
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -pedantic -g")
 
+# enable optimisation
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native")
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DXTENSOR_USE_XSIMD=ON")
 
 # load packages
 # -------------
 
 # add custom paths
 if(NOT "$ENV{INCLUDE_PATH}" STREQUAL "")
   string(REPLACE ":" ";" INCLUDE_LIST "$ENV{INCLUDE_PATH}")
   include_directories(${INCLUDE_LIST})
 endif()
 
 # load HDF5
 find_package(HDF5 COMPONENTS CXX REQUIRED)
 include_directories(${HDF5_INCLUDE_DIRS})
 set(PROJECT_LIBS ${PROJECT_LIBS} ${HDF5_LIBS} ${HDF5_LIBRARIES})
 
 # load eigen
 find_package(PkgConfig)
 pkg_check_modules(EIGEN3 REQUIRED eigen3)
 include_directories(SYSTEM ${EIGEN3_INCLUDE_DIRS})
 
 # create executable
 # -----------------
 
 add_executable(${PROJECT_NAME} main.cpp)
 
 target_link_libraries(${PROJECT_NAME} ${PROJECT_LIBS})
 
diff --git a/docs/examples/Dynamics/laminateElastic/noDamping/main.cpp b/docs/examples/Dynamics/laminateElastic/noDamping/main.cpp
index 63d696a..b432085 100644
--- a/docs/examples/Dynamics/laminateElastic/noDamping/main.cpp
+++ b/docs/examples/Dynamics/laminateElastic/noDamping/main.cpp
@@ -1,240 +1,239 @@
 
 #include <xGooseFEM/GooseFEM.h>
 #include <xElastoPlasticQPot/ElastoPlasticQPot.h>
 #include <HDF5pp.h>
 
 // -------------------------------------------------------------------------------------------------
 
 namespace GF = xGooseFEM;
 namespace QD = xGooseFEM::Element::Quad4;
 namespace GM = xElastoPlasticQPot::Cartesian2d;
 
 using T2 = xt::xtensor_fixed<double, xt::xshape<2,2>>;
 
 // -------------------------------------------------------------------------------------------------
 
 inline double sqdot(const xt::xtensor<double,1> &M, const xt::xtensor<double,1> &V)
 {
   double out = 0.;
 
   for ( size_t i = 0 ; i < M.size() ; ++i )
     out += M(i) * V(i) * V(i);
 
   return out;
 }
 
 // -------------------------------------------------------------------------------------------------
 
 class Geometry : public xGooseFEM::Dynamics::Geometry
 {
 private:
 
   // dimensions
   size_t m_nnode;
   size_t m_ndim;
   size_t m_nelem;
   size_t m_nne;
   size_t m_nip;
 
   // mesh
   xt::xtensor<size_t,2> m_conn;
   xt::xtensor<double,2> m_coor;
   xt::xtensor<double,2> m_u;
   xt::xtensor<double,2> m_v;
   xt::xtensor<double,2> m_a;
 
   // vector-definition: transform nodal vectors <-> DOF values
   GF::Vector m_vec;
 
   // mass matrix (diagonal)
   GF::MatrixDiagonal m_M;
 
   // numerical quadrature
   QD::Quadrature m_quad;
 
   // material definition
   GM::Matrix m_mat;
 
 public:
 
   // constructor
   Geometry(size_t nx);
 
   // apply update in macroscopic deformation gradient
   void add_dFbar(const T2 &dFbar);
 
   // compute total kinetic and potential energy
   double Ekin() const;
   double Epot() const;
 
   // return mesh
   xt::xtensor<size_t,2> conn() const { return m_conn; }
   xt::xtensor<double,2> coor() const { return m_coor; }
 
   // return nodal vectors
   xt::xtensor<double,2> u() const { return m_u; }
   xt::xtensor<double,2> v() const { return m_v; }
   xt::xtensor<double,2> a() const { return m_a; }
 
   // return DOF values
   xt::xtensor<double,1> dofs_v() const { return m_vec.asDofs(m_v); }
   xt::xtensor<double,1> dofs_a() const { return m_vec.asDofs(m_a); }
 
   // overwrite nodal vectors
   void set_u(const xt::xtensor<double,2> &nodevec) { m_u = nodevec; };
 
   // overwrite nodal vectors, reconstructed from DOF values
   void set_v(const xt::xtensor<double,1> &dofval) { m_v = m_vec.asNode(dofval); };
   void set_a(const xt::xtensor<double,1> &dofval) { m_a = m_vec.asNode(dofval); };
 
   // solve for DOF-accelerations
   xt::xtensor<double,1> solve_A();
 };
 
 // -------------------------------------------------------------------------------------------------
 
 inline Geometry::Geometry(size_t nx)
 {
   // get mesh
   GF::Mesh::Quad4::Regular mesh(nx,nx,1.);
 
   // vector-definition
   m_vec   = GF::Vector(mesh.conn(), mesh.dofsPeriodic());
 
   // quadrature
   m_quad  = QD::Quadrature(m_vec.asElement(mesh.coor()));
 
   // dimensions, connectivity, and coordinates
   m_nnode = mesh.nnode();
   m_ndim  = mesh.ndim();
   m_nelem = mesh.nelem();
   m_nne   = mesh.nne();
   m_nip   = m_quad.nip();
   m_conn  = mesh.conn();
   m_coor  = mesh.coor();
 
   // zero-initialize displacement, velocity, acceleration
   m_u = xt::zeros<double>({m_nnode, m_ndim});
   m_v = xt::zeros<double>({m_nnode, m_ndim});
   m_a = xt::zeros<double>({m_nnode, m_ndim});
 
   // material definition
   // - allocate
   m_mat = GM::Matrix({m_nelem, m_nip});
   // - phase indicators
   xt::xtensor<size_t,2> Ihard = xt::zeros<size_t>({m_nelem, m_nip});
   xt::xtensor<size_t,2> Isoft = xt::ones <size_t>({m_nelem, m_nip});
   // - set hard indicator
   for ( size_t e = 0 ; e < nx*nx/4 ; ++e )
     for ( size_t k = 0 ; k < m_nip ; ++k )
       Ihard(e,k) = 1;
   // - set soft indicator
   Isoft -= Ihard;
   // - set material definition
   m_mat.setElastic(Ihard, 100., 10.);
   m_mat.setElastic(Isoft, 100.,  1.);
   // - check that all points have been set (not strictly needed)
   m_mat.check();
 
   // mass matrix
   // - nodal quadrature
   QD::Quadrature q(m_vec.asElement(m_coor), QD::Nodal::xi(), QD::Nodal::w());
   // - set density
   xt::xtensor<double,2> rho = 1.0 * xt::ones<double>({m_nelem, m_nip});
   // - allocate mass matrix
   m_M = GF::MatrixDiagonal(m_conn, mesh.dofsPeriodic());
   // - compute (constant hereafter)
   m_M.assemble( q.int_N_scalar_NT_dV(rho) );
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline xt::xtensor<double,1> Geometry::solve_A()
 {
   xt::xtensor<double,4> Eps = m_quad.symGradN_vector( m_vec.asElement(m_u) );
   xt::xtensor<double,4> Sig = m_mat.Sig(Eps);
-  xt::xtensor<double,1> F = m_vec.assembleDofs( m_quad.int_gradN_dot_tensor2_dV(Sig) );
+  xt::xtensor<double,1> F   = m_vec.assembleDofs( m_quad.int_gradN_dot_tensor2_dV(Sig) );
 
   return m_M.solve( -F );
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline void Geometry::add_dFbar(const T2 &dFbar)
 {
   for ( size_t i = 0 ; i < m_nnode ; ++i )
     for ( size_t j = 0 ; j < m_ndim ; ++j )
       for ( size_t k = 0 ; k < m_ndim ; ++k )
         m_u(i,j) += dFbar(j,k) * ( m_coor(i,k) - m_coor(0,k) );
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline double Geometry::Ekin() const
 {
   xt::xtensor<double,1> V = m_vec.asDofs(m_v);
   xt::xtensor<double,1> M = m_M.asDiagonal();
 
   return 0.5 * sqdot(M,V);
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline double Geometry::Epot() const
 {
   xt::xtensor<double,4> Eps = m_quad.symGradN_vector( m_vec.asElement(m_u) );
   xt::xtensor<double,2> E   = m_mat.energy(Eps);
   xt::xtensor<double,2> dV  = m_quad.dV();
 
   return ( xt::sum( E*dV ) )[0];
 }
 
 // -------------------------------------------------------------------------------------------------
 
 int main()
 {
   // set simulation parameters
   double T     = 60.  ; // total time
   double dt    = 1.e-2; // time increment
   size_t nx    = 60   ; // number of elements in both directions
   double gamma = .05  ; // displacement step
 
   // define geometry
   Geometry geometry(nx);
 
   // define update in macroscopic deformation gradient
   T2 dFbar = xt::zeros<double>({2,2});
   dFbar(0,1) = gamma;
 
   // update displacement
   geometry.add_dFbar(dFbar);
 
   // output variables
   xt::xtensor<double,1> Epot = xt::zeros<double>({static_cast<size_t>(T/dt)});
   xt::xtensor<double,1> Ekin = xt::zeros<double>({static_cast<size_t>(T/dt)});
   xt::xtensor<double,1> t    = xt::zeros<double>({static_cast<size_t>(T/dt)});
 
   // loop over increments
   for ( size_t inc = 0 ; inc < static_cast<size_t>(Epot.size()) ; ++inc )
-  // for ( size_t inc = 0 ; inc < 2000 ; ++inc )
   {
     // - compute increment
     GF::Dynamics::velocityVerlet(geometry, dt);
 
     // - store output
     t   (inc) = static_cast<double>(inc) * dt;
     Ekin(inc) = geometry.Ekin();
     Epot(inc) = geometry.Epot();
   }
 
   // write to output file
   H5p::File f = H5p::File("example.hdf5","w");
   f.write("/global/Epot",Epot           );
   f.write("/global/Ekin",Ekin           );
   f.write("/global/t"   ,t              );
   f.write("/mesh/conn"  ,geometry.conn());
   f.write("/mesh/coor"  ,geometry.coor());
   f.write("/mesh/disp"  ,geometry.u()   );
 
   return 0;
 }
diff --git a/docs/examples/Dynamics/laminateElastic/noDamping/plot.py b/docs/examples/Dynamics/laminateElastic/noDamping/plot.py
index 07605d8..55a92f3 100644
--- a/docs/examples/Dynamics/laminateElastic/noDamping/plot.py
+++ b/docs/examples/Dynamics/laminateElastic/noDamping/plot.py
@@ -1,52 +1,52 @@
 
 import sys, os, re
 import h5py
 import numpy             as np
 import matplotlib.pyplot as plt
-import goosempl          as gplt
+import GooseMPL          as gplt
 import matplotlib        as mpl
 
 plt.style.use(['goose','goose-latex'])
 
 # --------------------------------------------------------------------------------------------------
 
 name = 'example.hdf5'
 f    = h5py.File(name,'r')
 
 fig, axes = plt.subplots(figsize=(20,10), ncols=2)
 
 # --------------------------------------------------------------------------------------------------
 
 ax = axes[0]
 
 t = f['/global/t'   ][:]
 E = f['/global/Ekin'][:] + f['/global/Epot'][:]; ax.plot(t,E,color='k',label=r'$E_\mathrm{pot} + E_\mathrm{kin}$')
 E = f['/global/Epot'][:]                       ; ax.plot(t,E,color='r',label=r'$E_\mathrm{pot}$')
 E = f['/global/Ekin'][:]                       ; ax.plot(t,E,color='b',label=r'$E_\mathrm{kin}$')
 
 ax.legend(loc='lower right')
 
 ax.set_xlabel(r'$t$')
 ax.set_ylabel(r'$E$')
 
 # --------------------------------------------------------------------------------------------------
 
 ax = axes[1]
 
 coor = f['/mesh/coor'][:]
 conn = f['/mesh/conn'][:]
 u    = f['/mesh/disp'][:]
 
 
 im = gplt.patch(coor=coor+u,conn=conn)
 
 Lx = np.max(coor[:,0])-np.min(coor[:,0])
 Ly = np.max(coor[:,1])-np.min(coor[:,1])
 
 ax.set_xlim([-0.1*Lx,1.3*Lx])
 ax.set_ylim([-0.1*Ly,1.3*Ly])
 
 ax.set_aspect('equal')
 
 # plt.savefig(re.sub(r'(.*)(\.hdf5)',r'\1.svg',name))
 plt.show()