diff --git a/include/GooseFEM/Matrix.hpp b/include/GooseFEM/Matrix.hpp
index adf511a..5a6e5b1 100644
--- a/include/GooseFEM/Matrix.hpp
+++ b/include/GooseFEM/Matrix.hpp
@@ -1,197 +1,196 @@
 /*
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 */
 
 #ifndef GOOSEFEM_MATRIX_HPP
 #define GOOSEFEM_MATRIX_HPP
 
 #include "Matrix.h"
 
 namespace GooseFEM {
 
 inline Matrix::Matrix(const xt::xtensor<size_t, 2>& conn, const xt::xtensor<size_t, 2>& dofs)
     : m_conn(conn), m_dofs(dofs)
 {
     m_nelem = m_conn.shape(0);
     m_nne = m_conn.shape(1);
     m_nnode = m_dofs.shape(0);
     m_ndim = m_dofs.shape(1);
     m_ndof = xt::amax(m_dofs)() + 1;
     m_T.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
     m_A.resize(m_ndof, m_ndof);
 
-    GOOSEFEM_ASSERT(xt::amax(m_conn)[0] + 1 == m_nnode);
+    GOOSEFEM_ASSERT(xt::amax(m_conn)() + 1 <= m_nnode);
     GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim);
 }
 
 inline size_t Matrix::nelem() const
 {
     return m_nelem;
 }
 
 inline size_t Matrix::nne() const
 {
     return m_nne;
 }
 
 inline size_t Matrix::nnode() const
 {
     return m_nnode;
 }
 
 inline size_t Matrix::ndim() const
 {
     return m_ndim;
 }
 
 inline size_t Matrix::ndof() const
 {
     return m_ndof;
 }
 
 inline xt::xtensor<size_t, 2> Matrix::dofs() const
 {
     return m_dofs;
 }
 
 inline void Matrix::assemble(const xt::xtensor<double, 3>& elemmat)
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim}));
 
     m_T.clear();
 
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 for (size_t n = 0; n < m_nne; ++n) {
                     for (size_t j = 0; j < m_ndim; ++j) {
                         m_T.push_back(Eigen::Triplet<double>(
                             m_dofs(m_conn(e, m), i),
                             m_dofs(m_conn(e, n), j),
                             elemmat(e, m * m_ndim + i, n * m_ndim + j)));
                     }
                 }
             }
         }
     }
 
     m_A.setFromTriplets(m_T.begin(), m_T.end());
     m_changed = true;
 }
 
 inline void Matrix::dot(const xt::xtensor<double, 2>& x, xt::xtensor<double, 2>& b) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(b, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(x, {m_nnode, m_ndim}));
     this->asNode(m_A * this->AsDofs(x), b);
 }
 
 inline void Matrix::dot(const xt::xtensor<double, 1>& x, xt::xtensor<double, 1>& b) const
 {
     GOOSEFEM_ASSERT(b.size() == m_ndof);
     GOOSEFEM_ASSERT(x.size() == m_ndof);
 
     Eigen::Map<Eigen::VectorXd>(b.data(), b.size()).noalias() =
         m_A * Eigen::Map<const Eigen::VectorXd>(x.data(), x.size());
 }
 
 inline xt::xtensor<double, 2> Matrix::Dot(const xt::xtensor<double, 2>& x) const
 {
     xt::xtensor<double, 2> b = xt::empty<double>({m_nnode, m_ndim});
     this->dot(x, b);
     return b;
 }
 
 inline xt::xtensor<double, 1> Matrix::Dot(const xt::xtensor<double, 1>& x) const
 {
     xt::xtensor<double, 1> b = xt::empty<double>({m_ndof});
     this->dot(x, b);
     return b;
 }
 
 inline Eigen::VectorXd Matrix::AsDofs(const xt::xtensor<double, 2>& nodevec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
-    Eigen::VectorXd dofval(m_ndof, 1);
+    Eigen::VectorXd dofval = Eigen::VectorXd::Zero(m_ndof, 1);
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             dofval(m_dofs(m, i)) = nodevec(m, i);
         }
     }
 
     return dofval;
 }
 
 inline void Matrix::asNode(const Eigen::VectorXd& dofval, xt::xtensor<double, 2>& nodevec) const
 {
     GOOSEFEM_ASSERT(static_cast<size_t>(dofval.size()) == m_ndof);
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             nodevec(m, i) = dofval(m_dofs(m, i));
         }
     }
 }
 
 template <class Solver>
 inline void MatrixSolver<Solver>::factorize(Matrix& matrix)
 {
     if (!matrix.m_changed && !m_factor) {
         return;
     }
     m_solver.compute(matrix.m_A);
     m_factor = false;
     matrix.m_changed = false;
 }
 
 template <class Solver>
 inline void MatrixSolver<Solver>::solve(
     Matrix& matrix, const xt::xtensor<double, 2>& b, xt::xtensor<double, 2>& x)
 {
     GOOSEFEM_ASSERT(xt::has_shape(b, {matrix.m_nnode, matrix.m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(x, {matrix.m_nnode, matrix.m_ndim}));
     this->factorize(matrix);
-    Eigen::VectorXd B = matrix.AsDofs(b);
-    Eigen::VectorXd X = m_solver.solve(B);
+    Eigen::VectorXd X = m_solver.solve(matrix.AsDofs(b));
     matrix.asNode(X, x);
 }
 
 template <class Solver>
 inline void MatrixSolver<Solver>::solve(
     Matrix& matrix, const xt::xtensor<double, 1>& b, xt::xtensor<double, 1>& x)
 {
     GOOSEFEM_ASSERT(b.size() == matrix.m_ndof);
     GOOSEFEM_ASSERT(x.size() == matrix.m_ndof);
     this->factorize(matrix);
     Eigen::Map<Eigen::VectorXd>(x.data(), x.size()).noalias() = m_solver.solve(
         Eigen::Map<const Eigen::VectorXd>(b.data(), matrix.m_ndof));
 }
 
 template <class Solver>
 inline xt::xtensor<double, 2>
 MatrixSolver<Solver>::Solve(Matrix& matrix, const xt::xtensor<double, 2>& b)
 {
     xt::xtensor<double, 2> x = xt::empty<double>({matrix.m_nnode, matrix.m_ndim});
     this->solve(matrix, b, x);
     return x;
 }
 
 template <class Solver>
 inline xt::xtensor<double, 1>
 MatrixSolver<Solver>::Solve(Matrix& matrix, const xt::xtensor<double, 1>& b)
 {
     xt::xtensor<double, 1> x = xt::empty<double>({matrix.m_ndof});
     this->solve(matrix, b, x);
     return x;
 }
 
 } // namespace GooseFEM
 
 #endif
diff --git a/include/GooseFEM/MatrixDiagonal.hpp b/include/GooseFEM/MatrixDiagonal.hpp
index 886d551..07eca8b 100644
--- a/include/GooseFEM/MatrixDiagonal.hpp
+++ b/include/GooseFEM/MatrixDiagonal.hpp
@@ -1,177 +1,178 @@
 /*
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 */
 
 #ifndef GOOSEFEM_MATRIXDIAGONAL_HPP
 #define GOOSEFEM_MATRIXDIAGONAL_HPP
 
 #include "MatrixDiagonal.h"
 
 namespace GooseFEM {
 
 inline MatrixDiagonal::MatrixDiagonal(
     const xt::xtensor<size_t, 2>& conn, const xt::xtensor<size_t, 2>& dofs)
     : m_conn(conn), m_dofs(dofs)
 {
     m_nelem = m_conn.shape(0);
     m_nne = m_conn.shape(1);
     m_nnode = m_dofs.shape(0);
     m_ndim = m_dofs.shape(1);
     m_ndof = xt::amax(m_dofs)() + 1;
     m_A = xt::empty<double>({m_ndof});
     m_inv = xt::empty<double>({m_ndof});
 
-    GOOSEFEM_ASSERT(xt::amax(m_conn)() + 1 == m_nnode);
+    GOOSEFEM_ASSERT(xt::amax(m_conn)() + 1 <= m_nnode);
     GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim);
 }
 
 inline size_t MatrixDiagonal::nelem() const
 {
     return m_nelem;
 }
 
 inline size_t MatrixDiagonal::nne() const
 {
     return m_nne;
 }
 
 inline size_t MatrixDiagonal::nnode() const
 {
     return m_nnode;
 }
 
 inline size_t MatrixDiagonal::ndim() const
 {
     return m_ndim;
 }
 
 inline size_t MatrixDiagonal::ndof() const
 {
     return m_ndof;
 }
 
 inline xt::xtensor<size_t, 2> MatrixDiagonal::dofs() const
 {
     return m_dofs;
 }
 
 inline void MatrixDiagonal::factorize()
 {
     if (!m_factor) {
         return;
     }
 
     #pragma omp parallel for
-    for (size_t d = 0; d < m_ndof; ++d)
+    for (size_t d = 0; d < m_ndof; ++d) {
         m_inv(d) = 1.0 / m_A(d);
+    }
 
     m_factor = false;
 }
 
 inline void MatrixDiagonal::set(const xt::xtensor<double, 1>& A)
 {
     GOOSEFEM_ASSERT(A.size() == m_ndof);
     std::copy(A.begin(), A.end(), m_A.begin());
     m_factor = true;
 }
 
 inline void MatrixDiagonal::assemble(const xt::xtensor<double, 3>& elemmat)
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim}));
     GOOSEFEM_ASSERT(Element::isDiagonal(elemmat));
 
     m_A.fill(0.0);
 
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 m_A(m_dofs(m_conn(e, m), i)) += elemmat(e, m * m_ndim + i, m * m_ndim + i);
             }
         }
     }
 
     m_factor = true;
 }
 
 inline void MatrixDiagonal::dot(const xt::xtensor<double, 2>& x, xt::xtensor<double, 2>& b) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(x, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(b, {m_nnode, m_ndim}));
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             b(m, i) = m_A(m_dofs(m, i)) * x(m, i);
         }
     }
 }
 
 inline void MatrixDiagonal::dot(const xt::xtensor<double, 1>& x, xt::xtensor<double, 1>& b) const
 {
     GOOSEFEM_ASSERT(x.size() == m_ndof);
     GOOSEFEM_ASSERT(b.size() == m_ndof);
 
     xt::noalias(b) = m_A * x;
 }
 
 inline void MatrixDiagonal::solve(const xt::xtensor<double, 2>& b, xt::xtensor<double, 2>& x)
 {
     GOOSEFEM_ASSERT(xt::has_shape(b, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(x, {m_nnode, m_ndim}));
 
     this->factorize();
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             x(m, i) = m_inv(m_dofs(m, i)) * b(m, i);
         }
     }
 }
 
 inline void MatrixDiagonal::solve(const xt::xtensor<double, 1>& b, xt::xtensor<double, 1>& x)
 {
     GOOSEFEM_ASSERT(b.size() == m_ndof);
     GOOSEFEM_ASSERT(x.size() == m_ndof);
     this->factorize();
     xt::noalias(x) = m_inv * b;
 }
 
 inline xt::xtensor<double, 1> MatrixDiagonal::AsDiagonal() const
 {
     return m_A;
 }
 
 inline xt::xtensor<double, 2> MatrixDiagonal::Dot(const xt::xtensor<double, 2>& x) const
 {
     xt::xtensor<double, 2> b = xt::empty<double>({m_nnode, m_ndim});
     this->dot(x, b);
     return b;
 }
 
 inline xt::xtensor<double, 1> MatrixDiagonal::Dot(const xt::xtensor<double, 1>& x) const
 {
     xt::xtensor<double, 1> b = xt::empty<double>({m_ndof});
     this->dot(x, b);
     return b;
 }
 
 inline xt::xtensor<double, 2> MatrixDiagonal::Solve(const xt::xtensor<double, 2>& b)
 {
     xt::xtensor<double, 2> x = xt::empty<double>({m_nnode, m_ndim});
     this->solve(b, x);
     return x;
 }
 
 inline xt::xtensor<double, 1> MatrixDiagonal::Solve(const xt::xtensor<double, 1>& b)
 {
     xt::xtensor<double, 1> x = xt::empty<double>({m_ndof});
     this->solve(b, x);
     return x;
 }
 
 } // namespace GooseFEM
 
 #endif
