diff --git a/python/Element.hpp b/python/Element.hpp
index bd9ba77..c22c2e6 100644
--- a/python/Element.hpp
+++ b/python/Element.hpp
@@ -1,35 +1,29 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
+#include <GooseFEM/GooseFEM.h>
 #include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
-#include "../include/GooseFEM/GooseFEM.h"
-
-// =================================================================================================
 
 namespace py = pybind11;
 
-// =================================================================================================
-
-void init_Element(py::module &m)
+void init_Element(py::module& m)
 {
 
-m.def("asElementVector",
-    &GooseFEM::Element::asElementVector,
-    "Convert nodal vector [nnode, ndim] to nodal vector stored per element [nelem, nne, ndim]",
-    py::arg("conn"),
-    py::arg("nodevec"));
-
-m.def("assembleElementVector",
-    &GooseFEM::Element::assembleNodeVector,
-    "Assemble nodal vector stored per element [nelem, nne, ndim] to nodal vector [nnode, ndim]",
-    py::arg("conn"),
-    py::arg("elemvec"));
-
+    m.def(
+        "asElementVector",
+        &GooseFEM::Element::asElementVector,
+        "Convert nodal vector [nnode, ndim] to nodal vector stored per element [nelem, nne, ndim]",
+        py::arg("conn"),
+        py::arg("nodevec"));
+
+    m.def(
+        "assembleElementVector",
+        &GooseFEM::Element::assembleNodeVector,
+        "Assemble nodal vector stored per element [nelem, nne, ndim] to nodal vector [nnode, ndim]",
+        py::arg("conn"),
+        py::arg("elemvec"));
 }
-
-// =================================================================================================
-
diff --git a/python/ElementHex8.hpp b/python/ElementHex8.hpp
index 4e0755a..8c64f97 100644
--- a/python/ElementHex8.hpp
+++ b/python/ElementHex8.hpp
@@ -1,144 +1,119 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
+#include <GooseFEM/GooseFEM.h>
 #include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
-#include "../include/GooseFEM/GooseFEM.h"
-
-// =================================================================================================
 
 namespace py = pybind11;
 
-// =================================================================================================
-
-void init_ElementHex8(py::module &m)
+void init_ElementHex8(py::module& m)
 {
 
-py::class_<GooseFEM::Element::Hex8::Quadrature>(m, "Quadrature")
-
-    .def(py::init<const xt::xtensor<double,3>&>(),
-        "Quadrature",
-        py::arg("x"))
-
-    .def(py::init<
-            const xt::xtensor<double,3>&,
-            const xt::xtensor<double,2>&,
-            const xt::xtensor<double,1>&>(),
-        "Quadrature",
-        py::arg("x"),
-        py::arg("xi"),
-        py::arg("w"))
-
-    .def("update_x",
-        &GooseFEM::Element::Hex8::Quadrature::update_x,
-        "Update the nodal positions")
-
-    .def("nelem",
-         &GooseFEM::Element::Hex8::Quadrature::nelem,
-         "Number of elements")
-
-    .def("nne",
-         &GooseFEM::Element::Hex8::Quadrature::nne,
-         "Number of nodes per element")
-
-    .def("ndim",
-         &GooseFEM::Element::Hex8::Quadrature::ndim,
-         "Number of dimensions")
-
-    .def("nip",
-         &GooseFEM::Element::Hex8::Quadrature::nip,
-         "Number of integration points")
-
-    .def("DV",
-        py::overload_cast<size_t>(&GooseFEM::Element::Hex8::Quadrature::DV, py::const_),
-        "Integration point volume (qtensor)")
-
-    .def("DV",
-        py::overload_cast<>(&GooseFEM::Element::Hex8::Quadrature::DV, py::const_),
-        "Integration point volume (qscalar)")
-
-    .def("GradN_vector",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::Element::Hex8::Quadrature::GradN_vector, py::const_),
-        "Dyadic product, returns 'qtensor'",
-        py::arg("elemvec"))
-
-    .def("GradN_vector_T",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::Element::Hex8::Quadrature::GradN_vector_T, py::const_),
-        "Dyadic product, returns 'qtensor'",
-        py::arg("elemvec"))
-
-    .def("SymGradN_vector",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::Element::Hex8::Quadrature::SymGradN_vector, py::const_),
-        "Dyadic product, returns 'qtensor'",
-        py::arg("elemvec"))
-
-    .def("Int_N_scalar_NT_dV",
-        py::overload_cast<const xt::xtensor<double,2>&>(
-            &GooseFEM::Element::Hex8::Quadrature::Int_N_scalar_NT_dV, py::const_),
-        "Integration, returns 'elemmat'",
-        py::arg("qscalar"))
-
-    .def("Int_gradN_dot_tensor2_dV",
-        py::overload_cast<const xt::xtensor<double,4>&>(
-            &GooseFEM::Element::Hex8::Quadrature::Int_gradN_dot_tensor2_dV, py::const_),
-        "Integration, returns 'elemvec'",
-        py::arg("qtensor"))
-
-    .def("Int_gradN_dot_tensor4_dot_gradNT_dV",
-        py::overload_cast<const xt::xtensor<double,6>&>(
-            &GooseFEM::Element::Hex8::Quadrature::Int_gradN_dot_tensor4_dot_gradNT_dV, py::const_),
-        "Integration, returns 'elemvec'",
-        py::arg("qtensor"))
-
-    .def("__repr__",
-        [](const GooseFEM::Element::Hex8::Quadrature&){
-            return "<GooseFEM.Element.Hex8.Quadrature>"; });
-
+    py::class_<GooseFEM::Element::Hex8::Quadrature>(m, "Quadrature")
+
+        .def(py::init<const xt::xtensor<double, 3>&>(), "Quadrature", py::arg("x"))
+
+        .def(
+            py::init<
+                const xt::xtensor<double, 3>&,
+                const xt::xtensor<double, 2>&,
+                const xt::xtensor<double, 1>&>(),
+            "Quadrature",
+            py::arg("x"),
+            py::arg("xi"),
+            py::arg("w"))
+
+        .def(
+            "update_x",
+            &GooseFEM::Element::Hex8::Quadrature::update_x,
+            "Update the nodal positions")
+
+        .def("nelem", &GooseFEM::Element::Hex8::Quadrature::nelem, "Number of elements")
+
+        .def("nne", &GooseFEM::Element::Hex8::Quadrature::nne, "Number of nodes per element")
+
+        .def("ndim", &GooseFEM::Element::Hex8::Quadrature::ndim, "Number of dimensions")
+
+        .def("nip", &GooseFEM::Element::Hex8::Quadrature::nip, "Number of integration points")
+
+        .def(
+            "DV",
+            py::overload_cast<size_t>(&GooseFEM::Element::Hex8::Quadrature::DV, py::const_),
+            "Integration point volume (qtensor)")
+
+        .def(
+            "DV",
+            py::overload_cast<>(&GooseFEM::Element::Hex8::Quadrature::DV, py::const_),
+            "Integration point volume (qscalar)")
+
+        .def(
+            "GradN_vector",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::Element::Hex8::Quadrature::GradN_vector, py::const_),
+            "Dyadic product, returns 'qtensor'",
+            py::arg("elemvec"))
+
+        .def(
+            "GradN_vector_T",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::Element::Hex8::Quadrature::GradN_vector_T, py::const_),
+            "Dyadic product, returns 'qtensor'",
+            py::arg("elemvec"))
+
+        .def(
+            "SymGradN_vector",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::Element::Hex8::Quadrature::SymGradN_vector, py::const_),
+            "Dyadic product, returns 'qtensor'",
+            py::arg("elemvec"))
+
+        .def(
+            "Int_N_scalar_NT_dV",
+            py::overload_cast<const xt::xtensor<double, 2>&>(
+                &GooseFEM::Element::Hex8::Quadrature::Int_N_scalar_NT_dV, py::const_),
+            "Integration, returns 'elemmat'",
+            py::arg("qscalar"))
+
+        .def(
+            "Int_gradN_dot_tensor2_dV",
+            py::overload_cast<const xt::xtensor<double, 4>&>(
+                &GooseFEM::Element::Hex8::Quadrature::Int_gradN_dot_tensor2_dV, py::const_),
+            "Integration, returns 'elemvec'",
+            py::arg("qtensor"))
+
+        .def(
+            "Int_gradN_dot_tensor4_dot_gradNT_dV",
+            py::overload_cast<const xt::xtensor<double, 6>&>(
+                &GooseFEM::Element::Hex8::Quadrature::Int_gradN_dot_tensor4_dot_gradNT_dV,
+                py::const_),
+            "Integration, returns 'elemvec'",
+            py::arg("qtensor"))
+
+        .def("__repr__", [](const GooseFEM::Element::Hex8::Quadrature&) {
+            return "<GooseFEM.Element.Hex8.Quadrature>";
+        });
 }
 
-// -------------------------------------------------------------------------------------------------
-
 void init_ElementHex8Gauss(py::module& m)
 {
 
-m.def("nip",
-    &GooseFEM::Element::Hex8::Gauss::nip,
-    "Return number of integration point");
-
-m.def("xi",
-    &GooseFEM::Element::Hex8::Gauss::xi,
-    "Return integration point coordinates");
+    m.def("nip", &GooseFEM::Element::Hex8::Gauss::nip, "Return number of integration point");
 
-m.def("w",
-    &GooseFEM::Element::Hex8::Gauss::w,
-    "Return integration point weights");
+    m.def("xi", &GooseFEM::Element::Hex8::Gauss::xi, "Return integration point coordinates");
 
+    m.def("w", &GooseFEM::Element::Hex8::Gauss::w, "Return integration point weights");
 }
 
-// -------------------------------------------------------------------------------------------------
-
 void init_ElementHex8Nodal(py::module& m)
 {
 
-m.def("nip",
-    &GooseFEM::Element::Hex8::Nodal::nip,
-    "Return number of integration point");
-
-m.def("xi",
-    &GooseFEM::Element::Hex8::Nodal::xi,
-    "Return integration point coordinates");
+    m.def("nip", &GooseFEM::Element::Hex8::Nodal::nip, "Return number of integration point");
 
-m.def("w",
-    &GooseFEM::Element::Hex8::Nodal::w,
-    "Return integration point weights");
+    m.def("xi", &GooseFEM::Element::Hex8::Nodal::xi, "Return integration point coordinates");
 
+    m.def("w", &GooseFEM::Element::Hex8::Nodal::w, "Return integration point weights");
 }
-
-// =================================================================================================
-
diff --git a/python/ElementQuad4.hpp b/python/ElementQuad4.hpp
index 7b8c9ba..9dfe9ac 100644
--- a/python/ElementQuad4.hpp
+++ b/python/ElementQuad4.hpp
@@ -1,144 +1,119 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
+#include <GooseFEM/GooseFEM.h>
 #include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
-#include "../include/GooseFEM/GooseFEM.h"
-
-// =================================================================================================
 
 namespace py = pybind11;
 
-// =================================================================================================
-
 void init_ElementQuad4(py::module& m)
 {
 
-py::class_<GooseFEM::Element::Quad4::Quadrature>(m, "Quadrature")
-
-    .def(py::init<const xt::xtensor<double,3>&>(),
-        "Quadrature",
-        py::arg("x"))
-
-    .def(py::init<
-            const xt::xtensor<double,3>&,
-            const xt::xtensor<double,2>&,
-            const xt::xtensor<double,1>&>(),
-        "Quadrature",
-        py::arg("x"),
-        py::arg("xi"),
-        py::arg("w"))
-
-    .def("update_x",
-        &GooseFEM::Element::Quad4::Quadrature::update_x,
-        "Update the nodal positions")
-
-    .def("nelem",
-        &GooseFEM::Element::Quad4::Quadrature::nelem,
-        "Number of elements")
-
-    .def("nne",
-        &GooseFEM::Element::Quad4::Quadrature::nne,
-        "Number of nodes per element")
-
-    .def("ndim",
-        &GooseFEM::Element::Quad4::Quadrature::ndim,
-        "Number of dimensions")
-
-    .def("nip",
-        &GooseFEM::Element::Quad4::Quadrature::nip,
-        "Number of integration points")
-
-    .def("DV",
-        py::overload_cast<size_t>(&GooseFEM::Element::Quad4::Quadrature::DV, py::const_),
-        "Integration point volume (qtensor)")
-
-    .def("DV",
-        py::overload_cast<>(&GooseFEM::Element::Quad4::Quadrature::DV, py::const_),
-        "Integration point volume (qscalar)")
-
-    .def("GradN_vector",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::Element::Quad4::Quadrature::GradN_vector, py::const_),
-        "Dyadic product, returns 'qtensor'",
-        py::arg("elemvec"))
-
-    .def("GradN_vector_T",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::Element::Quad4::Quadrature::GradN_vector_T, py::const_),
-        "Dyadic product, returns 'qtensor'",
-        py::arg("elemvec"))
-
-    .def("SymGradN_vector",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::Element::Quad4::Quadrature::SymGradN_vector, py::const_),
-        "Dyadic product, returns 'qtensor'",
-        py::arg("elemvec"))
-
-    .def("Int_N_scalar_NT_dV",
-        py::overload_cast<const xt::xtensor<double,2>&>(
-            &GooseFEM::Element::Quad4::Quadrature::Int_N_scalar_NT_dV, py::const_),
-        "Integration, returns 'elemmat'",
-        py::arg("qscalar"))
-
-    .def("Int_gradN_dot_tensor2_dV",
-        py::overload_cast<const xt::xtensor<double,4>&>(
-            &GooseFEM::Element::Quad4::Quadrature::Int_gradN_dot_tensor2_dV, py::const_),
-        "Integration, returns 'elemvec'",
-        py::arg("qtensor"))
-
-    .def("Int_gradN_dot_tensor4_dot_gradNT_dV",
-        py::overload_cast<const xt::xtensor<double,6>&>(
-            &GooseFEM::Element::Quad4::Quadrature::Int_gradN_dot_tensor4_dot_gradNT_dV, py::const_),
-        "Integration, returns 'elemvec'",
-        py::arg("qtensor"))
-
-    .def("__repr__",
-        [](const GooseFEM::Element::Quad4::Quadrature&){
-            return "<GooseFEM.Element.Quad4.Quadrature>"; });
-
+    py::class_<GooseFEM::Element::Quad4::Quadrature>(m, "Quadrature")
+
+        .def(py::init<const xt::xtensor<double, 3>&>(), "Quadrature", py::arg("x"))
+
+        .def(
+            py::init<
+                const xt::xtensor<double, 3>&,
+                const xt::xtensor<double, 2>&,
+                const xt::xtensor<double, 1>&>(),
+            "Quadrature",
+            py::arg("x"),
+            py::arg("xi"),
+            py::arg("w"))
+
+        .def(
+            "update_x",
+            &GooseFEM::Element::Quad4::Quadrature::update_x,
+            "Update the nodal positions")
+
+        .def("nelem", &GooseFEM::Element::Quad4::Quadrature::nelem, "Number of elements")
+
+        .def("nne", &GooseFEM::Element::Quad4::Quadrature::nne, "Number of nodes per element")
+
+        .def("ndim", &GooseFEM::Element::Quad4::Quadrature::ndim, "Number of dimensions")
+
+        .def("nip", &GooseFEM::Element::Quad4::Quadrature::nip, "Number of integration points")
+
+        .def(
+            "DV",
+            py::overload_cast<size_t>(&GooseFEM::Element::Quad4::Quadrature::DV, py::const_),
+            "Integration point volume (qtensor)")
+
+        .def(
+            "DV",
+            py::overload_cast<>(&GooseFEM::Element::Quad4::Quadrature::DV, py::const_),
+            "Integration point volume (qscalar)")
+
+        .def(
+            "GradN_vector",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::Element::Quad4::Quadrature::GradN_vector, py::const_),
+            "Dyadic product, returns 'qtensor'",
+            py::arg("elemvec"))
+
+        .def(
+            "GradN_vector_T",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::Element::Quad4::Quadrature::GradN_vector_T, py::const_),
+            "Dyadic product, returns 'qtensor'",
+            py::arg("elemvec"))
+
+        .def(
+            "SymGradN_vector",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::Element::Quad4::Quadrature::SymGradN_vector, py::const_),
+            "Dyadic product, returns 'qtensor'",
+            py::arg("elemvec"))
+
+        .def(
+            "Int_N_scalar_NT_dV",
+            py::overload_cast<const xt::xtensor<double, 2>&>(
+                &GooseFEM::Element::Quad4::Quadrature::Int_N_scalar_NT_dV, py::const_),
+            "Integration, returns 'elemmat'",
+            py::arg("qscalar"))
+
+        .def(
+            "Int_gradN_dot_tensor2_dV",
+            py::overload_cast<const xt::xtensor<double, 4>&>(
+                &GooseFEM::Element::Quad4::Quadrature::Int_gradN_dot_tensor2_dV, py::const_),
+            "Integration, returns 'elemvec'",
+            py::arg("qtensor"))
+
+        .def(
+            "Int_gradN_dot_tensor4_dot_gradNT_dV",
+            py::overload_cast<const xt::xtensor<double, 6>&>(
+                &GooseFEM::Element::Quad4::Quadrature::Int_gradN_dot_tensor4_dot_gradNT_dV,
+                py::const_),
+            "Integration, returns 'elemvec'",
+            py::arg("qtensor"))
+
+        .def("__repr__", [](const GooseFEM::Element::Quad4::Quadrature&) {
+            return "<GooseFEM.Element.Quad4.Quadrature>";
+        });
 }
 
-// -------------------------------------------------------------------------------------------------
-
 void init_ElementQuad4Gauss(py::module& m)
 {
 
-m.def("nip",
-    &GooseFEM::Element::Quad4::Gauss::nip,
-    "Return number of integration point");
-
-m.def("xi",
-    &GooseFEM::Element::Quad4::Gauss::xi,
-    "Return integration point coordinates");
+    m.def("nip", &GooseFEM::Element::Quad4::Gauss::nip, "Return number of integration point");
 
-m.def("w",
-    &GooseFEM::Element::Quad4::Gauss::w,
-    "Return integration point weights");
+    m.def("xi", &GooseFEM::Element::Quad4::Gauss::xi, "Return integration point coordinates");
 
+    m.def("w", &GooseFEM::Element::Quad4::Gauss::w, "Return integration point weights");
 }
 
-// -------------------------------------------------------------------------------------------------
-
 void init_ElementQuad4Nodal(py::module& m)
 {
 
-m.def("nip",
-    &GooseFEM::Element::Quad4::Nodal::nip,
-    "Return number of integration point");
-
-m.def("xi",
-    &GooseFEM::Element::Quad4::Nodal::xi,
-    "Return integration point coordinates");
+    m.def("nip", &GooseFEM::Element::Quad4::Nodal::nip, "Return number of integration point");
 
-m.def("w",
-    &GooseFEM::Element::Quad4::Nodal::w,
-    "Return integration point weights");
+    m.def("xi", &GooseFEM::Element::Quad4::Nodal::xi, "Return integration point coordinates");
 
+    m.def("w", &GooseFEM::Element::Quad4::Nodal::w, "Return integration point weights");
 }
-
-// =================================================================================================
-
diff --git a/python/ElementQuad4Axisymmetric.hpp b/python/ElementQuad4Axisymmetric.hpp
index 97e54fe..7faee06 100644
--- a/python/ElementQuad4Axisymmetric.hpp
+++ b/python/ElementQuad4Axisymmetric.hpp
@@ -1,106 +1,110 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
+#include <GooseFEM/GooseFEM.h>
 #include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
-#include "../include/GooseFEM/GooseFEM.h"
-
-// =================================================================================================
 
 namespace py = pybind11;
 
