diff --git a/docs/examples/CMakeLists.txt b/docs/examples/CMakeLists.txt
index 67c6041..898a4c5 100644
--- a/docs/examples/CMakeLists.txt
+++ b/docs/examples/CMakeLists.txt
@@ -1,81 +1,84 @@
 
-cmake_minimum_required(VERSION 3.0)
+cmake_minimum_required(VERSION 3.16)
 
 if(CMAKE_CURRENT_SOURCE_DIR STREQUAL CMAKE_SOURCE_DIR)
     project(GooseFEM-examples)
     find_package(GooseFEM REQUIRED CONFIG)
 endif()
 
 option(WARNINGS "Show build warnings" ON)
 option(SIMD "Enable xsimd" ON)
 option(DEBUG "Enable all assertions" OFF)
 
-find_package(HDF5 REQUIRED)
+set(HIGHFIVE_USE_BOOST 0)
+set(HIGHFIVE_USE_XTENSOR 1)
+find_package(HighFive REQUIRED)
+find_package(XDMFWrite_HighFive REQUIRED)
 
 add_library(libraries INTERFACE IMPORTED)
 
-target_link_libraries(libraries INTERFACE GooseFEM ${HDF5_LIBRARIES})
+target_link_libraries(libraries INTERFACE GooseFEM HighFive XDMFWrite_HighFive)
 
 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 (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)
 
 # 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)
 
 # 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)
 
 # 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)
 
 # 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)
 
 # 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)
 
 # 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)
 
 
diff --git a/docs/examples/dynamics/Elastic-VelocityVerlet/main.cpp b/docs/examples/dynamics/Elastic-VelocityVerlet/main.cpp
index fe3ff87..bd2d2d2 100644
--- a/docs/examples/dynamics/Elastic-VelocityVerlet/main.cpp
+++ b/docs/examples/dynamics/Elastic-VelocityVerlet/main.cpp
@@ -1,180 +1,179 @@
-
 #include <GMatElastoPlasticQPot/Cartesian2d.h>
 #include <GooseFEM/GooseFEM.h>
 #include <highfive/H5Easy.hpp>
 
 namespace GM = GMatElastoPlasticQPot::Cartesian2d;
 namespace GF = GooseFEM;
 namespace QD = GooseFEM::Element::Quad4;
 namespace H5 = H5Easy;
 
 inline double sqdot(const xt::xtensor<double, 1>& M, const xt::xtensor<double, 1>& V)
 {
     double ret = 0.;
 
     for (size_t i = 0; i < M.size(); ++i) {
         ret += M(i) * V(i) * V(i);
     }
 
     return ret;
 }
 
 int main()
 {
     // simulation parameters
 
     double T = 60.0;     // total time
     double dt = 1.0e-2;  // time increment
     size_t nx = 60;      // number of elements in both directions
     double gamma = 0.05; // displacement step
 
     // get mesh & quadrature
 
     GF::Mesh::Quad4::Regular mesh(nx, nx, 1.0);
 
     xt::xtensor<double, 2> coor = mesh.coor();
     xt::xtensor<size_t, 2> conn = mesh.conn();
     xt::xtensor<size_t, 2> dofs = mesh.dofsPeriodic();
 
     GF::Vector vector(conn, dofs);
 
     QD::Quadrature quad(vector.AsElement(coor));
 
     size_t ndim = mesh.ndim();
     size_t nne = mesh.nne();
     size_t nelem = mesh.nelem();
     size_t nip = quad.nip();
 
     xt::xtensor<double, 2> u = xt::zeros<double>(coor.shape());
     xt::xtensor<double, 2> v = xt::zeros<double>(coor.shape());
     xt::xtensor<double, 2> a = xt::zeros<double>(coor.shape());
     xt::xtensor<double, 2> v_n = xt::zeros<double>(coor.shape());
     xt::xtensor<double, 2> a_n = xt::zeros<double>(coor.shape());
     xt::xtensor<double, 2> fint = xt::zeros<double>(coor.shape());
     xt::xtensor<double, 2> fext = xt::zeros<double>(coor.shape());
     xt::xtensor<double, 2> fres = xt::zeros<double>(coor.shape());
 
     xt::xtensor<double, 3> ue = xt::zeros<double>({nelem, nne, ndim});
     xt::xtensor<double, 3> fe = xt::zeros<double>({nelem, nne, ndim});
 
     xt::xtensor<double, 4> Eps = xt::zeros<double>({nelem, nip, ndim, ndim});
     xt::xtensor<double, 4> Sig = xt::zeros<double>({nelem, nip, ndim, ndim});
 
     // material definition
 
     GM::Matrix material({nelem, nip});
 
     xt::xtensor<size_t, 2> Ihard = xt::zeros<size_t>({nelem, nip});
     xt::xtensor<size_t, 2> Isoft = xt::ones<size_t>({nelem, nip});
     xt::view(Ihard, xt::range(0, nx * nx / 4), xt::all()) = 1ul;
     Isoft -= Ihard;
 
     material.setElastic(Ihard, 100.0, 10.0);
     material.setElastic(Isoft, 100.0, 1.0);
 
     material.check();
 
     // mass matrix
 
     QD::Quadrature nodalQuad(vector.AsElement(coor), QD::Nodal::xi(), QD::Nodal::w());
 
     xt::xtensor<double, 2> rho = 1.0 * xt::ones<double>({nelem, nip});
 
     GF::MatrixDiagonal M(conn, dofs);
 
     M.assemble(nodalQuad.Int_N_scalar_NT_dV(rho));
 
     xt::xtensor<double, 1> mass = M.AsDiagonal();
 
     // update in macroscopic deformation gradient
 
     xt::xtensor<double, 2> dFbar = xt::zeros<double>({2, 2});
 
     dFbar(0, 1) = gamma;
 
     for (size_t j = 0; j < ndim; ++j) {
         for (size_t k = 0; k < ndim; ++k) {
             xt::view(u, xt::all(), j) += dFbar(j, k) * (xt::view(coor, xt::all(), k) - coor(0, k));
         }
     }
 
     // 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)});
 
     xt::xtensor<double, 2> dV = quad.DV();
 
     // loop over increments
 
     for (size_t inc = 0; inc < static_cast<size_t>(Epot.size()); ++inc) {
 
         // store history
 
         xt::noalias(v_n) = v;
         xt::noalias(a_n) = a;
 
         // new displacement
 
         xt::noalias(u) = u + dt * v + 0.5 * std::pow(dt, 2.0) * a;
 
         // compute strain/strain, and corresponding internal force
 
         vector.asElement(u, ue);
         quad.symGradN_vector(ue, Eps);
         material.stress(Eps, Sig);
         quad.int_gradN_dot_tensor2_dV(Sig, fe);
         vector.assembleNode(fe, fint);
 
         // estimate new velocity
 
         xt::noalias(v) = v_n + dt * a_n;
 
         // compute residual force & solve
 
         xt::noalias(fres) = fext - fint;
 
         M.solve(fres, a);
 
         // re-estimate new velocity
 
         xt::noalias(v) = v_n + 0.5 * dt * (a_n + a);
 
         // compute residual force & solve
 
         xt::noalias(fres) = fext - fint;
 
         M.solve(fres, a);
 
         // new velocity
 
         xt::noalias(v) = v_n + 0.5 * dt * (a_n + a);
 
         // compute residual force & solve
 
         xt::noalias(fres) = fext - fint;
 
         M.solve(fres, a);
 
         // store output variables
 
         xt::xtensor<double, 2> E = material.Energy(Eps);
         xt::xtensor<double, 1> V = vector.AsDofs(v);
 
         t(inc) = static_cast<double>(inc) * dt;
         Ekin(inc) = 0.5 * sqdot(mass, V);
         Epot(inc) = xt::sum(E * dV)[0];
     }
 
     // write output variables to file
     H5::File file("example.hdf5", H5::File::Overwrite);
     H5::dump(file, "/global/Epot", Epot);
     H5::dump(file, "/global/Ekin", Ekin);
     H5::dump(file, "/global/t", t);
     H5::dump(file, "/mesh/conn", conn);
     H5::dump(file, "/mesh/coor", coor);
     H5::dump(file, "/mesh/disp", u);
 
     return 0;
 }