diff --git a/include/GooseFEM/MatrixDiagonalPartitioned.hpp b/include/GooseFEM/MatrixDiagonalPartitioned.hpp
index 0dcf1e2..0eec5d1 100644
--- a/include/GooseFEM/MatrixDiagonalPartitioned.hpp
+++ b/include/GooseFEM/MatrixDiagonalPartitioned.hpp
@@ -1,396 +1,396 @@
 /*
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 */
 
 #ifndef GOOSEFEM_MATRIXDIAGONALPARTITIONED_HPP
 #define GOOSEFEM_MATRIXDIAGONALPARTITIONED_HPP
 
 #include "MatrixDiagonalPartitioned.h"
 #include "Mesh.h"
 
 namespace GooseFEM {
 
 inline MatrixDiagonalPartitioned::MatrixDiagonalPartitioned(
     const xt::xtensor<size_t, 2>& conn,
     const xt::xtensor<size_t, 2>& dofs,
     const xt::xtensor<size_t, 1>& iip)
     : m_conn(conn), m_dofs(dofs), m_iip(iip)
 {
     m_nelem = m_conn.shape(0);
     m_nne = m_conn.shape(1);
     m_nnode = m_dofs.shape(0);
     m_ndim = m_dofs.shape(1);
     m_iiu = xt::setdiff1d(dofs, iip);
     m_ndof = xt::amax(m_dofs)() + 1;
     m_nnp = m_iip.size();
     m_nnu = m_iiu.size();
     m_part = Mesh::Reorder({m_iiu, m_iip}).get(m_dofs);
     m_Auu = xt::empty<double>({m_nnu});
     m_App = xt::empty<double>({m_nnp});
     m_inv_uu = xt::empty<double>({m_nnu});
 
-    GOOSEFEM_ASSERT(xt::amax(m_conn)() + 1 == m_nnode);
+    GOOSEFEM_ASSERT(xt::amax(m_conn)() + 1 <= m_nnode);
     GOOSEFEM_ASSERT(xt::amax(m_iip)() <= xt::amax(m_dofs)());
     GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim);
 }
 
 inline size_t MatrixDiagonalPartitioned::nelem() const
 {
     return m_nelem;
 }
 
 inline size_t MatrixDiagonalPartitioned::nne() const
 {
     return m_nne;
 }
 
 inline size_t MatrixDiagonalPartitioned::nnode() const
 {
     return m_nnode;
 }
 
 inline size_t MatrixDiagonalPartitioned::ndim() const
 {
     return m_ndim;
 }
 
 inline size_t MatrixDiagonalPartitioned::ndof() const
 {
     return m_ndof;
 }
 
 inline size_t MatrixDiagonalPartitioned::nnu() const
 {
     return m_nnu;
 }
 
 inline size_t MatrixDiagonalPartitioned::nnp() const
 {
     return m_nnp;
 }
 
 inline xt::xtensor<size_t, 2> MatrixDiagonalPartitioned::dofs() const
 {
     return m_dofs;
 }
 
 inline xt::xtensor<size_t, 1> MatrixDiagonalPartitioned::iiu() const
 {
     return m_iiu;
 }
 
 inline xt::xtensor<size_t, 1> MatrixDiagonalPartitioned::iip() const
 {
     return m_iip;
 }
 
 inline void MatrixDiagonalPartitioned::factorize()
 {
     if (!m_factor) {
         return;
     }
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnu; ++d) {
         m_inv_uu(d) = 1.0 / m_Auu(d);
     }
 
     m_factor = false;
 }
 
 inline void MatrixDiagonalPartitioned::assemble(const xt::xtensor<double, 3>& elemmat)
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim}));
     GOOSEFEM_ASSERT(Element::isDiagonal(elemmat));
 
     m_Auu.fill(0.0);
     m_App.fill(0.0);
 
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
 
                 size_t d = m_part(m_conn(e, m), i);
 
                 if (d < m_nnu) {
                     m_Auu(d) += elemmat(e, m * m_ndim + i, m * m_ndim + i);
                 }
                 else {
                     m_App(d - m_nnu) += elemmat(e, m * m_ndim + i, m * m_ndim + i);
                 }
             }
         }
     }
 
     m_factor = true;
 }
 
 inline void
 MatrixDiagonalPartitioned::dot(const xt::xtensor<double, 2>& x, xt::xtensor<double, 2>& b) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(x, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(b, {m_nnode, m_ndim}));
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
 
             size_t d = m_part(m, i);
 
             if (d < m_nnu) {
                 b(m, i) = m_Auu(d) * x(m, i);
             }
             else {
                 b(m, i) = m_App(d - m_nnu) * x(m, i);
             }
         }
     }
 }
 
 inline void
 MatrixDiagonalPartitioned::dot(const xt::xtensor<double, 1>& x, xt::xtensor<double, 1>& b) const
 {
     GOOSEFEM_ASSERT(x.size() == m_ndof);
     GOOSEFEM_ASSERT(b.size() == m_ndof);
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnu; ++d) {
         b(m_iiu(d)) = m_Auu(d) * x(m_iiu(d));
     }
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnp; ++d) {
         b(m_iip(d)) = m_App(d) * x(m_iip(d));
     }
 }
 
 inline void MatrixDiagonalPartitioned::dot_u(
     const xt::xtensor<double, 1>& x_u,
     const xt::xtensor<double, 1>& x_p,
     xt::xtensor<double, 1>& b_u) const
 {
     UNUSED(x_p);
 
     GOOSEFEM_ASSERT(x_u.size() == m_nnu);
     GOOSEFEM_ASSERT(x_p.size() == m_nnp);
     GOOSEFEM_ASSERT(b_u.size() == m_nnu);
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnu; ++d) {
         b_u(d) = m_Auu(d) * x_u(d);
     }
 }
 
 inline void MatrixDiagonalPartitioned::dot_p(
     const xt::xtensor<double, 1>& x_u,
     const xt::xtensor<double, 1>& x_p,
     xt::xtensor<double, 1>& b_p) const
 {
     UNUSED(x_u);
 
     GOOSEFEM_ASSERT(x_u.size() == m_nnu);
     GOOSEFEM_ASSERT(x_p.size() == m_nnp);
     GOOSEFEM_ASSERT(b_p.size() == m_nnp);
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnp; ++d) {
         b_p(d) = m_App(d) * x_p(d);
     }
 }
 
 inline void
 MatrixDiagonalPartitioned::solve(const xt::xtensor<double, 2>& b, xt::xtensor<double, 2>& x)
 {
     GOOSEFEM_ASSERT(xt::has_shape(b, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(x, {m_nnode, m_ndim}));
 
     this->factorize();
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             if (m_part(m, i) < m_nnu) {
                 x(m, i) = m_inv_uu(m_part(m, i)) * b(m, i);
             }
         }
     }
 }
 
 inline void
 MatrixDiagonalPartitioned::solve(const xt::xtensor<double, 1>& b, xt::xtensor<double, 1>& x)
 {
     GOOSEFEM_ASSERT(b.size() == m_ndof);
     GOOSEFEM_ASSERT(x.size() == m_ndof);
 
     this->factorize();
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnu; ++d) {
         x(m_iiu(d)) = m_inv_uu(d) * b(m_iiu(d));
     }
 }
 
 inline void MatrixDiagonalPartitioned::solve_u(
     const xt::xtensor<double, 1>& b_u,
     const xt::xtensor<double, 1>& x_p,
     xt::xtensor<double, 1>& x_u)
 {
     UNUSED(x_p);
 
     GOOSEFEM_ASSERT(b_u.size() == m_nnu);
     GOOSEFEM_ASSERT(x_p.size() == m_nnp);
     GOOSEFEM_ASSERT(x_u.size() == m_nnu);
 
     this->factorize();
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnu; ++d) {
         x_u(d) = m_inv_uu(d) * b_u(d);
     }
 }
 
 inline void MatrixDiagonalPartitioned::reaction(
     const xt::xtensor<double, 2>& x, xt::xtensor<double, 2>& b) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(x, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(b, {m_nnode, m_ndim}));
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             if (m_part(m, i) >= m_nnu) {
                 b(m, i) = m_App(m_part(m, i) - m_nnu) * x(m, i);
             }
         }
     }
 }
 
 inline void MatrixDiagonalPartitioned::reaction(
     const xt::xtensor<double, 1>& x, xt::xtensor<double, 1>& b) const
 {
     GOOSEFEM_ASSERT(x.size() == m_ndof);
     GOOSEFEM_ASSERT(b.size() == m_ndof);
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnp; ++d) {
         b(m_iip(d)) = m_App(d) * x(m_iip(d));
     }
 }
 
 inline void MatrixDiagonalPartitioned::reaction_p(
     const xt::xtensor<double, 1>& x_u,
     const xt::xtensor<double, 1>& x_p,
     xt::xtensor<double, 1>& b_p) const
 {
     UNUSED(x_u);
 
     GOOSEFEM_ASSERT(x_u.size() == m_nnu);
     GOOSEFEM_ASSERT(x_p.size() == m_nnp);
     GOOSEFEM_ASSERT(b_p.size() == m_nnp);
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnp; ++d) {
         b_p(d) = m_App(d) * x_p(d);
     }
 }
 
 inline xt::xtensor<double, 1> MatrixDiagonalPartitioned::AsDiagonal() const
 {
-    xt::xtensor<double, 1> out = xt::zeros<double>({m_ndof});
+    xt::xtensor<double, 1> ret = xt::zeros<double>({m_ndof});
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnu; ++d) {
-        out(m_iiu(d)) = m_Auu(d);
+        ret(m_iiu(d)) = m_Auu(d);
     }
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnp; ++d) {
-        out(m_iip(d)) = m_App(d);
+        ret(m_iip(d)) = m_App(d);
     }
 