-// =================================================================================================
-
 void init_ElementQuad4Axisymmetric(py::module& m)
 {
 
-py::class_<GooseFEM::Element::Quad4::QuadratureAxisymmetric>(m, "QuadratureAxisymmetric")
-
-    .def(py::init<const xt::xtensor<double,3>&>(),
-        "QuadratureAxisymmetric",
-        py::arg("x"))
-
-    .def(py::init<
-            const xt::xtensor<double,3>&,
-            const xt::xtensor<double,2>&,
-            const xt::xtensor<double,1>&>(),
-        "QuadratureAxisymmetric",
-        py::arg("x"),
-        py::arg("xi"),
-        py::arg("w"))
-
-    .def("update_x",
-        &GooseFEM::Element::Quad4::QuadratureAxisymmetric::update_x,
-        "Update the nodal positions")
-
-    .def("nelem",
-        &GooseFEM::Element::Quad4::QuadratureAxisymmetric::nelem,
-        "Number of elements")
-
-    .def("nne",
-        &GooseFEM::Element::Quad4::QuadratureAxisymmetric::nne,
-        "Number of nodes per element")
-
-    .def("ndim",
-        &GooseFEM::Element::Quad4::QuadratureAxisymmetric::ndim,
-        "Number of dimensions")
-
-    .def("nip",
-        &GooseFEM::Element::Quad4::QuadratureAxisymmetric::nip,
-        "Number of integration points")
-
-    .def("DV",
-        py::overload_cast<size_t>(&GooseFEM::Element::Quad4::QuadratureAxisymmetric::DV, py::const_),
-        "Integration point volume (qtensor)")
-
-    .def("DV",
-        py::overload_cast<>(&GooseFEM::Element::Quad4::QuadratureAxisymmetric::DV, py::const_),
-        "Integration point volume (qscalar)")
-
-    .def("GradN_vector",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::Element::Quad4::QuadratureAxisymmetric::GradN_vector, py::const_),
-        "Dyadic product, returns 'qtensor'",
-        py::arg("elemvec"))
-
-    .def("GradN_vector_T",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::Element::Quad4::QuadratureAxisymmetric::GradN_vector_T, py::const_),
-        "Dyadic product, returns 'qtensor'",
-        py::arg("elemvec"))
-
-    .def("SymGradN_vector",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::Element::Quad4::QuadratureAxisymmetric::SymGradN_vector, py::const_),
-        "Dyadic product, returns 'qtensor'",
-        py::arg("elemvec"))
-
-    .def("Int_N_scalar_NT_dV",
-        py::overload_cast<const xt::xtensor<double,2>&>(
-            &GooseFEM::Element::Quad4::QuadratureAxisymmetric::Int_N_scalar_NT_dV, py::const_),
-        "Integration, returns 'elemmat'",
-        py::arg("qscalar"))
-
-    .def("Int_gradN_dot_tensor2_dV",
-        py::overload_cast<const xt::xtensor<double,4>&>(
-            &GooseFEM::Element::Quad4::QuadratureAxisymmetric::Int_gradN_dot_tensor2_dV, py::const_),
-        "Integration, returns 'elemvec'",
-        py::arg("qtensor"))
-
-    .def("Int_gradN_dot_tensor4_dot_gradNT_dV",
-        py::overload_cast<const xt::xtensor<double,6>&>(
-            &GooseFEM::Element::Quad4::QuadratureAxisymmetric::Int_gradN_dot_tensor4_dot_gradNT_dV, py::const_),
-        "Integration, returns 'elemvec'",
-        py::arg("qtensor"))
-
-    .def("__repr__",
-        [](const GooseFEM::Element::Quad4::QuadratureAxisymmetric&){
-            return "<GooseFEM.Element.Quad4.QuadratureAxisymmetric>"; });
-
+    py::class_<GooseFEM::Element::Quad4::QuadratureAxisymmetric>(m, "QuadratureAxisymmetric")
+
+        .def(py::init<const xt::xtensor<double, 3>&>(), "QuadratureAxisymmetric", py::arg("x"))
+
+        .def(
+            py::init<
+                const xt::xtensor<double, 3>&,
+                const xt::xtensor<double, 2>&,
+                const xt::xtensor<double, 1>&>(),
+            "QuadratureAxisymmetric",
+            py::arg("x"),
+            py::arg("xi"),
+            py::arg("w"))
+
+        .def(
+            "update_x",
+            &GooseFEM::Element::Quad4::QuadratureAxisymmetric::update_x,
+            "Update the nodal positions")
+
+        .def(
+            "nelem", &GooseFEM::Element::Quad4::QuadratureAxisymmetric::nelem, "Number of elements")
+
+        .def(
+            "nne",
+            &GooseFEM::Element::Quad4::QuadratureAxisymmetric::nne,
+            "Number of nodes per element")
+
+        .def(
+            "ndim", &GooseFEM::Element::Quad4::QuadratureAxisymmetric::ndim, "Number of dimensions")
+
+        .def(
+            "nip",
+            &GooseFEM::Element::Quad4::QuadratureAxisymmetric::nip,
+            "Number of integration points")
+
+        .def(
+            "DV",
+            py::overload_cast<size_t>(
+                &GooseFEM::Element::Quad4::QuadratureAxisymmetric::DV, py::const_),
+            "Integration point volume (qtensor)")
+
+        .def(
+            "DV",
+            py::overload_cast<>(&GooseFEM::Element::Quad4::QuadratureAxisymmetric::DV, py::const_),
+            "Integration point volume (qscalar)")
+
+        .def(
+            "GradN_vector",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::Element::Quad4::QuadratureAxisymmetric::GradN_vector, py::const_),
+            "Dyadic product, returns 'qtensor'",
+            py::arg("elemvec"))
+
+        .def(
+            "GradN_vector_T",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::Element::Quad4::QuadratureAxisymmetric::GradN_vector_T, py::const_),
+            "Dyadic product, returns 'qtensor'",
+            py::arg("elemvec"))
+
+        .def(
+            "SymGradN_vector",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::Element::Quad4::QuadratureAxisymmetric::SymGradN_vector, py::const_),
+            "Dyadic product, returns 'qtensor'",
+            py::arg("elemvec"))
+
+        .def(
+            "Int_N_scalar_NT_dV",
+            py::overload_cast<const xt::xtensor<double, 2>&>(
+                &GooseFEM::Element::Quad4::QuadratureAxisymmetric::Int_N_scalar_NT_dV, py::const_),
+            "Integration, returns 'elemmat'",
+            py::arg("qscalar"))
+
+        .def(
+            "Int_gradN_dot_tensor2_dV",
+            py::overload_cast<const xt::xtensor<double, 4>&>(
+                &GooseFEM::Element::Quad4::QuadratureAxisymmetric::Int_gradN_dot_tensor2_dV,
+                py::const_),
+            "Integration, returns 'elemvec'",
+            py::arg("qtensor"))
+
+        .def(
+            "Int_gradN_dot_tensor4_dot_gradNT_dV",
+            py::overload_cast<const xt::xtensor<double, 6>&>(
+                &GooseFEM::Element::Quad4::QuadratureAxisymmetric::
+                    Int_gradN_dot_tensor4_dot_gradNT_dV,
+                py::const_),
+            "Integration, returns 'elemvec'",
+            py::arg("qtensor"))
+
+        .def("__repr__", [](const GooseFEM::Element::Quad4::QuadratureAxisymmetric&) {
+            return "<GooseFEM.Element.Quad4.QuadratureAxisymmetric>";
+        });
 }
-
-// =================================================================================================
-
diff --git a/python/ElementQuad4Planar.hpp b/python/ElementQuad4Planar.hpp
index 05d6b0e..61a8a89 100644
--- a/python/ElementQuad4Planar.hpp
+++ b/python/ElementQuad4Planar.hpp
@@ -1,106 +1,100 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
+#include <GooseFEM/GooseFEM.h>
 #include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
-#include "../include/GooseFEM/GooseFEM.h"
-
-// =================================================================================================
 
 namespace py = pybind11;
 
-// =================================================================================================
-
 void init_ElementQuad4Planar(py::module& m)
 {
 
-py::class_<GooseFEM::Element::Quad4::QuadraturePlanar>(m, "QuadraturePlanar")
-
-    .def(py::init<const xt::xtensor<double,3>&>(),
-        "QuadraturePlanar",
-        py::arg("x"))
-
-    .def(py::init<
-            const xt::xtensor<double,3>&,
-            const xt::xtensor<double,2>&,
-            const xt::xtensor<double,1>&>(),
-        "QuadraturePlanar",
-        py::arg("x"),
-        py::arg("xi"),
-        py::arg("w"))
-
-    .def("update_x",
-        &GooseFEM::Element::Quad4::QuadraturePlanar::update_x,
-        "Update the nodal positions")
-
-    .def("nelem",
-        &GooseFEM::Element::Quad4::QuadraturePlanar::nelem,
-        "Number of elements")
-
-    .def("nne",
-        &GooseFEM::Element::Quad4::QuadraturePlanar::nne,
-        "Number of nodes per element")
-
-    .def("ndim",
-        &GooseFEM::Element::Quad4::QuadraturePlanar::ndim,
-        "Number of dimensions")
-
-    .def("nip",
-        &GooseFEM::Element::Quad4::QuadraturePlanar::nip,
-        "Number of integration points")
-
-    .def("DV",
-        py::overload_cast<size_t>(&GooseFEM::Element::Quad4::QuadraturePlanar::DV, py::const_),
-        "Integration point volume (qtensor)")
-
-    .def("DV",
-        py::overload_cast<>(&GooseFEM::Element::Quad4::QuadraturePlanar::DV, py::const_),
-        "Integration point volume (qscalar)")
-
-    .def("GradN_vector",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::Element::Quad4::QuadraturePlanar::GradN_vector, py::const_),
-        "Dyadic product, returns 'qtensor'",
-        py::arg("elemvec"))
-
-    .def("GradN_vector_T",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::Element::Quad4::QuadraturePlanar::GradN_vector_T, py::const_),
-        "Dyadic product, returns 'qtensor'",
-        py::arg("elemvec"))
-
-    .def("SymGradN_vector",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::Element::Quad4::QuadraturePlanar::SymGradN_vector, py::const_),
-        "Dyadic product, returns 'qtensor'",
-        py::arg("elemvec"))
-
-    .def("Int_N_scalar_NT_dV",
-        py::overload_cast<const xt::xtensor<double,2>&>(
-            &GooseFEM::Element::Quad4::QuadraturePlanar::Int_N_scalar_NT_dV, py::const_),
-        "Integration, returns 'elemmat'",
-        py::arg("qscalar"))
-
-    .def("Int_gradN_dot_tensor2_dV",
-        py::overload_cast<const xt::xtensor<double,4>&>(
-            &GooseFEM::Element::Quad4::QuadraturePlanar::Int_gradN_dot_tensor2_dV, py::const_),
-        "Integration, returns 'elemvec'",
-        py::arg("qtensor"))
-
-    .def("Int_gradN_dot_tensor4_dot_gradNT_dV",
-        py::overload_cast<const xt::xtensor<double,6>&>(
-            &GooseFEM::Element::Quad4::QuadraturePlanar::Int_gradN_dot_tensor4_dot_gradNT_dV, py::const_),
-        "Integration, returns 'elemvec'",
-        py::arg("qtensor"))
-
-    .def("__repr__",
-        [](const GooseFEM::Element::Quad4::QuadraturePlanar&){
-            return "<GooseFEM.Element.Quad4.QuadraturePlanar>"; });
-
+    py::class_<GooseFEM::Element::Quad4::QuadraturePlanar>(m, "QuadraturePlanar")
+
+        .def(py::init<const xt::xtensor<double, 3>&>(), "QuadraturePlanar", py::arg("x"))
+
+        .def(
+            py::init<
+                const xt::xtensor<double, 3>&,
+                const xt::xtensor<double, 2>&,
+                const xt::xtensor<double, 1>&>(),
+            "QuadraturePlanar",
+            py::arg("x"),
+            py::arg("xi"),
+            py::arg("w"))
+
+        .def(
+            "update_x",
+            &GooseFEM::Element::Quad4::QuadraturePlanar::update_x,
+            "Update the nodal positions")
+
+        .def("nelem", &GooseFEM::Element::Quad4::QuadraturePlanar::nelem, "Number of elements")
+
+        .def("nne", &GooseFEM::Element::Quad4::QuadraturePlanar::nne, "Number of nodes per element")
+
+        .def("ndim", &GooseFEM::Element::Quad4::QuadraturePlanar::ndim, "Number of dimensions")
+
+        .def(
+            "nip", &GooseFEM::Element::Quad4::QuadraturePlanar::nip, "Number of integration points")
+
+        .def(
+            "DV",
+            py::overload_cast<size_t>(&GooseFEM::Element::Quad4::QuadraturePlanar::DV, py::const_),
+            "Integration point volume (qtensor)")
+
+        .def(
+            "DV",
+            py::overload_cast<>(&GooseFEM::Element::Quad4::QuadraturePlanar::DV, py::const_),
+            "Integration point volume (qscalar)")
+
+        .def(
+            "GradN_vector",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::Element::Quad4::QuadraturePlanar::GradN_vector, py::const_),
+            "Dyadic product, returns 'qtensor'",
+            py::arg("elemvec"))
+
+        .def(
+            "GradN_vector_T",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::Element::Quad4::QuadraturePlanar::GradN_vector_T, py::const_),
+            "Dyadic product, returns 'qtensor'",
+            py::arg("elemvec"))
+
+        .def(
+            "SymGradN_vector",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::Element::Quad4::QuadraturePlanar::SymGradN_vector, py::const_),
+            "Dyadic product, returns 'qtensor'",
+            py::arg("elemvec"))
+
+        .def(
+            "Int_N_scalar_NT_dV",
+            py::overload_cast<const xt::xtensor<double, 2>&>(
+                &GooseFEM::Element::Quad4::QuadraturePlanar::Int_N_scalar_NT_dV, py::const_),
+            "Integration, returns 'elemmat'",
+            py::arg("qscalar"))
+
+        .def(
+            "Int_gradN_dot_tensor2_dV",
+            py::overload_cast<const xt::xtensor<double, 4>&>(
+                &GooseFEM::Element::Quad4::QuadraturePlanar::Int_gradN_dot_tensor2_dV, py::const_),
+            "Integration, returns 'elemvec'",
+            py::arg("qtensor"))
+
+        .def(
+            "Int_gradN_dot_tensor4_dot_gradNT_dV",
+            py::overload_cast<const xt::xtensor<double, 6>&>(
+                &GooseFEM::Element::Quad4::QuadraturePlanar::Int_gradN_dot_tensor4_dot_gradNT_dV,
+                py::const_),
+            "Integration, returns 'elemvec'",
+            py::arg("qtensor"))
+
+        .def("__repr__", [](const GooseFEM::Element::Quad4::QuadraturePlanar&) {
+            return "<GooseFEM.Element.Quad4.QuadraturePlanar>";
+        });
 }