diff --git a/docs/examples/dynamics/Elastic-Verlet/main.cpp b/docs/examples/dynamics/Elastic-Verlet/main.cpp
index 6f40094..8ba9399 100644
--- a/docs/examples/dynamics/Elastic-Verlet/main.cpp
+++ b/docs/examples/dynamics/Elastic-Verlet/main.cpp
@@ -1,160 +1,159 @@
-
 #include <GMatElastoPlasticQPot/Cartesian2d.h>
 #include <GooseFEM/GooseFEM.h>
 #include <highfive/H5Easy.hpp>
 
 namespace GM = GMatElastoPlasticQPot::Cartesian2d;
 namespace GF = GooseFEM;
 namespace QD = GooseFEM::Element::Quad4;
 namespace H5 = H5Easy;
 
 inline double sqdot(const xt::xtensor<double, 1>& M, const xt::xtensor<double, 1>& V)
 {
     double ret = 0.0;
 
     for (size_t i = 0; i < M.size(); ++i) {
         ret += M(i) * V(i) * V(i);
     }
 
     return ret;
 }
 
 int main()
 {
     // simulation parameters
 
     double T = 60.0;     // total time
     double dt = 1.0e-2;  // time increment
     size_t nx = 60;      // number of elements in both directions
     double gamma = 0.05; // displacement step
 
     // get mesh & quadrature
 
     GF::Mesh::Quad4::Regular mesh(nx, nx, 1.0);
 
     xt::xtensor<double, 2> coor = mesh.coor();
     xt::xtensor<size_t, 2> conn = mesh.conn();
     xt::xtensor<size_t, 2> dofs = mesh.dofsPeriodic();
 
     GF::Vector vector(conn, dofs);
 
     QD::Quadrature quad(vector.AsElement(coor));
 
     size_t ndim = mesh.ndim();
     size_t nne = mesh.nne();
     size_t nelem = mesh.nelem();
     size_t nip = quad.nip();
 
     xt::xtensor<double, 2> u = xt::zeros<double>(coor.shape());
     xt::xtensor<double, 2> v = xt::zeros<double>(coor.shape());
     xt::xtensor<double, 2> a = xt::zeros<double>(coor.shape());
     xt::xtensor<double, 2> v_n = xt::zeros<double>(coor.shape());
     xt::xtensor<double, 2> a_n = xt::zeros<double>(coor.shape());
     xt::xtensor<double, 2> fint = xt::zeros<double>(coor.shape());
     xt::xtensor<double, 2> fext = xt::zeros<double>(coor.shape());
     xt::xtensor<double, 2> fres = xt::zeros<double>(coor.shape());
 
     xt::xtensor<double, 3> ue = xt::zeros<double>({nelem, nne, ndim});
     xt::xtensor<double, 3> fe = xt::zeros<double>({nelem, nne, ndim});
 
     xt::xtensor<double, 4> Eps = xt::zeros<double>({nelem, nip, ndim, ndim});
     xt::xtensor<double, 4> Sig = xt::zeros<double>({nelem, nip, ndim, ndim});
 
     // material definition
 
     GM::Matrix material({nelem, nip});
 
     xt::xtensor<size_t, 2> Ihard = xt::zeros<size_t>({nelem, nip});
     xt::xtensor<size_t, 2> Isoft = xt::ones<size_t>({nelem, nip});
     xt::view(Ihard, xt::range(0, nx * nx / 4), xt::all()) = 1ul;
     Isoft -= Ihard;
 
     material.setElastic(Ihard, 100.0, 10.0);
     material.setElastic(Isoft, 100.0, 1.0);
 
     material.check();
 
     // mass matrix
 
     QD::Quadrature nodalQuad(vector.AsElement(coor), QD::Nodal::xi(), QD::Nodal::w());
 
     xt::xtensor<double, 2> rho = 1.0 * xt::ones<double>({nelem, nip});
 
     GF::MatrixDiagonal M(conn, dofs);
 
     M.assemble(nodalQuad.Int_N_scalar_NT_dV(rho));
 
     xt::xtensor<double, 1> mass = M.AsDiagonal();
 
     // update in macroscopic deformation gradient
 
     xt::xtensor<double, 2> dFbar = xt::zeros<double>({2, 2});
 
     dFbar(0, 1) = gamma;
 
     for (size_t j = 0; j < ndim; ++j) {
         for (size_t k = 0; k < ndim; ++k) {
             xt::view(u, xt::all(), j) += dFbar(j, k) * (xt::view(coor, xt::all(), k) - coor(0, k));
         }
     }
 
     // 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)});
 
     xt::xtensor<double, 2> dV = quad.DV();
 
     // loop over increments
 
     for (size_t inc = 0; inc < static_cast<size_t>(Epot.size()); ++inc) {
 
         // store history
 
         xt::noalias(v_n) = v;
         xt::noalias(a_n) = a;
 
         // new displacement
 
         xt::noalias(u) = u + dt * v + 0.5 * std::pow(dt, 2.0) * a;
 
         // compute strain/strain, and corresponding internal force
 
         vector.asElement(u, ue);
         quad.symGradN_vector(ue, Eps);
         material.stress(Eps, Sig);
         quad.int_gradN_dot_tensor2_dV(Sig, fe);
         vector.assembleNode(fe, fint);
 
         // compute residual force & solve
 
         xt::noalias(fres) = fext - fint;
 
         M.solve(fres, a);
 
         // new velocity
 
         xt::noalias(v) = v_n + 0.5 * dt * (a_n + a);
 
         // store output variables
 
         xt::xtensor<double, 2> E = material.Energy(Eps);
         xt::xtensor<double, 1> V = vector.AsDofs(v);
 
         t(inc) = static_cast<double>(inc) * dt;
         Ekin(inc) = 0.5 * sqdot(mass, V);
         Epot(inc) = xt::sum(E * dV)[0];
     }
 
     // write output variables to file
     H5::File file("example.hdf5", H5::File::Overwrite);
     H5::dump(file, "/global/Epot", Epot);
     H5::dump(file, "/global/Ekin", Ekin);
     H5::dump(file, "/global/t", t);
     H5::dump(file, "/mesh/conn", conn);
     H5::dump(file, "/mesh/coor", coor);
     H5::dump(file, "/mesh/disp", u);
 
     return 0;
 }
diff --git a/docs/examples/statics/Periodic_ElastoPlastic/main.cpp b/docs/examples/statics/Periodic_ElastoPlastic/main.cpp
index 627a7c8..2f53be2 100644
--- a/docs/examples/statics/Periodic_ElastoPlastic/main.cpp
+++ b/docs/examples/statics/Periodic_ElastoPlastic/main.cpp
@@ -1,232 +1,229 @@
 #include <GMatElastoPlastic/Cartesian3d.h>
 #include <GooseFEM/GooseFEM.h>
 #include <GooseFEM/MatrixPartitionedTyings.h>
 #include <GooseFEM/VectorPartitionedTyings.h>
 #include <GooseFEM/TyingsPeriodic.h>
-#include <GooseFEM/ParaView.h>
+#include <XDMFWrite_HighFive.hpp>
 #include <highfive/H5Easy.hpp>
 
 namespace GM = GMatElastoPlastic::Cartesian3d;
 namespace GF = GooseFEM;
-namespace PV = GooseFEM::ParaView::HDF5;
+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.elementMatrix();
 
     // 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 elem(vector.AsElement(coor));
     size_t nip = elem.nip();
     xt::xtensor<double, 4> dV = elem.AsTensor<2>(elem.dV());
 
     // 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();
 
     // 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> 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) {
             // strain
             vector.asElement(disp, ue);
             elem.symGradN_vector(ue, Eps);
 
             // stress & tangent
             mat.tangent(Eps, Sig, 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;
+                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;
         }
 
         // post-process
         // - compute strain and stress
         vector.asElement(disp, ue);
         elem.symGradN_vector(ue, Eps);
         mat.stress(Eps, Sig);
         // - 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});
         // - 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), PV::as3d(disp));
+        H5::dump(file, "/disp/" + std::to_string(inc), GF::as3d(disp));
         // - update ParaView meta-data
-        xdmf.push_back(PV::Increment(
-            PV::Connectivity(file, "/conn", mesh.getElementType()),
-            PV::Coordinates(file, "/coor"),
-            {
-                PV::Attribute(file, "/disp/" + std::to_string(inc), "Displacement", PV::AttributeType::Node),
-                PV::Attribute(file, "/sigeq/" + std::to_string(inc), "Eq. stress", PV::AttributeType::Cell),
-                PV::Attribute(file, "/epseq/" + std::to_string(inc), "Eq. strain", PV::AttributeType::Cell),
-            }));
+        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
-    xdmf.write("main.xdmf");
+    PV::write("main.xdmf", xdmf.get());
 
     return 0;
 }
diff --git a/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp b/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp
index 9588c44..c744575 100644
--- a/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp
+++ b/docs/examples/statics/Periodic_ElastoPlasticFiniteStrainSimo/main.cpp
@@ -1,240 +1,237 @@
 #include <Eigen/Eigen>
 #include <GMatElastoPlasticFiniteStrainSimo/Cartesian3d.h>
 #include <GooseFEM/GooseFEM.h>
-#include <GooseFEM/ParaView.h>
+#include <XDMFWrite_HighFive.hpp>
 #include <highfive/H5Easy.hpp>
 
 namespace GM = GMatElastoPlasticFiniteStrainSimo::Cartesian3d;
 namespace GF = GooseFEM;
-namespace PV = GooseFEM::ParaView::HDF5;
+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.elementMatrix();
 
     // 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();
 
     // 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);
 
             // 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;
+                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));
         }
 
         // 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);
         // - 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});
         // - 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), PV::as3d(disp));