-    return out;
+    return ret;
 }
 
 inline xt::xtensor<double, 2> MatrixDiagonalPartitioned::Dot(const xt::xtensor<double, 2>& x) const
 {
     xt::xtensor<double, 2> b = xt::empty<double>({m_nnode, m_ndim});
     this->dot(x, b);
     return b;
 }
 
 inline xt::xtensor<double, 1> MatrixDiagonalPartitioned::Dot(const xt::xtensor<double, 1>& x) const
 {
     xt::xtensor<double, 1> b = xt::empty<double>({m_ndof});
     this->dot(x, b);
     return b;
 }
 
 inline xt::xtensor<double, 1> MatrixDiagonalPartitioned::Dot_u(
     const xt::xtensor<double, 1>& x_u, const xt::xtensor<double, 1>& x_p) const
 {
     xt::xtensor<double, 1> b_u = xt::empty<double>({m_nnu});
     this->dot_u(x_u, x_p, b_u);
     return b_u;
 }
 
 inline xt::xtensor<double, 1> MatrixDiagonalPartitioned::Dot_p(
     const xt::xtensor<double, 1>& x_u, const xt::xtensor<double, 1>& x_p) const
 {
     xt::xtensor<double, 1> b_p = xt::empty<double>({m_nnp});
     this->dot_p(x_u, x_p, b_p);
     return b_p;
 }
 
 inline xt::xtensor<double, 2>
 MatrixDiagonalPartitioned::Solve(const xt::xtensor<double, 2>& b, const xt::xtensor<double, 2>& x)
 {
-    xt::xtensor<double, 2> out = x;
-    this->solve(b, out);
-    return out;
+    xt::xtensor<double, 2> ret = x;
+    this->solve(b, ret);
+    return ret;
 }
 
 inline xt::xtensor<double, 1>
 MatrixDiagonalPartitioned::Solve(const xt::xtensor<double, 1>& b, const xt::xtensor<double, 1>& x)
 {
-    xt::xtensor<double, 1> out = x;
-    this->solve(b, out);
-    return out;
+    xt::xtensor<double, 1> ret = x;
+    this->solve(b, ret);
+    return ret;
 }
 
 inline xt::xtensor<double, 1> MatrixDiagonalPartitioned::Solve_u(
     const xt::xtensor<double, 1>& b_u, const xt::xtensor<double, 1>& x_p)
 {
     xt::xtensor<double, 1> x_u = xt::empty<double>({m_nnu});
     this->solve_u(b_u, x_p, x_u);
     return x_u;
 }
 
 inline xt::xtensor<double, 2> MatrixDiagonalPartitioned::Reaction(
     const xt::xtensor<double, 2>& x, const xt::xtensor<double, 2>& b) const
 {
-    xt::xtensor<double, 2> out = b;
-    this->reaction(x, out);
-    return out;
+    xt::xtensor<double, 2> ret = b;
+    this->reaction(x, ret);
+    return ret;
 }
 
 inline xt::xtensor<double, 1> MatrixDiagonalPartitioned::Reaction(
     const xt::xtensor<double, 1>& x, const xt::xtensor<double, 1>& b) const
 {
-    xt::xtensor<double, 1> out = b;
-    this->reaction(x, out);
-    return out;
+    xt::xtensor<double, 1> ret = b;
+    this->reaction(x, ret);
+    return ret;
 }
 
 inline xt::xtensor<double, 1> MatrixDiagonalPartitioned::Reaction_p(
     const xt::xtensor<double, 1>& x_u, const xt::xtensor<double, 1>& x_p) const
 {
     xt::xtensor<double, 1> b_p = xt::empty<double>({m_nnp});
     this->reaction_p(x_u, x_p, b_p);
     return b_p;
 }
 
 } // namespace GooseFEM
 
 #endif
diff --git a/include/GooseFEM/MatrixPartitioned.hpp b/include/GooseFEM/MatrixPartitioned.hpp
index 170d886..c46c367 100644
--- a/include/GooseFEM/MatrixPartitioned.hpp
+++ b/include/GooseFEM/MatrixPartitioned.hpp
@@ -1,381 +1,381 @@
 /*
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 */
 
 #ifndef GOOSEFEM_MATRIXPARTITIONED_HPP
 #define GOOSEFEM_MATRIXPARTITIONED_HPP
 
 #include "MatrixPartitioned.h"
 #include "Mesh.h"
 
 namespace GooseFEM {
 
 inline MatrixPartitioned::MatrixPartitioned(
     const xt::xtensor<size_t, 2>& conn,
     const xt::xtensor<size_t, 2>& dofs,
     const xt::xtensor<size_t, 1>& iip)
     : m_conn(conn), m_dofs(dofs), m_iip(iip)
 {
     m_nelem = m_conn.shape(0);
     m_nne = m_conn.shape(1);
     m_nnode = m_dofs.shape(0);
     m_ndim = m_dofs.shape(1);
     m_iiu = xt::setdiff1d(dofs, iip);
     m_ndof = xt::amax(m_dofs)() + 1;
     m_nnp = m_iip.size();
     m_nnu = m_iiu.size();
     m_part = Mesh::Reorder({m_iiu, m_iip}).get(m_dofs);
     m_Tuu.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
     m_Tup.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
     m_Tpu.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
     m_Tpp.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
     m_Auu.resize(m_nnu, m_nnu);
     m_Aup.resize(m_nnu, m_nnp);
     m_Apu.resize(m_nnp, m_nnu);
     m_App.resize(m_nnp, m_nnp);
 
-    GOOSEFEM_ASSERT(xt::amax(m_conn)() + 1 == m_nnode);
+    GOOSEFEM_ASSERT(xt::amax(m_conn)() + 1 <= m_nnode);
     GOOSEFEM_ASSERT(xt::amax(m_iip)() <= xt::amax(m_dofs)());
     GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim);
 }
 
 inline size_t MatrixPartitioned::nelem() const
 {
     return m_nelem;
 }
 
 inline size_t MatrixPartitioned::nne() const
 {
     return m_nne;
 }
 
 inline size_t MatrixPartitioned::nnode() const
 {
     return m_nnode;
 }
 
 inline size_t MatrixPartitioned::ndim() const
 {
     return m_ndim;
 }
 
 inline size_t MatrixPartitioned::ndof() const
 {
     return m_ndof;
 }
 
 inline size_t MatrixPartitioned::nnu() const
 {
     return m_nnu;
 }
 
 inline size_t MatrixPartitioned::nnp() const
 {
     return m_nnp;
 }
 
 inline xt::xtensor<size_t, 2> MatrixPartitioned::dofs() const
 {
     return m_dofs;
 }
 
 inline xt::xtensor<size_t, 1> MatrixPartitioned::iiu() const
 {
     return m_iiu;
 }
 
 inline xt::xtensor<size_t, 1> MatrixPartitioned::iip() const
 {
     return m_iip;
 }
 
 inline void MatrixPartitioned::assemble(const xt::xtensor<double, 3>& elemmat)
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim}));
 
     m_Tuu.clear();
     m_Tup.clear();
     m_Tpu.clear();
     m_Tpp.clear();
 
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
 
                 size_t di = m_part(m_conn(e, m), i);
 
                 for (size_t n = 0; n < m_nne; ++n) {
                     for (size_t j = 0; j < m_ndim; ++j) {
 
                         size_t dj = m_part(m_conn(e, n), j);
 
                         if (di < m_nnu && dj < m_nnu) {
                             m_Tuu.push_back(Eigen::Triplet<double>(
                                 di, dj, elemmat(e, m * m_ndim + i, n * m_ndim + j)));
                         }
                         else if (di < m_nnu) {
                             m_Tup.push_back(Eigen::Triplet<double>(
                                 di, dj - m_nnu, elemmat(e, m * m_ndim + i, n * m_ndim + j)));
                         }
                         else if (dj < m_nnu) {
                             m_Tpu.push_back(Eigen::Triplet<double>(
                                 di - m_nnu, dj, elemmat(e, m * m_ndim + i, n * m_ndim + j)));
                         }
                         else {
                             m_Tpp.push_back(Eigen::Triplet<double>(
                                 di - m_nnu,
                                 dj - m_nnu,
                                 elemmat(e, m * m_ndim + i, n * m_ndim + j)));
                         }
                     }
                 }
             }
         }
     }
 
     m_Auu.setFromTriplets(m_Tuu.begin(), m_Tuu.end());
     m_Aup.setFromTriplets(m_Tup.begin(), m_Tup.end());
     m_Apu.setFromTriplets(m_Tpu.begin(), m_Tpu.end());
     m_App.setFromTriplets(m_Tpp.begin(), m_Tpp.end());
     m_changed = true;
 }
 
 inline void
 MatrixPartitioned::reaction(const xt::xtensor<double, 2>& x, xt::xtensor<double, 2>& b) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(x, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(b, {m_nnode, m_ndim}));
 
     Eigen::VectorXd X_u = this->AsDofs_u(x);
     Eigen::VectorXd X_p = this->AsDofs_p(x);
     Eigen::VectorXd B_p = m_Apu * X_u + m_App * X_p;
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             if (m_part(m, i) >= m_nnu) {
                 b(m, i) = B_p(m_part(m, i) - m_nnu);
             }
         }
     }
 }
 
 inline void
 MatrixPartitioned::reaction(const xt::xtensor<double, 1>& x, xt::xtensor<double, 1>& b) const
 {
     GOOSEFEM_ASSERT(x.size() == m_ndof);
     GOOSEFEM_ASSERT(b.size() == m_ndof);
 
     Eigen::VectorXd X_u = this->AsDofs_u(x);
     Eigen::VectorXd X_p = this->AsDofs_p(x);
     Eigen::VectorXd B_p = m_Apu * X_u + m_App * X_p;
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnp; ++d) {
         b(m_iip(d)) = B_p(d);
     }
 }
 
 inline void MatrixPartitioned::reaction_p(
     const xt::xtensor<double, 1>& x_u,
     const xt::xtensor<double, 1>& x_p,
     xt::xtensor<double, 1>& b_p) const
 {
     GOOSEFEM_ASSERT(x_u.size() == m_nnu);
     GOOSEFEM_ASSERT(x_p.size() == m_nnp);
     GOOSEFEM_ASSERT(b_p.size() == m_nnp);
 
     Eigen::Map<Eigen::VectorXd>(b_p.data(), b_p.size()).noalias() =
         m_Apu * Eigen::Map<const Eigen::VectorXd>(x_u.data(), x_u.size()) +
         m_App * Eigen::Map<const Eigen::VectorXd>(x_p.data(), x_p.size());
 }
 
 inline xt::xtensor<double, 2>
 MatrixPartitioned::Reaction(const xt::xtensor<double, 2>& x, const xt::xtensor<double, 2>& b) const
 {
-    xt::xtensor<double, 2> out = b;
-    this->reaction(x, out);
-    return out;
+    xt::xtensor<double, 2> ret = b;
+    this->reaction(x, ret);
+    return ret;
 }
 
 inline xt::xtensor<double, 1>
 MatrixPartitioned::Reaction(const xt::xtensor<double, 1>& x, const xt::xtensor<double, 1>& b) const
 {
-    xt::xtensor<double, 1> out = b;
-    this->reaction(x, out);
-    return out;
+    xt::xtensor<double, 1> ret = b;
+    this->reaction(x, ret);
+    return ret;
 }
 
 inline xt::xtensor<double, 1> MatrixPartitioned::Reaction_p(
     const xt::xtensor<double, 1>& x_u, const xt::xtensor<double, 1>& x_p) const
 {
     xt::xtensor<double, 1> b_p = xt::empty<double>({m_nnp});
     this->reaction_p(x_u, x_p, b_p);
     return b_p;
 }
 
 inline Eigen::VectorXd MatrixPartitioned::AsDofs_u(const xt::xtensor<double, 1>& dofval) const
 {
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
 
     Eigen::VectorXd dofval_u(m_nnu, 1);
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnu; ++d) {
         dofval_u(d) = dofval(m_iiu(d));
     }
 
     return dofval_u;
 }
 
 inline Eigen::VectorXd MatrixPartitioned::AsDofs_u(const xt::xtensor<double, 2>& nodevec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
-    Eigen::VectorXd dofval_u(m_nnu, 1);
+    Eigen::VectorXd dofval_u = Eigen::VectorXd::Zero(m_nnu, 1);
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             if (m_part(m, i) < m_nnu) {
                 dofval_u(m_part(m, i)) = nodevec(m, i);
             }
         }
     }
 
     return dofval_u;
 }
 
 inline Eigen::VectorXd MatrixPartitioned::AsDofs_p(const xt::xtensor<double, 1>& dofval) const
 {
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
 
     Eigen::VectorXd dofval_p(m_nnp, 1);
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnp; ++d) {
         dofval_p(d) = dofval(m_iip(d));
     }
 
     return dofval_p;
 }
 
 inline Eigen::VectorXd MatrixPartitioned::AsDofs_p(const xt::xtensor<double, 2>& nodevec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
-    Eigen::VectorXd dofval_p(m_nnp, 1);
+    Eigen::VectorXd dofval_p = Eigen::VectorXd::Zero(m_nnp, 1);
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             if (m_part(m, i) >= m_nnu) {
                 dofval_p(m_part(m, i) - m_nnu) = nodevec(m, i);
             }
         }
     }
 
     return dofval_p;
 }
 
 template <class Solver>
 inline void MatrixPartitionedSolver<Solver>::factorize(MatrixPartitioned& matrix)
 {
     if (!matrix.m_changed && !m_factor) {
         return;
     }
     m_solver.compute(matrix.m_Auu);
     m_factor = false;
     matrix.m_changed = false;
 }
 
 template <class Solver>
 inline void MatrixPartitionedSolver<Solver>::solve(
     MatrixPartitioned& matrix, const xt::xtensor<double, 2>& b, xt::xtensor<double, 2>& x)
 {
     GOOSEFEM_ASSERT(xt::has_shape(b, {matrix.m_nnode, matrix.m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(x, {matrix.m_nnode, matrix.m_ndim}));
 
     this->factorize(matrix);
     Eigen::VectorXd B_u = matrix.AsDofs_u(b);
     Eigen::VectorXd X_p = matrix.AsDofs_p(x);
     Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - matrix.m_Aup * X_p));
 
     #pragma omp parallel for
     for (size_t m = 0; m < matrix.m_nnode; ++m) {
         for (size_t i = 0; i < matrix.m_ndim; ++i) {
             if (matrix.m_part(m, i) < matrix.m_nnu) {
                 x(m, i) = X_u(matrix.m_part(m, i));
             }
         }
     }
 }
 
 template <class Solver>
 inline void MatrixPartitionedSolver<Solver>::solve(
     MatrixPartitioned& matrix, const xt::xtensor<double, 1>& b, xt::xtensor<double, 1>& x)
 {
     GOOSEFEM_ASSERT(b.size() == matrix.m_ndof);
     GOOSEFEM_ASSERT(x.size() == matrix.m_ndof);
 
     this->factorize(matrix);
     Eigen::VectorXd B_u = matrix.AsDofs_u(b);
     Eigen::VectorXd X_p = matrix.AsDofs_p(x);
     Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - matrix.m_Aup * X_p));
 
     #pragma omp parallel for
     for (size_t d = 0; d < matrix.m_nnu; ++d) {
         x(matrix.m_iiu(d)) = X_u(d);
     }
 }
 
 template <class Solver>
 inline void MatrixPartitionedSolver<Solver>::solve_u(
     MatrixPartitioned& matrix,
     const xt::xtensor<double, 1>& b_u,
     const xt::xtensor<double, 1>& x_p,
     xt::xtensor<double, 1>& x_u)
 {
     GOOSEFEM_ASSERT(b_u.size() == matrix.m_nnu);
     GOOSEFEM_ASSERT(x_p.size() == matrix.m_nnp);
     GOOSEFEM_ASSERT(x_u.size() == matrix.m_nnu);
 
     this->factorize(matrix);
 
     Eigen::Map<Eigen::VectorXd>(x_u.data(), x_u.size()).noalias() = m_solver.solve(Eigen::VectorXd(
         Eigen::Map<const Eigen::VectorXd>(b_u.data(), b_u.size()) -
         matrix.m_Aup * Eigen::Map<const Eigen::VectorXd>(x_p.data(), x_p.size())));
 }
 
 template <class Solver>
 inline xt::xtensor<double, 2> MatrixPartitionedSolver<Solver>::Solve(
     MatrixPartitioned& matrix, const xt::xtensor<double, 2>& b, const xt::xtensor<double, 2>& x)
 {
     xt::xtensor<double, 2> out = x;
     this->solve(matrix, b, out);
     return out;
 }
 
 template <class Solver>
 inline xt::xtensor<double, 1> MatrixPartitionedSolver<Solver>::Solve(
     MatrixPartitioned& matrix, const xt::xtensor<double, 1>& b, const xt::xtensor<double, 1>& x)
 {
     xt::xtensor<double, 1> out = x;
     this->solve(matrix, b, out);
     return out;
 }
 
 template <class Solver>
 inline xt::xtensor<double, 1> MatrixPartitionedSolver<Solver>::Solve_u(
     MatrixPartitioned& matrix, const xt::xtensor<double, 1>& b_u, const xt::xtensor<double, 1>& x_p)
 {
     xt::xtensor<double, 1> x_u = xt::empty<double>({matrix.m_nnu});
     this->solve_u(matrix, b_u, x_p, x_u);
     return x_u;
 }
 
 } // namespace GooseFEM
 
 #endif