-
-// =================================================================================================
-
diff --git a/python/Matrix.hpp b/python/Matrix.hpp
index ad5256d..335d0ff 100644
--- a/python/Matrix.hpp
+++ b/python/Matrix.hpp
@@ -1,93 +1,79 @@
 /* =================================================================================================
 
 (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 <GooseFEM/GooseFEM.h>
 #include <pybind11/eigen.h>
+#include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
-#include "../include/GooseFEM/GooseFEM.h"
 
 namespace py = pybind11;
 
 void init_Matrix(py::module& m)
 {
 
-py::class_<GooseFEM::Matrix>(m, "Matrix")
-
-    .def(py::init<
-            const xt::xtensor<size_t,2>&,
-            const xt::xtensor<size_t,2>&>(),
-        "Sparse matrix",
-        py::arg("conn"),
-        py::arg("dofs"))
-
-    .def("nelem",
-        &GooseFEM::Matrix::nelem,
-        "Number of element")
-
-    .def("nne",
-        &GooseFEM::Matrix::nne,
-        "Number of nodes per element")
-
-    .def("nnode",
-        &GooseFEM::Matrix::nnode,
-        "Number of nodes")
-
-    .def("ndim",
-        &GooseFEM::Matrix::ndim,
-        "Number of dimensions")
-
-    .def("ndof",
-        &GooseFEM::Matrix::ndof,
-        "Number of degrees-of-freedom")
-
-    .def("assemble",
-        &GooseFEM::Matrix::assemble,
-        "Assemble matrix from 'elemmat",
-        py::arg("elemmat"))
-
-    .def("dofs",
-        &GooseFEM::Matrix::dofs,
-        "Return degrees-of-freedom")
-
-    .def("Dot",
-        py::overload_cast<const xt::xtensor<double,1>&>(&GooseFEM::Matrix::Dot, py::const_),
-        "Dot",
-        py::arg("x"))
-
-    .def("Dot",
-        py::overload_cast<const xt::xtensor<double,2>&>(&GooseFEM::Matrix::Dot, py::const_),
-        "Dot",
-        py::arg("x"))
-
-    .def("__repr__",
-        [](const GooseFEM::Matrix&){
-            return "<GooseFEM.Matrix>"; });
-
-py::class_<GooseFEM::MatrixSolver<>>(m, "MatrixSolver")
-
-    .def(py::init<>(),
-        "Sparse matrix solver")
-
-    .def("Solve",
-        py::overload_cast<GooseFEM::Matrix&, const xt::xtensor<double,1>&>(
-            &GooseFEM::MatrixSolver<>::Solve),
-        "Solve",
-        py::arg("matrix"),
-        py::arg("b"))
-
-    .def("Solve",
-        py::overload_cast<GooseFEM::Matrix&, const xt::xtensor<double,2>&>(
-            &GooseFEM::MatrixSolver<>::Solve),
-        "Solve",
-        py::arg("matrix"),
-        py::arg("b"))
-
-    .def("__repr__",
-        [](const GooseFEM::MatrixSolver<>&){
-            return "<GooseFEM.MatrixSolver>"; });
+    py::class_<GooseFEM::Matrix>(m, "Matrix")
+
+        .def(
+            py::init<const xt::xtensor<size_t, 2>&, const xt::xtensor<size_t, 2>&>(),
+            "Sparse matrix",
+            py::arg("conn"),
+            py::arg("dofs"))
+
+        .def("nelem", &GooseFEM::Matrix::nelem, "Number of element")
+
+        .def("nne", &GooseFEM::Matrix::nne, "Number of nodes per element")
+
+        .def("nnode", &GooseFEM::Matrix::nnode, "Number of nodes")
+
+        .def("ndim", &GooseFEM::Matrix::ndim, "Number of dimensions")
+
+        .def("ndof", &GooseFEM::Matrix::ndof, "Number of degrees-of-freedom")
+
+        .def(
+            "assemble",
+            &GooseFEM::Matrix::assemble,
+            "Assemble matrix from 'elemmat",
+            py::arg("elemmat"))
+
+        .def("dofs", &GooseFEM::Matrix::dofs, "Return degrees-of-freedom")
+
+        .def(
+            "Dot",
+            py::overload_cast<const xt::xtensor<double, 1>&>(&GooseFEM::Matrix::Dot, py::const_),
+            "Dot",
+            py::arg("x"))
+
+        .def(
+            "Dot",
+            py::overload_cast<const xt::xtensor<double, 2>&>(&GooseFEM::Matrix::Dot, py::const_),
+            "Dot",
+            py::arg("x"))
+
+        .def("__repr__", [](const GooseFEM::Matrix&) { return "<GooseFEM.Matrix>"; });
+
+    py::class_<GooseFEM::MatrixSolver<>>(m, "MatrixSolver")
+
+        .def(py::init<>(), "Sparse matrix solver")
+
+        .def(
+            "Solve",
+            py::overload_cast<GooseFEM::Matrix&, const xt::xtensor<double, 1>&>(
+                &GooseFEM::MatrixSolver<>::Solve),
+            "Solve",
+            py::arg("matrix"),
+            py::arg("b"))
+
+        .def(
+            "Solve",
+            py::overload_cast<GooseFEM::Matrix&, const xt::xtensor<double, 2>&>(
+                &GooseFEM::MatrixSolver<>::Solve),
+            "Solve",
+            py::arg("matrix"),
+            py::arg("b"))
 
+        .def("__repr__", [](const GooseFEM::MatrixSolver<>&) { return "<GooseFEM.MatrixSolver>"; });
 }
diff --git a/python/MatrixDiagonal.hpp b/python/MatrixDiagonal.hpp
index 1e818d3..fbdf0b4 100644
--- a/python/MatrixDiagonal.hpp
+++ b/python/MatrixDiagonal.hpp
@@ -1,97 +1,78 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
+#include <GooseFEM/GooseFEM.h>
 #include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
-#include "../include/GooseFEM/GooseFEM.h"
-
-// =================================================================================================
 
 namespace py = pybind11;
 
-// =================================================================================================
-
 void init_MatrixDiagonal(py::module& m)
 {
 
-py::class_<GooseFEM::MatrixDiagonal>(m, "MatrixDiagonal")
-
-    .def(py::init<
-            const xt::xtensor<size_t,2>&,
-            const xt::xtensor<size_t,2>&>(),
-        "Diagonal matrix",
-        py::arg("conn"),
-        py::arg("dofs"))
-
-    .def("nelem",
-        &GooseFEM::MatrixDiagonal::nelem,
-        "Number of element")
-
-    .def("nne",
-        &GooseFEM::MatrixDiagonal::nne,
-        "Number of nodes per element")
-
-    .def("nnode",
-        &GooseFEM::MatrixDiagonal::nnode,
-        "Number of nodes")
-
-    .def("ndim",
-        &GooseFEM::MatrixDiagonal::ndim,
-        "Number of dimensions")
-
-    .def("ndof",
-        &GooseFEM::MatrixDiagonal::ndof,
-        "Number of degrees-of-freedom")
-
-    .def("set",
-        &GooseFEM::MatrixDiagonal::set,
-        "Set matrix components",
-        py::arg("A"))
-
-    .def("assemble",
-        &GooseFEM::MatrixDiagonal::assemble,
-        "Assemble matrix from 'elemmat",
-        py::arg("elemmat"))
-
-    .def("dofs",
-        &GooseFEM::MatrixDiagonal::dofs,
-        "Return degrees-of-freedom")
-
-    .def("AsDiagonal", &GooseFEM::MatrixDiagonal::AsDiagonal,
-        "Return as diagonal matrix (column)")
-
-    .def("Dot",
-        py::overload_cast<const xt::xtensor<double,1>&>(
-            &GooseFEM::MatrixDiagonal::Dot, py::const_),
-        "Dot product 'b_i = A_ij * x_j",
-        py::arg("x"))
-
-    .def("Dot",
-        py::overload_cast<const xt::xtensor<double,2>&>(
-            &GooseFEM::MatrixDiagonal::Dot, py::const_),
-        "Dot product 'b_i = A_ij * x_j",
-        py::arg("x"))
-
-    .def("Solve",
-        py::overload_cast<const xt::xtensor<double,1>&>(
-            &GooseFEM::MatrixDiagonal::Solve),
-        "Solve",
-        py::arg("b"))
-
-    .def("Solve",
-        py::overload_cast<const xt::xtensor<double,2>&>(
-            &GooseFEM::MatrixDiagonal::Solve),
-        "Solve",
-        py::arg("b"))
-
-    .def("__repr__",
-        [](const GooseFEM::MatrixDiagonal&){
-            return "<GooseFEM.MatrixDiagonal>"; });
+    py::class_<GooseFEM::MatrixDiagonal>(m, "MatrixDiagonal")
 
-}
+        .def(
+            py::init<const xt::xtensor<size_t, 2>&, const xt::xtensor<size_t, 2>&>(),
+            "Diagonal matrix",
+            py::arg("conn"),
+            py::arg("dofs"))
+
+        .def("nelem", &GooseFEM::MatrixDiagonal::nelem, "Number of element")
+
+        .def("nne", &GooseFEM::MatrixDiagonal::nne, "Number of nodes per element")
+
+        .def("nnode", &GooseFEM::MatrixDiagonal::nnode, "Number of nodes")
+
+        .def("ndim", &GooseFEM::MatrixDiagonal::ndim, "Number of dimensions")
+
+        .def("ndof", &GooseFEM::MatrixDiagonal::ndof, "Number of degrees-of-freedom")
 
-// =================================================================================================
+        .def("set", &GooseFEM::MatrixDiagonal::set, "Set matrix components", py::arg("A"))
 
+        .def(
+            "assemble",
+            &GooseFEM::MatrixDiagonal::assemble,
+            "Assemble matrix from 'elemmat",
+            py::arg("elemmat"))
+
+        .def("dofs", &GooseFEM::MatrixDiagonal::dofs, "Return degrees-of-freedom")
+
+        .def(
+            "AsDiagonal",
+            &GooseFEM::MatrixDiagonal::AsDiagonal,
+            "Return as diagonal matrix (column)")
+
+        .def(
+            "Dot",
+            py::overload_cast<const xt::xtensor<double, 1>&>(
+                &GooseFEM::MatrixDiagonal::Dot, py::const_),
+            "Dot product 'b_i = A_ij * x_j",
+            py::arg("x"))
+
+        .def(
+            "Dot",
+            py::overload_cast<const xt::xtensor<double, 2>&>(
+                &GooseFEM::MatrixDiagonal::Dot, py::const_),
+            "Dot product 'b_i = A_ij * x_j",
+            py::arg("x"))
+
+        .def(
+            "Solve",
+            py::overload_cast<const xt::xtensor<double, 1>&>(&GooseFEM::MatrixDiagonal::Solve),
+            "Solve",
+            py::arg("b"))
+
+        .def(
+            "Solve",
+            py::overload_cast<const xt::xtensor<double, 2>&>(&GooseFEM::MatrixDiagonal::Solve),
+            "Solve",
+            py::arg("b"))
+
+        .def("__repr__", [](const GooseFEM::MatrixDiagonal&) {
+            return "<GooseFEM.MatrixDiagonal>";
+        });
+}
diff --git a/python/MatrixDiagonalPartitioned.hpp b/python/MatrixDiagonalPartitioned.hpp
index 1b5311c..78a40d9 100644
--- a/python/MatrixDiagonalPartitioned.hpp
+++ b/python/MatrixDiagonalPartitioned.hpp
@@ -1,148 +1,142 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
+#include <GooseFEM/GooseFEM.h>
 #include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
-#include "../include/GooseFEM/GooseFEM.h"
-
-// =================================================================================================
 
 namespace py = pybind11;
 
-// =================================================================================================
-
 void init_MatrixDiagonalPartitioned(py::module& m)
 {
 
-py::class_<GooseFEM::MatrixDiagonalPartitioned>(m, "MatrixDiagonalPartitioned")
-
-    .def(py::init<
-            const xt::xtensor<size_t,2>&,
-            const xt::xtensor<size_t,2>&,
-            const xt::xtensor<size_t,1>&>(),
-        "Diagonal, partitioned, matrix",
-        py::arg("conn"),
-        py::arg("dofs"),
-        py::arg("iip"))
-
-    .def("nelem",
-        &GooseFEM::MatrixDiagonalPartitioned::nelem,
-        "Number of element")
-
-    .def("nne",
-        &GooseFEM::MatrixDiagonalPartitioned::nne,
-        "Number of nodes per element")
-
-    .def("nnode",
-        &GooseFEM::MatrixDiagonalPartitioned::nnode,
-        "Number of nodes")
-
-    .def("ndim",
-        &GooseFEM::MatrixDiagonalPartitioned::ndim,
-        "Number of dimensions")
-
-    .def("ndof",
-        &GooseFEM::MatrixDiagonalPartitioned::ndof,
-        "Number of degrees-of-freedom")
-
-    .def("nnu",
-        &GooseFEM::MatrixDiagonalPartitioned::nnu,
-        "Number of unknown degrees-of-freedom")
-
-    .def("nnp",
-        &GooseFEM::MatrixDiagonalPartitioned::nnp,
-        "Number of prescribed degrees-of-freedom")
-
-    .def("assemble",
-        &GooseFEM::MatrixDiagonalPartitioned::assemble,
-        "Assemble matrix from 'elemmat",
-        py::arg("elemmat"))
-
-    .def("dofs",
-        &GooseFEM::MatrixDiagonalPartitioned::dofs,
-        "Return degrees-of-freedom")
-
-    .def("iiu",
-        &GooseFEM::MatrixDiagonalPartitioned::iiu,
-        "Return unknown degrees-of-freedom")
-
-    .def("iip",
-        &GooseFEM::MatrixDiagonalPartitioned::iip,
-        "Return prescribed degrees-of-freedom")
-
-    .def("AsDiagonal", &GooseFEM::MatrixDiagonalPartitioned::AsDiagonal,
-        "Return as diagonal matrix (column)")
-
-    .def("Dot",
-        py::overload_cast<const xt::xtensor<double,1>&>(
-            &GooseFEM::MatrixDiagonalPartitioned::Dot, py::const_),
-        "Dot product 'b_i = A_ij * x_j",
-        py::arg("x"))
-
-    .def("Dot_u",
-        py::overload_cast<const xt::xtensor<double,1>&, const xt::xtensor<double,1>&>(
-            &GooseFEM::MatrixDiagonalPartitioned::Dot_u, py::const_),
-        "Dot product 'b_i = A_ij * x_j (b_u = A_uu * x_u + A_up * x_p == A_uu * x_u)",
-        py::arg("x_u"),
-        py::arg("x_p"))
-
-    .def("Dot_p",
-        py::overload_cast<const xt::xtensor<double,1>&, const xt::xtensor<double,1>&>(
-            &GooseFEM::MatrixDiagonalPartitioned::Dot_p, py::const_),
-        "Dot product 'b_i = A_ij * x_j (b_p = A_pu * x_u + A_pp * x_p == A_pp * x_p)",
-        py::arg("x_u"),
-        py::arg("x_p"))
-
-    .def("Solve",
-        py::overload_cast<const xt::xtensor<double,1>&, const xt::xtensor<double,1>&>(
-            &GooseFEM::MatrixDiagonalPartitioned::Solve),
-        "Solve",
-        py::arg("b"),
-        py::arg("x"))
-
-    .def("Solve",
-        py::overload_cast<const xt::xtensor<double,2>&, const xt::xtensor<double,2>&>(
-            &GooseFEM::MatrixDiagonalPartitioned::Solve),
-        "Solve",
-        py::arg("b"),
-        py::arg("x"))
-
-    .def("Solve_u",
-        py::overload_cast<const xt::xtensor<double,1>&, const xt::xtensor<double,1>&>(
-            &GooseFEM::MatrixDiagonalPartitioned::Solve_u),
-        "Solve_u",
-        py::arg("b_u"),
-        py::arg("x_p"))
-
-    .def("Reaction",
-        py::overload_cast<const xt::xtensor<double,1>&, const xt::xtensor<double,1>&>(
-            &GooseFEM::MatrixDiagonalPartitioned::Reaction, py::const_),
-        "Reaction",
-        py::arg("x"),
-        py::arg("b"))
-
-    .def("Reaction",
-        py::overload_cast<const xt::xtensor<double,2>&, const xt::xtensor<double,2>&>(
-            &GooseFEM::MatrixDiagonalPartitioned::Reaction, py::const_),
-        "Reaction",
-        py::arg("x"),
-        py::arg("b"))
-
-    .def("Reaction_p",
-        py::overload_cast<const xt::xtensor<double,1>&, const xt::xtensor<double,1>&>(
-            &GooseFEM::MatrixDiagonalPartitioned::Reaction_p, py::const_),
-        "Reaction_p",
-        py::arg("x_u"),
-        py::arg("x_p"))
-
-    .def("__repr__",
-        [](const GooseFEM::MatrixDiagonalPartitioned&){
-            return "<GooseFEM.MatrixDiagonalPartitioned>"; });
-
+    py::class_<GooseFEM::MatrixDiagonalPartitioned>(m, "MatrixDiagonalPartitioned")
+
+        .def(
+            py::init<
+                const xt::xtensor<size_t, 2>&,
+                const xt::xtensor<size_t, 2>&,
+                const xt::xtensor<size_t, 1>&>(),
+            "Diagonal, partitioned, matrix",
+            py::arg("conn"),
+            py::arg("dofs"),
+            py::arg("iip"))
+
+        .def("nelem", &GooseFEM::MatrixDiagonalPartitioned::nelem, "Number of element")
+
+        .def("nne", &GooseFEM::MatrixDiagonalPartitioned::nne, "Number of nodes per element")
+
+        .def("nnode", &GooseFEM::MatrixDiagonalPartitioned::nnode, "Number of nodes")
+
+        .def("ndim", &GooseFEM::MatrixDiagonalPartitioned::ndim, "Number of dimensions")
+
+        .def("ndof", &GooseFEM::MatrixDiagonalPartitioned::ndof, "Number of degrees-of-freedom")
+
+        .def(
+            "nnu",
+            &GooseFEM::MatrixDiagonalPartitioned::nnu,
+            "Number of unknown degrees-of-freedom")
+
+        .def(
+            "nnp",
+            &GooseFEM::MatrixDiagonalPartitioned::nnp,
+            "Number of prescribed degrees-of-freedom")
+
+        .def(
+            "assemble",
+            &GooseFEM::MatrixDiagonalPartitioned::assemble,
+            "Assemble matrix from 'elemmat",
+            py::arg("elemmat"))
+
+        .def("dofs", &GooseFEM::MatrixDiagonalPartitioned::dofs, "Return degrees-of-freedom")
+
+        .def("iiu", &GooseFEM::MatrixDiagonalPartitioned::iiu, "Return unknown degrees-of-freedom")
+
+        .def(
+            "iip",
+            &GooseFEM::MatrixDiagonalPartitioned::iip,
+            "Return prescribed degrees-of-freedom")
+
+        .def(
+            "AsDiagonal",
+            &GooseFEM::MatrixDiagonalPartitioned::AsDiagonal,
+            "Return as diagonal matrix (column)")
+
+        .def(
+            "Dot",
+            py::overload_cast<const xt::xtensor<double, 1>&>(
+                &GooseFEM::MatrixDiagonalPartitioned::Dot, py::const_),
+            "Dot product 'b_i = A_ij * x_j",
+            py::arg("x"))
+
+        .def(
+            "Dot_u",
+            py::overload_cast<const xt::xtensor<double, 1>&, const xt::xtensor<double, 1>&>(
+                &GooseFEM::MatrixDiagonalPartitioned::Dot_u, py::const_),
+            "Dot product 'b_i = A_ij * x_j (b_u = A_uu * x_u + A_up * x_p == A_uu * x_u)",
+            py::arg("x_u"),
+            py::arg("x_p"))
+
+        .def(
+            "Dot_p",
+            py::overload_cast<const xt::xtensor<double, 1>&, const xt::xtensor<double, 1>&>(
+                &GooseFEM::MatrixDiagonalPartitioned::Dot_p, py::const_),
+            "Dot product 'b_i = A_ij * x_j (b_p = A_pu * x_u + A_pp * x_p == A_pp * x_p)",
+            py::arg("x_u"),
+            py::arg("x_p"))
+
+        .def(
+            "Solve",
+            py::overload_cast<const xt::xtensor<double, 1>&, const xt::xtensor<double, 1>&>(
+                &GooseFEM::MatrixDiagonalPartitioned::Solve),
+            "Solve",
+            py::arg("b"),
+            py::arg("x"))
+
+        .def(
+            "Solve",
+            py::overload_cast<const xt::xtensor<double, 2>&, const xt::xtensor<double, 2>&>(
+                &GooseFEM::MatrixDiagonalPartitioned::Solve),
+            "Solve",
+            py::arg("b"),
+            py::arg("x"))
+
+        .def(
+            "Solve_u",
+            py::overload_cast<const xt::xtensor<double, 1>&, const xt::xtensor<double, 1>&>(
+                &GooseFEM::MatrixDiagonalPartitioned::Solve_u),
+            "Solve_u",
+            py::arg("b_u"),
+            py::arg("x_p"))
+
+        .def(
+            "Reaction",
+            py::overload_cast<const xt::xtensor<double, 1>&, const xt::xtensor<double, 1>&>(
+                &GooseFEM::MatrixDiagonalPartitioned::Reaction, py::const_),
+            "Reaction",
+            py::arg("x"),
+            py::arg("b"))
+
+        .def(
+            "Reaction",
+            py::overload_cast<const xt::xtensor<double, 2>&, const xt::xtensor<double, 2>&>(
+                &GooseFEM::MatrixDiagonalPartitioned::Reaction, py::const_),
+            "Reaction",
+            py::arg("x"),
+            py::arg("b"))
+
+        .def(
+            "Reaction_p",
+            py::overload_cast<const xt::xtensor<double, 1>&, const xt::xtensor<double, 1>&>(
+                &GooseFEM::MatrixDiagonalPartitioned::Reaction_p, py::const_),
+            "Reaction_p",
+            py::arg("x_u"),
+            py::arg("x_p"))
+
+        .def("__repr__", [](const GooseFEM::MatrixDiagonalPartitioned&) {
+            return "<GooseFEM.MatrixDiagonalPartitioned>";
+        });
 }
-
-// =================================================================================================
-
diff --git a/python/MatrixPartitioned.hpp b/python/MatrixPartitioned.hpp
index 9259e5c..f0fadf0 100644
--- a/python/MatrixPartitioned.hpp
+++ b/python/MatrixPartitioned.hpp
@@ -1,138 +1,124 @@
 /* =================================================================================================
 
 (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 <GooseFEM/GooseFEM.h>
 #include <pybind11/eigen.h>
+#include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
-#include "../include/GooseFEM/GooseFEM.h"
 
 namespace py = pybind11;
 
 void init_MatrixPartitioned(py::module& m)
 {
 
-py::class_<GooseFEM::MatrixPartitioned>(m, "MatrixPartitioned")
-
-    .def(py::init<
-            const xt::xtensor<size_t,2>&,
-            const xt::xtensor<size_t,2>&,
-            const xt::xtensor<size_t,1>&>(),
-        "Sparse, partitioned, matrix",
-        py::arg("conn"),
-        py::arg("dofs"),
-        py::arg("iip"))
-
-    .def("nelem",
-        &GooseFEM::MatrixPartitioned::nelem,
-        "Number of element")
-
-    .def("nne",
-        &GooseFEM::MatrixPartitioned::nne,
-        "Number of nodes per element")
-
-    .def("nnode",
-        &GooseFEM::MatrixPartitioned::nnode,
-        "Number of nodes")
-
-    .def("ndim",
-        &GooseFEM::MatrixPartitioned::ndim,
-        "Number of dimensions")
-
-    .def("ndof",
-        &GooseFEM::MatrixPartitioned::ndof,
-        "Number of degrees-of-freedom")
-
-    .def("nnu",
-        &GooseFEM::MatrixPartitioned::nnu,
-        "Number of unknown degrees-of-freedom")
-
-    .def("nnp",
-        &GooseFEM::MatrixPartitioned::nnp,
-        "Number of prescribed degrees-of-freedom")
-
-    .def("assemble",
-        &GooseFEM::MatrixPartitioned::assemble,
-        "Assemble matrix from 'elemmat",
-        py::arg("elemmat"))
-
-    .def("dofs",
-        &GooseFEM::MatrixPartitioned::dofs,
-        "Return degrees-of-freedom")
-
-    .def("iiu",
-        &GooseFEM::MatrixPartitioned::iiu,
-        "Return unknown degrees-of-freedom")
-
-    .def("iip",
-        &GooseFEM::MatrixPartitioned::iip,
-        "Return prescribed degrees-of-freedom")
-
-    .def("Reaction",
-        py::overload_cast<const xt::xtensor<double,1>&, const xt::xtensor<double,1>&>(
-            &GooseFEM::MatrixPartitioned::Reaction, py::const_),
-        "Reaction",
-        py::arg("x"),
-        py::arg("b"))
-
-    .def("Reaction",
-        py::overload_cast<const xt::xtensor<double,2>&, const xt::xtensor<double,2>&>(
-            &GooseFEM::MatrixPartitioned::Reaction, py::const_),
-        "Reaction",
-        py::arg("x"),
-        py::arg("b"))
-
-    .def("Reaction_p",
-        py::overload_cast<const xt::xtensor<double,1>&, const xt::xtensor<double,1>&>(
-            &GooseFEM::MatrixPartitioned::Reaction_p, py::const_),
-        "Reaction_p",
-        py::arg("x_u"),
-        py::arg("x_p"))
-
-    .def("__repr__",
-        [](const GooseFEM::MatrixPartitioned&){
-            return "<GooseFEM.MatrixPartitioned>"; });
-
-py::class_<GooseFEM::MatrixPartitionedSolver<>>(m, "MatrixPartitionedSolver")
-
-    .def(py::init<>(),
-        "Sparse, partitioned, matrix solver")
-
-    .def("Solve",
-        py::overload_cast<GooseFEM::MatrixPartitioned&,
-                          const xt::xtensor<double,1>&,
-                          const xt::xtensor<double,1>&>(
-            &GooseFEM::MatrixPartitionedSolver<>::Solve),
-        "Solve",
-        py::arg("matrix"),
-        py::arg("b"),
-        py::arg("x"))
-
-    .def("Solve",
-        py::overload_cast<GooseFEM::MatrixPartitioned&,
-                          const xt::xtensor<double,2>&,
-                          const xt::xtensor<double,2>&>(
-            &GooseFEM::MatrixPartitionedSolver<>::Solve),
-        "Solve",
-        py::arg("matrix"),
-        py::arg("b"),
-        py::arg("x"))
-
-    .def("Solve_u",
-        py::overload_cast<GooseFEM::MatrixPartitioned&,
-                          const xt::xtensor<double,1>&,
-                          const xt::xtensor<double,1>&>(
-            &GooseFEM::MatrixPartitionedSolver<>::Solve_u),
-        "Solve_u",
-        py::arg("matrix"),
-        py::arg("b_u"),
-        py::arg("x_p"))
-
-    .def("__repr__",
-        [](const GooseFEM::MatrixPartitionedSolver<>&){
-            return "<GooseFEM.MatrixPartitionedSolver>"; });
-
+    py::class_<GooseFEM::MatrixPartitioned>(m, "MatrixPartitioned")
+
+        .def(
+            py::init<
+                const xt::xtensor<size_t, 2>&,
+                const xt::xtensor<size_t, 2>&,
+                const xt::xtensor<size_t, 1>&>(),
+            "Sparse, partitioned, matrix",
+            py::arg("conn"),
+            py::arg("dofs"),
+            py::arg("iip"))
+
+        .def("nelem", &GooseFEM::MatrixPartitioned::nelem, "Number of element")
+
+        .def("nne", &GooseFEM::MatrixPartitioned::nne, "Number of nodes per element")
+
+        .def("nnode", &GooseFEM::MatrixPartitioned::nnode, "Number of nodes")
+
+        .def("ndim", &GooseFEM::MatrixPartitioned::ndim, "Number of dimensions")
+
+        .def("ndof", &GooseFEM::MatrixPartitioned::ndof, "Number of degrees-of-freedom")
+
+        .def("nnu", &GooseFEM::MatrixPartitioned::nnu, "Number of unknown degrees-of-freedom")
+
+        .def("nnp", &GooseFEM::MatrixPartitioned::nnp, "Number of prescribed degrees-of-freedom")
+
+        .def(
+            "assemble",
+            &GooseFEM::MatrixPartitioned::assemble,
+            "Assemble matrix from 'elemmat",
+            py::arg("elemmat"))
+
+        .def("dofs", &GooseFEM::MatrixPartitioned::dofs, "Return degrees-of-freedom")
+
+        .def("iiu", &GooseFEM::MatrixPartitioned::iiu, "Return unknown degrees-of-freedom")
+
+        .def("iip", &GooseFEM::MatrixPartitioned::iip, "Return prescribed degrees-of-freedom")
+
+        .def(
+            "Reaction",
+            py::overload_cast<const xt::xtensor<double, 1>&, const xt::xtensor<double, 1>&>(
+                &GooseFEM::MatrixPartitioned::Reaction, py::const_),
+            "Reaction",
+            py::arg("x"),
+            py::arg("b"))
+
+        .def(
+            "Reaction",
+            py::overload_cast<const xt::xtensor<double, 2>&, const xt::xtensor<double, 2>&>(
+                &GooseFEM::MatrixPartitioned::Reaction, py::const_),
+            "Reaction",
+            py::arg("x"),
+            py::arg("b"))
+
+        .def(
+            "Reaction_p",
+            py::overload_cast<const xt::xtensor<double, 1>&, const xt::xtensor<double, 1>&>(
+                &GooseFEM::MatrixPartitioned::Reaction_p, py::const_),
+            "Reaction_p",
+            py::arg("x_u"),
+            py::arg("x_p"))
+
+        .def("__repr__", [](const GooseFEM::MatrixPartitioned&) {
+            return "<GooseFEM.MatrixPartitioned>";
+        });
+
+    py::class_<GooseFEM::MatrixPartitionedSolver<>>(m, "MatrixPartitionedSolver")
+
+        .def(py::init<>(), "Sparse, partitioned, matrix solver")
+
+        .def(
+            "Solve",
+            py::overload_cast<
+                GooseFEM::MatrixPartitioned&,
+                const xt::xtensor<double, 1>&,
+                const xt::xtensor<double, 1>&>(&GooseFEM::MatrixPartitionedSolver<>::Solve),
+            "Solve",
+            py::arg("matrix"),
+            py::arg("b"),
+            py::arg("x"))
+
+        .def(
+            "Solve",
+            py::overload_cast<
+                GooseFEM::MatrixPartitioned&,
+                const xt::xtensor<double, 2>&,
+                const xt::xtensor<double, 2>&>(&GooseFEM::MatrixPartitionedSolver<>::Solve),
+            "Solve",
+            py::arg("matrix"),
+            py::arg("b"),
+            py::arg("x"))
+
+        .def(
+            "Solve_u",
+            py::overload_cast<
+                GooseFEM::MatrixPartitioned&,
+                const xt::xtensor<double, 1>&,
+                const xt::xtensor<double, 1>&>(&GooseFEM::MatrixPartitionedSolver<>::Solve_u),
+            "Solve_u",
+            py::arg("matrix"),
+            py::arg("b_u"),
+            py::arg("x_p"))
+
+        .def("__repr__", [](const GooseFEM::MatrixPartitionedSolver<>&) {
+            return "<GooseFEM.MatrixPartitionedSolver>";
+        });
 }
diff --git a/python/MatrixPartitionedTyings.hpp b/python/MatrixPartitionedTyings.hpp
index 105f4a3..91860b6 100644
--- a/python/MatrixPartitionedTyings.hpp
+++ b/python/MatrixPartitionedTyings.hpp
@@ -1,148 +1,112 @@
 /* =================================================================================================
 
 (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 <GooseFEM/GooseFEM.h>
 #include <pybind11/eigen.h>
+#include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
-#include "../include/GooseFEM/GooseFEM.h"
-
-// =================================================================================================
 
 namespace py = pybind11;
 
-// =================================================================================================
-
 void init_MatrixPartitionedTyings(py::module& m)
 {
 
-py::class_<GooseFEM::MatrixPartitionedTyings>(m, "MatrixPartitionedTyings")
-
-    .def(py::init<
-            const xt::xtensor<size_t,2>&,
-            const xt::xtensor<size_t,2>&,
-            const Eigen::SparseMatrix<double>&,
-            const Eigen::SparseMatrix<double>&>(),
-        "Sparse, partitioned, matrix",
-        py::arg("conn"),
-        py::arg("dofs"),
-        py::arg("Cdu"),
-        py::arg("Cdp"))
-
-    .def("nelem",
-        &GooseFEM::MatrixPartitionedTyings::nelem,
-        "Number of element")
-
-    .def("nne",
-        &GooseFEM::MatrixPartitionedTyings::nne,
-        "Number of nodes per element")
-
-    .def("nnode",
-        &GooseFEM::MatrixPartitionedTyings::nnode,
-        "Number of nodes")
-
-    .def("ndim",
-        &GooseFEM::MatrixPartitionedTyings::ndim,
-        "Number of dimensions")
-
-    .def("ndof",
-        &GooseFEM::MatrixPartitionedTyings::ndof,
-        "Number of degrees-of-freedom")
-
-    .def("nnu",
-        &GooseFEM::MatrixPartitionedTyings::nnu,
-        "Number of unknown degrees-of-freedom")
-
-    .def("nnp",
-        &GooseFEM::MatrixPartitionedTyings::nnp,
-        "Number of prescribed degrees-of-freedom")
-
-    .def("nni",
-        &GooseFEM::MatrixPartitionedTyings::nni,
-        "Number of independent degrees-of-freedom")
-
-    .def("nnd",
-        &GooseFEM::MatrixPartitionedTyings::nnd,
-        "Number of dependent degrees-of-freedom")
-
-    .def("assemble",
-        &GooseFEM::MatrixPartitionedTyings::assemble,
-        "Assemble matrix from 'elemmat",
-        py::arg("elemmat"))
-
-    .def("dofs",
-        &GooseFEM::MatrixPartitionedTyings::dofs,
-        "Return degrees-of-freedom")
-
-    .def("iiu",
-        &GooseFEM::MatrixPartitionedTyings::iiu,
-        "Return unknown degrees-of-freedom")
-
-    .def("iip",
-        &GooseFEM::MatrixPartitionedTyings::iip,
-        "Return prescribed degrees-of-freedom")
-
-    .def("iii",
-        &GooseFEM::MatrixPartitionedTyings::iii,
-        "Return independent degrees-of-freedom")
-
-    .def("iid",
-        &GooseFEM::MatrixPartitionedTyings::iid,
-        "Return dependent degrees-of-freedom")
-
-    .def("__repr__",
-        [](const GooseFEM::MatrixPartitionedTyings&){
-            return "<GooseFEM.MatrixPartitionedTyings>"; });
-
-
-py::class_<GooseFEM::MatrixPartitionedTyingsSolver<>>(m, "MatrixPartitionedTyingsSolver")
-
-    .def(py::init<>(),
-        "Sparse, partitioned, matrix solver")
-
-    .def("Solve",
-        py::overload_cast<
-            GooseFEM::MatrixPartitionedTyings&,
-            const xt::xtensor<double,1>&,
-            const xt::xtensor<double,1>&>(
-                &GooseFEM::MatrixPartitionedTyingsSolver<>::Solve),
-        "Solve",
-        py::arg("matrix"),
-        py::arg("b"),
-        py::arg("x"))
-
-    .def("Solve",
-        py::overload_cast<
-            GooseFEM::MatrixPartitionedTyings&,
-            const xt::xtensor<double,2>&,
-            const xt::xtensor<double,2>&>(
-                &GooseFEM::MatrixPartitionedTyingsSolver<>::Solve),
-        "Solve",
-        py::arg("matrix"),
-        py::arg("b"),
-        py::arg("x"))
-
-    .def("Solve_u",
-        py::overload_cast<
-            GooseFEM::MatrixPartitionedTyings&,
-            const xt::xtensor<double,1>&,
-            const xt::xtensor<double,1>&,
-            const xt::xtensor<double,1>&>(
-                &GooseFEM::MatrixPartitionedTyingsSolver<>::Solve_u),
-        "Solve_u",
-        py::arg("matrix"),
-        py::arg("b_u"),
-        py::arg("b_d"),
-        py::arg("x_p"))
-
-    .def("__repr__",
-        [](const GooseFEM::MatrixPartitionedTyingsSolver<>&){
-            return "<GooseFEM.MatrixPartitionedTyingsSolver>"; });
+    py::class_<GooseFEM::MatrixPartitionedTyings>(m, "MatrixPartitionedTyings")
 
-}
+        .def(
+            py::init<
+                const xt::xtensor<size_t, 2>&,
+                const xt::xtensor<size_t, 2>&,
+                const Eigen::SparseMatrix<double>&,
+                const Eigen::SparseMatrix<double>&>(),
+            "Sparse, partitioned, matrix",
+            py::arg("conn"),
+            py::arg("dofs"),
+            py::arg("Cdu"),
+            py::arg("Cdp"))
+
+        .def("nelem", &GooseFEM::MatrixPartitionedTyings::nelem, "Number of element")
+
+        .def("nne", &GooseFEM::MatrixPartitionedTyings::nne, "Number of nodes per element")
+
+        .def("nnode", &GooseFEM::MatrixPartitionedTyings::nnode, "Number of nodes")
+
+        .def("ndim", &GooseFEM::MatrixPartitionedTyings::ndim, "Number of dimensions")
+
+        .def("ndof", &GooseFEM::MatrixPartitionedTyings::ndof, "Number of DOFs")
+
+        .def("nnu", &GooseFEM::MatrixPartitionedTyings::nnu, "Number of unknown DOFs")
+
+        .def("nnp", &GooseFEM::MatrixPartitionedTyings::nnp, "Number of prescribed DOFs")
+
+        .def("nni", &GooseFEM::MatrixPartitionedTyings::nni, "Number of independent DOFs")
 
-// =================================================================================================
+        .def("nnd", &GooseFEM::MatrixPartitionedTyings::nnd, "Number of dependent DOFs")
 
+        .def(
+            "assemble",
+            &GooseFEM::MatrixPartitionedTyings::assemble,
+            "Assemble matrix from 'elemmat",
+            py::arg("elemmat"))
+
+        .def("dofs", &GooseFEM::MatrixPartitionedTyings::dofs, "Degrees-of-freedom")
+
+        .def("iiu", &GooseFEM::MatrixPartitionedTyings::iiu, "Unknown DOFs")
+
+        .def("iip", &GooseFEM::MatrixPartitionedTyings::iip, "Prescribed DOFs")
+
+        .def("iii", &GooseFEM::MatrixPartitionedTyings::iii, "Independent DOFs")
+
+        .def("iid", &GooseFEM::MatrixPartitionedTyings::iid, "Dependent DOFs")
+
+        .def("__repr__", [](const GooseFEM::MatrixPartitionedTyings&) {
+            return "<GooseFEM.MatrixPartitionedTyings>";
+        });
+
+    py::class_<GooseFEM::MatrixPartitionedTyingsSolver<>>(m, "MatrixPartitionedTyingsSolver")
+
+        .def(py::init<>(), "Sparse, partitioned, matrix solver")
+
+        .def(
+            "Solve",
+            py::overload_cast<
+                GooseFEM::MatrixPartitionedTyings&,
+                const xt::xtensor<double, 1>&,
+                const xt::xtensor<double, 1>&>(&GooseFEM::MatrixPartitionedTyingsSolver<>::Solve),
+            "Solve",
+            py::arg("matrix"),
+            py::arg("b"),
+            py::arg("x"))
+
+        .def(
+            "Solve",
+            py::overload_cast<
+                GooseFEM::MatrixPartitionedTyings&,
+                const xt::xtensor<double, 2>&,
+                const xt::xtensor<double, 2>&>(&GooseFEM::MatrixPartitionedTyingsSolver<>::Solve),
+            "Solve",
+            py::arg("matrix"),
+            py::arg("b"),
+            py::arg("x"))
+
+        .def(
+            "Solve_u",
+            py::overload_cast<
+                GooseFEM::MatrixPartitionedTyings&,
+                const xt::xtensor<double, 1>&,
+                const xt::xtensor<double, 1>&,
+                const xt::xtensor<double, 1>&>(&GooseFEM::MatrixPartitionedTyingsSolver<>::Solve_u),
+            "Solve_u",
+            py::arg("matrix"),
+            py::arg("b_u"),
+            py::arg("b_d"),
+            py::arg("x_p"))
+
+        .def("__repr__", [](const GooseFEM::MatrixPartitionedTyingsSolver<>&) {
+            return "<GooseFEM.MatrixPartitionedTyingsSolver>";
+        });
+}
diff --git a/python/Mesh.hpp b/python/Mesh.hpp
index 0e2d9c5..6692020 100644
--- a/python/Mesh.hpp
+++ b/python/Mesh.hpp
@@ -1,117 +1,113 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
+#include <GooseFEM/GooseFEM.h>
 #include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
-#include "../include/GooseFEM/GooseFEM.h"
-
-// =================================================================================================
 
 namespace py = pybind11;
 
-// =================================================================================================
-
 void init_Mesh(py::module& m)
 {
 
-// -------------------------------------------------------------------------------------------------
-
-py::enum_<GooseFEM::Mesh::ElementType>(m, "ElementType", "ElementType")
-    .value("Tri3", GooseFEM::Mesh::ElementType::Tri3)
-    .value("Quad4", GooseFEM::Mesh::ElementType::Quad4)
-    .value("Hex8", GooseFEM::Mesh::ElementType::Hex8)
-    .export_values();
-
-// -------------------------------------------------------------------------------------------------
-
-py::class_<GooseFEM::Mesh::Renumber>(m, "Renumber")
-
-    .def(py::init<const xt::xarray<size_t>&>(),
-        "Class to renumber to the lowest possible indices. Use ``Renumber(...).get()`` to get the renumbered result",
-        py::arg("dofs"))
-
-    .def("get",
-        &GooseFEM::Mesh::Renumber::get,
-        "Get result of renumbering")
-
-    .def("index",
-        &GooseFEM::Mesh::Renumber::index,
-        "Get index list to apply renumbering. Apply renumbering using ``index[dofs]``")
-
-    .def("__repr__",
-        [](const GooseFEM::Mesh::Renumber&){ return "<GooseFEM.Mesh.Renumber>"; });
-
-// -------------------------------------------------------------------------------------------------
-
-py::class_<GooseFEM::Mesh::Reorder>(m, "Reorder")
-
-    .def(py::init([](xt::xtensor<size_t,1>& a)
-        { return new GooseFEM::Mesh::Reorder({a}); }))
-
-    .def(py::init([](xt::xtensor<size_t,1>& a, xt::xtensor<size_t,1>& b)
-        { return new GooseFEM::Mesh::Reorder({a,b}); }))
-
-    .def(py::init([](xt::xtensor<size_t,1>& a, xt::xtensor<size_t,1>& b, xt::xtensor<size_t,1>& c)
-        { return new GooseFEM::Mesh::Reorder({a,b,c}); }))
-
-    .def(py::init([](xt::xtensor<size_t,1>& a, xt::xtensor<size_t,1>& b, xt::xtensor<size_t,1>& c, xt::xtensor<size_t,1>& d)
-        { return new GooseFEM::Mesh::Reorder({a,b,c,d}); }))
-
-    .def("get",
-        &GooseFEM::Mesh::Reorder::get,
-        "Reorder matrix (e.g. ``dofs``)")
-
-    .def("index",
-        &GooseFEM::Mesh::Reorder::index,
-        "Get index list to apply renumbering. Apply renumbering using ``index[dofs]``")
-
-    .def("__repr__",
-        [](const GooseFEM::Mesh::Reorder &){ return "<GooseFEM.Mesh.Reorder>"; });
-
-// -------------------------------------------------------------------------------------------------
-
-m.def("dofs",
-    &GooseFEM::Mesh::dofs,
-    "List with DOF-numbers (in sequential order)",
-    py::arg("nnode"),
-    py::arg("ndim"));
-
-m.def("renumber",
-    &GooseFEM::Mesh::renumber,
-    "Renumber to lowest possible indices. Use ``GooseFEM.Mesh.Renumber`` for advanced functionality",
-    py::arg("dofs"));
-
-m.def("coordination",
-    &GooseFEM::Mesh::coordination,
-    "Coordination number of each node (number of elements connected to each node)",
-    py::arg("conn"));
-
-m.def("elem2node",
-    &GooseFEM::Mesh::elem2node,
-    "Element-numbers connected to each node",
-    py::arg("conn"));
-
-m.def("edgesize",
-    py::overload_cast<const xt::xtensor<double,2>&, const xt::xtensor<size_t,2>&>(
-        &GooseFEM::Mesh::edgesize),
-    "Get the edge size of all elements",
-    py::arg("coor"),
-    py::arg("conn"));
-
-m.def("edgesize",
-    py::overload_cast<
-        const xt::xtensor<double,2>&,
-        const xt::xtensor<size_t,2>&,
-        GooseFEM::Mesh::ElementType>(&GooseFEM::Mesh::edgesize),
-    "Get the edge size of all elements",
-    py::arg("coor"),
-    py::arg("conn"),
-    py::arg("type"));
-
+    py::enum_<GooseFEM::Mesh::ElementType>(m, "ElementType", "ElementType")
+        .value("Tri3", GooseFEM::Mesh::ElementType::Tri3)
+        .value("Quad4", GooseFEM::Mesh::ElementType::Quad4)
+        .value("Hex8", GooseFEM::Mesh::ElementType::Hex8)
+        .export_values();
+
+    py::class_<GooseFEM::Mesh::Renumber>(m, "Renumber")
+
+        .def(
+            py::init<const xt::xarray<size_t>&>(),
+            "Class to renumber to the lowest possible indices. Use ``Renumber(...).get()`` to get "
+            "the renumbered result",
+            py::arg("dofs"))
+
+        .def("get", &GooseFEM::Mesh::Renumber::get, "Get result of renumbering")
+
+        .def(
+            "index",
+            &GooseFEM::Mesh::Renumber::index,
+            "Get index list to apply renumbering. Apply renumbering using ``index[dofs]``")
+
+        .def(
+            "__repr__", [](const GooseFEM::Mesh::Renumber&) { return "<GooseFEM.Mesh.Renumber>"; });
+
+    py::class_<GooseFEM::Mesh::Reorder>(m, "Reorder")
+
+        .def(py::init([](xt::xtensor<size_t, 1>& a) { return new GooseFEM::Mesh::Reorder({a}); }))
+
+        .def(py::init([](xt::xtensor<size_t, 1>& a, xt::xtensor<size_t, 1>& b) {
+            return new GooseFEM::Mesh::Reorder({a, b});
+        }))
+
+        .def(py::init(
+            [](xt::xtensor<size_t, 1>& a, xt::xtensor<size_t, 1>& b, xt::xtensor<size_t, 1>& c) {
+                return new GooseFEM::Mesh::Reorder({a, b, c});
+            }))
+
+        .def(py::init([](xt::xtensor<size_t, 1>& a,
+                         xt::xtensor<size_t, 1>& b,
+                         xt::xtensor<size_t, 1>& c,
+                         xt::xtensor<size_t, 1>& d) {
+            return new GooseFEM::Mesh::Reorder({a, b, c, d});
+        }))
+
+        .def("get", &GooseFEM::Mesh::Reorder::get, "Reorder matrix (e.g. ``dofs``)")
+
+        .def(
+            "index",
+            &GooseFEM::Mesh::Reorder::index,
+            "Get index list to apply renumbering. Apply renumbering using ``index[dofs]``")
+
+        .def("__repr__", [](const GooseFEM::Mesh::Reorder&) { return "<GooseFEM.Mesh.Reorder>"; });
+
+    m.def(
+        "dofs",
+        &GooseFEM::Mesh::dofs,
+        "List with DOF-numbers (in sequential order)",
+        py::arg("nnode"),
+        py::arg("ndim"));
+
+    m.def(
+        "renumber",
+        &GooseFEM::Mesh::renumber,
+        "Renumber to lowest possible indices. Use ``GooseFEM.Mesh.Renumber`` for advanced "
+        "functionality",
+        py::arg("dofs"));
+
+    m.def(
+        "coordination",
+        &GooseFEM::Mesh::coordination,
+        "Coordination number of each node (number of elements connected to each node)",
+        py::arg("conn"));
+
+    m.def(
+        "elem2node",
+        &GooseFEM::Mesh::elem2node,
+        "Element-numbers connected to each node",
+        py::arg("conn"));
+
+    m.def(
+        "edgesize",
+        py::overload_cast<const xt::xtensor<double, 2>&, const xt::xtensor<size_t, 2>&>(
+            &GooseFEM::Mesh::edgesize),
+        "Get the edge size of all elements",
+        py::arg("coor"),
+        py::arg("conn"));
+
+    m.def(
+        "edgesize",
+        py::overload_cast<
+            const xt::xtensor<double, 2>&,
+            const xt::xtensor<size_t, 2>&,
+            GooseFEM::Mesh::ElementType>(&GooseFEM::Mesh::edgesize),
+        "Get the edge size of all elements",
+        py::arg("coor"),
+        py::arg("conn"),
+        py::arg("type"));
 }
-
-// =================================================================================================
-
diff --git a/python/MeshHex8.hpp b/python/MeshHex8.hpp
index 721f739..323b0f1 100644
--- a/python/MeshHex8.hpp
+++ b/python/MeshHex8.hpp
@@ -1,776 +1,612 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
+#include <GooseFEM/GooseFEM.h>
 #include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
-#include "../include/GooseFEM/GooseFEM.h"
-
-// =================================================================================================
 
 namespace py = pybind11;
 
-// =================================================================================================
-
 void init_MeshHex8(py::module& m)
 {
 
-py::class_<GooseFEM::Mesh::Hex8::Regular>(m, "Regular")
-
-    .def(py::init<size_t,size_t,size_t,double>(),
-        "Mesh with nx*ny*nz 'pixels' and edge size h",
-        py::arg("nx"),
-        py::arg("ny"),
-        py::arg("nz"),
-        py::arg("h")=1.)
+    py::class_<GooseFEM::Mesh::Hex8::Regular>(m, "Regular")
 
-    .def("nelem",
-        &GooseFEM::Mesh::Hex8::Regular::nelem)
+        .def(
+            py::init<size_t, size_t, size_t, double>(),
+            "Mesh with nx*ny*nz 'pixels' and edge size h",
+            py::arg("nx"),
+            py::arg("ny"),
+            py::arg("nz"),
+            py::arg("h") = 1.)
 
-    .def("nnode",
-        &GooseFEM::Mesh::Hex8::Regular::nnode)
+        .def("nelem", &GooseFEM::Mesh::Hex8::Regular::nelem)
 
-    .def("nne",
-        &GooseFEM::Mesh::Hex8::Regular::nne)
+        .def("nnode", &GooseFEM::Mesh::Hex8::Regular::nnode)
 
-    .def("ndim",
-        &GooseFEM::Mesh::Hex8::Regular::ndim)
+        .def("nne", &GooseFEM::Mesh::Hex8::Regular::nne)
 
-    .def("coor",
-        &GooseFEM::Mesh::Hex8::Regular::coor)
+        .def("ndim", &GooseFEM::Mesh::Hex8::Regular::ndim)
 
-    .def("conn",
-        &GooseFEM::Mesh::Hex8::Regular::conn)
+        .def("coor", &GooseFEM::Mesh::Hex8::Regular::coor)
 
-    .def("getElementType",
-        &GooseFEM::Mesh::Hex8::Regular::getElementType)
+        .def("conn", &GooseFEM::Mesh::Hex8::Regular::conn)
 
-    .def("nodesFront",
-        &GooseFEM::Mesh::Hex8::Regular::nodesFront)
+        .def("getElementType", &GooseFEM::Mesh::Hex8::Regular::getElementType)
 
-    .def("nodesBack",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBack)
+        .def("nodesFront", &GooseFEM::Mesh::Hex8::Regular::nodesFront)
 
-    .def("nodesLeft",
-        &GooseFEM::Mesh::Hex8::Regular::nodesLeft)
+        .def("nodesBack", &GooseFEM::Mesh::Hex8::Regular::nodesBack)
 
-    .def("nodesRight",
-        &GooseFEM::Mesh::Hex8::Regular::nodesRight)
+        .def("nodesLeft", &GooseFEM::Mesh::Hex8::Regular::nodesLeft)
 
-    .def("nodesBottom",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBottom)
+        .def("nodesRight", &GooseFEM::Mesh::Hex8::Regular::nodesRight)
 
-    .def("nodesTop",
-        &GooseFEM::Mesh::Hex8::Regular::nodesTop)
+        .def("nodesBottom", &GooseFEM::Mesh::Hex8::Regular::nodesBottom)
 
-    .def("nodesFrontFace",
-        &GooseFEM::Mesh::Hex8::Regular::nodesFrontFace)
+        .def("nodesTop", &GooseFEM::Mesh::Hex8::Regular::nodesTop)
 
-    .def("nodesBackFace",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBackFace)
+        .def("nodesFrontFace", &GooseFEM::Mesh::Hex8::Regular::nodesFrontFace)
 
-    .def("nodesLeftFace",
-        &GooseFEM::Mesh::Hex8::Regular::nodesLeftFace)
+        .def("nodesBackFace", &GooseFEM::Mesh::Hex8::Regular::nodesBackFace)
 
-    .def("nodesRightFace",
-        &GooseFEM::Mesh::Hex8::Regular::nodesRightFace)
+        .def("nodesLeftFace", &GooseFEM::Mesh::Hex8::Regular::nodesLeftFace)
 
-    .def("nodesBottomFace",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBottomFace)
+        .def("nodesRightFace", &GooseFEM::Mesh::Hex8::Regular::nodesRightFace)
 
-    .def("nodesTopFace",
-        &GooseFEM::Mesh::Hex8::Regular::nodesTopFace)
+        .def("nodesBottomFace", &GooseFEM::Mesh::Hex8::Regular::nodesBottomFace)
 
-    .def("nodesFrontBottomEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomEdge)
+        .def("nodesTopFace", &GooseFEM::Mesh::Hex8::Regular::nodesTopFace)
 
-    .def("nodesFrontTopEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopEdge)
+        .def("nodesFrontBottomEdge", &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomEdge)
 
-    .def("nodesFrontLeftEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftEdge)
+        .def("nodesFrontTopEdge", &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopEdge)
 
-    .def("nodesFrontRightEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightEdge)
+        .def("nodesFrontLeftEdge", &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftEdge)
 
-    .def("nodesBackBottomEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomEdge)
+        .def("nodesFrontRightEdge", &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightEdge)
 
-    .def("nodesBackTopEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBackTopEdge)
+        .def("nodesBackBottomEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomEdge)
 
-    .def("nodesBackLeftEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftEdge)
+        .def("nodesBackTopEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBackTopEdge)
 
-    .def("nodesBackRightEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBackRightEdge)
+        .def("nodesBackLeftEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftEdge)
 
-    .def("nodesBottomLeftEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftEdge)
+        .def("nodesBackRightEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBackRightEdge)
 
-    .def("nodesBottomRightEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightEdge)
+        .def("nodesBottomLeftEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftEdge)
 
-    .def("nodesTopLeftEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftEdge)
+        .def("nodesBottomRightEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightEdge)
 
-    .def("nodesTopRightEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesTopRightEdge)
+        .def("nodesTopLeftEdge", &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftEdge)
 
-    .def("nodesBottomFrontEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontEdge)
+        .def("nodesTopRightEdge", &GooseFEM::Mesh::Hex8::Regular::nodesTopRightEdge)
 
-    .def("nodesBottomBackEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackEdge)
+        .def("nodesBottomFrontEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontEdge)
 
-    .def("nodesTopFrontEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontEdge)
+        .def("nodesBottomBackEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackEdge)
 
-    .def("nodesTopBackEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesTopBackEdge)
+        .def("nodesTopFrontEdge", &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontEdge)
 
-    .def("nodesLeftBottomEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomEdge)
+        .def("nodesTopBackEdge", &GooseFEM::Mesh::Hex8::Regular::nodesTopBackEdge)
 
-    .def("nodesLeftFrontEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontEdge)
+        .def("nodesLeftBottomEdge", &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomEdge)
 
-    .def("nodesLeftBackEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackEdge)
+        .def("nodesLeftFrontEdge", &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontEdge)
 
-    .def("nodesLeftTopEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopEdge)
+        .def("nodesLeftBackEdge", &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackEdge)
 
-    .def("nodesRightBottomEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomEdge)
+        .def("nodesLeftTopEdge", &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopEdge)
 
-    .def("nodesRightTopEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesRightTopEdge)
+        .def("nodesRightBottomEdge", &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomEdge)
 
-    .def("nodesRightFrontEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontEdge)
+        .def("nodesRightTopEdge", &GooseFEM::Mesh::Hex8::Regular::nodesRightTopEdge)
 
-    .def("nodesRightBackEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesRightBackEdge)
+        .def("nodesRightFrontEdge", &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontEdge)
 
-    .def("nodesFrontBottomOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomOpenEdge)
+        .def("nodesRightBackEdge", &GooseFEM::Mesh::Hex8::Regular::nodesRightBackEdge)
 
-    .def("nodesFrontTopOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopOpenEdge)
+        .def("nodesFrontBottomOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomOpenEdge)
 
-    .def("nodesFrontLeftOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftOpenEdge)
+        .def("nodesFrontTopOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopOpenEdge)
 
-    .def("nodesFrontRightOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightOpenEdge)
+        .def("nodesFrontLeftOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftOpenEdge)
 
-    .def("nodesBackBottomOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomOpenEdge)
+        .def("nodesFrontRightOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightOpenEdge)
 
-    .def("nodesBackTopOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBackTopOpenEdge)
+        .def("nodesBackBottomOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomOpenEdge)
 
-    .def("nodesBackLeftOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftOpenEdge)
+        .def("nodesBackTopOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBackTopOpenEdge)
 
-    .def("nodesBackRightOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBackRightOpenEdge)
+        .def("nodesBackLeftOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftOpenEdge)
 
-    .def("nodesBottomLeftOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftOpenEdge)
+        .def("nodesBackRightOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBackRightOpenEdge)
 
-    .def("nodesBottomRightOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightOpenEdge)
+        .def("nodesBottomLeftOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftOpenEdge)
 
-    .def("nodesTopLeftOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftOpenEdge)
+        .def("nodesBottomRightOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightOpenEdge)
 
-    .def("nodesTopRightOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesTopRightOpenEdge)
+        .def("nodesTopLeftOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftOpenEdge)
 
-    .def("nodesBottomFrontOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontOpenEdge)
+        .def("nodesTopRightOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesTopRightOpenEdge)
 
-    .def("nodesBottomBackOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackOpenEdge)
+        .def("nodesBottomFrontOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontOpenEdge)
 
-    .def("nodesTopFrontOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontOpenEdge)
+        .def("nodesBottomBackOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackOpenEdge)
 
-    .def("nodesTopBackOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesTopBackOpenEdge)
+        .def("nodesTopFrontOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontOpenEdge)
 
-    .def("nodesLeftBottomOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomOpenEdge)
+        .def("nodesTopBackOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesTopBackOpenEdge)
 
-    .def("nodesLeftFrontOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontOpenEdge)
+        .def("nodesLeftBottomOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomOpenEdge)
 
-    .def("nodesLeftBackOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackOpenEdge)
+        .def("nodesLeftFrontOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontOpenEdge)
 
-    .def("nodesLeftTopOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopOpenEdge)
+        .def("nodesLeftBackOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackOpenEdge)
 
-    .def("nodesRightBottomOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomOpenEdge)
+        .def("nodesLeftTopOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopOpenEdge)
 
-    .def("nodesRightTopOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesRightTopOpenEdge)
+        .def("nodesRightBottomOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomOpenEdge)
 
-    .def("nodesRightFrontOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontOpenEdge)
+        .def("nodesRightTopOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesRightTopOpenEdge)
 
-    .def("nodesRightBackOpenEdge",
-        &GooseFEM::Mesh::Hex8::Regular::nodesRightBackOpenEdge)
+        .def("nodesRightFrontOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontOpenEdge)
 
-    .def("nodesFrontBottomLeftCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomLeftCorner)
+        .def("nodesRightBackOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesRightBackOpenEdge)
 
-    .def("nodesFrontBottomRightCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomRightCorner)
+        .def(
+            "nodesFrontBottomLeftCorner",
+            &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomLeftCorner)
 
-    .def("nodesFrontTopLeftCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopLeftCorner)
+        .def(
+            "nodesFrontBottomRightCorner",
+            &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomRightCorner)
 
-    .def("nodesFrontTopRightCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopRightCorner)
+        .def("nodesFrontTopLeftCorner", &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopLeftCorner)
 
-    .def("nodesBackBottomLeftCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomLeftCorner)
+        .def("nodesFrontTopRightCorner", &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopRightCorner)
 
-    .def("nodesBackBottomRightCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomRightCorner)
+        .def("nodesBackBottomLeftCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomLeftCorner)
 
-    .def("nodesBackTopLeftCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBackTopLeftCorner)
+        .def(
+            "nodesBackBottomRightCorner",
+            &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomRightCorner)
 
-    .def("nodesBackTopRightCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBackTopRightCorner)
+        .def("nodesBackTopLeftCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBackTopLeftCorner)
 
-    .def("nodesFrontLeftBottomCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftBottomCorner)
+        .def("nodesBackTopRightCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBackTopRightCorner)
 
-    .def("nodesBottomFrontLeftCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontLeftCorner)
+        .def(
+            "nodesFrontLeftBottomCorner",
+            &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftBottomCorner)
 
-    .def("nodesBottomLeftFrontCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftFrontCorner)
+        .def(
+            "nodesBottomFrontLeftCorner",
+            &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontLeftCorner)
 
-    .def("nodesLeftFrontBottomCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontBottomCorner)
+        .def(
+            "nodesBottomLeftFrontCorner",
+            &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftFrontCorner)
 
-    .def("nodesLeftBottomFrontCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomFrontCorner)
+        .def(
+            "nodesLeftFrontBottomCorner",
+            &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontBottomCorner)
 
-    .def("nodesFrontRightBottomCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightBottomCorner)
+        .def(
+            "nodesLeftBottomFrontCorner",
+            &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomFrontCorner)
 
-    .def("nodesBottomFrontRightCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontRightCorner)
+        .def(
+            "nodesFrontRightBottomCorner",
+            &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightBottomCorner)
 
-    .def("nodesBottomRightFrontCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightFrontCorner)
+        .def(
+            "nodesBottomFrontRightCorner",
+            &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontRightCorner)
 
-    .def("nodesRightFrontBottomCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontBottomCorner)
+        .def(
+            "nodesBottomRightFrontCorner",
+            &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightFrontCorner)
 
-    .def("nodesRightBottomFrontCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomFrontCorner)
+        .def(
+            "nodesRightFrontBottomCorner",
+            &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontBottomCorner)
 
-    .def("nodesFrontLeftTopCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftTopCorner)
+        .def(
+            "nodesRightBottomFrontCorner",
+            &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomFrontCorner)
 
-    .def("nodesTopFrontLeftCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontLeftCorner)
+        .def("nodesFrontLeftTopCorner", &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftTopCorner)
 
-    .def("nodesTopLeftFrontCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftFrontCorner)
+        .def("nodesTopFrontLeftCorner", &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontLeftCorner)
 
-    .def("nodesLeftFrontTopCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontTopCorner)
+        .def("nodesTopLeftFrontCorner", &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftFrontCorner)
 
-    .def("nodesLeftTopFrontCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopFrontCorner)
+        .def("nodesLeftFrontTopCorner", &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontTopCorner)
 
-    .def("nodesFrontRightTopCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightTopCorner)
+        .def("nodesLeftTopFrontCorner", &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopFrontCorner)
 
-    .def("nodesTopFrontRightCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontRightCorner)
+        .def("nodesFrontRightTopCorner", &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightTopCorner)
 
-    .def("nodesTopRightFrontCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesTopRightFrontCorner)
+        .def("nodesTopFrontRightCorner", &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontRightCorner)
 
-    .def("nodesRightFrontTopCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontTopCorner)
+        .def("nodesTopRightFrontCorner", &GooseFEM::Mesh::Hex8::Regular::nodesTopRightFrontCorner)
 
-    .def("nodesRightTopFrontCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesRightTopFrontCorner)
+        .def("nodesRightFrontTopCorner", &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontTopCorner)
 
-    .def("nodesBackLeftBottomCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftBottomCorner)
+        .def("nodesRightTopFrontCorner", &GooseFEM::Mesh::Hex8::Regular::nodesRightTopFrontCorner)
 
-    .def("nodesBottomBackLeftCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackLeftCorner)
+        .def("nodesBackLeftBottomCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftBottomCorner)
 
-    .def("nodesBottomLeftBackCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftBackCorner)
+        .def("nodesBottomBackLeftCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackLeftCorner)
 
-    .def("nodesLeftBackBottomCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackBottomCorner)
+        .def("nodesBottomLeftBackCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftBackCorner)
 
-    .def("nodesLeftBottomBackCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomBackCorner)
+        .def("nodesLeftBackBottomCorner", &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackBottomCorner)
 
-    .def("nodesBackRightBottomCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBackRightBottomCorner)
+        .def("nodesLeftBottomBackCorner", &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomBackCorner)
 
-    .def("nodesBottomBackRightCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackRightCorner)
+        .def(
+            "nodesBackRightBottomCorner",
+            &GooseFEM::Mesh::Hex8::Regular::nodesBackRightBottomCorner)
 
-    .def("nodesBottomRightBackCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightBackCorner)
+        .def(
+            "nodesBottomBackRightCorner",
+            &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackRightCorner)
 
-    .def("nodesRightBackBottomCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesRightBackBottomCorner)
+        .def(
+            "nodesBottomRightBackCorner",
+            &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightBackCorner)
 
-    .def("nodesRightBottomBackCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomBackCorner)
+        .def(
+            "nodesRightBackBottomCorner",
+            &GooseFEM::Mesh::Hex8::Regular::nodesRightBackBottomCorner)
 
-    .def("nodesBackLeftTopCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftTopCorner)
+        .def(
+            "nodesRightBottomBackCorner",
+            &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomBackCorner)
 
-    .def("nodesTopBackLeftCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesTopBackLeftCorner)
+        .def("nodesBackLeftTopCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftTopCorner)
 
-    .def("nodesTopLeftBackCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftBackCorner)
+        .def("nodesTopBackLeftCorner", &GooseFEM::Mesh::Hex8::Regular::nodesTopBackLeftCorner)
 
-    .def("nodesLeftBackTopCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackTopCorner)
+        .def("nodesTopLeftBackCorner", &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftBackCorner)
 
-    .def("nodesLeftTopBackCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopBackCorner)
+        .def("nodesLeftBackTopCorner", &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackTopCorner)
 
-    .def("nodesBackRightTopCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesBackRightTopCorner)
+        .def("nodesLeftTopBackCorner", &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopBackCorner)
 
-    .def("nodesTopBackRightCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesTopBackRightCorner)
+        .def("nodesBackRightTopCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBackRightTopCorner)
 
-    .def("nodesTopRightBackCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesTopRightBackCorner)
+        .def("nodesTopBackRightCorner", &GooseFEM::Mesh::Hex8::Regular::nodesTopBackRightCorner)
 
-    .def("nodesRightBackTopCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesRightBackTopCorner)
+        .def("nodesTopRightBackCorner", &GooseFEM::Mesh::Hex8::Regular::nodesTopRightBackCorner)
 
-    .def("nodesRightTopBackCorner",
-        &GooseFEM::Mesh::Hex8::Regular::nodesRightTopBackCorner)
+        .def("nodesRightBackTopCorner", &GooseFEM::Mesh::Hex8::Regular::nodesRightBackTopCorner)
 
-    .def("nodesPeriodic",
-        &GooseFEM::Mesh::Hex8::Regular::nodesPeriodic)
+        .def("nodesRightTopBackCorner", &GooseFEM::Mesh::Hex8::Regular::nodesRightTopBackCorner)
 
-    .def("nodesOrigin",
-        &GooseFEM::Mesh::Hex8::Regular::nodesOrigin)
+        .def("nodesPeriodic", &GooseFEM::Mesh::Hex8::Regular::nodesPeriodic)
 
-    .def("dofs",
-        &GooseFEM::Mesh::Hex8::Regular::dofs)
+        .def("nodesOrigin", &GooseFEM::Mesh::Hex8::Regular::nodesOrigin)
 
-    .def("dofsPeriodic",
-        &GooseFEM::Mesh::Hex8::Regular::dofsPeriodic)
+        .def("dofs", &GooseFEM::Mesh::Hex8::Regular::dofs)
 
-    .def("__repr__",
-        [](const GooseFEM::Mesh::Hex8::Regular &){ return "<GooseFEM.Mesh.Hex8.Regular>"; });
+        .def("dofsPeriodic", &GooseFEM::Mesh::Hex8::Regular::dofsPeriodic)
 
-// -------------------------------------------------------------------------------------------------
+        .def("__repr__", [](const GooseFEM::Mesh::Hex8::Regular&) {
+            return "<GooseFEM.Mesh.Hex8.Regular>";
+        });
 
-py::class_<GooseFEM::Mesh::Hex8::FineLayer>(m, "FineLayer")
+    py::class_<GooseFEM::Mesh::Hex8::FineLayer>(m, "FineLayer")
 
-    .def(py::init<size_t,size_t,size_t,double,size_t>(),
-        "Mesh with nx*ny*nz 'pixels' and edge size h",
-        py::arg("nx"),
-        py::arg("ny"),
-        py::arg("nz"),
-        py::arg("h")=1.0,
-        py::arg("nfine")=1)
+        .def(
+            py::init<size_t, size_t, size_t, double, size_t>(),
+            "Mesh with nx*ny*nz 'pixels' and edge size h",
+            py::arg("nx"),
+            py::arg("ny"),
+            py::arg("nz"),
+            py::arg("h") = 1.0,
+            py::arg("nfine") = 1)
 
-    .def("nelem",
-        &GooseFEM::Mesh::Hex8::FineLayer::nelem)
+        .def("nelem", &GooseFEM::Mesh::Hex8::FineLayer::nelem)
 
-    .def("nnode",
-        &GooseFEM::Mesh::Hex8::FineLayer::nnode)
+        .def("nnode", &GooseFEM::Mesh::Hex8::FineLayer::nnode)
 
-    .def("nne",
-        &GooseFEM::Mesh::Hex8::FineLayer::nne)
+        .def("nne", &GooseFEM::Mesh::Hex8::FineLayer::nne)
 
-    .def("ndim",
-        &GooseFEM::Mesh::Hex8::FineLayer::ndim)
+        .def("ndim", &GooseFEM::Mesh::Hex8::FineLayer::ndim)
 
-    .def("nelx",
-        &GooseFEM::Mesh::Hex8::FineLayer::nelx)
+        .def("nelx", &GooseFEM::Mesh::Hex8::FineLayer::nelx)
 
-    .def("nely",
-        &GooseFEM::Mesh::Hex8::FineLayer::nely)
+        .def("nely", &GooseFEM::Mesh::Hex8::FineLayer::nely)
 
-    .def("nelz",
-        &GooseFEM::Mesh::Hex8::FineLayer::nelz)
+        .def("nelz", &GooseFEM::Mesh::Hex8::FineLayer::nelz)
 
-    .def("coor",
-        &GooseFEM::Mesh::Hex8::FineLayer::coor)
+        .def("coor", &GooseFEM::Mesh::Hex8::FineLayer::coor)
 
-    .def("conn",
-        &GooseFEM::Mesh::Hex8::FineLayer::conn)
+        .def("conn", &GooseFEM::Mesh::Hex8::FineLayer::conn)
 
-    .def("getElementType",
-        &GooseFEM::Mesh::Hex8::FineLayer::getElementType)
+        .def("getElementType", &GooseFEM::Mesh::Hex8::FineLayer::getElementType)
 
-    .def("elementsMiddleLayer",
-        &GooseFEM::Mesh::Hex8::FineLayer::elementsMiddleLayer)
+        .def("elementsMiddleLayer", &GooseFEM::Mesh::Hex8::FineLayer::elementsMiddleLayer)
 
-    .def("nodesFront",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesFront)
+        .def("nodesFront", &GooseFEM::Mesh::Hex8::FineLayer::nodesFront)
 
-    .def("nodesBack",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBack)
+        .def("nodesBack", &GooseFEM::Mesh::Hex8::FineLayer::nodesBack)
 
-    .def("nodesLeft",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesLeft)
+        .def("nodesLeft", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeft)
 
-    .def("nodesRight",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesRight)
+        .def("nodesRight", &GooseFEM::Mesh::Hex8::FineLayer::nodesRight)
 
-    .def("nodesBottom",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBottom)
+        .def("nodesBottom", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottom)
 
-    .def("nodesTop",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesTop)
+        .def("nodesTop", &GooseFEM::Mesh::Hex8::FineLayer::nodesTop)
 
-    .def("nodesFrontFace",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontFace)
+        .def("nodesFrontFace", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontFace)
 
-    .def("nodesBackFace",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBackFace)
+        .def("nodesBackFace", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackFace)
 
-    .def("nodesLeftFace",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFace)
+        .def("nodesLeftFace", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFace)
 
-    .def("nodesRightFace",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFace)
+        .def("nodesRightFace", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFace)
 
-    .def("nodesBottomFace",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFace)
+        .def("nodesBottomFace", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFace)
 
-    .def("nodesTopFace",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFace)
+        .def("nodesTopFace", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFace)
 
-    .def("nodesFrontBottomEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomEdge)
+        .def("nodesFrontBottomEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomEdge)
 
-    .def("nodesFrontTopEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopEdge)
+        .def("nodesFrontTopEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopEdge)
 
-    .def("nodesFrontLeftEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftEdge)
+        .def("nodesFrontLeftEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftEdge)
 
-    .def("nodesFrontRightEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightEdge)
+        .def("nodesFrontRightEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightEdge)
 
-    .def("nodesBackBottomEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomEdge)
+        .def("nodesBackBottomEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomEdge)
 
-    .def("nodesBackTopEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopEdge)
+        .def("nodesBackTopEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopEdge)
 
-    .def("nodesBackLeftEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftEdge)
+        .def("nodesBackLeftEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftEdge)
 
-    .def("nodesBackRightEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightEdge)
+        .def("nodesBackRightEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightEdge)
 
-    .def("nodesBottomLeftEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftEdge)
+        .def("nodesBottomLeftEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftEdge)
 
-    .def("nodesBottomRightEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightEdge)
+        .def("nodesBottomRightEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightEdge)
 
-    .def("nodesTopLeftEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftEdge)
+        .def("nodesTopLeftEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftEdge)
 
-    .def("nodesTopRightEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightEdge)
+        .def("nodesTopRightEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightEdge)
 
-    .def("nodesBottomFrontEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontEdge)
+        .def("nodesBottomFrontEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontEdge)
 
-    .def("nodesBottomBackEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackEdge)
+        .def("nodesBottomBackEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackEdge)
 
-    .def("nodesTopFrontEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontEdge)
+        .def("nodesTopFrontEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontEdge)
 
-    .def("nodesTopBackEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackEdge)
+        .def("nodesTopBackEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackEdge)
 
-    .def("nodesLeftBottomEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomEdge)
+        .def("nodesLeftBottomEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomEdge)
 
-    .def("nodesLeftFrontEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontEdge)
+        .def("nodesLeftFrontEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontEdge)
 
-    .def("nodesLeftBackEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackEdge)
+        .def("nodesLeftBackEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackEdge)
 
-    .def("nodesLeftTopEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopEdge)
+        .def("nodesLeftTopEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopEdge)
 
-    .def("nodesRightBottomEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomEdge)
+        .def("nodesRightBottomEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomEdge)
 
-    .def("nodesRightTopEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopEdge)
+        .def("nodesRightTopEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopEdge)
 
-    .def("nodesRightFrontEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontEdge)
+        .def("nodesRightFrontEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontEdge)
 
-    .def("nodesRightBackEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackEdge)
+        .def("nodesRightBackEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackEdge)
 
-    .def("nodesFrontBottomOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomOpenEdge)
+        .def("nodesFrontBottomOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomOpenEdge)
 
-    .def("nodesFrontTopOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopOpenEdge)
+        .def("nodesFrontTopOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopOpenEdge)
 
-    .def("nodesFrontLeftOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftOpenEdge)
+        .def("nodesFrontLeftOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftOpenEdge)
 
-    .def("nodesFrontRightOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightOpenEdge)
+        .def("nodesFrontRightOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightOpenEdge)
 
-    .def("nodesBackBottomOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomOpenEdge)
+        .def("nodesBackBottomOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomOpenEdge)
 
-    .def("nodesBackTopOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopOpenEdge)
+        .def("nodesBackTopOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopOpenEdge)
 
-    .def("nodesBackLeftOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftOpenEdge)
+        .def("nodesBackLeftOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftOpenEdge)
 
-    .def("nodesBackRightOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightOpenEdge)
+        .def("nodesBackRightOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightOpenEdge)
 
-    .def("nodesBottomLeftOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftOpenEdge)
+        .def("nodesBottomLeftOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftOpenEdge)
 
-    .def("nodesBottomRightOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightOpenEdge)
+        .def("nodesBottomRightOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightOpenEdge)
 
-    .def("nodesTopLeftOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftOpenEdge)
+        .def("nodesTopLeftOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftOpenEdge)
 
-    .def("nodesTopRightOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightOpenEdge)
+        .def("nodesTopRightOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightOpenEdge)
 
-    .def("nodesBottomFrontOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontOpenEdge)
+        .def("nodesBottomFrontOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontOpenEdge)
 
-    .def("nodesBottomBackOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackOpenEdge)
+        .def("nodesBottomBackOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackOpenEdge)
 
-    .def("nodesTopFrontOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontOpenEdge)
+        .def("nodesTopFrontOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontOpenEdge)
 
-    .def("nodesTopBackOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackOpenEdge)
+        .def("nodesTopBackOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackOpenEdge)
 
-    .def("nodesLeftBottomOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomOpenEdge)
+        .def("nodesLeftBottomOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomOpenEdge)
 
-    .def("nodesLeftFrontOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontOpenEdge)
+        .def("nodesLeftFrontOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontOpenEdge)
 
-    .def("nodesLeftBackOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackOpenEdge)
+        .def("nodesLeftBackOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackOpenEdge)
 
-    .def("nodesLeftTopOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopOpenEdge)
+        .def("nodesLeftTopOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopOpenEdge)
 
-    .def("nodesRightBottomOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomOpenEdge)
+        .def("nodesRightBottomOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomOpenEdge)
 
-    .def("nodesRightTopOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopOpenEdge)
+        .def("nodesRightTopOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopOpenEdge)
 
-    .def("nodesRightFrontOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontOpenEdge)
+        .def("nodesRightFrontOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontOpenEdge)
 
-    .def("nodesRightBackOpenEdge",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackOpenEdge)
+        .def("nodesRightBackOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackOpenEdge)
 
-    .def("nodesFrontBottomLeftCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomLeftCorner)
+        .def(
+            "nodesFrontBottomLeftCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomLeftCorner)
 
-    .def("nodesFrontBottomRightCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomRightCorner)
+        .def(
+            "nodesFrontBottomRightCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomRightCorner)
 
-    .def("nodesFrontTopLeftCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopLeftCorner)
+        .def("nodesFrontTopLeftCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopLeftCorner)
 
-    .def("nodesFrontTopRightCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopRightCorner)
+        .def("nodesFrontTopRightCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopRightCorner)
 
-    .def("nodesBackBottomLeftCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomLeftCorner)
+        .def(
+            "nodesBackBottomLeftCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomLeftCorner)
 
-    .def("nodesBackBottomRightCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomRightCorner)
+        .def(
+            "nodesBackBottomRightCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomRightCorner)
 
-    .def("nodesBackTopLeftCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopLeftCorner)
+        .def("nodesBackTopLeftCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopLeftCorner)
 
-    .def("nodesBackTopRightCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopRightCorner)
+        .def("nodesBackTopRightCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopRightCorner)
 
-    .def("nodesFrontLeftBottomCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftBottomCorner)
+        .def(
+            "nodesFrontLeftBottomCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftBottomCorner)
 
-    .def("nodesBottomFrontLeftCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontLeftCorner)
+        .def(
+            "nodesBottomFrontLeftCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontLeftCorner)
 
-    .def("nodesBottomLeftFrontCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftFrontCorner)
+        .def(
+            "nodesBottomLeftFrontCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftFrontCorner)
 
-    .def("nodesLeftFrontBottomCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontBottomCorner)
+        .def(
+            "nodesLeftFrontBottomCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontBottomCorner)
 
-    .def("nodesLeftBottomFrontCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomFrontCorner)
+        .def(
+            "nodesLeftBottomFrontCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomFrontCorner)
 
-    .def("nodesFrontRightBottomCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightBottomCorner)
+        .def(
+            "nodesFrontRightBottomCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightBottomCorner)
 
-    .def("nodesBottomFrontRightCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontRightCorner)
+        .def(
+            "nodesBottomFrontRightCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontRightCorner)
 
-    .def("nodesBottomRightFrontCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightFrontCorner)
+        .def(
+            "nodesBottomRightFrontCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightFrontCorner)
 
-    .def("nodesRightFrontBottomCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontBottomCorner)
+        .def(
+            "nodesRightFrontBottomCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontBottomCorner)
 
-    .def("nodesRightBottomFrontCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomFrontCorner)
+        .def(
+            "nodesRightBottomFrontCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomFrontCorner)
 
-    .def("nodesFrontLeftTopCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftTopCorner)
+        .def("nodesFrontLeftTopCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftTopCorner)
 
-    .def("nodesTopFrontLeftCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontLeftCorner)
+        .def("nodesTopFrontLeftCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontLeftCorner)
 
-    .def("nodesTopLeftFrontCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftFrontCorner)
+        .def("nodesTopLeftFrontCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftFrontCorner)
 
-    .def("nodesLeftFrontTopCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontTopCorner)
+        .def("nodesLeftFrontTopCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontTopCorner)
 
-    .def("nodesLeftTopFrontCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopFrontCorner)
+        .def("nodesLeftTopFrontCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopFrontCorner)
 
-    .def("nodesFrontRightTopCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightTopCorner)
+        .def("nodesFrontRightTopCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightTopCorner)
 
-    .def("nodesTopFrontRightCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontRightCorner)
+        .def("nodesTopFrontRightCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontRightCorner)
 
-    .def("nodesTopRightFrontCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightFrontCorner)
+        .def("nodesTopRightFrontCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightFrontCorner)
 
-    .def("nodesRightFrontTopCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontTopCorner)
+        .def("nodesRightFrontTopCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontTopCorner)
 
-    .def("nodesRightTopFrontCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopFrontCorner)
+        .def("nodesRightTopFrontCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopFrontCorner)
 
-    .def("nodesBackLeftBottomCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftBottomCorner)
+        .def(
+            "nodesBackLeftBottomCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftBottomCorner)
 
-    .def("nodesBottomBackLeftCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackLeftCorner)
+        .def(
+            "nodesBottomBackLeftCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackLeftCorner)
 
-    .def("nodesBottomLeftBackCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftBackCorner)
+        .def(
+            "nodesBottomLeftBackCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftBackCorner)
 
-    .def("nodesLeftBackBottomCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackBottomCorner)
+        .def(
+            "nodesLeftBackBottomCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackBottomCorner)
 
-    .def("nodesLeftBottomBackCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomBackCorner)
+        .def(
+            "nodesLeftBottomBackCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomBackCorner)
 
-    .def("nodesBackRightBottomCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightBottomCorner)
+        .def(
+            "nodesBackRightBottomCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightBottomCorner)
 
-    .def("nodesBottomBackRightCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackRightCorner)
+        .def(
+            "nodesBottomBackRightCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackRightCorner)
 
-    .def("nodesBottomRightBackCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightBackCorner)
+        .def(
+            "nodesBottomRightBackCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightBackCorner)
 
-    .def("nodesRightBackBottomCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackBottomCorner)
+        .def(
+            "nodesRightBackBottomCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackBottomCorner)
 
-    .def("nodesRightBottomBackCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomBackCorner)
+        .def(
+            "nodesRightBottomBackCorner",
+            &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomBackCorner)
 
-    .def("nodesBackLeftTopCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftTopCorner)
+        .def("nodesBackLeftTopCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftTopCorner)
 
-    .def("nodesTopBackLeftCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackLeftCorner)
+        .def("nodesTopBackLeftCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackLeftCorner)
 
-    .def("nodesTopLeftBackCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftBackCorner)
+        .def("nodesTopLeftBackCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftBackCorner)
 
-    .def("nodesLeftBackTopCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackTopCorner)
+        .def("nodesLeftBackTopCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackTopCorner)
 
-    .def("nodesLeftTopBackCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopBackCorner)
+        .def("nodesLeftTopBackCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopBackCorner)
 
-    .def("nodesBackRightTopCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightTopCorner)
+        .def("nodesBackRightTopCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightTopCorner)
 
-    .def("nodesTopBackRightCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackRightCorner)
+        .def("nodesTopBackRightCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackRightCorner)
 
-    .def("nodesTopRightBackCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightBackCorner)
+        .def("nodesTopRightBackCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightBackCorner)
 
-    .def("nodesRightBackTopCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackTopCorner)
+        .def("nodesRightBackTopCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackTopCorner)
 
-    .def("nodesRightTopBackCorner",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopBackCorner)
+        .def("nodesRightTopBackCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopBackCorner)
 
-    .def("nodesPeriodic",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesPeriodic)
+        .def("nodesPeriodic", &GooseFEM::Mesh::Hex8::FineLayer::nodesPeriodic)
 
-    .def("nodesOrigin",
-        &GooseFEM::Mesh::Hex8::FineLayer::nodesOrigin)
+        .def("nodesOrigin", &GooseFEM::Mesh::Hex8::FineLayer::nodesOrigin)
 
-    .def("dofs",
-        &GooseFEM::Mesh::Hex8::FineLayer::dofs)
+        .def("dofs", &GooseFEM::Mesh::Hex8::FineLayer::dofs)
 
-    .def("dofsPeriodic",
-        &GooseFEM::Mesh::Hex8::FineLayer::dofsPeriodic)
-
-    .def("__repr__",
-        [](const GooseFEM::Mesh::Hex8::FineLayer &){ return "<GooseFEM.Mesh.Hex8.FineLayer>"; });
+        .def("dofsPeriodic", &GooseFEM::Mesh::Hex8::FineLayer::dofsPeriodic)
 
+        .def("__repr__", [](const GooseFEM::Mesh::Hex8::FineLayer&) {
+            return "<GooseFEM.Mesh.Hex8.FineLayer>";
+        });
 }
-
-// =================================================================================================
-
diff --git a/python/MeshQuad4.hpp b/python/MeshQuad4.hpp
index 3ab9c4d..2a038c1 100644
--- a/python/MeshQuad4.hpp
+++ b/python/MeshQuad4.hpp
@@ -1,327 +1,255 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
+#include <GooseFEM/GooseFEM.h>
 #include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
-#include "../include/GooseFEM/GooseFEM.h"
-
-// =================================================================================================
 
 namespace py = pybind11;
 
-// =================================================================================================
-
 void init_MeshQuad4(py::module& m)
 {
 
-py::class_<GooseFEM::Mesh::Quad4::Regular>(m, "Regular")
-
-    .def(py::init<size_t,size_t,double>(),
-        "Regular mesh: 'nx' pixels in horizontal direction, 'ny' in vertical direction, edge size 'h'",
-        py::arg("nx"),
-        py::arg("ny"),
-        py::arg("h")=1.)
+    py::class_<GooseFEM::Mesh::Quad4::Regular>(m, "Regular")
 
-    .def("coor",
-        &GooseFEM::Mesh::Quad4::Regular::coor)
+        .def(
+            py::init<size_t, size_t, double>(),
+            "Regular mesh: 'nx' pixels in horizontal direction, 'ny' in vertical direction, edge "
+            "size 'h'",
+            py::arg("nx"),
+            py::arg("ny"),
+            py::arg("h") = 1.)
 
-    .def("conn",
-        &GooseFEM::Mesh::Quad4::Regular::conn)
+        .def("coor", &GooseFEM::Mesh::Quad4::Regular::coor)
 
-    .def("nelem",
-        &GooseFEM::Mesh::Quad4::Regular::nelem)
+        .def("conn", &GooseFEM::Mesh::Quad4::Regular::conn)
 
-    .def("nnode",
-        &GooseFEM::Mesh::Quad4::Regular::nnode)
+        .def("nelem", &GooseFEM::Mesh::Quad4::Regular::nelem)
 
-    .def("nne",
-        &GooseFEM::Mesh::Quad4::Regular::nne)
+        .def("nnode", &GooseFEM::Mesh::Quad4::Regular::nnode)
 
-    .def("ndim",
-        &GooseFEM::Mesh::Quad4::Regular::ndim)
+        .def("nne", &GooseFEM::Mesh::Quad4::Regular::nne)
 
-    .def("nelx",
-        &GooseFEM::Mesh::Quad4::Regular::nelx)
+        .def("ndim", &GooseFEM::Mesh::Quad4::Regular::ndim)
 
-    .def("nely",
-        &GooseFEM::Mesh::Quad4::Regular::nely)
+        .def("nelx", &GooseFEM::Mesh::Quad4::Regular::nelx)
 
-    .def("h",
-        &GooseFEM::Mesh::Quad4::Regular::h)
+        .def("nely", &GooseFEM::Mesh::Quad4::Regular::nely)
 
-    .def("getElementType",
-        &GooseFEM::Mesh::Quad4::Regular::getElementType)
+        .def("h", &GooseFEM::Mesh::Quad4::Regular::h)
 
-    .def("nodesBottomEdge",
-        &GooseFEM::Mesh::Quad4::Regular::nodesBottomEdge)
+        .def("getElementType", &GooseFEM::Mesh::Quad4::Regular::getElementType)
 
-    .def("nodesTopEdge",
-        &GooseFEM::Mesh::Quad4::Regular::nodesTopEdge)
+        .def("nodesBottomEdge", &GooseFEM::Mesh::Quad4::Regular::nodesBottomEdge)
 
-    .def("nodesLeftEdge",
-        &GooseFEM::Mesh::Quad4::Regular::nodesLeftEdge)
+        .def("nodesTopEdge", &GooseFEM::Mesh::Quad4::Regular::nodesTopEdge)
 
-    .def("nodesRightEdge",
-        &GooseFEM::Mesh::Quad4::Regular::nodesRightEdge)
+        .def("nodesLeftEdge", &GooseFEM::Mesh::Quad4::Regular::nodesLeftEdge)
 
-    .def("nodesBottomOpenEdge",
-        &GooseFEM::Mesh::Quad4::Regular::nodesBottomOpenEdge)
+        .def("nodesRightEdge", &GooseFEM::Mesh::Quad4::Regular::nodesRightEdge)
 
-    .def("nodesTopOpenEdge",
-        &GooseFEM::Mesh::Quad4::Regular::nodesTopOpenEdge)
+        .def("nodesBottomOpenEdge", &GooseFEM::Mesh::Quad4::Regular::nodesBottomOpenEdge)
 
-    .def("nodesLeftOpenEdge",
-        &GooseFEM::Mesh::Quad4::Regular::nodesLeftOpenEdge)
+        .def("nodesTopOpenEdge", &GooseFEM::Mesh::Quad4::Regular::nodesTopOpenEdge)
 
-    .def("nodesRightOpenEdge",
-        &GooseFEM::Mesh::Quad4::Regular::nodesRightOpenEdge)
+        .def("nodesLeftOpenEdge", &GooseFEM::Mesh::Quad4::Regular::nodesLeftOpenEdge)
 
-    .def("nodesBottomLeftCorner",
-        &GooseFEM::Mesh::Quad4::Regular::nodesBottomLeftCorner)
+        .def("nodesRightOpenEdge", &GooseFEM::Mesh::Quad4::Regular::nodesRightOpenEdge)
 
-    .def("nodesBottomRightCorner",
-        &GooseFEM::Mesh::Quad4::Regular::nodesBottomRightCorner)
+        .def("nodesBottomLeftCorner", &GooseFEM::Mesh::Quad4::Regular::nodesBottomLeftCorner)
 
-    .def("nodesTopLeftCorner",
-        &GooseFEM::Mesh::Quad4::Regular::nodesTopLeftCorner)
+        .def("nodesBottomRightCorner", &GooseFEM::Mesh::Quad4::Regular::nodesBottomRightCorner)
 
-    .def("nodesTopRightCorner",
-        &GooseFEM::Mesh::Quad4::Regular::nodesTopRightCorner)
+        .def("nodesTopLeftCorner", &GooseFEM::Mesh::Quad4::Regular::nodesTopLeftCorner)
 
-    .def("nodesLeftBottomCorner",
-        &GooseFEM::Mesh::Quad4::Regular::nodesLeftBottomCorner)
+        .def("nodesTopRightCorner", &GooseFEM::Mesh::Quad4::Regular::nodesTopRightCorner)
 
-    .def("nodesLeftTopCorner",
-        &GooseFEM::Mesh::Quad4::Regular::nodesLeftTopCorner)
+        .def("nodesLeftBottomCorner", &GooseFEM::Mesh::Quad4::Regular::nodesLeftBottomCorner)
 
-    .def("nodesRightBottomCorner",
-        &GooseFEM::Mesh::Quad4::Regular::nodesRightBottomCorner)
+        .def("nodesLeftTopCorner", &GooseFEM::Mesh::Quad4::Regular::nodesLeftTopCorner)
 
-    .def("nodesRightTopCorner",
-        &GooseFEM::Mesh::Quad4::Regular::nodesRightTopCorner)
+        .def("nodesRightBottomCorner", &GooseFEM::Mesh::Quad4::Regular::nodesRightBottomCorner)
 
-    .def("dofs",
-        &GooseFEM::Mesh::Quad4::Regular::dofs)
+        .def("nodesRightTopCorner", &GooseFEM::Mesh::Quad4::Regular::nodesRightTopCorner)
 
-    .def("nodesPeriodic",
-        &GooseFEM::Mesh::Quad4::Regular::nodesPeriodic)
+        .def("dofs", &GooseFEM::Mesh::Quad4::Regular::dofs)
 
-    .def("nodesOrigin",
-        &GooseFEM::Mesh::Quad4::Regular::nodesOrigin)
+        .def("nodesPeriodic", &GooseFEM::Mesh::Quad4::Regular::nodesPeriodic)
 
-    .def("dofsPeriodic",
-        &GooseFEM::Mesh::Quad4::Regular::dofsPeriodic)
+        .def("nodesOrigin", &GooseFEM::Mesh::Quad4::Regular::nodesOrigin)
 
-    .def("elementMatrix",
-        &GooseFEM::Mesh::Quad4::Regular::elementMatrix)
+        .def("dofsPeriodic", &GooseFEM::Mesh::Quad4::Regular::dofsPeriodic)
 
-    .def("__repr__",
-        [](const GooseFEM::Mesh::Quad4::Regular&){ return "<GooseFEM.Mesh.Quad4.Regular>"; });
+        .def("elementMatrix", &GooseFEM::Mesh::Quad4::Regular::elementMatrix)
 
-// -------------------------------------------------------------------------------------------------
+        .def("__repr__", [](const GooseFEM::Mesh::Quad4::Regular&) {
+            return "<GooseFEM.Mesh.Quad4.Regular>";
+        });
 
-py::class_<GooseFEM::Mesh::Quad4::FineLayer>(m, "FineLayer")
+    py::class_<GooseFEM::Mesh::Quad4::FineLayer>(m, "FineLayer")
 
-    .def(py::init<size_t,size_t,double,size_t>(),
-        "FineLayer mesh: 'nx' pixels in horizontal direction (length 'Lx'), idem in vertical direction",
-        py::arg("nx"),
-        py::arg("ny"),
-        py::arg("h")=1.,
-        py::arg("nfine")=1)
+        .def(
+            py::init<size_t, size_t, double, size_t>(),
+            "FineLayer mesh: 'nx' pixels in horizontal direction (length 'Lx'), idem in vertical "
+            "direction",
+            py::arg("nx"),
+            py::arg("ny"),
+            py::arg("h") = 1.,
+            py::arg("nfine") = 1)
 
-    .def("coor",
-        &GooseFEM::Mesh::Quad4::FineLayer::coor)
+        .def("coor", &GooseFEM::Mesh::Quad4::FineLayer::coor)
 
-    .def("conn",
-        &GooseFEM::Mesh::Quad4::FineLayer::conn)
+        .def("conn", &GooseFEM::Mesh::Quad4::FineLayer::conn)
 
-    .def("nelem",
-        &GooseFEM::Mesh::Quad4::FineLayer::nelem)
+        .def("nelem", &GooseFEM::Mesh::Quad4::FineLayer::nelem)
 
-    .def("nnode",
-        &GooseFEM::Mesh::Quad4::FineLayer::nnode)
+        .def("nnode", &GooseFEM::Mesh::Quad4::FineLayer::nnode)
 
-    .def("nne",
-        &GooseFEM::Mesh::Quad4::FineLayer::nne)
+        .def("nne", &GooseFEM::Mesh::Quad4::FineLayer::nne)
 
-    .def("ndim",
-        &GooseFEM::Mesh::Quad4::FineLayer::ndim)
+        .def("ndim", &GooseFEM::Mesh::Quad4::FineLayer::ndim)
 
-    .def("nelx",
-        &GooseFEM::Mesh::Quad4::FineLayer::nelx)
+        .def("nelx", &GooseFEM::Mesh::Quad4::FineLayer::nelx)
 
-    .def("nely",
-        &GooseFEM::Mesh::Quad4::FineLayer::nely)
+        .def("nely", &GooseFEM::Mesh::Quad4::FineLayer::nely)
 
-    .def("h",
-        &GooseFEM::Mesh::Quad4::FineLayer::h)
+        .def("h", &GooseFEM::Mesh::Quad4::FineLayer::h)
 
-    .def("getElementType",
-        &GooseFEM::Mesh::Quad4::FineLayer::getElementType)
+        .def("getElementType", &GooseFEM::Mesh::Quad4::FineLayer::getElementType)
 
-    .def("elementsMiddleLayer",
-        &GooseFEM::Mesh::Quad4::FineLayer::elementsMiddleLayer)
+        .def("elementsMiddleLayer", &GooseFEM::Mesh::Quad4::FineLayer::elementsMiddleLayer)
 
-    .def("nodesBottomEdge",
-        &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomEdge)
+        .def("nodesBottomEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomEdge)
 
-    .def("nodesTopEdge",
-        &GooseFEM::Mesh::Quad4::FineLayer::nodesTopEdge)
+        .def("nodesTopEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesTopEdge)
 
-    .def("nodesLeftEdge",
-        &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftEdge)
+        .def("nodesLeftEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftEdge)
 
-    .def("nodesRightEdge",
-        &GooseFEM::Mesh::Quad4::FineLayer::nodesRightEdge)
+        .def("nodesRightEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesRightEdge)
 
-    .def("nodesBottomOpenEdge",
-        &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomOpenEdge)
+        .def("nodesBottomOpenEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomOpenEdge)
 
-    .def("nodesTopOpenEdge",
-        &GooseFEM::Mesh::Quad4::FineLayer::nodesTopOpenEdge)
+        .def("nodesTopOpenEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesTopOpenEdge)
 
-    .def("nodesLeftOpenEdge",
-        &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftOpenEdge)
+        .def("nodesLeftOpenEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftOpenEdge)
 
-    .def("nodesRightOpenEdge",
-        &GooseFEM::Mesh::Quad4::FineLayer::nodesRightOpenEdge)
+        .def("nodesRightOpenEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesRightOpenEdge)
 
-    .def("nodesBottomLeftCorner",
-        &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomLeftCorner)
+        .def("nodesBottomLeftCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomLeftCorner)
 
-    .def("nodesBottomRightCorner",
-        &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomRightCorner)
+        .def("nodesBottomRightCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomRightCorner)
 
-    .def("nodesTopLeftCorner",
-        &GooseFEM::Mesh::Quad4::FineLayer::nodesTopLeftCorner)
+        .def("nodesTopLeftCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesTopLeftCorner)
 
-    .def("nodesTopRightCorner",
-        &GooseFEM::Mesh::Quad4::FineLayer::nodesTopRightCorner)
+        .def("nodesTopRightCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesTopRightCorner)
 
-    .def("nodesLeftBottomCorner",
-        &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftBottomCorner)
+        .def("nodesLeftBottomCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftBottomCorner)
 
-    .def("nodesLeftTopCorner",
-        &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftTopCorner)
+        .def("nodesLeftTopCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftTopCorner)
 
-    .def("nodesRightBottomCorner",
-        &GooseFEM::Mesh::Quad4::FineLayer::nodesRightBottomCorner)
+        .def("nodesRightBottomCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesRightBottomCorner)
 
-    .def("nodesRightTopCorner",
-        &GooseFEM::Mesh::Quad4::FineLayer::nodesRightTopCorner)
+        .def("nodesRightTopCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesRightTopCorner)
 
-    .def("dofs",
-        &GooseFEM::Mesh::Quad4::FineLayer::dofs)
+        .def("dofs", &GooseFEM::Mesh::Quad4::FineLayer::dofs)
 
-    .def("nodesPeriodic",
-        &GooseFEM::Mesh::Quad4::FineLayer::nodesPeriodic)
+        .def("nodesPeriodic", &GooseFEM::Mesh::Quad4::FineLayer::nodesPeriodic)
 
-    .def("nodesOrigin",
-        &GooseFEM::Mesh::Quad4::FineLayer::nodesOrigin)
+        .def("nodesOrigin", &GooseFEM::Mesh::Quad4::FineLayer::nodesOrigin)
 
-    .def("dofsPeriodic",
-        &GooseFEM::Mesh::Quad4::FineLayer::dofsPeriodic)
-
-    .def("__repr__",
-        [](const GooseFEM::Mesh::Quad4::FineLayer&){ return "<GooseFEM.Mesh.Quad4.FineLayer>"; }
-    );
+        .def("dofsPeriodic", &GooseFEM::Mesh::Quad4::FineLayer::dofsPeriodic)
 
+        .def("__repr__", [](const GooseFEM::Mesh::Quad4::FineLayer&) {
+            return "<GooseFEM.Mesh.Quad4.FineLayer>";
+        });
 }
 
-// =================================================================================================
-
 void init_MeshQuad4Map(py::module& m)
 {
 
-// -------------------------------------------------------------------------------------------------
-
-py::class_<GooseFEM::Mesh::Quad4::Map::RefineRegular>(m, "RefineRegular")
+    py::class_<GooseFEM::Mesh::Quad4::Map::RefineRegular>(m, "RefineRegular")
 
-    .def(py::init<const GooseFEM::Mesh::Quad4::Regular &, size_t, size_t>(),
-        "Refine a regular mesh",
-        py::arg("mesh"),
-        py::arg("nx"),
-        py::arg("ny"))
+        .def(
+            py::init<const GooseFEM::Mesh::Quad4::Regular&, size_t, size_t>(),
+            "Refine a regular mesh",
+            py::arg("mesh"),
+            py::arg("nx"),
+            py::arg("ny"))
 
-    .def("getCoarseMesh",
-        &GooseFEM::Mesh::Quad4::Map::RefineRegular::getCoarseMesh)
+        .def("getCoarseMesh", &GooseFEM::Mesh::Quad4::Map::RefineRegular::getCoarseMesh)
 
-    .def("getFineMesh",
-        &GooseFEM::Mesh::Quad4::Map::RefineRegular::getFineMesh)
+        .def("getFineMesh", &GooseFEM::Mesh::Quad4::Map::RefineRegular::getFineMesh)
 
-    .def("getMap",
-        &GooseFEM::Mesh::Quad4::Map::RefineRegular::getMap)
+        .def("getMap", &GooseFEM::Mesh::Quad4::Map::RefineRegular::getMap)
 
-    .def("mapToCoarse",
-        py::overload_cast<const xt::xtensor<double,1>&>(
-            &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToCoarse, py::const_))
+        .def(
+            "mapToCoarse",
+            py::overload_cast<const xt::xtensor<double, 1>&>(
+                &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToCoarse, py::const_))
 
-    .def("mapToCoarse",
-        py::overload_cast<const xt::xtensor<double,2>&>(
-            &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToCoarse, py::const_))
+        .def(
+            "mapToCoarse",
+            py::overload_cast<const xt::xtensor<double, 2>&>(
+                &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToCoarse, py::const_))
 
-    .def("mapToCoarse",
-        py::overload_cast<const xt::xtensor<double,4>&>(
-            &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToCoarse, py::const_))
+        .def(
+            "mapToCoarse",
+            py::overload_cast<const xt::xtensor<double, 4>&>(
+                &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToCoarse, py::const_))
 
-    .def("mapToFine",
-        py::overload_cast<const xt::xtensor<double,1>&>(
-            &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine, py::const_))
+        .def(
+            "mapToFine",
+            py::overload_cast<const xt::xtensor<double, 1>&>(
+                &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine, py::const_))
 
-    .def("mapToFine",
-        py::overload_cast<const xt::xtensor<double,2>&>(
-            &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine, py::const_))
+        .def(
+            "mapToFine",
+            py::overload_cast<const xt::xtensor<double, 2>&>(
+                &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine, py::const_))
 
-    .def("mapToFine",
-        py::overload_cast<const xt::xtensor<double,4>&>(
-            &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine, py::const_))
+        .def(
+            "mapToFine",
+            py::overload_cast<const xt::xtensor<double, 4>&>(
+                &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine, py::const_))
 
-    .def("__repr__",
-        [](const GooseFEM::Mesh::Quad4::Map::RefineRegular &){
-            return "<GooseFEM.Mesh.Quad4.Map.RefineRegular>"; });
+        .def("__repr__", [](const GooseFEM::Mesh::Quad4::Map::RefineRegular&) {
+            return "<GooseFEM.Mesh.Quad4.Map.RefineRegular>";
+        });
 
-// -------------------------------------------------------------------------------------------------
+    py::class_<GooseFEM::Mesh::Quad4::Map::FineLayer2Regular>(m, "FineLayer2Regular")
 
-py::class_<GooseFEM::Mesh::Quad4::Map::FineLayer2Regular>(m, "FineLayer2Regular")
+        .def(
+            py::init<const GooseFEM::Mesh::Quad4::FineLayer&>(),
+            "Map a FineLayer mesh to a Regular mesh",
+            py::arg("mesh"))
 
-    .def(py::init<const GooseFEM::Mesh::Quad4::FineLayer &>(),
-        "Map a FineLayer mesh to a Regular mesh",
-        py::arg("mesh"))
+        .def("getRegularMesh", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getRegularMesh)
 
-    .def("getRegularMesh",
-        &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getRegularMesh)
+        .def("getFineLayerMesh", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getFineLayerMesh)
 
-    .def("getFineLayerMesh",
-        &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getFineLayerMesh)
+        .def("getMap", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getMap)
 
-    .def("getMap",
-        &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getMap)
+        .def("getMapFraction", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getMapFraction)
 
-    .def("getMapFraction",
-        &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getMapFraction)
+        .def(
+            "mapToRegular",
+            py::overload_cast<const xt::xtensor<double, 1>&>(
+                &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular, py::const_))
 
-    .def("mapToRegular",
-        py::overload_cast<const xt::xtensor<double,1>&>(
-            &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular, py::const_))
+        .def(
+            "mapToRegular",
+            py::overload_cast<const xt::xtensor<double, 2>&>(
+                &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular, py::const_))
 
-    .def("mapToRegular",
-        py::overload_cast<const xt::xtensor<double,2>&>(
-            &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular, py::const_))
-
-    .def("mapToRegular",
-        py::overload_cast<const xt::xtensor<double,4>&>(
-            &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular, py::const_))
-
-    .def("__repr__",
-        [](const GooseFEM::Mesh::Quad4::Map::FineLayer2Regular &){
-            return "<GooseFEM.Mesh.Quad4.Map.FineLayer2Regular>"; });
-
-// -------------------------------------------------------------------------------------------------
+        .def(
+            "mapToRegular",
+            py::overload_cast<const xt::xtensor<double, 4>&>(
+                &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular, py::const_))
 
+        .def("__repr__", [](const GooseFEM::Mesh::Quad4::Map::FineLayer2Regular&) {
+            return "<GooseFEM.Mesh.Quad4.Map.FineLayer2Regular>";
+        });
 }
-
-// =================================================================================================
-
diff --git a/python/MeshTri3.hpp b/python/MeshTri3.hpp
index b77fe6f..79a33f1 100644
--- a/python/MeshTri3.hpp
+++ b/python/MeshTri3.hpp
@@ -1,145 +1,112 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
+#include <GooseFEM/GooseFEM.h>
 #include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
-#include "../include/GooseFEM/GooseFEM.h"
-
-// =================================================================================================
 
 namespace py = pybind11;
 
-// =================================================================================================
-
 void init_MeshTri3(py::module& m)
 {
 
-py::class_<GooseFEM::Mesh::Tri3::Regular>(m, "Regular")
-
-    .def(
-        py::init<size_t,size_t,double>(),
-        "Regular mesh: 'nx' pixels in horizontal direction, 'ny' in vertical direction, edge size 'h'",
-        py::arg("nx"),
-        py::arg("ny"),
-        py::arg("h")=1.)
+    py::class_<GooseFEM::Mesh::Tri3::Regular>(m, "Regular")
 
-    .def("coor",
-        &GooseFEM::Mesh::Tri3::Regular::coor)
+        .def(
+            py::init<size_t, size_t, double>(),
+            "Regular mesh: 'nx' pixels in horizontal direction, 'ny' in vertical direction, edge "
+            "size 'h'",
+            py::arg("nx"),
+            py::arg("ny"),
+            py::arg("h") = 1.)
 
-    .def("conn",
-        &GooseFEM::Mesh::Tri3::Regular::conn)
+        .def("coor", &GooseFEM::Mesh::Tri3::Regular::coor)
 
-    .def("nelem",
-        &GooseFEM::Mesh::Tri3::Regular::nelem)
+        .def("conn", &GooseFEM::Mesh::Tri3::Regular::conn)
 
-    .def("nnode",
-        &GooseFEM::Mesh::Tri3::Regular::nnode)
+        .def("nelem", &GooseFEM::Mesh::Tri3::Regular::nelem)
 
-    .def("nne",
-        &GooseFEM::Mesh::Tri3::Regular::nne)
+        .def("nnode", &GooseFEM::Mesh::Tri3::Regular::nnode)
 
-    .def("ndim",
-        &GooseFEM::Mesh::Tri3::Regular::ndim)
+        .def("nne", &GooseFEM::Mesh::Tri3::Regular::nne)
 
-    .def("getElementType",
-        &GooseFEM::Mesh::Tri3::Regular::getElementType)
+        .def("ndim", &GooseFEM::Mesh::Tri3::Regular::ndim)
 
-    .def("nodesBottomEdge",
-        &GooseFEM::Mesh::Tri3::Regular::nodesBottomEdge)
+        .def("getElementType", &GooseFEM::Mesh::Tri3::Regular::getElementType)
 
-    .def("nodesTopEdge",
-        &GooseFEM::Mesh::Tri3::Regular::nodesTopEdge)
+        .def("nodesBottomEdge", &GooseFEM::Mesh::Tri3::Regular::nodesBottomEdge)
 
-    .def("nodesLeftEdge",
-        &GooseFEM::Mesh::Tri3::Regular::nodesLeftEdge)
+        .def("nodesTopEdge", &GooseFEM::Mesh::Tri3::Regular::nodesTopEdge)
 
-    .def("nodesRightEdge",
-        &GooseFEM::Mesh::Tri3::Regular::nodesRightEdge)
+        .def("nodesLeftEdge", &GooseFEM::Mesh::Tri3::Regular::nodesLeftEdge)
 
-    .def("nodesBottomOpenEdge",
-        &GooseFEM::Mesh::Tri3::Regular::nodesBottomOpenEdge)
+        .def("nodesRightEdge", &GooseFEM::Mesh::Tri3::Regular::nodesRightEdge)
 
-    .def("nodesTopOpenEdge",
-        &GooseFEM::Mesh::Tri3::Regular::nodesTopOpenEdge)
+        .def("nodesBottomOpenEdge", &GooseFEM::Mesh::Tri3::Regular::nodesBottomOpenEdge)
 
-    .def("nodesLeftOpenEdge",
-        &GooseFEM::Mesh::Tri3::Regular::nodesLeftOpenEdge)
+        .def("nodesTopOpenEdge", &GooseFEM::Mesh::Tri3::Regular::nodesTopOpenEdge)
 
-    .def("nodesRightOpenEdge",
-        &GooseFEM::Mesh::Tri3::Regular::nodesRightOpenEdge)
+        .def("nodesLeftOpenEdge", &GooseFEM::Mesh::Tri3::Regular::nodesLeftOpenEdge)
 
-    .def("nodesBottomLeftCorner",
-        &GooseFEM::Mesh::Tri3::Regular::nodesBottomLeftCorner)
+        .def("nodesRightOpenEdge", &GooseFEM::Mesh::Tri3::Regular::nodesRightOpenEdge)
 
-    .def("nodesBottomRightCorner",
-        &GooseFEM::Mesh::Tri3::Regular::nodesBottomRightCorner)
+        .def("nodesBottomLeftCorner", &GooseFEM::Mesh::Tri3::Regular::nodesBottomLeftCorner)
 
-    .def("nodesTopLeftCorner",
-        &GooseFEM::Mesh::Tri3::Regular::nodesTopLeftCorner)
+        .def("nodesBottomRightCorner", &GooseFEM::Mesh::Tri3::Regular::nodesBottomRightCorner)
 
-    .def("nodesTopRightCorner",
-        &GooseFEM::Mesh::Tri3::Regular::nodesTopRightCorner)
+        .def("nodesTopLeftCorner", &GooseFEM::Mesh::Tri3::Regular::nodesTopLeftCorner)
 
-    .def("nodesLeftBottomCorner",
-        &GooseFEM::Mesh::Tri3::Regular::nodesLeftBottomCorner)
+        .def("nodesTopRightCorner", &GooseFEM::Mesh::Tri3::Regular::nodesTopRightCorner)
 
-    .def("nodesLeftTopCorner",
-        &GooseFEM::Mesh::Tri3::Regular::nodesLeftTopCorner)
+        .def("nodesLeftBottomCorner", &GooseFEM::Mesh::Tri3::Regular::nodesLeftBottomCorner)
 
-    .def("nodesRightBottomCorner",
-        &GooseFEM::Mesh::Tri3::Regular::nodesRightBottomCorner)
+        .def("nodesLeftTopCorner", &GooseFEM::Mesh::Tri3::Regular::nodesLeftTopCorner)
 
-    .def("nodesRightTopCorner",
-        &GooseFEM::Mesh::Tri3::Regular::nodesRightTopCorner)
+        .def("nodesRightBottomCorner", &GooseFEM::Mesh::Tri3::Regular::nodesRightBottomCorner)
 
-    .def("nodesPeriodic",
-        &GooseFEM::Mesh::Tri3::Regular::nodesPeriodic)
+        .def("nodesRightTopCorner", &GooseFEM::Mesh::Tri3::Regular::nodesRightTopCorner)
 
-    .def("nodesOrigin",
-        &GooseFEM::Mesh::Tri3::Regular::nodesOrigin)
+        .def("nodesPeriodic", &GooseFEM::Mesh::Tri3::Regular::nodesPeriodic)
 
-    .def("dofs",
-        &GooseFEM::Mesh::Tri3::Regular::dofs)
+        .def("nodesOrigin", &GooseFEM::Mesh::Tri3::Regular::nodesOrigin)
 
-    .def("dofsPeriodic",
-        &GooseFEM::Mesh::Tri3::Regular::dofsPeriodic)
+        .def("dofs", &GooseFEM::Mesh::Tri3::Regular::dofs)
 
-    .def("__repr__",
-        [](const GooseFEM::Mesh::Tri3::Regular&){ return "<GooseFEM.Mesh.Tri3.Regular>"; });
+        .def("dofsPeriodic", &GooseFEM::Mesh::Tri3::Regular::dofsPeriodic)
 
-// -------------------------------------------------------------------------------------------------
+        .def("__repr__", [](const GooseFEM::Mesh::Tri3::Regular&) {
+            return "<GooseFEM.Mesh.Tri3.Regular>";
+        });
 
-m.def("getOrientation",
-    &GooseFEM::Mesh::Tri3::getOrientation,
-    "Get the orientation of each element",
-    py::arg("coor"),
-    py::arg("conn"));
+    m.def(
+        "getOrientation",
+        &GooseFEM::Mesh::Tri3::getOrientation,
+        "Get the orientation of each element",
+        py::arg("coor"),
+        py::arg("conn"));
 
-m.def("setOrientation",
-    py::overload_cast<
-        const xt::xtensor<double,2>&,
-        const xt::xtensor<size_t,2>&,
-        int>(&GooseFEM::Mesh::Tri3::setOrientation),
-    "Set the orientation of each element",
-    py::arg("coor"),
-    py::arg("conn"),
-    py::arg("orientation"));
+    m.def(
+        "setOrientation",
+        py::overload_cast<const xt::xtensor<double, 2>&, const xt::xtensor<size_t, 2>&, int>(
+            &GooseFEM::Mesh::Tri3::setOrientation),
+        "Set the orientation of each element",
+        py::arg("coor"),
+        py::arg("conn"),
+        py::arg("orientation"));
 
-m.def("setOrientation",
-    py::overload_cast<
-        const xt::xtensor<double,2>&,
-        const xt::xtensor<size_t,2>&,
-        const xt::xtensor<int,1>&,
-        int>(&GooseFEM::Mesh::Tri3::setOrientation),
-    "Set the orientation of each element",
-    py::arg("coor"),
-    py::arg("conn"),
-    py::arg("val"),
-    py::arg("orientation"));
+    m.def(
+        "setOrientation",
+        py::overload_cast<
+            const xt::xtensor<double, 2>&,
+            const xt::xtensor<size_t, 2>&,
+            const xt::xtensor<int, 1>&,
+            int>(&GooseFEM::Mesh::Tri3::setOrientation),
+        "Set the orientation of each element",
+        py::arg("coor"),
+        py::arg("conn"),
+        py::arg("val"),
+        py::arg("orientation"));
 }
-
-// =================================================================================================
-
diff --git a/python/ParaView.hpp b/python/ParaView.hpp
index 345e7f7..f3c948c 100644
--- a/python/ParaView.hpp
+++ b/python/ParaView.hpp
@@ -1,220 +1,191 @@
 /* =================================================================================================
 
 (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/GooseFEM.h"
 #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)
+    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("xdmf",
-        &GooseFEM::ParaView::HDF5::TimeSeries::xdmf,
-        py::arg("indent")=4)
+        .def(
+            py::init<
+                const GooseFEM::ParaView::HDF5::Connectivity&,
+                const GooseFEM::ParaView::HDF5::Coordinates&>(),
+            "Increment",
+            py::arg("conn"),
+            py::arg("coor"))
 
-    .def("write",
-        &GooseFEM::ParaView::HDF5::TimeSeries::write,
-        py::arg("fname"),
-        py::arg("indent")=4)
+        .def(
+            "push_back",
+            py::overload_cast<const GooseFEM::ParaView::HDF5::Connectivity&>(
+                &GooseFEM::ParaView::HDF5::Increment::push_back))
 
-    .def("__repr__",
-        [](const GooseFEM::ParaView::HDF5::TimeSeries &){
-            return "<GooseFEM.ParaView.HDF5.TimeSeries>"; });
+        .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/Vector.hpp b/python/Vector.hpp
index 8b95be7..a2e09fc 100644
--- a/python/Vector.hpp
+++ b/python/Vector.hpp
@@ -1,102 +1,96 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
+#include <GooseFEM/GooseFEM.h>
 #include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
-#include "../include/GooseFEM/GooseFEM.h"
-
-// =================================================================================================
 
 namespace py = pybind11;
 
-// =================================================================================================
-
 void init_Vector(py::module& m)
 {
 
-py::class_<GooseFEM::Vector>(m, "Vector")
-
-    .def(py::init<const xt::xtensor<size_t,2> &, const xt::xtensor<size_t,2> &>(),
-        "Switch between dofval/nodevec/elemvec",
-        py::arg("conn"),
-        py::arg("dofs"))
-
-    .def("nelem",
-        &GooseFEM::Vector::nelem,
-        "Number of element")
-
-    .def("nne",
-        &GooseFEM::Vector::nne,
-        "Number of nodes per element")
-
-    .def("nnode",
-        &GooseFEM::Vector::nnode,
-        "Number of nodes")
-
-    .def("ndim",
-        &GooseFEM::Vector::ndim,
-        "Number of dimensions")
-
-    .def("ndof",
-        &GooseFEM::Vector::ndof,
-        "Number of degrees-of-freedom")
-
-    .def("dofs",
-        &GooseFEM::Vector::dofs,
-        "Return degrees-of-freedom")
-
-    .def("AsDofs",
-        py::overload_cast<const xt::xtensor<double,2>&>(&GooseFEM::Vector::AsDofs, py::const_),
-        "Set 'dofval",
-        py::arg("nodevec"))
-
-    .def("AsDofs",
-        py::overload_cast<const xt::xtensor<double,3>&>(&GooseFEM::Vector::AsDofs, py::const_),
-        "Set 'dofval",
-        py::arg("elemvec"))
-
-    .def("AsNode",
-        py::overload_cast<const xt::xtensor<double,1>&>(&GooseFEM::Vector::AsNode, py::const_),
-        "Set 'nodevec",
-        py::arg("dofval"))
-
-    .def("AsNode",
-        py::overload_cast<const xt::xtensor<double,3>&>(&GooseFEM::Vector::AsNode, py::const_),
-        "Set 'nodevec",
-        py::arg("elemvec"))
-
-    .def("AsElement",
-        py::overload_cast<const xt::xtensor<double,1>&>(&GooseFEM::Vector::AsElement, py::const_),
-        "Set 'elemvec",
-        py::arg("dofval"))
-
-    .def("AsElement",
-        py::overload_cast<const xt::xtensor<double,2>&>(&GooseFEM::Vector::AsElement, py::const_),
-        "Set 'elemvec",
-        py::arg("nodevec"))
-
-    .def("AssembleDofs",
-        py::overload_cast<const xt::xtensor<double,2>&>(&GooseFEM::Vector::AssembleDofs, py::const_),
-        "Assemble 'dofval'",
-        py::arg("nodevec"))
-
-    .def("AssembleDofs",
-        py::overload_cast<const xt::xtensor<double,3>&>(&GooseFEM::Vector::AssembleDofs, py::const_),
-        "Assemble 'dofval'",
-        py::arg("elemvec"))
-
-    .def("AssembleNode",
-        py::overload_cast<const xt::xtensor<double,3>&>(&GooseFEM::Vector::AssembleNode, py::const_),
-        "Assemble 'nodevec'",
-        py::arg("elemvec"))
-
-    .def("__repr__",
-        [](const GooseFEM::Vector&){ return "<GooseFEM.Vector>"; });
-
+    py::class_<GooseFEM::Vector>(m, "Vector")
+
+        .def(
+            py::init<const xt::xtensor<size_t, 2>&, const xt::xtensor<size_t, 2>&>(),
+            "Switch between dofval/nodevec/elemvec",
+            py::arg("conn"),
+            py::arg("dofs"))
+
+        .def("nelem", &GooseFEM::Vector::nelem, "Number of element")
+
+        .def("nne", &GooseFEM::Vector::nne, "Number of nodes per element")
+
+        .def("nnode", &GooseFEM::Vector::nnode, "Number of nodes")
+
+        .def("ndim", &GooseFEM::Vector::ndim, "Number of dimensions")
+
+        .def("ndof", &GooseFEM::Vector::ndof, "Number of degrees-of-freedom")
+
+        .def("dofs", &GooseFEM::Vector::dofs, "Return degrees-of-freedom")
+
+        .def(
+            "AsDofs",
+            py::overload_cast<const xt::xtensor<double, 2>&>(&GooseFEM::Vector::AsDofs, py::const_),
+            "Set 'dofval",
+            py::arg("nodevec"))
+
+        .def(
+            "AsDofs",
+            py::overload_cast<const xt::xtensor<double, 3>&>(&GooseFEM::Vector::AsDofs, py::const_),
+            "Set 'dofval",
+            py::arg("elemvec"))
+
+        .def(
+            "AsNode",
+            py::overload_cast<const xt::xtensor<double, 1>&>(&GooseFEM::Vector::AsNode, py::const_),
+            "Set 'nodevec",
+            py::arg("dofval"))
+
+        .def(
+            "AsNode",
+            py::overload_cast<const xt::xtensor<double, 3>&>(&GooseFEM::Vector::AsNode, py::const_),
+            "Set 'nodevec",
+            py::arg("elemvec"))
+
+        .def(
+            "AsElement",
+            py::overload_cast<const xt::xtensor<double, 1>&>(
+                &GooseFEM::Vector::AsElement, py::const_),
+            "Set 'elemvec",
+            py::arg("dofval"))
+
+        .def(
+            "AsElement",
+            py::overload_cast<const xt::xtensor<double, 2>&>(
+                &GooseFEM::Vector::AsElement, py::const_),
+            "Set 'elemvec",
+            py::arg("nodevec"))
+
+        .def(
+            "AssembleDofs",
+            py::overload_cast<const xt::xtensor<double, 2>&>(
+                &GooseFEM::Vector::AssembleDofs, py::const_),
+            "Assemble 'dofval'",
+            py::arg("nodevec"))
+
+        .def(
+            "AssembleDofs",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::Vector::AssembleDofs, py::const_),
+            "Assemble 'dofval'",
+            py::arg("elemvec"))
+
+        .def(
+            "AssembleNode",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::Vector::AssembleNode, py::const_),
+            "Assemble 'nodevec'",
+            py::arg("elemvec"))
+
+        .def("__repr__", [](const GooseFEM::Vector&) { return "<GooseFEM.Vector>"; });
 }
-
-// =================================================================================================
-
diff --git a/python/VectorPartitioned.hpp b/python/VectorPartitioned.hpp
index 3792030..c94c288 100644
--- a/python/VectorPartitioned.hpp
+++ b/python/VectorPartitioned.hpp
@@ -1,214 +1,200 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
+#include <GooseFEM/GooseFEM.h>
 #include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
-#include "../include/GooseFEM/GooseFEM.h"
-
-// =================================================================================================
 
 namespace py = pybind11;
 
-// =================================================================================================
-
 void init_VectorPartitioned(py::module& m)
 {
 
-py::class_<GooseFEM::VectorPartitioned>(m, "VectorPartitioned")
-
-    .def(py::init<
-            const xt::xtensor<size_t,2>&,
-            const xt::xtensor<size_t,2>&,
-            const xt::xtensor<size_t,1>&>(),
-        "Switch between dofval/nodevec/elemvec",
-        py::arg("conn"),
-        py::arg("dofs"),
-        py::arg("iip"))
-
-    .def("nelem",
-        &GooseFEM::VectorPartitioned::nelem,
-        "Number of element")
-
-    .def("nne",
-        &GooseFEM::VectorPartitioned::nne,
-        "Number of nodes per element")
-
-    .def("nnode",
-        &GooseFEM::VectorPartitioned::nnode,
-        "Number of nodes")
-
-    .def("ndim",
-        &GooseFEM::VectorPartitioned::ndim,
-        "Number of dimensions")
-
-    .def("ndof",
-        &GooseFEM::VectorPartitioned::ndof,
-        "Number of degrees-of-freedom")
-
-    .def("nnu",
-        &GooseFEM::VectorPartitioned::nnu,
-        "Number of unknown degrees-of-freedom")
-
-    .def("nnp",
-        &GooseFEM::VectorPartitioned::nnp,
-        "Number of prescribed degrees-of-freedom")
-
-    .def("dofs",
-        &GooseFEM::VectorPartitioned::dofs,
-        "Return degrees-of-freedom")
-
-    .def("iiu",
-        &GooseFEM::VectorPartitioned::iiu,
-        "Return unknown degrees-of-freedom")
-
-    .def("iip",
-        &GooseFEM::VectorPartitioned::iip,
-        "Return prescribed degrees-of-freedom")
-
-
-    .def("AsDofs",
-        py::overload_cast<const xt::xtensor<double,1>&,const xt::xtensor<double,1>&>(
-            &GooseFEM::VectorPartitioned::AsDofs, py::const_),
-        "Set 'dofval",
-        py::arg("dofval_u"),
-        py::arg("dofval_p"))
-
-    .def("AsDofs",
-        py::overload_cast<const xt::xtensor<double,2>&>(
-            &GooseFEM::VectorPartitioned::AsDofs, py::const_),
-        "Set 'dofval",
-        py::arg("nodevec"))
-
-    .def("AsDofs",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::VectorPartitioned::AsDofs, py::const_),
-        "Set 'dofval",
-        py::arg("elemvec"))
-
-
-    .def("AsDofs_u",
-        py::overload_cast<const xt::xtensor<double,2>&>(
-            &GooseFEM::VectorPartitioned::AsDofs_u, py::const_),
-        "Set 'dofval",
-        py::arg("nodevec"))
-
-    .def("AsDofs_u",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::VectorPartitioned::AsDofs_u, py::const_),
-        "Set 'dofval",
-        py::arg("elemvec"))
-
-    .def("AsDofs_p",
-        py::overload_cast<const xt::xtensor<double,2>&>(
-            &GooseFEM::VectorPartitioned::AsDofs_p, py::const_),
-        "Set 'dofval",
-        py::arg("nodevec"))
-
-    .def("AsDofs_p",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::VectorPartitioned::AsDofs_p, py::const_),
-        "Set 'dofval",
-        py::arg("elemvec"))
-
-    .def("AsNode",
-        py::overload_cast<const xt::xtensor<double,1>&,const xt::xtensor<double,1>&>(
-            &GooseFEM::VectorPartitioned::AsNode, py::const_),
-        "Set 'nodevec",
-        py::arg("dofval_u"),
-        py::arg("dofval_p"))
-
-    .def("AsNode",
-        py::overload_cast<const xt::xtensor<double,1>&>(
-            &GooseFEM::VectorPartitioned::AsNode, py::const_),
-        "Set 'nodevec",
-        py::arg("dofval"))
-
-    .def("AsNode",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::VectorPartitioned::AsNode, py::const_),
-        "Set 'nodevec",
-        py::arg("elemvec"))
-
-    .def("AsElement",
-        py::overload_cast<const xt::xtensor<double,1>&,const xt::xtensor<double,1>&>(
-            &GooseFEM::VectorPartitioned::AsElement, py::const_),
-        "Set 'elemvec",
-        py::arg("dofval_u"),
-        py::arg("dofval_p"))
-
-    .def("AsElement",
-        py::overload_cast<const xt::xtensor<double,1>&>(
-            &GooseFEM::VectorPartitioned::AsElement, py::const_),
-        "Set 'elemvec",
-        py::arg("dofval"))
-
-    .def("AsElement",
-        py::overload_cast<const xt::xtensor<double,2>&>(
-            &GooseFEM::VectorPartitioned::AsElement, py::const_),
-        "Set 'elemvec",
-        py::arg("nodevec"))
-
-    .def("AssembleDofs",
-        py::overload_cast<const xt::xtensor<double,2>&>(
-            &GooseFEM::VectorPartitioned::AssembleDofs, py::const_),
-        "Assemble 'dofval'",
-        py::arg("nodevec"))
-
-    .def("AssembleDofs",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::VectorPartitioned::AssembleDofs, py::const_),
-        "Assemble 'dofval'",
-        py::arg("elemvec"))
-
-    .def("AssembleDofs_u",
-        py::overload_cast<const xt::xtensor<double,2>&>(
-            &GooseFEM::VectorPartitioned::AssembleDofs_u, py::const_),
-        "Assemble 'dofval'",
-        py::arg("nodevec"))
-
-    .def("AssembleDofs_u",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::VectorPartitioned::AssembleDofs_u, py::const_),
-        "Assemble 'dofval'",
-        py::arg("elemvec"))
-
-    .def("AssembleDofs_p",
-        py::overload_cast<const xt::xtensor<double,2>&>(
-            &GooseFEM::VectorPartitioned::AssembleDofs_p, py::const_),
-        "Assemble 'dofval'",
-        py::arg("nodevec"))
-
-    .def("AssembleDofs_p",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::VectorPartitioned::AssembleDofs_p, py::const_),
-        "Assemble 'dofval'",
-        py::arg("elemvec"))
-
-    .def("AssembleNode",
-        py::overload_cast<const xt::xtensor<double,3>&>(
-            &GooseFEM::VectorPartitioned::AssembleNode, py::const_),
-        "Assemble 'nodevec'",
-        py::arg("elemvec"))
-
-    .def("Copy",
-        &GooseFEM::VectorPartitioned::Copy,
-        "Copy")
-
-    .def("Copy_u",
-        &GooseFEM::VectorPartitioned::Copy_u,
-        "Copy iiu")
-
-    .def("Copy_p",
-        &GooseFEM::VectorPartitioned::Copy_p,
-        "Copy iip")
-
-    .def("__repr__", [](
-        const GooseFEM::VectorPartitioned&){ return "<GooseFEM.VectorPartitioned>"; });
-
+    py::class_<GooseFEM::VectorPartitioned>(m, "VectorPartitioned")
+
+        .def(
+            py::init<
+                const xt::xtensor<size_t, 2>&,
+                const xt::xtensor<size_t, 2>&,
+                const xt::xtensor<size_t, 1>&>(),
+            "Switch between dofval/nodevec/elemvec",
+            py::arg("conn"),
+            py::arg("dofs"),
+            py::arg("iip"))
+
+        .def("nelem", &GooseFEM::VectorPartitioned::nelem, "Number of element")
+
+        .def("nne", &GooseFEM::VectorPartitioned::nne, "Number of nodes per element")
+
+        .def("nnode", &GooseFEM::VectorPartitioned::nnode, "Number of nodes")
+
+        .def("ndim", &GooseFEM::VectorPartitioned::ndim, "Number of dimensions")
+
+        .def("ndof", &GooseFEM::VectorPartitioned::ndof, "Number of degrees-of-freedom")
+
+        .def("nnu", &GooseFEM::VectorPartitioned::nnu, "Number of unknown degrees-of-freedom")
+
+        .def("nnp", &GooseFEM::VectorPartitioned::nnp, "Number of prescribed degrees-of-freedom")
+
+        .def("dofs", &GooseFEM::VectorPartitioned::dofs, "Return degrees-of-freedom")
+
+        .def("iiu", &GooseFEM::VectorPartitioned::iiu, "Return unknown degrees-of-freedom")
+
+        .def("iip", &GooseFEM::VectorPartitioned::iip, "Return prescribed degrees-of-freedom")
+
+        .def(
+            "AsDofs",
+            py::overload_cast<const xt::xtensor<double, 1>&, const xt::xtensor<double, 1>&>(
+                &GooseFEM::VectorPartitioned::AsDofs, py::const_),
+            "Set 'dofval",
+            py::arg("dofval_u"),
+            py::arg("dofval_p"))
+
+        .def(
+            "AsDofs",
+            py::overload_cast<const xt::xtensor<double, 2>&>(
+                &GooseFEM::VectorPartitioned::AsDofs, py::const_),
+            "Set 'dofval",
+            py::arg("nodevec"))
+
+        .def(
+            "AsDofs",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::VectorPartitioned::AsDofs, py::const_),
+            "Set 'dofval",
+            py::arg("elemvec"))
+
+        .def(
+            "AsDofs_u",
+            py::overload_cast<const xt::xtensor<double, 2>&>(
+                &GooseFEM::VectorPartitioned::AsDofs_u, py::const_),
+            "Set 'dofval",
+            py::arg("nodevec"))
+
+        .def(
+            "AsDofs_u",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::VectorPartitioned::AsDofs_u, py::const_),
+            "Set 'dofval",
+            py::arg("elemvec"))
+
+        .def(
+            "AsDofs_p",
+            py::overload_cast<const xt::xtensor<double, 2>&>(
+                &GooseFEM::VectorPartitioned::AsDofs_p, py::const_),
+            "Set 'dofval",
+            py::arg("nodevec"))
+
+        .def(
+            "AsDofs_p",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::VectorPartitioned::AsDofs_p, py::const_),
+            "Set 'dofval",
+            py::arg("elemvec"))
+
+        .def(
+            "AsNode",
+            py::overload_cast<const xt::xtensor<double, 1>&, const xt::xtensor<double, 1>&>(
+                &GooseFEM::VectorPartitioned::AsNode, py::const_),
+            "Set 'nodevec",
+            py::arg("dofval_u"),
+            py::arg("dofval_p"))
+
+        .def(
+            "AsNode",
+            py::overload_cast<const xt::xtensor<double, 1>&>(
+                &GooseFEM::VectorPartitioned::AsNode, py::const_),
+            "Set 'nodevec",
+            py::arg("dofval"))
+
+        .def(
+            "AsNode",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::VectorPartitioned::AsNode, py::const_),
+            "Set 'nodevec",
+            py::arg("elemvec"))
+
+        .def(
+            "AsElement",
+            py::overload_cast<const xt::xtensor<double, 1>&, const xt::xtensor<double, 1>&>(
+                &GooseFEM::VectorPartitioned::AsElement, py::const_),
+            "Set 'elemvec",
+            py::arg("dofval_u"),
+            py::arg("dofval_p"))
+
+        .def(
+            "AsElement",
+            py::overload_cast<const xt::xtensor<double, 1>&>(
+                &GooseFEM::VectorPartitioned::AsElement, py::const_),
+            "Set 'elemvec",
+            py::arg("dofval"))
+
+        .def(
+            "AsElement",
+            py::overload_cast<const xt::xtensor<double, 2>&>(
+                &GooseFEM::VectorPartitioned::AsElement, py::const_),
+            "Set 'elemvec",
+            py::arg("nodevec"))
+
+        .def(
+            "AssembleDofs",
+            py::overload_cast<const xt::xtensor<double, 2>&>(
+                &GooseFEM::VectorPartitioned::AssembleDofs, py::const_),
+            "Assemble 'dofval'",
+            py::arg("nodevec"))
+
+        .def(
+            "AssembleDofs",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::VectorPartitioned::AssembleDofs, py::const_),
+            "Assemble 'dofval'",
+            py::arg("elemvec"))
+
+        .def(
+            "AssembleDofs_u",
+            py::overload_cast<const xt::xtensor<double, 2>&>(
+                &GooseFEM::VectorPartitioned::AssembleDofs_u, py::const_),
+            "Assemble 'dofval'",
+            py::arg("nodevec"))
+
+        .def(
+            "AssembleDofs_u",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::VectorPartitioned::AssembleDofs_u, py::const_),
+            "Assemble 'dofval'",
+            py::arg("elemvec"))
+
+        .def(
+            "AssembleDofs_p",
+            py::overload_cast<const xt::xtensor<double, 2>&>(
+                &GooseFEM::VectorPartitioned::AssembleDofs_p, py::const_),
+            "Assemble 'dofval'",
+            py::arg("nodevec"))
+
+        .def(
+            "AssembleDofs_p",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::VectorPartitioned::AssembleDofs_p, py::const_),
+            "Assemble 'dofval'",
+            py::arg("elemvec"))
+
+        .def(
+            "AssembleNode",
+            py::overload_cast<const xt::xtensor<double, 3>&>(
+                &GooseFEM::VectorPartitioned::AssembleNode, py::const_),
+            "Assemble 'nodevec'",
+            py::arg("elemvec"))
+
+        .def("Copy", &GooseFEM::VectorPartitioned::Copy, "Copy")
+
+        .def("Copy_u", &GooseFEM::VectorPartitioned::Copy_u, "Copy iiu")
+
+        .def("Copy_p", &GooseFEM::VectorPartitioned::Copy_p, "Copy iip")
+
+        .def("__repr__", [](const GooseFEM::VectorPartitioned&) {
+            return "<GooseFEM.VectorPartitioned>";
+        });
 }
-
-// =================================================================================================
-
diff --git a/python/VectorPartitionedTyings.hpp b/python/VectorPartitionedTyings.hpp
index e9e940f..6f9a336 100644
--- a/python/VectorPartitionedTyings.hpp
+++ b/python/VectorPartitionedTyings.hpp
@@ -1,122 +1,90 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
+#include <GooseFEM/GooseFEM.h>
 #include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
-#include "../include/GooseFEM/GooseFEM.h"
-
-// =================================================================================================
 
 namespace py = pybind11;
 
-// =================================================================================================
-
 void init_VectorPartitionedTyings(py::module& m)
 {
 
-py::class_<GooseFEM::VectorPartitionedTyings>(m, "VectorPartitionedTyings")
-
-    .def(py::init<
-            const xt::xtensor<size_t,2>&,
-            const xt::xtensor<size_t,2>&,
-            const Eigen::SparseMatrix<double>&,
-            const Eigen::SparseMatrix<double>&,
-            const Eigen::SparseMatrix<double>&>(),
-        "Switch between dofval/nodevec/elemvec",
-        py::arg("conn"),
-        py::arg("dofs"),
-        py::arg("Cdu"),
-        py::arg("Cdp"),
-        py::arg("Cdi"))
-
-    .def("nelem",
-        &GooseFEM::VectorPartitionedTyings::nelem,
-        "Number of element")
-
-    .def("nne",
-        &GooseFEM::VectorPartitionedTyings::nne,
-        "Number of nodes per element")
-
-    .def("nnode",
-        &GooseFEM::VectorPartitionedTyings::nnode,
-        "Number of nodes")
-
-    .def("ndim",
-        &GooseFEM::VectorPartitionedTyings::ndim,
-        "Number of dimensions")
-
-    .def("ndof",
-        &GooseFEM::VectorPartitionedTyings::ndof,
-        "Number of degrees-of-freedom")
-
-    .def("nnu",
-        &GooseFEM::VectorPartitionedTyings::nnu,
-        "Number of unknown degrees-of-freedom")
-
-    .def("nnp",
-        &GooseFEM::VectorPartitionedTyings::nnp,
-        "Number of prescribed degrees-of-freedom")
-
-    .def("nni",
-        &GooseFEM::VectorPartitionedTyings::nni,
-        "Number of independent degrees-of-freedom")
-
-    .def("nnd",
-        &GooseFEM::VectorPartitionedTyings::nnd,
-        "Number of dependent degrees-of-freedom")
-
-    .def("dofs",
-        &GooseFEM::VectorPartitionedTyings::dofs,
-        "Return degrees-of-freedom")
-
-    .def("iiu",
-        &GooseFEM::VectorPartitionedTyings::iiu,
-        "Return unknown degrees-of-freedom")
-
-    .def("iip",
-        &GooseFEM::VectorPartitionedTyings::iip,
-        "Return prescribed degrees-of-freedom")
-
-    .def("iii",
-        &GooseFEM::VectorPartitionedTyings::iii,
-        "Return independent degrees-of-freedom")
-
-    .def("iid",
-        &GooseFEM::VectorPartitionedTyings::iid,
-        "Return dependent degrees-of-freedom")
-
-    .def("AsDofs_i",
-        &GooseFEM::VectorPartitionedTyings::AsDofs_i,
-        "Set 'dofval",
-        py::arg("nodevec"))
-
-    .def("AsNode",
-        &GooseFEM::VectorPartitionedTyings::AsNode,
-        "Set 'nodevec",
-        py::arg("dofval"))
-
-    .def("AsElement",
-        &GooseFEM::VectorPartitionedTyings::AsElement,
-        "Set 'elemvec",
-        py::arg("nodevec"))
-
-    .def("AssembleDofs",
-        &GooseFEM::VectorPartitionedTyings::AssembleDofs,
-        "Assemble 'dofval'",
-        py::arg("elemvec"))
-
-    .def("AssembleNode",
-        &GooseFEM::VectorPartitionedTyings::AssembleNode,
-        "Assemble 'nodevec'",
-        py::arg("elemvec"))
-
-    .def("__repr__", [](
-        const GooseFEM::VectorPartitionedTyings&){ return "<GooseFEM.VectorPartitionedTyings>"; });
+    py::class_<GooseFEM::VectorPartitionedTyings>(m, "VectorPartitionedTyings")
 
-}
+        .def(
+            py::init<
+                const xt::xtensor<size_t, 2>&,
+                const xt::xtensor<size_t, 2>&,
+                const Eigen::SparseMatrix<double>&,
+                const Eigen::SparseMatrix<double>&,
+                const Eigen::SparseMatrix<double>&>(),
+            "Switch between dofval/nodevec/elemvec",
+            py::arg("conn"),
+            py::arg("dofs"),
+            py::arg("Cdu"),
+            py::arg("Cdp"),
+            py::arg("Cdi"))
+
+        .def("nelem", &GooseFEM::VectorPartitionedTyings::nelem, "Number of element")
+
+        .def("nne", &GooseFEM::VectorPartitionedTyings::nne, "Number of nodes per element")
+
+        .def("nnode", &GooseFEM::VectorPartitionedTyings::nnode, "Number of nodes")
+
+        .def("ndim", &GooseFEM::VectorPartitionedTyings::ndim, "Number of dimensions")
+
+        .def("ndof", &GooseFEM::VectorPartitionedTyings::ndof, "Number of DOFs")
+
+        .def("nnu", &GooseFEM::VectorPartitionedTyings::nnu, "Number of unknown DOFs")
+
+        .def("nnp", &GooseFEM::VectorPartitionedTyings::nnp, "Number of prescribed DOFs")
+
+        .def("nni", &GooseFEM::VectorPartitionedTyings::nni, "Number of independent DOFs")
 
-// =================================================================================================
+        .def("nnd", &GooseFEM::VectorPartitionedTyings::nnd, "Number of dependent DOFs")
 
+        .def("dofs", &GooseFEM::VectorPartitionedTyings::dofs, "DOFs")
+
+        .def("iiu", &GooseFEM::VectorPartitionedTyings::iiu, "Unknown DOFs")
+
+        .def("iip", &GooseFEM::VectorPartitionedTyings::iip, "Prescribed DOFs")
+
+        .def("iii", &GooseFEM::VectorPartitionedTyings::iii, "Independent DOFs")
+
+        .def("iid", &GooseFEM::VectorPartitionedTyings::iid, "Dependent DOFs")
+
+        .def(
+            "AsDofs_i",
+            &GooseFEM::VectorPartitionedTyings::AsDofs_i,
+            "Set 'dofval",
+            py::arg("nodevec"))
+
+        .def(
+            "AsNode", &GooseFEM::VectorPartitionedTyings::AsNode, "Set 'nodevec", py::arg("dofval"))
+
+        .def(
+            "AsElement",
+            &GooseFEM::VectorPartitionedTyings::AsElement,
+            "Set 'elemvec",
+            py::arg("nodevec"))
+
+        .def(
+            "AssembleDofs",
+            &GooseFEM::VectorPartitionedTyings::AssembleDofs,
+            "Assemble 'dofval'",
+            py::arg("elemvec"))
+
+        .def(
+            "AssembleNode",
+            &GooseFEM::VectorPartitionedTyings::AssembleNode,
+            "Assemble 'nodevec'",
+            py::arg("elemvec"))
+
+        .def("__repr__", [](const GooseFEM::VectorPartitionedTyings&) {
+            return "<GooseFEM.VectorPartitionedTyings>";
+        });
+}
diff --git a/python/main.cpp b/python/main.cpp
index 1bb9f14..831b116 100644
--- a/python/main.cpp
+++ b/python/main.cpp
@@ -1,140 +1,136 @@
 /* =================================================================================================
 
 (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 "../include/GooseFEM/GooseFEM.h"
+#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);
 
-// =================================================================================================
-
 }