+        H5::dump(file, "/disp/" + std::to_string(inc), GF::as3d(disp));
         // - update ParaView meta-data
-        xdmf.push_back(PV::Increment(
-            PV::Connectivity(file, "/conn", mesh.getElementType()),
-            PV::Coordinates(file, "/coor"),
-            {
-                PV::Attribute(file, "/disp/" + std::to_string(inc), "Displacement", PV::AttributeType::Node),
-                PV::Attribute(file, "/sigeq/" + std::to_string(inc), "Eq. stress", PV::AttributeType::Cell),
-                PV::Attribute(file, "/epseq/" + std::to_string(inc), "Eq. strain", PV::AttributeType::Cell),
-            }));
+        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
-    xdmf.write("main.xdmf");
+    PV::write("main.xdmf", xdmf.get());
 
     return 0;
 }
diff --git a/docs/examples/statics/Periodic_LinearElastic/main.cpp b/docs/examples/statics/Periodic_LinearElastic/main.cpp
index d0e9891..0a46b7b 100644
--- a/docs/examples/statics/Periodic_LinearElastic/main.cpp
+++ b/docs/examples/statics/Periodic_LinearElastic/main.cpp
@@ -1,168 +1,166 @@
 #include <Eigen/Eigen>
 #include <GMatElastic/Cartesian3d.h>
 #include <GooseFEM/GooseFEM.h>
-#include <GooseFEM/ParaView.h>
+#include <XDMFWrite_HighFive.hpp>
 #include <highfive/H5Easy.hpp>
 
 namespace GM = GMatElastic::Cartesian3d;
 namespace GF = GooseFEM;
-namespace PV = GooseFEM::ParaView::HDF5;
+namespace PV = XDMFWrite_HighFive;
 namespace H5 = H5Easy;
 
 int main()
 {
     // mesh
     // ----
 
     // define mesh
     GF::Mesh::Quad4::Regular mesh(5 * 10, 5 * 10);
 
     // 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.elementMatrix();
 
     // 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> 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});
 
     // element/material definition
     // ---------------------------
 
     // FEM quadrature
     GF::Element::Quad4::QuadraturePlanar elem(vector.AsElement(coor));
     size_t nip = elem.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();
 
     // some artificial material definition
     xt::xtensor<size_t, 1> ehard = xt::ravel(xt::view(elmat, xt::range(0, 2 * 10), xt::range(0, 2 * 10)));
     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.setElastic(Isoft, 10.0, 1.0);
     mat.setElastic(Ihard, 10.0, 10.0);
 
     // solve
     // -----
 
     // allocate tensors
     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;
 
     // strain
     vector.asElement(disp, ue);
     elem.symGradN_vector(ue, Eps);
 
     // stress & tangent
     mat.tangent(Eps, Sig, 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;
 
     // set fixed displacements
     disp(control_nodes(0), 1) = 0.1;
 
     // solve
     Solver.solve(K, fres, disp);
 
     // post-process
     // - 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);
     // - integration point volume
     xt::xtensor<double, 4> dV = elem.AsTensor<2>(elem.dV());
     // - compute strain and stress
     vector.asElement(disp, ue);
     elem.symGradN_vector(ue, Eps);
     mat.stress(Eps, Sig);
     // - 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: macroscopic response
     H5::dump(file, "/macroscopic/sigeq", GM::Sigeq(Sigbar));
     H5::dump(file, "/macroscopic/epseq", GM::Epseq(Epsbar));
     // - write to output-file: element quantities
     H5::dump(file, "/sigeq", GM::Sigeq(Sigelem));
     H5::dump(file, "/epseq", GM::Epseq(Epselem));
-    H5::dump(file, "/disp", PV::as3d(disp));
+    H5::dump(file, "/disp", GF::as3d(disp));
     // - update ParaView meta-data
-    xdmf.push_back(PV::Increment(
-        PV::Connectivity(file, "/conn", mesh.getElementType()),
-        PV::Coordinates(file, "/coor"),
-        {
-            PV::Attribute(file, "/disp", "Displacement", PV::AttributeType::Node),
-            PV::Attribute(file, "/sigeq", "Eq. stress", PV::AttributeType::Cell),
-            PV::Attribute(file, "/epseq", "Eq. strain", PV::AttributeType::Cell),
-        }));
+    xdmf.push_back({
+        PV::Topology(file, "/conn", mesh.getElementType()),
+        PV::Geometry(file, "/coor"),
+        PV::Attribute(file, "/disp", PV::AttributeCenter::Node, "Displacement"),
+        PV::Attribute(file, "/sigeq", PV::AttributeCenter::Cell, "Eq. stress"),
+        PV::Attribute(file, "/epseq", PV::AttributeCenter::Cell, "Eq. strain")});
 
     // write ParaView meta-data
-    xdmf.write("main.xdmf");
+    PV::write("main.xdmf", xdmf.get());
 
     return 0;
 }
diff --git a/docs/examples/statics/Periodic_NonLinearElastic/main.cpp b/docs/examples/statics/Periodic_NonLinearElastic/main.cpp
index f529ef3..a3b6fe7 100644
--- a/docs/examples/statics/Periodic_NonLinearElastic/main.cpp
+++ b/docs/examples/statics/Periodic_NonLinearElastic/main.cpp
@@ -1,214 +1,213 @@
 #include <Eigen/Eigen>
 #include <GMatNonLinearElastic/Cartesian3d.h>
 #include <GooseFEM/GooseFEM.h>
-#include <GooseFEM/ParaView.h>
+#include <XDMFWrite_HighFive.hpp>
 #include <highfive/H5Easy.hpp>
 
 namespace GM = GMatNonLinearElastic::Cartesian3d;
 namespace GF = GooseFEM;
-namespace PV = GooseFEM::ParaView::HDF5;
+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.elementMatrix();
 
     // 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 elem(vector.AsElement(coor));
     size_t nip = elem.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();
 
     // 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.setNonLinearElastic(Isoft, 10.0, 0.1, 0.1, 2.0);
     mat.setNonLinearElastic(Ihard, 10.0, 1.0, 0.1, 2.0);
 
     // solve
     // -----
 
     // allocate tensors
     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;
 
     // iterate
     for (size_t iter = 0;; ++iter) {
         // strain
         vector.asElement(disp, ue);
         elem.symGradN_vector(ue, Eps);
 
         // stress & tangent
         mat.tangent(Eps, Sig, 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 << "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) = 0.1;
         }
 
         // solve
         Solver.solve(K, fres, du);
 
         // add displacement update
         disp += du;
     }
 
     // post-process
     // - 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);
     // - integration point volume
     xt::xtensor<double, 4> dV = elem.AsTensor<2>(elem.dV());
     // - compute strain and stress
     vector.asElement(disp, ue);
     elem.symGradN_vector(ue, Eps);
     mat.stress(Eps, Sig);
     // - 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: macroscopic response
     H5::dump(file, "/macroscopic/sigeq", GM::Sigeq(Sigbar));
     H5::dump(file, "/macroscopic/epseq", GM::Epseq(Epsbar));
     // - write to output-file: element quantities
     H5::dump(file, "/sigeq", GM::Sigeq(Sigelem));
     H5::dump(file, "/epseq", GM::Epseq(Epselem));
-    H5::dump(file, "/disp", PV::as3d(disp));
+    H5::dump(file, "/disp", GF::as3d(disp));
     // - update ParaView meta-data
-    xdmf.push_back(PV::Increment(
-        PV::Connectivity(file, "/conn", mesh.getElementType()),
-        PV::Coordinates(file, "/coor"),
-        {
-            PV::Attribute(file, "/disp", "Displacement", PV::AttributeType::Node),
-            PV::Attribute(file, "/sigeq", "Eq. stress", PV::AttributeType::Cell),
-            PV::Attribute(file, "/epseq", "Eq. strain", PV::AttributeType::Cell),
-        }));
+    xdmf.push_back({
+        PV::Topology(file, "/conn", mesh.getElementType()),
+        PV::Geometry(file, "/coor"),
+        PV::Attribute(file, "/disp", PV::AttributeCenter::Node, "Displacement"),
+        PV::Attribute(file, "/sigeq", PV::AttributeCenter::Cell, "Eq. stress"),
+        PV::Attribute(file, "/epseq", PV::AttributeCenter::Cell, "Eq. strain")});
+
 
     // write ParaView meta-data
-    xdmf.write("main.xdmf");
+    PV::write("main.xdmf", xdmf.get());
 
     return 0;
 }