diff --git a/include/GooseFEM/MatrixPartitionedTyings.hpp b/include/GooseFEM/MatrixPartitionedTyings.hpp
index cb8f829..8c2ed70 100644
--- a/include/GooseFEM/MatrixPartitionedTyings.hpp
+++ b/include/GooseFEM/MatrixPartitionedTyings.hpp
@@ -1,454 +1,454 @@
 /*
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 */
 
 #ifndef GOOSEFEM_MATRIXPARTITIONEDTYINGS_HPP
 #define GOOSEFEM_MATRIXPARTITIONEDTYINGS_HPP
 
 #include "MatrixPartitionedTyings.h"
 
 namespace GooseFEM {
 
 inline MatrixPartitionedTyings::MatrixPartitionedTyings(
     const xt::xtensor<size_t, 2>& conn,
     const xt::xtensor<size_t, 2>& dofs,
     const Eigen::SparseMatrix<double>& Cdu,
     const Eigen::SparseMatrix<double>& Cdp)
     : m_conn(conn), m_dofs(dofs), m_Cdu(Cdu), m_Cdp(Cdp)
 {
     GOOSEFEM_ASSERT(Cdu.rows() == Cdp.rows());
 
     m_nnu = static_cast<size_t>(m_Cdu.cols());
     m_nnp = static_cast<size_t>(m_Cdp.cols());
     m_nnd = static_cast<size_t>(m_Cdp.rows());
     m_nni = m_nnu + m_nnp;
     m_ndof = m_nni + m_nnd;
     m_iiu = xt::arange<size_t>(m_nnu);
     m_iip = xt::arange<size_t>(m_nnu, m_nnu + m_nnp);
     m_iid = xt::arange<size_t>(m_nni, m_nni + m_nnd);
     m_nelem = m_conn.shape(0);
     m_nne = m_conn.shape(1);
     m_nnode = m_dofs.shape(0);
     m_ndim = m_dofs.shape(1);
     m_Cud = m_Cdu.transpose();
     m_Cpd = m_Cdp.transpose();
     m_Tuu.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
     m_Tup.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
     m_Tpu.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
     m_Tpp.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
     m_Tud.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
     m_Tpd.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
     m_Tdu.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
     m_Tdp.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
     m_Tdd.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
     m_Auu.resize(m_nnu, m_nnu);
     m_Aup.resize(m_nnu, m_nnp);
     m_Apu.resize(m_nnp, m_nnu);
     m_App.resize(m_nnp, m_nnp);
     m_Aud.resize(m_nnu, m_nnd);
     m_Apd.resize(m_nnp, m_nnd);
     m_Adu.resize(m_nnd, m_nnu);
     m_Adp.resize(m_nnd, m_nnp);
     m_Add.resize(m_nnd, m_nnd);
 
     GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim);
     GOOSEFEM_ASSERT(m_ndof == xt::amax(m_dofs)() + 1);
 }
 
 inline size_t MatrixPartitionedTyings::nelem() const
 {
     return m_nelem;
 }
 
 inline size_t MatrixPartitionedTyings::nne() const
 {
     return m_nne;
 }
 
 inline size_t MatrixPartitionedTyings::nnode() const
 {
     return m_nnode;
 }
 
 inline size_t MatrixPartitionedTyings::ndim() const
 {
     return m_ndim;
 }
 
 inline size_t MatrixPartitionedTyings::ndof() const
 {
     return m_ndof;
 }
 
 inline size_t MatrixPartitionedTyings::nnu() const
 {
     return m_nnu;
 }
 
 inline size_t MatrixPartitionedTyings::nnp() const
 {
     return m_nnp;
 }
 
 inline size_t MatrixPartitionedTyings::nni() const
 {
     return m_nni;
 }
 
 inline size_t MatrixPartitionedTyings::nnd() const
 {
     return m_nnd;
 }
 
 inline xt::xtensor<size_t, 2> MatrixPartitionedTyings::dofs() const
 {
     return m_dofs;
 }
 
 inline xt::xtensor<size_t, 1> MatrixPartitionedTyings::iiu() const
 {
     return m_iiu;
 }
 
 inline xt::xtensor<size_t, 1> MatrixPartitionedTyings::iip() const
 {
     return m_iip;
 }
 
 inline xt::xtensor<size_t, 1> MatrixPartitionedTyings::iii() const
 {
     return xt::arange<size_t>(m_nni);
 }
 
 inline xt::xtensor<size_t, 1> MatrixPartitionedTyings::iid() const
 {
     return m_iid;
 }
 
 inline void MatrixPartitionedTyings::assemble(const xt::xtensor<double, 3>& elemmat)
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim}));
 
     m_Tuu.clear();
     m_Tup.clear();
     m_Tpu.clear();
     m_Tpp.clear();
     m_Tud.clear();
     m_Tpd.clear();
     m_Tdu.clear();
     m_Tdp.clear();
     m_Tdd.clear();
 
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
 
                 size_t di = m_dofs(m_conn(e, m), i);
 
                 for (size_t n = 0; n < m_nne; ++n) {
                     for (size_t j = 0; j < m_ndim; ++j) {
 
                         size_t dj = m_dofs(m_conn(e, n), j);
 
                         if (di < m_nnu && dj < m_nnu) {
                             m_Tuu.push_back(Eigen::Triplet<double>(
                                 di, dj, elemmat(e, m * m_ndim + i, n * m_ndim + j)));
                         }
                         else if (di < m_nnu && dj < m_nni) {
                             m_Tup.push_back(Eigen::Triplet<double>(
                                 di, dj - m_nnu, elemmat(e, m * m_ndim + i, n * m_ndim + j)));
                         }
                         else if (di < m_nnu) {
                             m_Tud.push_back(Eigen::Triplet<double>(
                                 di, dj - m_nni, elemmat(e, m * m_ndim + i, n * m_ndim + j)));
                         }
                         else if (di < m_nni && dj < m_nnu) {
                             m_Tpu.push_back(Eigen::Triplet<double>(
                                 di - m_nnu, dj, elemmat(e, m * m_ndim + i, n * m_ndim + j)));
                         }
                         else if (di < m_nni && dj < m_nni) {
                             m_Tpp.push_back(Eigen::Triplet<double>(
                                 di - m_nnu,
                                 dj - m_nnu,
                                 elemmat(e, m * m_ndim + i, n * m_ndim + j)));
                         }
                         else if (di < m_nni) {
                             m_Tpd.push_back(Eigen::Triplet<double>(
                                 di - m_nnu,
                                 dj - m_nni,
                                 elemmat(e, m * m_ndim + i, n * m_ndim + j)));
                         }
                         else if (dj < m_nnu) {
                             m_Tdu.push_back(Eigen::Triplet<double>(
                                 di - m_nni, dj, elemmat(e, m * m_ndim + i, n * m_ndim + j)));
                         }
                         else if (dj < m_nni) {
                             m_Tdp.push_back(Eigen::Triplet<double>(
                                 di - m_nni,
                                 dj - m_nnu,
                                 elemmat(e, m * m_ndim + i, n * m_ndim + j)));
                         }
                         else {
                             m_Tdd.push_back(Eigen::Triplet<double>(
                                 di - m_nni,
                                 dj - m_nni,
                                 elemmat(e, m * m_ndim + i, n * m_ndim + j)));
                         }
                     }
                 }
             }
         }
     }
 
     m_Auu.setFromTriplets(m_Tuu.begin(), m_Tuu.end());
     m_Aup.setFromTriplets(m_Tup.begin(), m_Tup.end());
     m_Apu.setFromTriplets(m_Tpu.begin(), m_Tpu.end());
     m_App.setFromTriplets(m_Tpp.begin(), m_Tpp.end());
     m_Aud.setFromTriplets(m_Tud.begin(), m_Tud.end());
     m_Apd.setFromTriplets(m_Tpd.begin(), m_Tpd.end());
     m_Adu.setFromTriplets(m_Tdu.begin(), m_Tdu.end());
     m_Adp.setFromTriplets(m_Tdp.begin(), m_Tdp.end());
     m_Add.setFromTriplets(m_Tdd.begin(), m_Tdd.end());
     m_changed = true;
 }
 
 inline Eigen::VectorXd MatrixPartitionedTyings::AsDofs_u(const xt::xtensor<double, 1>& dofval) const
 {
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
 
     Eigen::VectorXd dofval_u(m_nnu, 1);
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnu; ++d) {
         dofval_u(d) = dofval(m_iiu(d));
     }
 
     return dofval_u;
 }
 
 inline Eigen::VectorXd
 MatrixPartitionedTyings::AsDofs_u(const xt::xtensor<double, 2>& nodevec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
-    Eigen::VectorXd dofval_u(m_nnu, 1);
+    Eigen::VectorXd dofval_u = Eigen::VectorXd::Zero(m_nnu, 1);
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             if (m_dofs(m, i) < m_nnu) {
                 dofval_u(m_dofs(m, i)) = nodevec(m, i);
             }
         }
     }
 
     return dofval_u;
 }
 
 inline Eigen::VectorXd MatrixPartitionedTyings::AsDofs_p(const xt::xtensor<double, 1>& dofval) const
 {
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
 
     Eigen::VectorXd dofval_p(m_nnp, 1);
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnp; ++d) {
         dofval_p(d) = dofval(m_iip(d));
     }
 
     return dofval_p;
 }
 
 inline Eigen::VectorXd
 MatrixPartitionedTyings::AsDofs_p(const xt::xtensor<double, 2>& nodevec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
-    Eigen::VectorXd dofval_p(m_nnp, 1);
+    Eigen::VectorXd dofval_p = Eigen::VectorXd::Zero(m_nnp, 1);
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             if (m_dofs(m, i) >= m_nnu && m_dofs(m, i) < m_nni) {
                 dofval_p(m_dofs(m, i) - m_nnu) = nodevec(m, i);
             }
         }
     }
 
     return dofval_p;
 }
 
 inline Eigen::VectorXd MatrixPartitionedTyings::AsDofs_d(const xt::xtensor<double, 1>& dofval) const
 {
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
 
     Eigen::VectorXd dofval_d(m_nnd, 1);
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnd; ++d) {
         dofval_d(d) = dofval(m_iip(d));
     }
 
     return dofval_d;
 }
 
 inline Eigen::VectorXd
 MatrixPartitionedTyings::AsDofs_d(const xt::xtensor<double, 2>& nodevec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
-    Eigen::VectorXd dofval_d(m_nnd, 1);
+    Eigen::VectorXd dofval_d = Eigen::VectorXd::Zero(m_nnd, 1);
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             if (m_dofs(m, i) >= m_nni) {
                 dofval_d(m_dofs(m, i) - m_nni) = nodevec(m, i);
             }
         }
     }
 
     return dofval_d;
 }
 
 template <class Solver>
 inline void MatrixPartitionedTyingsSolver<Solver>::factorize(MatrixPartitionedTyings& matrix)
 {
     if (!matrix.m_changed && !m_factor) {
         return;
     }
 
     matrix.m_ACuu = matrix.m_Auu + matrix.m_Aud * matrix.m_Cdu + matrix.m_Cud * matrix.m_Adu +
                     matrix.m_Cud * matrix.m_Add * matrix.m_Cdu;
 
     matrix.m_ACup = matrix.m_Aup + matrix.m_Aud * matrix.m_Cdp + matrix.m_Cud * matrix.m_Adp +
                     matrix.m_Cud * matrix.m_Add * matrix.m_Cdp;
 
     // matrix.m_ACpu = matrix.m_Apu + matrix.m_Apd * matrix.m_Cdu + matrix.m_Cpd * matrix.m_Adu
     //     + matrix.m_Cpd * matrix.m_Add * matrix.m_Cdu;
 
     // matrix.m_ACpp = matrix.m_App + matrix.m_Apd * matrix.m_Cdp + matrix.m_Cpd * matrix.m_Adp
     //     + matrix.m_Cpd * matrix.m_Add * matrix.m_Cdp;
 
     m_solver.compute(matrix.m_ACuu);
     m_factor = false;
     matrix.m_changed = false;
 }
 
 template <class Solver>
 inline void MatrixPartitionedTyingsSolver<Solver>::solve(
     MatrixPartitionedTyings& matrix, const xt::xtensor<double, 2>& b, xt::xtensor<double, 2>& x)
 {
     GOOSEFEM_ASSERT(xt::has_shape(b, {matrix.m_nnode, matrix.m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(x, {matrix.m_nnode, matrix.m_ndim}));
 
     this->factorize(matrix);
 
     Eigen::VectorXd B_u = matrix.AsDofs_u(b);
     Eigen::VectorXd B_d = matrix.AsDofs_d(b);
     Eigen::VectorXd X_p = matrix.AsDofs_p(x);
 
     B_u += matrix.m_Cud * B_d;
 
     Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - matrix.m_ACup * X_p));
     Eigen::VectorXd X_d = matrix.m_Cdu * X_u + matrix.m_Cdp * X_p;
 
     #pragma omp parallel for
     for (size_t m = 0; m < matrix.m_nnode; ++m) {
         for (size_t i = 0; i < matrix.m_ndim; ++i) {
             if (matrix.m_dofs(m, i) < matrix.m_nnu) {
                 x(m, i) = X_u(matrix.m_dofs(m, i));
             }
             else if (matrix.m_dofs(m, i) >= matrix.m_nni) {
                 x(m, i) = X_d(matrix.m_dofs(m, i) - matrix.m_nni);
             }
         }
     }
 }
 
 template <class Solver>
 inline void MatrixPartitionedTyingsSolver<Solver>::solve(
     MatrixPartitionedTyings& matrix, const xt::xtensor<double, 1>& b, xt::xtensor<double, 1>& x)
 {
     GOOSEFEM_ASSERT(b.size() == matrix.m_ndof);
     GOOSEFEM_ASSERT(x.size() == matrix.m_ndof);
 
     this->factorize(matrix);
 
     Eigen::VectorXd B_u = matrix.AsDofs_u(b);
     Eigen::VectorXd B_d = matrix.AsDofs_d(b);
     Eigen::VectorXd X_p = matrix.AsDofs_p(x);
 
     Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - matrix.m_ACup * X_p));
     Eigen::VectorXd X_d = matrix.m_Cdu * X_u + matrix.m_Cdp * X_p;
 
     #pragma omp parallel for
     for (size_t d = 0; d < matrix.m_nnu; ++d) {
         x(matrix.m_iiu(d)) = X_u(d);
     }
 
     #pragma omp parallel for
     for (size_t d = 0; d < matrix.m_nnd; ++d) {
         x(matrix.m_iid(d)) = X_d(d);
     }
 }
 
 template <class Solver>
 inline void MatrixPartitionedTyingsSolver<Solver>::solve_u(
     MatrixPartitionedTyings& matrix,
     const xt::xtensor<double, 1>& b_u,
     const xt::xtensor<double, 1>& b_d,
     const xt::xtensor<double, 1>& x_p,
     xt::xtensor<double, 1>& x_u)
 {
     GOOSEFEM_ASSERT(b_u.size() == matrix.m_nnu);
     GOOSEFEM_ASSERT(b_d.size() == matrix.m_nnd);
     GOOSEFEM_ASSERT(x_p.size() == matrix.m_nnp);
     GOOSEFEM_ASSERT(x_u.size() == matrix.m_nnu);
 
     this->factorize(matrix);
 
     Eigen::Map<Eigen::VectorXd>(x_u.data(), x_u.size()).noalias() = m_solver.solve(Eigen::VectorXd(
         Eigen::Map<const Eigen::VectorXd>(b_u.data(), b_u.size()) -
         matrix.m_ACup * Eigen::Map<const Eigen::VectorXd>(x_p.data(), x_p.size())));
 }
 
 template <class Solver>
 inline xt::xtensor<double, 2> MatrixPartitionedTyingsSolver<Solver>::Solve(
     MatrixPartitionedTyings& matrix,
     const xt::xtensor<double, 2>& b,
     const xt::xtensor<double, 2>& x)
 {
-    xt::xtensor<double, 2> out = x;
-    this->solve(matrix, b, out);
-    return out;
+    xt::xtensor<double, 2> ret = x;
+    this->solve(matrix, b, ret);
+    return ret;
 }
 
 template <class Solver>
 inline xt::xtensor<double, 1> MatrixPartitionedTyingsSolver<Solver>::Solve(
     MatrixPartitionedTyings& matrix,
     const xt::xtensor<double, 1>& b,
     const xt::xtensor<double, 1>& x)
 {
-    xt::xtensor<double, 1> out = x;
-    this->solve(matrix, b, out);
-    return out;
+    xt::xtensor<double, 1> ret = x;
+    this->solve(matrix, b, ret);
+    return ret;
 }
 
 template <class Solver>
 inline xt::xtensor<double, 1> MatrixPartitionedTyingsSolver<Solver>::Solve_u(
     MatrixPartitionedTyings& matrix,
     const xt::xtensor<double, 1>& b_u,
     const xt::xtensor<double, 1>& b_d,
     const xt::xtensor<double, 1>& x_p)
 {
     xt::xtensor<double, 1> x_u = xt::empty<double>({matrix.m_nnu});
     this->solve_u(matrix, b_u, b_d, x_p, x_u);
     return x_u;
 }
 
 } // namespace GooseFEM
 
 #endif
diff --git a/include/GooseFEM/Vector.hpp b/include/GooseFEM/Vector.hpp
index 9ef5eb2..a4813df 100644
--- a/include/GooseFEM/Vector.hpp
+++ b/include/GooseFEM/Vector.hpp
@@ -1,265 +1,271 @@
 /*
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 */
 
 #ifndef GOOSEFEM_VECTOR_HPP
 #define GOOSEFEM_VECTOR_HPP
 
 #include "Vector.h"
 
 namespace GooseFEM {
 
 inline Vector::Vector(const xt::xtensor<size_t, 2>& conn, const xt::xtensor<size_t, 2>& dofs)
     : m_conn(conn), m_dofs(dofs)
 {
     m_nelem = m_conn.shape(0);
     m_nne = m_conn.shape(1);
     m_nnode = m_dofs.shape(0);
     m_ndim = m_dofs.shape(1);
     m_ndof = xt::amax(m_dofs)() + 1;
 
-    GOOSEFEM_ASSERT(xt::amax(m_conn)[0] + 1 == m_nnode);
+    GOOSEFEM_ASSERT(xt::amax(m_conn)() + 1 <= m_nnode);
     GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim);
 }
 
 inline size_t Vector::nelem() const
 {
     return m_nelem;
 }
 
 inline size_t Vector::nne() const
 {
     return m_nne;
 }
 
 inline size_t Vector::nnode() const
 {
     return m_nnode;
 }
 
 inline size_t Vector::ndim() const
 {
     return m_ndim;
 }
 
 inline size_t Vector::ndof() const
 {
     return m_ndof;
 }
 
 inline xt::xtensor<size_t, 2> Vector::dofs() const
 {
     return m_dofs;
 }
 
 inline void
 Vector::copy(const xt::xtensor<double, 2>& nodevec_src, xt::xtensor<double, 2>& nodevec_dest) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec_src, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(nodevec_dest, {m_nnode, m_ndim}));
 
     xt::noalias(nodevec_dest) = nodevec_src;
 }
 
 inline void
 Vector::asDofs(const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 1>& dofval) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
 