diff --git a/include/GooseFEM/Allocate.h b/include/GooseFEM/Allocate.h
index 0f80207..6333321 100644
--- a/include/GooseFEM/Allocate.h
+++ b/include/GooseFEM/Allocate.h
@@ -1,36 +1,41 @@
 /*
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 */
 
 #ifndef GOOSEFEM_ALLOCATE_H
 #define GOOSEFEM_ALLOCATE_H
 
 #include "config.h"
 
 namespace GooseFEM {
 
-// "Broadcast"
+// "Broadcast" a scalar stored in an array (e.g. [m, n]) to the same scalar of all tensor components
+// of a tensor of certain rank (e.g. for rank 2 [m, n, i, j])
 
 template <size_t dim, size_t rank>
 inline void asTensor(const xt::xtensor<double, dim>& arg, xt::xtensor<double, dim + rank>& ret);
 
 template <size_t dim, size_t rank>
 inline xt::xtensor<double, dim + rank>
 AsTensor(const xt::xtensor<double, dim>& arg, const std::array<size_t, rank>& shape);
 
 template <size_t dim, size_t rank>
 inline xt::xtensor<double, dim + rank> AsTensor(const xt::xtensor<double, dim>& arg, size_t n);
 
 template <class T>
 inline xt::xarray<double> AsTensor(size_t rank, const T& arg, const std::vector<size_t>& shape);
 
 template <class T>
 inline xt::xarray<double> AsTensor(size_t rank, const T& arg, size_t n);
 
+// Zero-pad columns to a matrix until is that shape [m, 3]
+
+inline xt::xtensor<double, 2> as3d(const xt::xtensor<double, 2>& data);
+
 } // namespace GooseFEM
 
 #include "Allocate.hpp"
 
 #endif
diff --git a/include/GooseFEM/Allocate.hpp b/include/GooseFEM/Allocate.hpp
index 88d4a7a..3543239 100644
--- a/include/GooseFEM/Allocate.hpp
+++ b/include/GooseFEM/Allocate.hpp
@@ -1,93 +1,114 @@
 /*
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 */
 
 #ifndef GOOSEFEM_ALLOCATE_HPP
 #define GOOSEFEM_ALLOCATE_HPP
 
 #include "Allocate.h"
 
 namespace GooseFEM {
 
 namespace detail {
 
 template <class T, class S>
 inline bool has_shape_begin(const T& t, const S& s)
 {
     return s.dimension() >= t.dimension() &&
            std::equal(t.shape().cbegin(), t.shape().cend(), s.shape().begin());
 }
 
 } // namespace detail
 
 template <size_t dim, size_t rank>
 inline void asTensor(const xt::xtensor<double, dim>& arg, xt::xtensor<double, dim + rank>& ret)
 {
     using strides_type = typename xt::xtensor<double, dim>::strides_type::value_type;
     GOOSEFEM_ASSERT(detail::has_shape_begin(arg, ret));
     std::array<strides_type, dim + rank> ret_strides;
     std::copy(arg.strides().begin(), arg.strides().end(), ret_strides.begin());
     std::fill(ret_strides.begin() + dim, ret_strides.end(), 0);
     ret = xt::strided_view(arg, ret.shape(), std::move(ret_strides), 0ul, xt::layout_type::dynamic);
 }
 
 template <size_t dim, size_t rank>
 inline xt::xtensor<double, dim + rank>
 AsTensor(const xt::xtensor<double, dim>& arg, const std::array<size_t, rank>& shape)
 {
     std::array<size_t, dim + rank> ret_shape;
     std::copy(arg.shape().begin(), arg.shape().end(), ret_shape.begin());
     std::copy(shape.begin(), shape.end(), ret_shape.begin() + dim);
     xt::xtensor<double, dim + rank> ret = xt::empty<double>(ret_shape);
     GooseFEM::asTensor<dim, rank>(arg, ret);
     return ret;
 }
 
 template <size_t dim, size_t rank>
 inline xt::xtensor<double, dim + rank> AsTensor(const xt::xtensor<double, dim>& arg, size_t n)
 {
     std::array<size_t, dim + rank> ret_shape;
     std::copy(arg.shape().begin(), arg.shape().end(), ret_shape.begin());
     std::fill(ret_shape.begin() + dim, ret_shape.end(), n);
     xt::xtensor<double, dim + rank> ret = xt::empty<double>(ret_shape);
     GooseFEM::asTensor<dim, rank>(arg, ret);
     return ret;
 }
 
 template <class T>
 inline xt::xarray<double> AsTensor(size_t rank, const T& arg, const std::vector<size_t>& shape)
 {
     GOOSEFEM_ASSERT(rank == shape.size());
     size_t dim = arg.dimension();
     std::vector<size_t> ret_shape(dim + rank);
     xt::dynamic_shape<ptrdiff_t> ret_strides(dim + rank);
     std::copy(arg.shape().begin(), arg.shape().end(), ret_shape.begin());
     std::copy(arg.strides().begin(), arg.strides().end(), ret_strides.begin());
     std::copy(shape.begin(), shape.end(), ret_shape.begin() + dim);
     std::fill(ret_strides.begin() + dim, ret_strides.end(), 0);
     xt::xarray<double> ret = xt::empty<double>(ret_shape);
     ret = xt::strided_view(arg, ret.shape(), std::move(ret_strides), 0ul, xt::layout_type::dynamic);
     return ret;
 }
 
 template <class T>
 inline xt::xarray<double> AsTensor(size_t rank, const T& arg, size_t n)
 {
     size_t dim = arg.dimension();
     using strides_type = typename T::strides_type::value_type;
     std::vector<size_t> ret_shape(dim + rank);
     xt::svector<strides_type> ret_strides(dim + rank);
     std::copy(arg.shape().begin(), arg.shape().end(), ret_shape.begin());
     std::copy(arg.strides().begin(), arg.strides().end(), ret_strides.begin());
     std::fill(ret_shape.begin() + dim, ret_shape.end(), n);
     std::fill(ret_strides.begin() + dim, ret_strides.end(), 0);
     xt::xarray<double> ret = xt::empty<double>(ret_shape);
     ret = xt::strided_view(arg, ret.shape(), std::move(ret_strides), 0ul, xt::layout_type::dynamic);
     return ret;
 }
 
+inline xt::xtensor<double, 2> as3d(const xt::xtensor<double, 2>& data)
+{
+    GOOSEFEM_ASSERT(data.shape(1) > 0 && data.shape(1) < 4)
+
+    if (data.shape(1) == 3ul) {
+        return data;
+    }
+
+    xt::xtensor<double, 2> ret = xt::zeros<double>(std::array<size_t, 2>{data.shape(0), 3ul});
+
+    if (data.shape(1) == 2ul) {
+        xt::view(ret, xt::all(), xt::keep(0, 1)) = data;
+    }
+
+    if (data.shape(1) == 1ul) {
+        xt::view(ret, xt::all(), xt::keep(0)) = data;
+    }
+
+    return ret;
+}
+
 } // namespace GooseFEM
 
 #endif