+    dofval.fill(0.0);
+
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             dofval(m_dofs(m, i)) = nodevec(m, i);
         }
     }
 }
 
 inline void
 Vector::asDofs(const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 1>& dofval) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
 
+    dofval.fill(0.0);
+
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 dofval(m_dofs(m_conn(e, m), i)) = elemvec(e, m, i);
             }
         }
     }
 }
 
 inline void
 Vector::asNode(const xt::xtensor<double, 1>& dofval, xt::xtensor<double, 2>& nodevec) const
 {
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             nodevec(m, i) = dofval(m_dofs(m, i));
         }
     }
 }
 
 inline void
 Vector::asNode(const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 2>& nodevec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
+    nodevec.fill(0.0);
+
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 nodevec(m_conn(e, m), i) = elemvec(e, m, i);
             }
         }
     }
 }
 
 inline void
 Vector::asElement(const xt::xtensor<double, 1>& dofval, xt::xtensor<double, 3>& elemvec) const
 {
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
 
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 elemvec(e, m, i) = dofval(m_dofs(m_conn(e, m), i));
             }
         }
     }
 }
 
 inline void
 Vector::asElement(const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 3>& elemvec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
 
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 elemvec(e, m, i) = nodevec(m_conn(e, m), i);
             }
         }
     }
 }
 
 inline void
 Vector::assembleDofs(const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 1>& dofval) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
 
     dofval.fill(0.0);
 
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             dofval(m_dofs(m, i)) += nodevec(m, i);
         }
     }
 }
 
 inline void
 Vector::assembleDofs(const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 1>& dofval) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
 
     dofval.fill(0.0);
 
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 dofval(m_dofs(m_conn(e, m), i)) += elemvec(e, m, i);
             }
         }
     }
 }
 
 inline void
 Vector::assembleNode(const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 2>& nodevec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
     xt::xtensor<double, 1> dofval = this->AssembleDofs(elemvec);
     this->asNode(dofval, nodevec);
 }
 
 inline xt::xtensor<double, 1> Vector::AsDofs(const xt::xtensor<double, 2>& nodevec) const
 {
     xt::xtensor<double, 1> dofval = xt::empty<double>({m_ndof});
     this->asDofs(nodevec, dofval);
     return dofval;
 }
 
 inline xt::xtensor<double, 1> Vector::AsDofs(const xt::xtensor<double, 3>& elemvec) const
 {
     xt::xtensor<double, 1> dofval = xt::empty<double>({m_ndof});
     this->asDofs(elemvec, dofval);
     return dofval;
 }
 
 inline xt::xtensor<double, 2> Vector::AsNode(const xt::xtensor<double, 1>& dofval) const
 {
     xt::xtensor<double, 2> nodevec = xt::empty<double>({m_nnode, m_ndim});
     this->asNode(dofval, nodevec);
     return nodevec;
 }
 
 inline xt::xtensor<double, 2> Vector::AsNode(const xt::xtensor<double, 3>& elemvec) const
 {
     xt::xtensor<double, 2> nodevec = xt::empty<double>({m_nnode, m_ndim});
     this->asNode(elemvec, nodevec);
     return nodevec;
 }
 
 inline xt::xtensor<double, 3> Vector::AsElement(const xt::xtensor<double, 1>& dofval) const
 {
     xt::xtensor<double, 3> elemvec = xt::empty<double>({m_nelem, m_nne, m_ndim});
     this->asElement(dofval, elemvec);
     return elemvec;
 }
 
 inline xt::xtensor<double, 3> Vector::AsElement(const xt::xtensor<double, 2>& nodevec) const
 {
     xt::xtensor<double, 3> elemvec = xt::empty<double>({m_nelem, m_nne, m_ndim});
     this->asElement(nodevec, elemvec);
     return elemvec;
 }
 
 inline xt::xtensor<double, 1> Vector::AssembleDofs(const xt::xtensor<double, 2>& nodevec) const
 {
     xt::xtensor<double, 1> dofval = xt::empty<double>({m_ndof});
     this->assembleDofs(nodevec, dofval);
     return dofval;
 }
 
 inline xt::xtensor<double, 1> Vector::AssembleDofs(const xt::xtensor<double, 3>& elemvec) const
 {
     xt::xtensor<double, 1> dofval = xt::empty<double>({m_ndof});
     this->assembleDofs(elemvec, dofval);
     return dofval;
 }
 
 inline xt::xtensor<double, 2> Vector::AssembleNode(const xt::xtensor<double, 3>& elemvec) const
 {
     xt::xtensor<double, 2> nodevec = xt::empty<double>({m_nnode, m_ndim});
     this->assembleNode(elemvec, nodevec);
     return nodevec;
 }
 
 } // namespace GooseFEM
 
 #endif
diff --git a/include/GooseFEM/VectorPartitioned.hpp b/include/GooseFEM/VectorPartitioned.hpp
index 3e6de9c..b64be1c 100644
--- a/include/GooseFEM/VectorPartitioned.hpp
+++ b/include/GooseFEM/VectorPartitioned.hpp
@@ -1,689 +1,704 @@
 /*
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 */
 
 #ifndef GOOSEFEM_VECTORPARTITIONED_HPP
 #define GOOSEFEM_VECTORPARTITIONED_HPP
 
 #include "Mesh.h"
 #include "VectorPartitioned.h"
 
 namespace GooseFEM {
 
 inline VectorPartitioned::VectorPartitioned(
     const xt::xtensor<size_t, 2>& conn,
     const xt::xtensor<size_t, 2>& dofs,
     const xt::xtensor<size_t, 1>& iip)
     : m_conn(conn), m_dofs(dofs), m_iip(iip)
 {
     m_nelem = m_conn.shape(0);
     m_nne = m_conn.shape(1);
     m_nnode = m_dofs.shape(0);
     m_ndim = m_dofs.shape(1);
     m_iiu = xt::setdiff1d(dofs, iip);
     m_ndof = xt::amax(m_dofs)() + 1;
     m_nnp = m_iip.size();
     m_nnu = m_iiu.size();
     m_part = Mesh::Reorder({m_iiu, m_iip}).get(m_dofs);
 
-    GOOSEFEM_ASSERT(xt::amax(m_conn)[0] + 1 == m_nnode);
+    GOOSEFEM_ASSERT(xt::amax(m_conn)() + 1 <= m_nnode);
     GOOSEFEM_ASSERT(xt::amax(m_iip)() <= xt::amax(m_dofs)());
     GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim);
+    GOOSEFEM_ASSERT(m_ndof == m_nnu + m_nnp);
 }
 
 inline size_t VectorPartitioned::nelem() const
 {
     return m_nelem;
 }
 
 inline size_t VectorPartitioned::nne() const
 {
     return m_nne;
 }
 
 inline size_t VectorPartitioned::nnode() const
 {
     return m_nnode;
 }
 
 inline size_t VectorPartitioned::ndim() const
 {
     return m_ndim;
 }
 
 inline size_t VectorPartitioned::ndof() const
 {
     return m_ndof;
 }
 
 inline size_t VectorPartitioned::nnu() const
 {
     return m_nnu;
 }
 
 inline size_t VectorPartitioned::nnp() const
 {
     return m_nnp;
 }
 
 inline xt::xtensor<size_t, 2> VectorPartitioned::dofs() const
 {
     return m_dofs;
 }
 
 inline xt::xtensor<size_t, 1> VectorPartitioned::iiu() const
 {
     return m_iiu;
 }
 
 inline xt::xtensor<size_t, 1> VectorPartitioned::iip() const
 {
     return m_iip;
 }
 
 inline void VectorPartitioned::copy(
     const xt::xtensor<double, 2>& nodevec_src, xt::xtensor<double, 2>& nodevec_dest) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec_src, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(nodevec_dest, {m_nnode, m_ndim}));
 
     xt::noalias(nodevec_dest) = nodevec_src;
 }
 
 inline void VectorPartitioned::copy_u(
     const xt::xtensor<double, 2>& nodevec_src, xt::xtensor<double, 2>& nodevec_dest) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec_src, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(nodevec_dest, {m_nnode, m_ndim}));
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             if (m_part(m, i) < m_nnu) {
                 nodevec_dest(m, i) = nodevec_src(m, i);
             }
         }
     }
 }
 
 inline void VectorPartitioned::copy_p(
     const xt::xtensor<double, 2>& nodevec_src, xt::xtensor<double, 2>& nodevec_dest) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec_src, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(nodevec_dest, {m_nnode, m_ndim}));
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             if (m_part(m, i) >= m_nnu) {
                 nodevec_dest(m, i) = nodevec_src(m, i);
             }
         }
     }
 }
 
 inline void VectorPartitioned::asDofs(
     const xt::xtensor<double, 1>& dofval_u,
     const xt::xtensor<double, 1>& dofval_p,
     xt::xtensor<double, 1>& dofval) const
 {
     GOOSEFEM_ASSERT(dofval_u.size() == m_nnu);
     GOOSEFEM_ASSERT(dofval_p.size() == m_nnp);
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnu; ++d) {
         dofval(m_iiu(d)) = dofval_u(d);
     }
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnp; ++d) {
         dofval(m_iip(d)) = dofval_p(d);
     }
 }
 
 inline void VectorPartitioned::asDofs(
     const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 1>& dofval) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
 
+    dofval.fill(0.0);
+
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             dofval(m_dofs(m, i)) = nodevec(m, i);
         }
     }
 }
 
 inline void VectorPartitioned::asDofs_u(
     const xt::xtensor<double, 1>& dofval, xt::xtensor<double, 1>& dofval_u) const
 {
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
     GOOSEFEM_ASSERT(dofval_u.size() == m_nnu);
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnu; ++d) {
         dofval_u(d) = dofval(m_iiu(d));
     }
 }
 
 inline void VectorPartitioned::asDofs_u(
     const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 1>& dofval_u) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(dofval_u.size() == m_nnu);
 