diff --git a/include/GooseFEM/ParaView.h b/include/GooseFEM/ParaView.h
deleted file mode 100644
index e8f6646..0000000
--- a/include/GooseFEM/ParaView.h
+++ /dev/null
@@ -1,266 +0,0 @@
-/*
-
-(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
-
-*/
-
-#ifndef GOOSEFEM_PARAVIEW_H
-#define GOOSEFEM_PARAVIEW_H
-
-// See: http://xdmf.org/index.php/XDMF_Model_and_Format
-
-#include "config.h"
-
-#include <fstream>
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-#include <highfive/H5Easy.hpp>
-#endif
-
-namespace GooseFEM {
-namespace ParaView {
-namespace HDF5 {
-
-inline std::string join(const std::vector<std::string>& lines, const std::string& sep = "\n");
-inline std::string indent(size_t n);
-xt::xtensor<double, 2> as3d(const xt::xtensor<double, 2>& data);
-
-enum class ElementType {
-    Triangle,
-    Quadrilateral,
-    Hexahedron };
-
-enum class AttributeType {
-    Cell,
-    Node };
-
-inline ElementType convert(GooseFEM::Mesh::ElementType type);
-inline std::string to_string(ElementType type);
-inline std::string to_string(AttributeType type);
-
-class Connectivity {
-public:
-    Connectivity() = default;
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-    Connectivity(
-        const H5Easy::File& data, const std::string& dataset, GooseFEM::Mesh::ElementType type);
-#endif
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-    Connectivity(
-        const std::string& fname, const std::string& dataset, GooseFEM::Mesh::ElementType type);
-#endif
-
-    Connectivity(
-        const std::string& fname,
-        const std::string& dataset,
-        GooseFEM::Mesh::ElementType type,
-        const std::vector<size_t>& shape);
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-    Connectivity(const H5Easy::File& data, const std::string& dataset, ElementType type);
-#endif
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-    Connectivity(const std::string& fname, const std::string& dataset, ElementType type);
-#endif
-
-    Connectivity(
-        const std::string& fname,
-        const std::string& dataset,
-        ElementType type,
-        const std::vector<size_t>& shape);
-
-    size_t nelem() const;
-    size_t nne() const;
-    std::vector<size_t> shape() const;
-    std::string fname() const;
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-    void checkShape();
-#endif
-
-    std::vector<std::string> xdmf(size_t indent = 4) const;
-
-private:
-#ifndef GOOSEFEM_NO_HIGHFIVE
-    void readShape(const H5Easy::File& data);
-#endif
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-    void init(const H5Easy::File& data, const std::string& dataset, ElementType type);
-#endif
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-    void init(const std::string& fname, const std::string& dataset, ElementType type);
-#endif
-
-    void init(
-        const std::string& fname,
-        const std::string& dataset,
-        ElementType type,
-        const std::vector<size_t>& shape);
-
-    ElementType m_type;
-    std::string m_fname;
-    std::string m_dataset;
-    std::vector<size_t> m_shape;
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-    bool m_verified = false; // if true: shape read from file, not need to check
-#endif
-};
-
-class Coordinates {
-public:
-    Coordinates() = default;
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-    Coordinates(const H5Easy::File& data, const std::string& dataset);
-#endif
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-    Coordinates(const std::string& fname, const std::string& dataset);
-#endif
-
-    Coordinates(
-        const std::string& fname, const std::string& dataset, const std::vector<size_t>& shape);
-
-    size_t nnode() const;
-    size_t ndim() const;
-    std::vector<size_t> shape() const;
-    std::string fname() const;
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-    void checkShape();
-#endif
-
-    std::vector<std::string> xdmf(size_t indent = 4) const;
-
-private:
-#ifndef GOOSEFEM_NO_HIGHFIVE
-    void readShape(const H5Easy::File& data);
-#endif
-
-    std::string m_fname;
-    std::string m_dataset;
-    std::vector<size_t> m_shape;
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-    bool m_verified = false; // if true: shape read from file, not need to check
-#endif
-};
-
-class Attribute {
-public:
-    Attribute() = default;
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-    Attribute(
-        const H5Easy::File& data,
-        const std::string& dataset,
-        const std::string& name,
-        AttributeType type);
-#endif
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-    Attribute(
-        const std::string& fname,
-        const std::string& dataset,
-        const std::string& name,
-        AttributeType type);
-#endif
-
-    Attribute(
-        const std::string& fname,
-        const std::string& dataset,
-        const std::string& name,
-        AttributeType type,
-        const std::vector<size_t>& shape);
-
-    std::vector<size_t> shape() const;
-    std::string fname() const;
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-    void checkShape();
-#endif
-
-    std::vector<std::string> xdmf(size_t indent = 4) const;
-
-private:
-#ifndef GOOSEFEM_NO_HIGHFIVE
-    void readShape(const H5Easy::File& data);
-#endif
-
-    AttributeType m_type;
-    std::string m_fname;
-    std::string m_dataset;
-    std::string m_name;
-    std::vector<size_t> m_shape;
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-    bool m_verified = false; // if true: shape read from file, not need to check
-#endif
-};
-
-class Mesh {
-public:
-    Mesh() = default;
-    Mesh(const Connectivity& conn, const Coordinates& coor);
-
-    void push_back(const Attribute& data);
-
-    std::vector<std::string> xdmf(size_t indent = 4) const;
-
-    void write(const std::string& fname, size_t n_indent = 4) const;
-
-private:
-    Connectivity m_conn;
-    Coordinates m_coor;
-    std::vector<Attribute> m_attr;
-};
-
-class Increment {
-public:
-    Increment() = default;
-
-    Increment(const Connectivity& conn, const Coordinates& coor);
-
-    Increment(
-        const Connectivity& conn, const Coordinates& coor, const std::vector<Attribute>& attr);
-
-    void push_back(const Connectivity& data);
-    void push_back(const Coordinates& data);
-    void push_back(const Attribute& data);
-
-    std::vector<std::string> xdmf(size_t indent = 4) const;
-
-private:
-    std::vector<Connectivity> m_conn;
-    std::vector<Coordinates> m_coor;
-    std::vector<Attribute> m_attr;
-};
-
-class TimeSeries {
-public:
-    TimeSeries() = default;
-    TimeSeries(const std::vector<Increment>& data);
-
-    void push_back(const Increment& data);
-
-    std::vector<std::string> xdmf(size_t indent = 4) const;
-
-    void write(const std::string& fname, size_t n_indent = 4) const;
-
-private:
-    std::vector<Increment> m_data;
-};
-
-} // namespace HDF5
-} // namespace ParaView
-} // namespace GooseFEM
-
-#include "ParaView.hpp"
-
-#endif
diff --git a/include/GooseFEM/ParaView.hpp b/include/GooseFEM/ParaView.hpp
deleted file mode 100644
index e201dcc..0000000
--- a/include/GooseFEM/ParaView.hpp
+++ /dev/null
@@ -1,671 +0,0 @@
-/*
-
-(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
-
-*/
-
-#ifndef GOOSEFEM_PARAVIEW_HPP
-#define GOOSEFEM_PARAVIEW_HPP
-
-#include "ParaView.h"
-
-namespace GooseFEM {
-namespace ParaView {
-namespace HDF5 {
-
-inline std::string join(const std::vector<std::string>& lines, const std::string& sep)
-{
-    if (lines.size() == 1) {
-        return lines[0];
-    }
-
-    std::string ret = "";
-
-    for (auto line : lines) {
-        if (ret.size() == 0) {
-            ret += line;
-            continue;
-        }
-
-        if (line[0] == sep[0]) {
-            ret += line;
-        }
-        else if (ret[ret.size() - 1] == sep[0]) {
-            ret += line;
-        }
-        else {
-            ret += sep + line;
-        }
-    }
-
-    return ret;
-}
-
-inline std::string indent(size_t n)
-{
-    std::string ret = "";
-
-    for (size_t i = 0; i < n; ++i) {
-        ret += " ";
-    }
-
-    return ret;
-}
-
-xt::xtensor<double, 2> as3d(const xt::xtensor<double, 2>& data)
-{
-    GOOSEFEM_ASSERT(data.shape(1) > 0 && data.shape(1) < 4)
-
-    if (data.shape(1) == 3ul) {
-        return data;
-    }
-
-    xt::xtensor<double, 2> ret = xt::zeros<double>(std::array<size_t, 2>{data.shape(0), 3ul});
-
-    if (data.shape(1) == 2ul) {
-        xt::view(ret, xt::all(), xt::keep(0, 1)) = data;
-    }
-
-    if (data.shape(1) == 1ul) {
-        xt::view(ret, xt::all(), xt::keep(0)) = data;
-    }
-
-    return ret;
-}
-
-inline ElementType convert(GooseFEM::Mesh::ElementType type)
-{
-    if (type == GooseFEM::Mesh::ElementType::Tri3) {
-        return ElementType::Triangle;
-    }
-
-    if (type == GooseFEM::Mesh::ElementType::Quad4) {
-        return ElementType::Quadrilateral;
-    }
-
-    if (type == GooseFEM::Mesh::ElementType::Hex8) {
-        return ElementType::Hexahedron;
-    }
-
-    throw std::runtime_error("Unknown GooseFEM::Mesh::ElementType");
-}
-
-inline std::string to_string(ElementType type)
-{
-    if (type == ElementType::Triangle) {
-        return "Triangle";
-    }
-
-    if (type == ElementType::Quadrilateral) {
-        return "Quadrilateral";
-    }
-
-    if (type == ElementType::Hexahedron) {
-        return "Hexahedron";
-    }
-
-    throw std::runtime_error("Unknown GooseFEM::ParaView::HDF5::ElementType");
-}
-
-inline std::string to_string(AttributeType type)
-{
-    if (type == AttributeType::Cell) {
-        return "Cell";
-    }
-
-    if (type == AttributeType::Node) {
-        return "Node";
-    }
-
-    throw std::runtime_error("Unknown GooseFEM::ParaView::HDF5::AttributeType");
-}
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-inline Connectivity::Connectivity(
-    const H5Easy::File& data, const std::string& dataset, GooseFEM::Mesh::ElementType type)
-{
-    init(data, dataset, convert(type));
-}
-#endif
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-inline Connectivity::Connectivity(
-    const std::string& fname, const std::string& dataset, GooseFEM::Mesh::ElementType type)
-{
-    init(fname, dataset, convert(type));
-}
-#endif
-
-inline Connectivity::Connectivity(
-    const std::string& fname,
-    const std::string& dataset,
-    GooseFEM::Mesh::ElementType type,
-    const std::vector<size_t>& shape)
-{
-    init(fname, dataset, convert(type), shape);
-}
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-inline Connectivity::Connectivity(
-    const H5Easy::File& data, const std::string& dataset, ElementType type)
-{
-    init(data, dataset, type);
-}
-#endif
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-inline Connectivity::Connectivity(
-    const std::string& fname, const std::string& dataset, ElementType type)
-{
-    init(fname, dataset, type);
-}
-#endif
-
-inline Connectivity::Connectivity(
-    const std::string& fname,
-    const std::string& dataset,
-    ElementType type,
-    const std::vector<size_t>& shape)
-{
-    init(fname, dataset, type, shape);
-}
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-inline void
-Connectivity::init(const H5Easy::File& data, const std::string& dataset, ElementType type)
-{
-    m_type = type;
-    m_dataset = dataset;
-
-    this->readShape(data);
-}
-#endif
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-inline void
-Connectivity::init(const std::string& fname, const std::string& dataset, ElementType type)
-{
-    m_type = type;
-    m_fname = fname;
-    m_dataset = dataset;
-
-    H5Easy::File data(m_fname, H5Easy::File::ReadOnly);
-
-    this->readShape(data);
-}
-#endif
-
-inline void Connectivity::init(
-    const std::string& fname,
-    const std::string& dataset,
-    ElementType type,
-    const std::vector<size_t>& shape)
-{
-    m_type = type;
-    m_fname = fname;
-    m_dataset = dataset;
-    m_shape = shape;
-
-    GOOSEFEM_ASSERT(m_shape.size() == 2);
-
-#ifdef GOOSEFEM_ENABLE_ASSERT
-    if (m_type == ElementType::Triangle) {
-        GOOSEFEM_ASSERT(m_shape[1] == 3);
-    }
-    else if (m_type == ElementType::Quadrilateral) {
-        GOOSEFEM_ASSERT(m_shape[1] == 4);
-    }
-    else if (m_type == ElementType::Hexahedron) {
-        GOOSEFEM_ASSERT(m_shape[1] == 8);
-    }
-#endif
-}
-
-inline size_t Connectivity::nelem() const
-{
-    return m_shape[0];
-}
-
-inline size_t Connectivity::nne() const
-{
-    return m_shape[1];
-}
-
-inline std::vector<size_t> Connectivity::shape() const
-{
-    return m_shape;
-}
-
-inline std::string Connectivity::fname() const
-{
-    return m_fname;
-}
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-inline void Connectivity::checkShape()
-{
-    if (m_verified)
-        return;
-
-    H5Easy::File data(m_fname, H5Easy::File::ReadOnly);
-
-    GOOSEFEM_CHECK(m_shape == H5Easy::getShape(data, m_dataset));
-
-    m_verified = true;
-}
-#endif
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-inline void Connectivity::readShape(const H5Easy::File& data)
-{
-    if (m_fname.size() == 0)
-        m_fname = data.getName();
-
-    m_shape = H5Easy::getShape(data, m_dataset);
-
-    GOOSEFEM_ASSERT(m_shape.size() == 2);
-
-#ifdef GOOSEFEM_ENABLE_ASSERT
-    if (m_type == ElementType::Triangle) {
-        GOOSEFEM_ASSERT(m_shape[1] == 3);
-    }
-    else if (m_type == ElementType::Quadrilateral) {
-        GOOSEFEM_ASSERT(m_shape[1] == 4);
-    }
-    else if (m_type == ElementType::Hexahedron) {
-        GOOSEFEM_ASSERT(m_shape[1] == 8);
-    }
-#endif
-
-    m_verified = true;
-}
-#endif
-
-inline std::vector<std::string> Connectivity::xdmf(size_t n_indent) const
-{
-    std::vector<std::string> ret;
-
-    ret.push_back(
-        "<Topology NumberOfElements=\"" + std::to_string(m_shape[0]) + "\" TopologyType=\"" +
-        to_string(m_type) + "\">");
-
-    ret.push_back(
-        indent(n_indent) + "<DataItem Dimensions=\"" + std::to_string(m_shape[0]) + " " +
-        std::to_string(m_shape[1]) + "\" Format=\"HDF\">" + m_fname + ":" + m_dataset +
-        "</DataItem>");
-
-    ret.push_back("</Topology>");
-
-    return ret;
-}
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-inline Coordinates::Coordinates(const H5Easy::File& data, const std::string& dataset)
-    : m_dataset(dataset)
-{
-    this->readShape(data);
-}
-#endif
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-inline Coordinates::Coordinates(const std::string& fname, const std::string& dataset)
-    : m_fname(fname), m_dataset(dataset)
-{
-    H5Easy::File data(m_fname, H5Easy::File::ReadOnly);
-
-    this->readShape(data);
-}
-#endif
-
-inline Coordinates::Coordinates(
-    const std::string& fname, const std::string& dataset, const std::vector<size_t>& shape)
-    : m_fname(fname), m_dataset(dataset), m_shape(shape)
-{
-    GOOSEFEM_ASSERT(m_shape.size() == 2);
-}
-
-inline size_t Coordinates::nnode() const
-{
-    return m_shape[0];
-}
-
-inline size_t Coordinates::ndim() const
-{
-    return m_shape[1];
-}
-
-inline std::vector<size_t> Coordinates::shape() const
-{
-    return m_shape;
-}
-
-inline std::string Coordinates::fname() const
-{
-    return m_fname;
-}
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-inline void Coordinates::checkShape()
-{
-    if (m_verified)
-        return;
-
-    H5Easy::File data(m_fname, H5Easy::File::ReadOnly);
-
-    GOOSEFEM_CHECK(m_shape == H5Easy::getShape(data, m_dataset));
-
-    m_verified = true;
-}
-#endif
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-inline void Coordinates::readShape(const H5Easy::File& data)
-{
-    if (m_fname.size() == 0)
-        m_fname = data.getName();
-
-    m_shape = H5Easy::getShape(data, m_dataset);
-
-    GOOSEFEM_ASSERT(m_shape.size() == 2);
-
-    m_verified = true;
-}
-#endif
-
-inline std::vector<std::string> Coordinates::xdmf(size_t n_indent) const
-{
-    std::vector<std::string> ret;
-
-    if (m_shape[1] == 1) {
-        ret.push_back("<Geometry GeometryType=\"X\">");
-    }
-    else if (m_shape[1] == 2) {
-        ret.push_back("<Geometry GeometryType=\"XY\">");
-    }
-    else if (m_shape[1] == 3) {
-        ret.push_back("<Geometry GeometryType=\"XYZ\">");
-    }
-
-    ret.push_back(
-        indent(n_indent) + "<DataItem Dimensions=\"" + std::to_string(m_shape[0]) + " " +
-        std::to_string(m_shape[1]) + "\" Format=\"HDF\">" + m_fname + ":" + m_dataset +
-        "</DataItem>");
-
-    ret.push_back("</Geometry>)");
-
-    return ret;
-}
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-inline Attribute::Attribute(
-    const H5Easy::File& data,
-    const std::string& dataset,
-    const std::string& name,
-    AttributeType type)
-    : m_type(type), m_dataset(dataset), m_name(name)
-{
-    this->readShape(data);
-}
-#endif
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-inline Attribute::Attribute(
-    const std::string& fname,
-    const std::string& dataset,
-    const std::string& name,
-    AttributeType type)
-    : m_type(type), m_fname(fname), m_dataset(dataset), m_name(name)
-{
-    H5Easy::File data(m_fname, H5Easy::File::ReadOnly);
-
-    this->readShape(data);
-}
-#endif
-
-inline Attribute::Attribute(
-    const std::string& fname,
-    const std::string& dataset,
-    const std::string& name,
-    AttributeType type,
-    const std::vector<size_t>& shape)
-    : m_type(type), m_fname(fname), m_dataset(dataset), m_name(name), m_shape(shape)
-{
-    GOOSEFEM_ASSERT(m_shape.size() > 0);
-}
-
-inline std::vector<size_t> Attribute::shape() const
-{
-    return m_shape;
-}
-
-inline std::string Attribute::fname() const
-{
-    return m_fname;
-}
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-inline void Attribute::checkShape()
-{
-    if (m_verified)
-        return;
-
-    H5Easy::File data(m_fname, H5Easy::File::ReadOnly);
-
-    GOOSEFEM_CHECK(m_shape == H5Easy::getShape(data, m_dataset));
-
-    m_verified = true;
-}
-#endif
-
-#ifndef GOOSEFEM_NO_HIGHFIVE
-inline void Attribute::readShape(const H5Easy::File& data)
-{
-    if (m_fname.size() == 0) {
-        m_fname = data.getName();
-    }
-
-    m_shape = H5Easy::getShape(data, m_dataset);
-
-    GOOSEFEM_ASSERT(m_shape.size() > 0);
-
-    m_verified = true;
-}
-#endif
-
-inline std::vector<std::string> Attribute::xdmf(size_t n_indent) const
-{
-    GOOSEFEM_ASSERT(m_shape.size() > 0);
-    GOOSEFEM_ASSERT(m_shape.size() < 3);
-
-    std::vector<std::string> ret;
-
-    if (m_shape.size() == 1) {
-        ret.push_back(
-            "<Attribute AttributeType=\"Scalar\" Center=\"" + to_string(m_type) + "\" Name=\"" +
-            m_name + "\">");
-
-        ret.push_back(
-            indent(n_indent) + "<DataItem Dimensions=\"" + std::to_string(m_shape[0]) +
-            "\" Format=\"HDF\">" + m_fname + ":" + m_dataset + "</DataItem>");
-    }
-    else if (m_shape.size() == 2) {
-        ret.push_back(
-            "<Attribute AttributeType=\"Vector\" Center=\"" + to_string(m_type) + "\" Name=\"" +
-            m_name + "\">");
-
-        ret.push_back(
-            indent(n_indent) + "<DataItem Dimensions=\"" + std::to_string(m_shape[0]) + " " +
-            std::to_string(m_shape[1]) + "\" Format=\"HDF\">" + m_fname + ":" + m_dataset +
-            "</DataItem>");
-    }
-
-    ret.push_back("</Attribute>)");
-
-    return ret;
-}
-
-inline Mesh::Mesh(const Connectivity& conn, const Coordinates& coor) : m_conn(conn), m_coor(coor)
-{
-}
-
-inline void Mesh::push_back(const Attribute& data)
-{
-    m_attr.push_back(data);
-}
-
-inline std::vector<std::string> Mesh::xdmf(size_t n_indent) const
-{
-    std::vector<std::string> ret;
-
-    {
-        std::vector<std::string> lines = m_conn.xdmf(n_indent);
-
-        for (auto& line : lines) {
-            ret.push_back(line);
-        }
-    }
-
-    {
-        std::vector<std::string> lines = m_coor.xdmf(n_indent);
-
-        for (auto& line : lines) {
-            ret.push_back(line);
-        }
-    }
-
-    for (auto& i : m_attr) {
-        std::vector<std::string> lines = i.xdmf(n_indent);
-
-        for (auto& line : lines) {
-            ret.push_back(line);
-        }
-    }
-
-    return ret;
-}
-
-inline void Mesh::write(const std::string& fname, size_t n_indent) const
-{
-    TimeSeries({Increment(m_conn, m_coor, m_attr)}).write(fname, n_indent);
-}
-
-inline Increment::Increment(const Connectivity& conn, const Coordinates& coor)
-{
-    m_conn.push_back(conn);
-    m_coor.push_back(coor);
-}
-
-inline Increment::Increment(
-    const Connectivity& conn, const Coordinates& coor, const std::vector<Attribute>& attr)
-    : m_attr(attr)
-{
-    m_conn.push_back(conn);
-    m_coor.push_back(coor);
-}
-
-inline void Increment::push_back(const Connectivity& data)
-{
-    m_conn.push_back(data);
-}
-
-inline void Increment::push_back(const Coordinates& data)
-{
-    m_coor.push_back(data);
-}
-
-inline void Increment::push_back(const Attribute& data)
-{
-    m_attr.push_back(data);
-}
-
-inline std::vector<std::string> Increment::xdmf(size_t n_indent) const
-{
-    std::vector<std::string> ret;
-
-    for (auto& i : m_conn) {
-        std::vector<std::string> lines = i.xdmf(n_indent);
-
-        for (auto& line : lines) {
-            ret.push_back(line);
-        }
-    }
-
-    for (auto& i : m_coor) {
-        std::vector<std::string> lines = i.xdmf(n_indent);
-
-        for (auto& line : lines) {
-            ret.push_back(line);
-        }
-    }
-
-    for (auto& i : m_attr) {
-        std::vector<std::string> lines = i.xdmf(n_indent);
-
-        for (auto& line : lines) {
-            ret.push_back(line);
-        }
-    }
-
-    return ret;
-}
-
-inline TimeSeries::TimeSeries(const std::vector<Increment>& data) : m_data(data)
-{
-}
-
-inline void TimeSeries::push_back(const Increment& data)
-{
-    m_data.push_back(data);
-}
-
-inline std::vector<std::string> TimeSeries::xdmf(size_t n_indent) const
-{
-    std::vector<std::string> ret;
-
-    for (size_t inc = 0; inc < m_data.size(); ++inc) {
-        ret.push_back("<Grid Name=\"Increment " + std::to_string(inc) + "\">");
-
-        ret.push_back(indent(n_indent) + "<Time Value=\"" + std::to_string(inc) + "\"/>");
-
-        std::vector<std::string> lines = m_data[inc].xdmf(n_indent);
-
-        for (auto& line : lines) {
-            ret.push_back(indent(n_indent) + line);
-        }
-
-        ret.push_back("</Grid>)");
-    }
-
-    return ret;
-}
-
-inline void TimeSeries::write(const std::string& fname, size_t n_indent) const
-{
-    std::ofstream myfile;
-
-    myfile.open(fname);
-
-    myfile << "<Xdmf Version=\"2.0\">" << std::endl;
-    myfile << indent(n_indent) + "<Domain>" << std::endl;
-    myfile << indent(n_indent * 2) +
-                  "<Grid CollectionType=\"Temporal\" GridType=\"Collection\" Name=\"TimeSeries\">"
-           << std::endl;
-
-    std::vector<std::string> lines = this->xdmf(n_indent);
-
-    for (auto& line : lines) {
-        myfile << indent(n_indent * 3) + line << std::endl;
-    }
-
-    myfile << indent(n_indent * 2) + "</Grid>" << std::endl;
-    myfile << indent(n_indent) + "</Domain>" << std::endl;
-    myfile << "</Xdmf>" << std::endl;
-
-    myfile.close();
-}
-
-} // namespace HDF5
-} // namespace ParaView
-} // namespace GooseFEM
-
-#endif
diff --git a/python/ParaView.hpp b/python/ParaView.hpp
deleted file mode 100644
index f3c948c..0000000
--- a/python/ParaView.hpp
+++ /dev/null
@@ -1,191 +0,0 @@
-/* =================================================================================================
-
-(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
-
-================================================================================================= */
-
-#include <pybind11/pybind11.h>
-#include <pyxtensor/pyxtensor.hpp>
-#define GOOSEFEM_NO_HIGHFIVE
-#include "../include/GooseFEM/ParaView.h"
-#include <GooseFEM/GooseFEM.h>
-
-namespace py = pybind11;
-
-void init_ParaViewHDF5(py::module& m)
-{
-
-    m.def("as3d", &GooseFEM::ParaView::HDF5::as3d);
-
-    py::enum_<GooseFEM::ParaView::HDF5::ElementType>(m, "ElementType", "ElementType")
-        .value("Triangle", GooseFEM::ParaView::HDF5::ElementType::Triangle)
-        .value("Quadrilateral", GooseFEM::ParaView::HDF5::ElementType::Quadrilateral)
-        .value("Hexahedron", GooseFEM::ParaView::HDF5::ElementType::Hexahedron)
-        .export_values();
-
-    py::enum_<GooseFEM::ParaView::HDF5::AttributeType>(m, "AttributeType", "AttributeType")
-        .value("Cell", GooseFEM::ParaView::HDF5::AttributeType::Cell)
-        .value("Node", GooseFEM::ParaView::HDF5::AttributeType::Node)
-        .export_values();
-
-    py::class_<GooseFEM::ParaView::HDF5::Connectivity>(m, "Connectivity")
-
-        .def(
-            py::init<
-                const std::string&,
-                const std::string&,
-                GooseFEM::ParaView::HDF5::ElementType,
-                const std::vector<size_t>&>(),
-            "Connectivity",
-            py::arg("fname"),
-            py::arg("dataset"),
-            py::arg("ElementType"),
-            py::arg("shape"))
-
-        .def(
-            py::init<
-                const std::string&,
-                const std::string&,
-                GooseFEM::Mesh::ElementType,
-                const std::vector<size_t>&>(),
-            "Connectivity",
-            py::arg("fname"),
-            py::arg("dataset"),
-            py::arg("ElementType"),
-            py::arg("shape"))
-
-        .def("nelem", &GooseFEM::ParaView::HDF5::Connectivity::nelem)
-
-        .def("nne", &GooseFEM::ParaView::HDF5::Connectivity::nne)
-
-        .def("shape", &GooseFEM::ParaView::HDF5::Connectivity::shape)
-
-        .def("fname", &GooseFEM::ParaView::HDF5::Connectivity::fname)
-
-        .def("xdmf", &GooseFEM::ParaView::HDF5::Connectivity::xdmf, py::arg("indent") = 4)
-
-        .def("__repr__", [](const GooseFEM::ParaView::HDF5::Connectivity&) {
-            return "<GooseFEM.ParaView.HDF5.Connectivity>";
-        });
-
-    py::class_<GooseFEM::ParaView::HDF5::Coordinates>(m, "Coordinates")
-
-        .def(
-            py::init<const std::string&, const std::string&, const std::vector<size_t>&>(),
-            "Coordinates",
-            py::arg("fname"),
-            py::arg("dataset"),
-            py::arg("shape"))
-
-        .def("nnode", &GooseFEM::ParaView::HDF5::Coordinates::nnode)
-        .def("ndim", &GooseFEM::ParaView::HDF5::Coordinates::ndim)
-        .def("shape", &GooseFEM::ParaView::HDF5::Coordinates::shape)
-        .def("fname", &GooseFEM::ParaView::HDF5::Coordinates::fname)
-        .def("xdmf", &GooseFEM::ParaView::HDF5::Coordinates::xdmf, py::arg("indent") = 4)
-
-        .def("__repr__", [](const GooseFEM::ParaView::HDF5::Coordinates&) {
-            return "<GooseFEM.ParaView.HDF5.Coordinates>";
-        });
-
-    py::class_<GooseFEM::ParaView::HDF5::Attribute>(m, "Attribute")
-
-        .def(
-            py::init<
-                const std::string&,
-                const std::string&,
-                const std::string&,
-                GooseFEM::ParaView::HDF5::AttributeType,
-                const std::vector<size_t>&>(),
-            "Attribute",
-            py::arg("fname"),
-            py::arg("dataset"),
-            py::arg("name"),
-            py::arg("AttributeType"),
-            py::arg("shape"))
-
-        .def("shape", &GooseFEM::ParaView::HDF5::Attribute::shape)
-
-        .def("fname", &GooseFEM::ParaView::HDF5::Attribute::fname)
-
-        .def("xdmf", &GooseFEM::ParaView::HDF5::Attribute::xdmf, py::arg("indent") = 4)
-
-        .def("__repr__", [](const GooseFEM::ParaView::HDF5::Attribute&) {
-            return "<GooseFEM.ParaView.HDF5.Attribute>";
-        });
-
-    py::class_<GooseFEM::ParaView::HDF5::Mesh>(m, "Mesh")
-
-        .def(
-            py::init<
-                const GooseFEM::ParaView::HDF5::Connectivity&,
-                const GooseFEM::ParaView::HDF5::Coordinates&>(),
-            "Mesh",
-            py::arg("conn"),
-            py::arg("coor"))
-
-        .def(
-            "push_back",
-            py::overload_cast<const GooseFEM::ParaView::HDF5::Attribute&>(
-                &GooseFEM::ParaView::HDF5::Mesh::push_back))
-
-        .def("xdmf", &GooseFEM::ParaView::HDF5::Mesh::xdmf, py::arg("indent") = 4)
-
-        .def(
-            "write",
-            &GooseFEM::ParaView::HDF5::Mesh::write,
-            py::arg("fname"),
-            py::arg("indent") = 4)
-
-        .def("__repr__", [](const GooseFEM::ParaView::HDF5::Mesh&) {
-            return "<GooseFEM.ParaView.HDF5.Mesh>";
-        });
-
-    py::class_<GooseFEM::ParaView::HDF5::Increment>(m, "Increment")
-
-        .def(
-            py::init<
-                const GooseFEM::ParaView::HDF5::Connectivity&,
-                const GooseFEM::ParaView::HDF5::Coordinates&>(),
-            "Increment",
-            py::arg("conn"),
-            py::arg("coor"))
-
-        .def(
-            "push_back",
-            py::overload_cast<const GooseFEM::ParaView::HDF5::Connectivity&>(
-                &GooseFEM::ParaView::HDF5::Increment::push_back))
-
-        .def(
-            "push_back",
-            py::overload_cast<const GooseFEM::ParaView::HDF5::Coordinates&>(
-                &GooseFEM::ParaView::HDF5::Increment::push_back))
-
-        .def(
-            "push_back",
-            py::overload_cast<const GooseFEM::ParaView::HDF5::Attribute&>(
-                &GooseFEM::ParaView::HDF5::Increment::push_back))
-
-        .def("xdmf", &GooseFEM::ParaView::HDF5::Increment::xdmf, py::arg("indent") = 4)
-
-        .def("__repr__", [](const GooseFEM::ParaView::HDF5::Increment&) {
-            return "<GooseFEM.ParaView.HDF5.Increment>";
-        });
-
-    py::class_<GooseFEM::ParaView::HDF5::TimeSeries>(m, "TimeSeries")
-
-        .def(py::init<>(), "TimeSeries")
-
-        .def("push_back", &GooseFEM::ParaView::HDF5::TimeSeries::push_back)
-
-        .def("xdmf", &GooseFEM::ParaView::HDF5::TimeSeries::xdmf, py::arg("indent") = 4)
-
-        .def(
-            "write",
-            &GooseFEM::ParaView::HDF5::TimeSeries::write,
-            py::arg("fname"),
-            py::arg("indent") = 4)
-
-        .def("__repr__", [](const GooseFEM::ParaView::HDF5::TimeSeries&) {
-            return "<GooseFEM.ParaView.HDF5.TimeSeries>";
-        });
-}
diff --git a/python/main.cpp b/python/main.cpp
index 831b116..3a64638 100644
--- a/python/main.cpp
+++ b/python/main.cpp
@@ -1,136 +1,125 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
 #include <Eigen/Eigen>
 
 #include <pybind11/pybind11.h>
 #include <pybind11/eigen.h>
 
 #include <pyxtensor/pyxtensor.hpp>
 
 #define GOOSEFEM_ENABLE_ASSERT
 #include <GooseFEM/GooseFEM.h>
 
 namespace py = pybind11;
 
 #include "Vector.hpp"
 #include "VectorPartitioned.hpp"
 #include "VectorPartitionedTyings.hpp"
 #include "Matrix.hpp"
 #include "MatrixPartitioned.hpp"
 #include "MatrixPartitionedTyings.hpp"
 #include "MatrixDiagonal.hpp"
 #include "MatrixDiagonalPartitioned.hpp"
 #include "Element.hpp"
 #include "ElementQuad4.hpp"
 #include "ElementQuad4Planar.hpp"
 #include "ElementQuad4Axisymmetric.hpp"
 #include "ElementHex8.hpp"
 #include "Mesh.hpp"
 #include "MeshTri3.hpp"
 #include "MeshQuad4.hpp"
 #include "MeshHex8.hpp"