+    dofval_u.fill(0.0);
+
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             if (m_part(m, i) < m_nnu) {
                 dofval_u(m_part(m, i)) = nodevec(m, i);
             }
         }
     }
 }
 
 inline void VectorPartitioned::asDofs_p(
     const xt::xtensor<double, 1>& dofval, xt::xtensor<double, 1>& dofval_p) const
 {
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
     GOOSEFEM_ASSERT(dofval_p.size() == m_nnp);
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnp; ++d) {
         dofval_p(d) = dofval(m_iip(d));
     }
 }
 
 inline void VectorPartitioned::asDofs_p(
     const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 1>& dofval_p) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(dofval_p.size() == m_nnp);
 
+    dofval_p.fill(0.0);
+
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             if (m_part(m, i) >= m_nnu) {
                 dofval_p(m_part(m, i) - m_nnu) = nodevec(m, i);
             }
         }
     }
 }
 
 inline void VectorPartitioned::asDofs(
     const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 1>& dofval) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
 
+    dofval.fill(0.0);
+
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 dofval(m_dofs(m_conn(e, m), i)) = elemvec(e, m, i);
             }
         }
     }
 }
 
 inline void VectorPartitioned::asDofs_u(
     const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 1>& dofval_u) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
     GOOSEFEM_ASSERT(dofval_u.size() == m_nnu);
 
+    dofval_u.fill(0.0);
+
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 if (m_part(m_conn(e, m), i) < m_nnu) {
                     dofval_u(m_part(m_conn(e, m), i)) = elemvec(e, m, i);
                 }
             }
         }
     }
 }
 
 inline void VectorPartitioned::asDofs_p(
     const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 1>& dofval_p) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
     GOOSEFEM_ASSERT(dofval_p.size() == m_nnp);
 
+    dofval_p.fill(0.0);
+
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 if (m_part(m_conn(e, m), i) >= m_nnu) {
                     dofval_p(m_part(m_conn(e, m), i) - m_nnu) = elemvec(e, m, i);
                 }
             }
         }
     }
 }
 
 inline void VectorPartitioned::asNode(
     const xt::xtensor<double, 1>& dofval, xt::xtensor<double, 2>& nodevec) const
 {
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             nodevec(m, i) = dofval(m_dofs(m, i));
         }
     }
 }
 
 inline void VectorPartitioned::asNode(
     const xt::xtensor<double, 1>& dofval_u,
     const xt::xtensor<double, 1>& dofval_p,
     xt::xtensor<double, 2>& nodevec) const
 {
     GOOSEFEM_ASSERT(dofval_u.size() == m_nnu);
     GOOSEFEM_ASSERT(dofval_p.size() == m_nnp);
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             if (m_part(m, i) < m_nnu) {
                 nodevec(m, i) = dofval_u(m_part(m, i));
             }
             else {
                 nodevec(m, i) = dofval_p(m_part(m, i) - m_nnu);
             }
         }
     }
 }
 
 inline void VectorPartitioned::asNode(
     const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 2>& nodevec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
+    nodevec.fill(0.0);
+
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 nodevec(m_conn(e, m), i) = elemvec(e, m, i);
             }
         }
     }
 }
 
 inline void VectorPartitioned::asElement(
     const xt::xtensor<double, 1>& dofval, xt::xtensor<double, 3>& elemvec) const
 {
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
 
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 elemvec(e, m, i) = dofval(m_dofs(m_conn(e, m), i));
             }
         }
     }
 }
 
 inline void VectorPartitioned::asElement(
     const xt::xtensor<double, 1>& dofval_u,
     const xt::xtensor<double, 1>& dofval_p,
     xt::xtensor<double, 3>& elemvec) const
 {
     GOOSEFEM_ASSERT(dofval_u.size() == m_nnu);
     GOOSEFEM_ASSERT(dofval_p.size() == m_nnp);
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
 
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 if (m_part(m_conn(e, m), i) < m_nnu) {
                     elemvec(e, m, i) = dofval_u(m_part(m_conn(e, m), i));
                 }
                 else {
                     elemvec(e, m, i) = dofval_p(m_part(m_conn(e, m), i) - m_nnu);
                 }
             }
         }
     }
 }
 
 inline void VectorPartitioned::asElement(
     const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 3>& elemvec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
 
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 elemvec(e, m, i) = nodevec(m_conn(e, m), i);
             }
         }
     }
 }
 
 inline void VectorPartitioned::assembleDofs(
     const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 1>& dofval) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
 
     dofval.fill(0.0);
 
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             dofval(m_dofs(m, i)) += nodevec(m, i);
         }
     }
 }
 
 inline void VectorPartitioned::assembleDofs_u(
     const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 1>& dofval_u) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(dofval_u.size() == m_nnu);
 
     dofval_u.fill(0.0);
 
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             if (m_part(m, i) < m_nnu) {
                 dofval_u(m_part(m, i)) += nodevec(m, i);
             }
         }
     }
 }
 
 inline void VectorPartitioned::assembleDofs_p(
     const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 1>& dofval_p) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(dofval_p.size() == m_nnp);
 
     dofval_p.fill(0.0);
 
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             if (m_part(m, i) >= m_nnu) {
                 dofval_p(m_part(m, i) - m_nnu) += nodevec(m, i);
             }
         }
     }
 }
 
 inline void VectorPartitioned::assembleDofs(
     const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 1>& dofval) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
 
     dofval.fill(0.0);
 
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 dofval(m_dofs(m_conn(e, m), i)) += elemvec(e, m, i);
             }
         }
     }
 }
 
 inline void VectorPartitioned::assembleDofs_u(
     const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 1>& dofval_u) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
     GOOSEFEM_ASSERT(dofval_u.size() == m_nnu);
 
     dofval_u.fill(0.0);
 
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 if (m_part(m_conn(e, m), i) < m_nnu) {
                     dofval_u(m_part(m_conn(e, m), i)) += elemvec(e, m, i);
                 }
             }
         }
     }
 }
 
 inline void VectorPartitioned::assembleDofs_p(
     const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 1>& dofval_p) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
     GOOSEFEM_ASSERT(dofval_p.size() == m_nnp);
 
     dofval_p.fill(0.0);
 
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 if (m_part(m_conn(e, m), i) >= m_nnu) {
                     dofval_p(m_part(m_conn(e, m), i) - m_nnu) += elemvec(e, m, i);
                 }
             }
         }
     }
 }
 
 inline void VectorPartitioned::assembleNode(
     const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 2>& nodevec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
     xt::xtensor<double, 1> dofval = this->AssembleDofs(elemvec);
     this->asNode(dofval, nodevec);
 }
 
 inline xt::xtensor<double, 1> VectorPartitioned::AsDofs(
     const xt::xtensor<double, 1>& dofval_u, const xt::xtensor<double, 1>& dofval_p) const
 {
     xt::xtensor<double, 1> dofval = xt::empty<double>({m_ndof});
     this->asDofs(dofval_u, dofval_p, dofval);
     return dofval;
 }
 
 inline xt::xtensor<double, 1> VectorPartitioned::AsDofs(const xt::xtensor<double, 2>& nodevec) const
 {
     xt::xtensor<double, 1> dofval = xt::empty<double>({m_ndof});
     this->asDofs(nodevec, dofval);
     return dofval;
 }
 
 inline xt::xtensor<double, 1>
 VectorPartitioned::AsDofs_u(const xt::xtensor<double, 1>& dofval) const
 {
     xt::xtensor<double, 1> dofval_u = xt::empty<double>({m_nnu});
     this->asDofs_u(dofval, dofval_u);
     return dofval_u;
 }
 
 inline xt::xtensor<double, 1>
 VectorPartitioned::AsDofs_u(const xt::xtensor<double, 2>& nodevec) const
 {
     xt::xtensor<double, 1> dofval_u = xt::empty<double>({m_nnu});
     this->asDofs_u(nodevec, dofval_u);
     return dofval_u;
 }
 
 inline xt::xtensor<double, 1>
 VectorPartitioned::AsDofs_p(const xt::xtensor<double, 1>& dofval) const
 {
     xt::xtensor<double, 1> dofval_p = xt::empty<double>({m_nnp});
     this->asDofs_p(dofval, dofval_p);
     return dofval_p;
 }
 
 inline xt::xtensor<double, 1>
 VectorPartitioned::AsDofs_p(const xt::xtensor<double, 2>& nodevec) const
 {
     xt::xtensor<double, 1> dofval_p = xt::empty<double>({m_nnp});
     this->asDofs_p(nodevec, dofval_p);
     return dofval_p;
 }
 
 inline xt::xtensor<double, 1> VectorPartitioned::AsDofs(const xt::xtensor<double, 3>& elemvec) const
 {
     xt::xtensor<double, 1> dofval = xt::empty<double>({m_ndof});
     this->asDofs(elemvec, dofval);
     return dofval;
 }
 
 inline xt::xtensor<double, 1>
 VectorPartitioned::AsDofs_u(const xt::xtensor<double, 3>& elemvec) const
 {
     xt::xtensor<double, 1> dofval_u = xt::empty<double>({m_nnu});
     this->asDofs_u(elemvec, dofval_u);
     return dofval_u;
 }
 
 inline xt::xtensor<double, 1>
 VectorPartitioned::AsDofs_p(const xt::xtensor<double, 3>& elemvec) const
 {
     xt::xtensor<double, 1> dofval_p = xt::empty<double>({m_nnp});
     this->asDofs_p(elemvec, dofval_p);
     return dofval_p;
 }
 
 inline xt::xtensor<double, 2> VectorPartitioned::AsNode(const xt::xtensor<double, 1>& dofval) const
 {
     xt::xtensor<double, 2> nodevec = xt::empty<double>({m_nnode, m_ndim});
     this->asNode(dofval, nodevec);
     return nodevec;
 }
 
 inline xt::xtensor<double, 2> VectorPartitioned::AsNode(
     const xt::xtensor<double, 1>& dofval_u, const xt::xtensor<double, 1>& dofval_p) const
 {
     xt::xtensor<double, 2> nodevec = xt::empty<double>({m_nnode, m_ndim});
     this->asNode(dofval_u, dofval_p, nodevec);
     return nodevec;
 }
 
 inline xt::xtensor<double, 2> VectorPartitioned::AsNode(const xt::xtensor<double, 3>& elemvec) const
 {
     xt::xtensor<double, 2> nodevec = xt::empty<double>({m_nnode, m_ndim});
     this->asNode(elemvec, nodevec);
     return nodevec;
 }
 
 inline xt::xtensor<double, 3>
 VectorPartitioned::AsElement(const xt::xtensor<double, 1>& dofval) const
 {
     xt::xtensor<double, 3> elemvec = xt::empty<double>({m_nelem, m_nne, m_ndim});
     this->asElement(dofval, elemvec);
     return elemvec;
 }
 
 inline xt::xtensor<double, 3> VectorPartitioned::AsElement(
     const xt::xtensor<double, 1>& dofval_u, const xt::xtensor<double, 1>& dofval_p) const
 {
     xt::xtensor<double, 3> elemvec = xt::empty<double>({m_nelem, m_nne, m_ndim});
     this->asElement(dofval_u, dofval_p, elemvec);
     return elemvec;
 }
 
 inline xt::xtensor<double, 3>
 VectorPartitioned::AsElement(const xt::xtensor<double, 2>& nodevec) const
 {
     xt::xtensor<double, 3> elemvec = xt::empty<double>({m_nelem, m_nne, m_ndim});
     this->asElement(nodevec, elemvec);
     return elemvec;
 }
 
 inline xt::xtensor<double, 1>
 VectorPartitioned::AssembleDofs(const xt::xtensor<double, 2>& nodevec) const
 {
     xt::xtensor<double, 1> dofval = xt::empty<double>({m_ndof});
     this->assembleDofs(nodevec, dofval);
     return dofval;
 }
 
 inline xt::xtensor<double, 1>
 VectorPartitioned::AssembleDofs_u(const xt::xtensor<double, 2>& nodevec) const
 {
     xt::xtensor<double, 1> dofval_u = xt::empty<double>({m_nnu});
     this->assembleDofs_u(nodevec, dofval_u);
     return dofval_u;
 }
 
 inline xt::xtensor<double, 1>
 VectorPartitioned::AssembleDofs_p(const xt::xtensor<double, 2>& nodevec) const
 {
     xt::xtensor<double, 1> dofval_p = xt::empty<double>({m_nnp});
     this->assembleDofs_p(nodevec, dofval_p);
     return dofval_p;
 }
 
 inline xt::xtensor<double, 1>
 VectorPartitioned::AssembleDofs(const xt::xtensor<double, 3>& elemvec) const
 {
     xt::xtensor<double, 1> dofval = xt::empty<double>({m_ndof});
     this->assembleDofs(elemvec, dofval);
     return dofval;
 }
 
 inline xt::xtensor<double, 1>
 VectorPartitioned::AssembleDofs_u(const xt::xtensor<double, 3>& elemvec) const
 {
     xt::xtensor<double, 1> dofval_u = xt::empty<double>({m_nnu});
     this->assembleDofs_u(elemvec, dofval_u);
     return dofval_u;
 }
 
 inline xt::xtensor<double, 1>
 VectorPartitioned::AssembleDofs_p(const xt::xtensor<double, 3>& elemvec) const
 {
     xt::xtensor<double, 1> dofval_p = xt::empty<double>({m_nnp});
     this->assembleDofs_p(elemvec, dofval_p);
     return dofval_p;
 }
 
 inline xt::xtensor<double, 2>
 VectorPartitioned::AssembleNode(const xt::xtensor<double, 3>& elemvec) const
 {
     xt::xtensor<double, 2> nodevec = xt::empty<double>({m_nnode, m_ndim});
     this->assembleNode(elemvec, nodevec);
     return nodevec;
 }
 
 inline xt::xtensor<double, 2> VectorPartitioned::Copy(
     const xt::xtensor<double, 2>& nodevec_src, const xt::xtensor<double, 2>& nodevec_dest) const
 {
     xt::xtensor<double, 2> out = nodevec_dest;
     this->copy(nodevec_src, out);
     return out;
 }
 
 inline xt::xtensor<double, 2> VectorPartitioned::Copy_u(
     const xt::xtensor<double, 2>& nodevec_src, const xt::xtensor<double, 2>& nodevec_dest) const
 {
     xt::xtensor<double, 2> out = nodevec_dest;
     this->copy_u(nodevec_src, out);
     return out;
 }
 
 inline xt::xtensor<double, 2> VectorPartitioned::Copy_p(
     const xt::xtensor<double, 2>& nodevec_src, const xt::xtensor<double, 2>& nodevec_dest) const
 {
     xt::xtensor<double, 2> out = nodevec_dest;
     this->copy_p(nodevec_src, out);
     return out;
 }
 
 } // namespace GooseFEM
 
 #endif
diff --git a/include/GooseFEM/VectorPartitionedTyings.hpp b/include/GooseFEM/VectorPartitionedTyings.hpp
index 09ce044..dde0412 100644
--- a/include/GooseFEM/VectorPartitionedTyings.hpp
+++ b/include/GooseFEM/VectorPartitionedTyings.hpp
@@ -1,275 +1,278 @@
 /*
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 */
 
 #ifndef GOOSEFEM_VECTORPARTITIONEDTYINGS_HPP
 #define GOOSEFEM_VECTORPARTITIONEDTYINGS_HPP
 
 #include "VectorPartitionedTyings.h"
 
 namespace GooseFEM {
 
 inline VectorPartitionedTyings::VectorPartitionedTyings(
     const xt::xtensor<size_t, 2>& conn,
     const xt::xtensor<size_t, 2>& dofs,
     const Eigen::SparseMatrix<double>& Cdu,
     const Eigen::SparseMatrix<double>& Cdp,
     const Eigen::SparseMatrix<double>& Cdi)
     : m_conn(conn), m_dofs(dofs), m_Cdu(Cdu), m_Cdp(Cdp), m_Cdi(Cdi)
 {
     GOOSEFEM_ASSERT(Cdu.rows() == Cdp.rows());
     GOOSEFEM_ASSERT(Cdi.rows() == Cdp.rows());
 
     m_nnu = static_cast<size_t>(m_Cdu.cols());
     m_nnp = static_cast<size_t>(m_Cdp.cols());
     m_nnd = static_cast<size_t>(m_Cdp.rows());
     m_nni = m_nnu + m_nnp;
     m_ndof = m_nni + m_nnd;
     m_iiu = xt::arange<size_t>(m_nnu);
     m_iip = xt::arange<size_t>(m_nnu, m_nnu + m_nnp);
     m_iid = xt::arange<size_t>(m_nni, m_nni + m_nnd);
     m_nelem = m_conn.shape(0);
     m_nne = m_conn.shape(1);
     m_nnode = m_dofs.shape(0);
     m_ndim = m_dofs.shape(1);
     m_Cud = m_Cdu.transpose();
     m_Cpd = m_Cdp.transpose();
     m_Cid = m_Cdi.transpose();
 
     GOOSEFEM_ASSERT(static_cast<size_t>(m_Cdi.cols()) == m_nni);
     GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim);
     GOOSEFEM_ASSERT(m_ndof == xt::amax(m_dofs)() + 1);
 }
 
 inline size_t VectorPartitionedTyings::nelem() const
 {
     return m_nelem;
 }
 
 inline size_t VectorPartitionedTyings::nne() const
 {
     return m_nne;
 }
 
 inline size_t VectorPartitionedTyings::nnode() const
 {
     return m_nnode;
 }
 
 inline size_t VectorPartitionedTyings::ndim() const
 {
     return m_ndim;
 }
 
 inline size_t VectorPartitionedTyings::ndof() const
 {
     return m_ndof;
 }
 
 inline size_t VectorPartitionedTyings::nnu() const
 {
     return m_nnu;
 }
 
 inline size_t VectorPartitionedTyings::nnp() const
 {
     return m_nnp;
 }
 
 inline size_t VectorPartitionedTyings::nni() const
 {
     return m_nni;
 }
 
 inline size_t VectorPartitionedTyings::nnd() const
 {
     return m_nnd;
 }
 
 inline xt::xtensor<size_t, 2> VectorPartitionedTyings::dofs() const
 {
     return m_dofs;
 }
 
 inline xt::xtensor<size_t, 1> VectorPartitionedTyings::iiu() const
 {
     return m_iiu;
 }
 
 inline xt::xtensor<size_t, 1> VectorPartitionedTyings::iip() const
 {
     return m_iip;
 }
 
 inline xt::xtensor<size_t, 1> VectorPartitionedTyings::iii() const
 {
     return xt::arange<size_t>(m_nni);
 }
 
 inline xt::xtensor<size_t, 1> VectorPartitionedTyings::iid() const
 {
     return m_iid;
 }
 
 inline void VectorPartitionedTyings::copy_p(
     const xt::xtensor<double, 1>& dofval_src, xt::xtensor<double, 1>& dofval_dest) const
 {
     GOOSEFEM_ASSERT(dofval_src.size() == m_ndof || dofval_src.size() == m_nni);
     GOOSEFEM_ASSERT(dofval_dest.size() == m_ndof || dofval_dest.size() == m_nni);
 
     #pragma omp parallel for
     for (size_t i = m_nnu; i < m_nni; ++i) {
         dofval_dest(i) = dofval_src(i);
     }
 }
 
 inline void VectorPartitionedTyings::asDofs_i(
     const xt::xtensor<double, 2>& nodevec,
     xt::xtensor<double, 1>& dofval_i,
     bool apply_tyings) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(dofval_i.size() == m_nni);
 
+    dofval_i.fill(0.0);
+
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             if (m_dofs(m, i) < m_nni) {
                 dofval_i(m_dofs(m, i)) = nodevec(m, i);
             }
         }
     }
 
-    if (!apply_tyings)
+    if (!apply_tyings) {
         return;
+    }
 
     Eigen::VectorXd Dofval_d = this->Eigen_asDofs_d(nodevec);
     Eigen::VectorXd Dofval_i = m_Cid * Dofval_d;
 
     #pragma omp parallel for
     for (size_t i = 0; i < m_nni; ++i) {
         dofval_i(i) += Dofval_i(i);
     }
 }
 
 inline void VectorPartitionedTyings::asNode(
     const xt::xtensor<double, 1>& dofval, xt::xtensor<double, 2>& nodevec) const
 {
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             nodevec(m, i) = dofval(m_dofs(m, i));
         }
     }
 }
 
 inline void VectorPartitionedTyings::asElement(
     const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 3>& elemvec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
 
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 elemvec(e, m, i) = nodevec(m_conn(e, m), i);
             }
         }
     }
 }
 
 inline void VectorPartitionedTyings::assembleDofs(
     const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 1>& dofval) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
 
     dofval.fill(0.0);
 
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 dofval(m_dofs(m_conn(e, m), i)) += elemvec(e, m, i);
             }
         }
     }
 }
 
 inline void VectorPartitionedTyings::assembleNode(
     const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 2>& nodevec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
     xt::xtensor<double, 1> dofval = this->AssembleDofs(elemvec);
     this->asNode(dofval, nodevec);
 }
 
 inline xt::xtensor<double, 1>
 VectorPartitionedTyings::AsDofs_i(const xt::xtensor<double, 2>& nodevec) const
 {
     xt::xtensor<double, 1> dofval = xt::empty<double>({m_nni});
     this->asDofs_i(nodevec, dofval);
     return dofval;
 }
 
 inline xt::xtensor<double, 2>
 VectorPartitionedTyings::AsNode(const xt::xtensor<double, 1>& dofval) const
 {
     xt::xtensor<double, 2> nodevec = xt::empty<double>({m_nnode, m_ndim});
     this->asNode(dofval, nodevec);
     return nodevec;
 }
 
 inline xt::xtensor<double, 3>
 VectorPartitionedTyings::AsElement(const xt::xtensor<double, 2>& nodevec) const
 {
     xt::xtensor<double, 3> elemvec = xt::empty<double>({m_nelem, m_nne, m_ndim});
     this->asElement(nodevec, elemvec);
     return elemvec;
 }
 
 inline xt::xtensor<double, 1>
 VectorPartitionedTyings::AssembleDofs(const xt::xtensor<double, 3>& elemvec) const
 {
     xt::xtensor<double, 1> dofval = xt::empty<double>({m_ndof});
     this->assembleDofs(elemvec, dofval);
     return dofval;
 }
 
 inline xt::xtensor<double, 2>
 VectorPartitionedTyings::AssembleNode(const xt::xtensor<double, 3>& elemvec) const
 {
     xt::xtensor<double, 2> nodevec = xt::empty<double>({m_nnode, m_ndim});
     this->assembleNode(elemvec, nodevec);
     return nodevec;
 }
 
 inline Eigen::VectorXd
 VectorPartitionedTyings::Eigen_asDofs_d(const xt::xtensor<double, 2>& nodevec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
     Eigen::VectorXd dofval_d(m_nnd, 1);
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             if (m_dofs(m, i) >= m_nni) {
                 dofval_d(m_dofs(m, i) - m_nni) = nodevec(m, i);
             }
         }
     }
 
     return dofval_d;
 }
 
 } // namespace GooseFEM
 
 #endif