-#include "ParaView.hpp"
 
 PYBIND11_MODULE(GooseFEM, m) {
 
 // --------
 // GooseFEM
 // --------
 
 m.doc() = "Some simple finite element meshes and operations";
 
 init_Vector(m);
 init_VectorPartitioned(m);
 init_VectorPartitionedTyings(m);
 init_Matrix(m);
 init_MatrixPartitioned(m);
 init_MatrixPartitionedTyings(m);
 init_MatrixDiagonal(m);
 init_MatrixDiagonalPartitioned(m);
 
 // ----------------
 // GooseFEM.Element
 // ----------------
 
 py::module mElement = m.def_submodule("Element", "Generic element routines");
 
 init_Element(mElement);
 
 // ----------------------
 // GooseFEM.Element.Quad4
 // ----------------------
 
 py::module mElementQuad4 = mElement.def_submodule("Quad4", "Linear quadrilateral elements (2D)");
 py::module mElementQuad4Gauss = mElementQuad4.def_submodule("Gauss", "Gauss quadrature");
 py::module mElementQuad4Nodal = mElementQuad4.def_submodule("Nodal", "Nodal quadrature");
 
 init_ElementQuad4(mElementQuad4);
 init_ElementQuad4Planar(mElementQuad4);
 init_ElementQuad4Axisymmetric(mElementQuad4);
 init_ElementQuad4Gauss(mElementQuad4Gauss);
 init_ElementQuad4Nodal(mElementQuad4Nodal);
 
 // ---------------------
 // GooseFEM.Element.Hex8
 // ---------------------
 
 py::module mElementHex8 = mElement.def_submodule("Hex8", "Linear hexahedron (brick) elements (3D)");
 py::module mElementHex8Gauss = mElementHex8.def_submodule("Gauss", "Gauss quadrature");
 py::module mElementHex8Nodal = mElementHex8.def_submodule("Nodal", "Nodal quadrature");
 
 init_ElementHex8(mElementHex8);
 init_ElementHex8Gauss(mElementHex8Gauss);
 init_ElementHex8Nodal(mElementHex8Nodal);
 
 // -------------
 // GooseFEM.Mesh
 // -------------
 
 py::module mMesh = m.def_submodule("Mesh", "Generic mesh routines");
 
 init_Mesh(mMesh);
 
 // ------------------
 // GooseFEM.Mesh.Tri3
 // ------------------
 
 py::module mMeshTri3 = mMesh.def_submodule("Tri3", "Linear triangular elements (2D)");
 
 init_MeshTri3(mMeshTri3);
 
 // -------------------
 // GooseFEM.Mesh.Quad4
 // -------------------
 
 py::module mMeshQuad4 = mMesh.def_submodule("Quad4", "Linear quadrilateral elements (2D)");
 
 init_MeshQuad4(mMeshQuad4);
 
 py::module mMeshQuad4Map = mMeshQuad4.def_submodule("Map", "Map mesh objects");
 
 init_MeshQuad4Map(mMeshQuad4Map);
 
 // ------------------
 // GooseFEM.Mesh.Hex8
 // ------------------
 
 py::module mMeshHex8 = mMesh.def_submodule("Hex8", "Linear hexahedron (brick) elements (3D)");
 
 init_MeshHex8(mMeshHex8);
 
-// -----------------
-// GooseFEM.ParaView
-// -----------------
-
-py::module mParaView = m.def_submodule("ParaView", "ParaView output files");
-
-py::module mParaViewHDF5 = mParaView.def_submodule("HDF5", "ParaView/HDF5 support using XDMF files");
-
-init_ParaViewHDF5(mParaViewHDF5);
-
 }