diff --git a/include/GooseFEM/Matrix.h b/include/GooseFEM/Matrix.h
index 1128efa..6798c0c 100644
--- a/include/GooseFEM/Matrix.h
+++ b/include/GooseFEM/Matrix.h
@@ -1,121 +1,121 @@
 /*
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 */
 
 #ifndef GOOSEFEM_MATRIX_H
 #define GOOSEFEM_MATRIX_H
 
 #include "config.h"
 
 #include <Eigen/Eigen>
 #include <Eigen/Sparse>
 #include <Eigen/SparseCholesky>
 
 namespace GooseFEM {
 
 // forward declaration
 template <class> class MatrixSolver;
 
 class Matrix {
 public:
     // Constructors
     Matrix() = default;
     Matrix(const xt::xtensor<size_t, 2>& conn, const xt::xtensor<size_t, 2>& dofs);
 
     // Dimensions
     size_t nelem() const; // number of elements
     size_t nne() const;   // number of nodes per element
     size_t nnode() const; // number of nodes
     size_t ndim() const;  // number of dimensions
     size_t ndof() const;  // number of DOFs
 
     // DOF lists
     xt::xtensor<size_t, 2> dofs() const; // DOFs
 
     // Assemble from matrices stored per element [nelem, nne*ndim, nne*ndim]
     void assemble(const xt::xtensor<double, 3>& elemmat);
 
     // Overwrite with a dense (sub-) matrix
     void set(
         const xt::xtensor<size_t, 1>& rows,
         const xt::xtensor<size_t, 1>& cols,
         const xt::xtensor<double, 2>& matrix);
 
     // Add a dense (sub-) matrix to the current matrix
     void add(
         const xt::xtensor<size_t, 1>& rows,
         const xt::xtensor<size_t, 1>& cols,
         const xt::xtensor<double, 2>& matrix);
 
     // Return as dense matrix
     void todense(xt::xtensor<double, 2>& ret) const;
 
     // Dot-product:
     // b_i = A_ij * x_j
     void dot(const xt::xtensor<double, 2>& x, xt::xtensor<double, 2>& b) const;
     void dot(const xt::xtensor<double, 1>& x, xt::xtensor<double, 1>& b) const;
 
     // Auto-allocation of the functions above
     xt::xtensor<double, 2> Todense() const;
     xt::xtensor<double, 2> Dot(const xt::xtensor<double, 2>& x) const;
     xt::xtensor<double, 1> Dot(const xt::xtensor<double, 1>& x) const;
 
 private:
     // The matrix
     Eigen::SparseMatrix<double> m_A;
 
     // Matrix entries
     std::vector<Eigen::Triplet<double>> m_T;
 
     // Signal changes to data
     bool m_changed = true;
 
     // Bookkeeping
     xt::xtensor<size_t, 2> m_conn; // connectivity [nelem, nne]
     xt::xtensor<size_t, 2> m_dofs; // DOF-numbers per node [nnode, ndim]
 
     // Dimensions
     size_t m_nelem; // number of elements
     size_t m_nne;   // number of nodes per element
     size_t m_nnode; // number of nodes
     size_t m_ndim;  // number of dimensions
     size_t m_ndof;  // number of DOFs
 
     // grant access to solver class
     template <class> friend class MatrixSolver;
 
     // Convert arrays (Eigen version of Vector, which contains public functions)
     Eigen::VectorXd AsDofs(const xt::xtensor<double, 2>& nodevec) const;
 
     void asNode(const Eigen::VectorXd& dofval, xt::xtensor<double, 2>& nodevec) const;
 };
 
 
 template <class Solver = Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>>>
 class MatrixSolver {
 public:
     // Constructors
     MatrixSolver() = default;
 
     // Solve
     // x_u = A_uu \ (b_u - A_up * x_p)
     void solve(Matrix& matrix, const xt::xtensor<double, 2>& b, xt::xtensor<double, 2>& x);
     void solve(Matrix& matrix, const xt::xtensor<double, 1>& b, xt::xtensor<double, 1>& x);
 
     // Auto-allocation of the functions above
     xt::xtensor<double, 2> Solve(Matrix& matrix, const xt::xtensor<double, 2>& b);
     xt::xtensor<double, 1> Solve(Matrix& matrix, const xt::xtensor<double, 1>& b);
 
 private:
     Solver m_solver; // solver
-    bool m_factor = false; // signal to force factorization
+    bool m_factor = true; // signal to force factorization
     void factorize(Matrix& matrix); // compute inverse (evaluated by "solve")
 };
 
 } // namespace GooseFEM
 
 #include "Matrix.hpp"
 
 #endif
diff --git a/include/GooseFEM/MatrixDiagonal.h b/include/GooseFEM/MatrixDiagonal.h
index baaddde..ba276bf 100644
--- a/include/GooseFEM/MatrixDiagonal.h
+++ b/include/GooseFEM/MatrixDiagonal.h
@@ -1,83 +1,83 @@
 /*
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 */
 
 #ifndef GOOSEFEM_MATRIXDIAGONAL_H
 #define GOOSEFEM_MATRIXDIAGONAL_H
 
 #include "config.h"
 
 namespace GooseFEM {
 
 class MatrixDiagonal {
 public:
     // Constructors
     MatrixDiagonal() = default;
     MatrixDiagonal(const xt::xtensor<size_t, 2>& conn, const xt::xtensor<size_t, 2>& dofs);
 
     // Dimensions
     size_t nelem() const; // number of elements
     size_t nne() const;   // number of nodes per element
     size_t nnode() const; // number of nodes
     size_t ndim() const;  // number of dimensions
     size_t ndof() const;  // number of DOFs
 
     // DOF lists
     xt::xtensor<size_t, 2> dofs() const; // DOFs
 
     // Set matrix components
     void set(const xt::xtensor<double, 1>& A);
 
     // assemble from matrices stored per element [nelem, nne*ndim, nne*ndim]
     // WARNING: ignores any off-diagonal terms
     void assemble(const xt::xtensor<double, 3>& elemmat);
 
     // Dot-product:
     // b_i = A_ij * x_j
     void dot(const xt::xtensor<double, 2>& x, xt::xtensor<double, 2>& b) const;
     void dot(const xt::xtensor<double, 1>& x, xt::xtensor<double, 1>& b) const;
 
     // Solve:
     // x = A \ b
     void solve(const xt::xtensor<double, 2>& b, xt::xtensor<double, 2>& x);
     void solve(const xt::xtensor<double, 1>& b, xt::xtensor<double, 1>& x);
 
     // Return matrix as diagonal matrix (column)
     xt::xtensor<double, 1> AsDiagonal() const;
 
     // Auto-allocation of the functions above
     xt::xtensor<double, 2> Dot(const xt::xtensor<double, 2>& x) const;
     xt::xtensor<double, 1> Dot(const xt::xtensor<double, 1>& x) const;
     xt::xtensor<double, 2> Solve(const xt::xtensor<double, 2>& b);
     xt::xtensor<double, 1> Solve(const xt::xtensor<double, 1>& b);
 
 private:
     // The diagonal matrix, and its inverse (re-used to solve different RHS)
     xt::xtensor<double, 1> m_A;
     xt::xtensor<double, 1> m_inv;
 
     // Signal changes to data compare to the last inverse
-    bool m_factor = false;
+    bool m_factor = true;
 
     // Bookkeeping
     xt::xtensor<size_t, 2> m_conn; // connectivity [nelem, nne]
     xt::xtensor<size_t, 2> m_dofs; // DOF-numbers per node [nnode, ndim]
 
     // Dimensions
     size_t m_nelem; // number of elements
     size_t m_nne;   // number of nodes per element
     size_t m_nnode; // number of nodes
     size_t m_ndim;  // number of dimensions
     size_t m_ndof;  // number of DOFs
 
     // Compute inverse (automatically evaluated by "solve")
     void factorize();
 };
 
 } // namespace GooseFEM
 
 #include "MatrixDiagonal.hpp"
 
 #endif
diff --git a/include/GooseFEM/MatrixDiagonalPartitioned.h b/include/GooseFEM/MatrixDiagonalPartitioned.h
index 3b9cde3..041ad9e 100644
--- a/include/GooseFEM/MatrixDiagonalPartitioned.h
+++ b/include/GooseFEM/MatrixDiagonalPartitioned.h
@@ -1,142 +1,142 @@
 /*
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 */
 
 #ifndef GOOSEFEM_MATRIXDIAGONALPARTITIONED_H
 #define GOOSEFEM_MATRIXDIAGONALPARTITIONED_H
 
 #include "config.h"
 
 namespace GooseFEM {
 
 class MatrixDiagonalPartitioned {
 public:
     // Constructors
     MatrixDiagonalPartitioned() = default;
 
     MatrixDiagonalPartitioned(
         const xt::xtensor<size_t, 2>& conn,
         const xt::xtensor<size_t, 2>& dofs,
         const xt::xtensor<size_t, 1>& iip);
 
     // Dimensions
     size_t nelem() const; // number of elements
     size_t nne() const;   // number of nodes per element
     size_t nnode() const; // number of nodes
     size_t ndim() const;  // number of dimensions
     size_t ndof() const;  // number of DOFs
     size_t nnu() const;   // number of unknown DOFs
     size_t nnp() const;   // number of prescribed DOFs
 
     // DOF lists
     xt::xtensor<size_t, 2> dofs() const; // DOFs
     xt::xtensor<size_t, 1> iiu() const;  // unknown DOFs
     xt::xtensor<size_t, 1> iip() const;  // prescribed DOFs
 
     // Assemble from matrices stored per element [nelem, nne*ndim, nne*ndim]
     // WARNING: ignores any off-diagonal terms
     void assemble(const xt::xtensor<double, 3>& elemmat);
 
     // Dot-product:
     // b_i = A_ij * x_j
     void dot(const xt::xtensor<double, 2>& x, xt::xtensor<double, 2>& b) const;
     void dot(const xt::xtensor<double, 1>& x, xt::xtensor<double, 1>& b) const;
 
     void dot_u(
         const xt::xtensor<double, 1>& x_u,
         const xt::xtensor<double, 1>& x_p,
         xt::xtensor<double, 1>& b_u) const;
 
     void dot_p(
         const xt::xtensor<double, 1>& x_u,
         const xt::xtensor<double, 1>& x_p,
         xt::xtensor<double, 1>& b_p) const;
 
     // Solve:
     // x_u = A_uu \ ( b_u - A_up * x_p ) = A_uu \ b_u
     void solve(const xt::xtensor<double, 2>& b, xt::xtensor<double, 2>& x); // modified with "x_u"
 
     void solve(const xt::xtensor<double, 1>& b, xt::xtensor<double, 1>& x); // modified with "x_u"
 
     void solve_u(
         const xt::xtensor<double, 1>& b_u,
         const xt::xtensor<double, 1>& x_p,
         xt::xtensor<double, 1>& x_u);
 
     // Get right-hand-size for corresponding to the prescribed DOFs:
     // b_p = A_pu * x_u + A_pp * x_p = A_pp * x_p
     void reaction(
         const xt::xtensor<double, 2>& x, xt::xtensor<double, 2>& b) const; // modified with "b_p"
 
     void reaction(
         const xt::xtensor<double, 1>& x, xt::xtensor<double, 1>& b) const; // modified with "b_p"
 
     void reaction_p(
         const xt::xtensor<double, 1>& x_u,
         const xt::xtensor<double, 1>& x_p,
         xt::xtensor<double, 1>& b_p) const;
 
     // Return matrix as diagonal matrix (column)
     xt::xtensor<double, 1> AsDiagonal() const;
 
     // Auto-allocation of the functions above
     xt::xtensor<double, 2> Dot(const xt::xtensor<double, 2>& x) const;
     xt::xtensor<double, 1> Dot(const xt::xtensor<double, 1>& x) const;
 
     xt::xtensor<double, 1> Dot_u(
         const xt::xtensor<double, 1>& x_u, const xt::xtensor<double, 1>& x_p) const;
 
     xt::xtensor<double, 1> Dot_p(
         const xt::xtensor<double, 1>& x_u, const xt::xtensor<double, 1>& x_p) const;
 
     xt::xtensor<double, 2> Solve(const xt::xtensor<double, 2>& b, const xt::xtensor<double, 2>& x);
     xt::xtensor<double, 1> Solve(const xt::xtensor<double, 1>& b, const xt::xtensor<double, 1>& x);
 
     xt::xtensor<double, 1> Solve_u(
         const xt::xtensor<double, 1>& b_u, const xt::xtensor<double, 1>& x_p);
 
     xt::xtensor<double, 2> Reaction(
         const xt::xtensor<double, 2>& x, const xt::xtensor<double, 2>& b) const;
 
     xt::xtensor<double, 1> Reaction(
         const xt::xtensor<double, 1>& x, const xt::xtensor<double, 1>& b) const;
 
     xt::xtensor<double, 1> Reaction_p(
         const xt::xtensor<double, 1>& x_u, const xt::xtensor<double, 1>& x_p) const;
 
 private:
     // The diagonal matrix, and its inverse (re-used to solve different RHS)
     xt::xtensor<double, 1> m_Auu;
     xt::xtensor<double, 1> m_App;
     xt::xtensor<double, 1> m_inv_uu;
 
     // Signal changes to data compare to the last inverse
-    bool m_factor = false;
+    bool m_factor = true;
 
     // Bookkeeping
     xt::xtensor<size_t, 2> m_conn; // connectivity                      [nelem, nne ]
     xt::xtensor<size_t, 2> m_dofs; // DOF-numbers per node              [nnode, ndim]
     xt::xtensor<size_t, 2> m_part; // DOF-numbers per node, renumbered  [nnode, ndim]
     xt::xtensor<size_t, 1> m_iiu;  // DOF-numbers that are unknown      [nnu]
     xt::xtensor<size_t, 1> m_iip;  // DOF-numbers that are prescribed   [nnp]
 
     // Dimensions
     size_t m_nelem; // number of elements
     size_t m_nne;   // number of nodes per element
     size_t m_nnode; // number of nodes
     size_t m_ndim;  // number of dimensions
     size_t m_ndof;  // number of DOFs
     size_t m_nnu;   // number of unknown DOFs
     size_t m_nnp;   // number of prescribed DOFs
 
     // Compute inverse (automatically evaluated by "solve")
     void factorize();
 };
 
 } // namespace GooseFEM
 
 #include "MatrixDiagonalPartitioned.hpp"
 
 #endif
diff --git a/include/GooseFEM/MatrixPartitioned.h b/include/GooseFEM/MatrixPartitioned.h
index ee7747f..4d4ca9b 100644
--- a/include/GooseFEM/MatrixPartitioned.h
+++ b/include/GooseFEM/MatrixPartitioned.h
@@ -1,185 +1,185 @@
 /*
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 */
 
 #ifndef GOOSEFEM_MATRIXPARTITIONED_H
 #define GOOSEFEM_MATRIXPARTITIONED_H
 
 #include "config.h"
 
 #include <Eigen/Eigen>
 #include <Eigen/Sparse>
 #include <Eigen/SparseCholesky>
 
 namespace GooseFEM {
 
 // forward declaration
 template <class> class MatrixPartitionedSolver;
 
 class MatrixPartitioned {
 public:
     // Constructors
     MatrixPartitioned() = default;
 
     MatrixPartitioned(
         const xt::xtensor<size_t, 2>& conn,
         const xt::xtensor<size_t, 2>& dofs,
         const xt::xtensor<size_t, 1>& iip);
 
     // Dimensions
     size_t nelem() const; // number of elements
     size_t nne() const;   // number of nodes per element
     size_t nnode() const; // number of nodes
     size_t ndim() const;  // number of dimensions
     size_t ndof() const;  // number of DOFs
     size_t nnu() const;   // number of unknown DOFs
     size_t nnp() const;   // number of prescribed DOFs
 
     // DOF lists
     xt::xtensor<size_t, 2> dofs() const; // DOFs
     xt::xtensor<size_t, 1> iiu() const;  // unknown DOFs
     xt::xtensor<size_t, 1> iip() const;  // prescribed DOFs
 
     // Assemble from matrices stored per element [nelem, nne*ndim, nne*ndim]
     void assemble(const xt::xtensor<double, 3>& elemmat);
 
     // Overwrite with a dense (sub-) matrix
     void set(
         const xt::xtensor<size_t, 1>& rows,
         const xt::xtensor<size_t, 1>& cols,
         const xt::xtensor<double, 2>& matrix);
 
     // Add a dense (sub-) matrix to the current matrix
     void add(
         const xt::xtensor<size_t, 1>& rows,
         const xt::xtensor<size_t, 1>& cols,
         const xt::xtensor<double, 2>& matrix);
 
     // Return as dense matrix
     void todense(xt::xtensor<double, 2>& ret) const;
 
     // Dot-product:
     // b_i = A_ij * x_j
     void dot(const xt::xtensor<double, 2>& x, xt::xtensor<double, 2>& b) const;
     void dot(const xt::xtensor<double, 1>& x, xt::xtensor<double, 1>& b) const;
 
     // Get right-hand-size for corresponding to the prescribed DOFs:
     // b_p = A_pu * x_u + A_pp * x_p = A_pp * x_p
     void reaction(
         const xt::xtensor<double, 2>& x,
         xt::xtensor<double, 2>& b) const; // modified with "b_p"
 
     void reaction(
         const xt::xtensor<double, 1>& x,
         xt::xtensor<double, 1>& b) const; // modified with "b_p"
 
     void reaction_p(
         const xt::xtensor<double, 1>& x_u,
         const xt::xtensor<double, 1>& x_p,
         xt::xtensor<double, 1>& b_p) const;
 
     // Auto-allocation of the functions above
     xt::xtensor<double, 2> Todense() const;
     xt::xtensor<double, 2> Dot(const xt::xtensor<double, 2>& x) const;
     xt::xtensor<double, 1> Dot(const xt::xtensor<double, 1>& x) const;
     xt::xtensor<double, 2> Reaction(
         const xt::xtensor<double, 2>& x, const xt::xtensor<double, 2>& b) const;
     xt::xtensor<double, 1> Reaction(
         const xt::xtensor<double, 1>& x, const xt::xtensor<double, 1>& b) const;
     xt::xtensor<double, 1> Reaction_p(
         const xt::xtensor<double, 1>& x_u, const xt::xtensor<double, 1>& x_p) const;
 
 private:
     // The matrix
     Eigen::SparseMatrix<double> m_Auu;
     Eigen::SparseMatrix<double> m_Aup;
     Eigen::SparseMatrix<double> m_Apu;
     Eigen::SparseMatrix<double> m_App;
 
     // Matrix entries
     std::vector<Eigen::Triplet<double>> m_Tuu;
     std::vector<Eigen::Triplet<double>> m_Tup;
     std::vector<Eigen::Triplet<double>> m_Tpu;
     std::vector<Eigen::Triplet<double>> m_Tpp;
 
     // Signal changes to data compare to the last inverse
     bool m_changed = true;
 
     // Bookkeeping
     xt::xtensor<size_t, 2> m_conn; // connectivity                      [nelem, nne ]
     xt::xtensor<size_t, 2> m_dofs; // DOF-numbers per node              [nnode, ndim]
     xt::xtensor<size_t, 2> m_part; // DOF-numbers per node, renumbered  [nnode, ndim]
     xt::xtensor<size_t, 1> m_iiu;  // unknown    DOFs                   [nnu]
     xt::xtensor<size_t, 1> m_iip;  // prescribed DOFs                   [nnp]
 
     // Dimensions
     size_t m_nelem; // number of elements
     size_t m_nne;   // number of nodes per element
     size_t m_nnode; // number of nodes
     size_t m_ndim;  // number of dimensions
     size_t m_ndof;  // number of DOFs
     size_t m_nnu;   // number of unknown DOFs
     size_t m_nnp;   // number of prescribed DOFs
 
     // grant access to solver class
     template <class> friend class MatrixPartitionedSolver;
 
     // Convert arrays (Eigen version of VectorPartitioned, which contains public functions)
     Eigen::VectorXd AsDofs_u(const xt::xtensor<double, 1>& dofval) const;
     Eigen::VectorXd AsDofs_u(const xt::xtensor<double, 2>& nodevec) const;
     Eigen::VectorXd AsDofs_p(const xt::xtensor<double, 1>& dofval) const;
     Eigen::VectorXd AsDofs_p(const xt::xtensor<double, 2>& nodevec) const;
 };
 
 template <class Solver = Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>>>
 class MatrixPartitionedSolver {
 public:
     // Constructors
     MatrixPartitionedSolver() = default;
 
     // Solve:
     // x_u = A_uu \ ( b_u - A_up * x_p )
     void solve(
         MatrixPartitioned& matrix,
         const xt::xtensor<double, 2>& b,
         xt::xtensor<double, 2>& x); // modified with "x_u"
 
     void solve(
         MatrixPartitioned& matrix,
         const xt::xtensor<double, 1>& b,
         xt::xtensor<double, 1>& x); // modified with "x_u"
 
     void solve_u(
         MatrixPartitioned& matrix,
         const xt::xtensor<double, 1>& b_u,
         const xt::xtensor<double, 1>& x_p,
         xt::xtensor<double, 1>& x_u);
 
     // Auto-allocation of the functions above
     xt::xtensor<double, 2> Solve(
         MatrixPartitioned& matrix,
         const xt::xtensor<double, 2>& b,
         const xt::xtensor<double, 2>& x);
     xt::xtensor<double, 1> Solve(
         MatrixPartitioned& matrix,
         const xt::xtensor<double, 1>& b,
         const xt::xtensor<double, 1>& x);
 
     xt::xtensor<double, 1> Solve_u(
         MatrixPartitioned& matrix,
         const xt::xtensor<double, 1>& b_u,
         const xt::xtensor<double, 1>& x_p);
 
 private:
     Solver m_solver; // solver
-    bool m_factor = false; // signal to force factorization
+    bool m_factor = true; // signal to force factorization
     void factorize(MatrixPartitioned& matrix); // compute inverse (evaluated by "solve")
 };
 
 } // namespace GooseFEM
 
 #include "MatrixPartitioned.hpp"
 
 #endif
diff --git a/include/GooseFEM/MatrixPartitionedTyings.h b/include/GooseFEM/MatrixPartitionedTyings.h
index 070cb41..edafddf 100644
--- a/include/GooseFEM/MatrixPartitionedTyings.h
+++ b/include/GooseFEM/MatrixPartitionedTyings.h
@@ -1,177 +1,177 @@
 /*
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 */
 
 #ifndef GOOSEFEM_MATRIXPARTITIONEDTYINGS_H
 #define GOOSEFEM_MATRIXPARTITIONEDTYINGS_H
 
 #include "config.h"
 
 #include <Eigen/Eigen>
 #include <Eigen/Sparse>
 #include <Eigen/SparseCholesky>
 
 namespace GooseFEM {
 
 // forward declaration
 template <class> class MatrixPartitionedTyingsSolver;
 
 class MatrixPartitionedTyings {
 public:
     // Constructors
     MatrixPartitionedTyings() = default;
 
     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);
 
     // Dimensions
     size_t nelem() const; // number of elements
     size_t nne() const;   // number of nodes per element
     size_t nnode() const; // number of nodes
     size_t ndim() const;  // number of dimensions
     size_t ndof() const;  // number of DOFs
     size_t nnu() const;   // number of independent, unknown DOFs
     size_t nnp() const;   // number of independent, prescribed DOFs
     size_t nni() const;   // number of independent DOFs
     size_t nnd() const;   // number of dependent DOFs
 
     // DOF lists
     xt::xtensor<size_t, 2> dofs() const; // DOFs
     xt::xtensor<size_t, 1> iiu() const;  // independent, unknown DOFs
     xt::xtensor<size_t, 1> iip() const;  // independent, prescribed DOFs
     xt::xtensor<size_t, 1> iii() const;  // independent DOFs
     xt::xtensor<size_t, 1> iid() const;  // dependent DOFs
 
     // Assemble from matrices stored per element [nelem, nne*ndim, nne*ndim]
     void assemble(const xt::xtensor<double, 3>& elemmat);
 
 private:
     // The matrix
     Eigen::SparseMatrix<double> m_Auu;
     Eigen::SparseMatrix<double> m_Aup;
     Eigen::SparseMatrix<double> m_Apu;
     Eigen::SparseMatrix<double> m_App;
     Eigen::SparseMatrix<double> m_Aud;
     Eigen::SparseMatrix<double> m_Apd;
     Eigen::SparseMatrix<double> m_Adu;
     Eigen::SparseMatrix<double> m_Adp;
     Eigen::SparseMatrix<double> m_Add;
 
     // The matrix for which the tyings have been applied
     Eigen::SparseMatrix<double> m_ACuu;
     Eigen::SparseMatrix<double> m_ACup;
     Eigen::SparseMatrix<double> m_ACpu;
     Eigen::SparseMatrix<double> m_ACpp;
 
     // Matrix entries
     std::vector<Eigen::Triplet<double>> m_Tuu;
     std::vector<Eigen::Triplet<double>> m_Tup;
     std::vector<Eigen::Triplet<double>> m_Tpu;
     std::vector<Eigen::Triplet<double>> m_Tpp;
     std::vector<Eigen::Triplet<double>> m_Tud;
     std::vector<Eigen::Triplet<double>> m_Tpd;
     std::vector<Eigen::Triplet<double>> m_Tdu;
     std::vector<Eigen::Triplet<double>> m_Tdp;
     std::vector<Eigen::Triplet<double>> m_Tdd;
 
     // Signal changes to data
     bool m_changed = true;
 
     // Bookkeeping
     xt::xtensor<size_t, 2> m_conn; // connectivity          [nelem, nne ]
     xt::xtensor<size_t, 2> m_dofs; // DOF-numbers per node  [nnode, ndim]
     xt::xtensor<size_t, 1> m_iiu;  // unknown     DOFs      [nnu]
     xt::xtensor<size_t, 1> m_iip;  // prescribed  DOFs      [nnp]
     xt::xtensor<size_t, 1> m_iid;  // dependent   DOFs      [nnd]
 
     // Dimensions
     size_t m_nelem; // number of elements
     size_t m_nne;   // number of nodes per element
     size_t m_nnode; // number of nodes
     size_t m_ndim;  // number of dimensions
     size_t m_ndof;  // number of DOFs
     size_t m_nnu;   // number of independent, unknown DOFs
     size_t m_nnp;   // number of independent, prescribed DOFs
     size_t m_nni;   // number of independent DOFs
     size_t m_nnd;   // number of dependent DOFs
 
     // Tyings
     Eigen::SparseMatrix<double> m_Cdu;
     Eigen::SparseMatrix<double> m_Cdp;
     Eigen::SparseMatrix<double> m_Cud;
     Eigen::SparseMatrix<double> m_Cpd;
 
     // grant access to solver class
     template <class> friend class MatrixPartitionedTyingsSolver;
 
     // Convert arrays (Eigen version of VectorPartitioned, which contains public functions)
     Eigen::VectorXd AsDofs_u(const xt::xtensor<double, 1>& dofval) const;
     Eigen::VectorXd AsDofs_u(const xt::xtensor<double, 2>& nodevec) const;
     Eigen::VectorXd AsDofs_p(const xt::xtensor<double, 1>& dofval) const;
     Eigen::VectorXd AsDofs_p(const xt::xtensor<double, 2>& nodevec) const;
     Eigen::VectorXd AsDofs_d(const xt::xtensor<double, 1>& dofval) const;
     Eigen::VectorXd AsDofs_d(const xt::xtensor<double, 2>& nodevec) const;
 };
 
 template <class Solver = Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>>>
 class MatrixPartitionedTyingsSolver {
 public:
     // Constructors
     MatrixPartitionedTyingsSolver() = default;
 
     // Solve:
     // A' = A_ii + K_id * C_di + C_di^T * K_di + C_di^T * K_dd * C_di
     // b' = b_i + C_di^T * b_d
     // x_u = A'_uu \ ( b'_u - A'_up * x_p )
     // x_i = [x_u, x_p]
     // x_d = C_di * x_i
     void solve(
         MatrixPartitionedTyings& matrix,
         const xt::xtensor<double, 2>& b,
         xt::xtensor<double, 2>& x); // updates x_u and x_d
 
     void solve(
         MatrixPartitionedTyings& matrix,
         const xt::xtensor<double, 1>& b,
         xt::xtensor<double, 1>& x); // updates x_u and x_d
 
     void 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);
 
     // Auto-allocation of the functions above
     xt::xtensor<double, 2> Solve(
         MatrixPartitionedTyings& matrix,
         const xt::xtensor<double, 2>& b,
         const xt::xtensor<double, 2>& x);
 
     xt::xtensor<double, 1> Solve(
         MatrixPartitionedTyings& matrix,
         const xt::xtensor<double, 1>& b,
         const xt::xtensor<double, 1>& x);
 
     xt::xtensor<double, 1> 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);
 
 private:
     Solver m_solver; // solver
-    bool m_factor = false; // signal to force factorization
+    bool m_factor = true; // signal to force factorization
     void factorize(MatrixPartitionedTyings& matrix); // compute inverse (evaluated by "solve")
 };
 
 } // namespace GooseFEM
 
 #include "MatrixPartitionedTyings.hpp"
 
 #endif
diff --git a/test/basic/MatrixParitioned.cpp b/test/basic/MatrixParitioned.cpp
index 1af23e7..7a5736a 100644
--- a/test/basic/MatrixParitioned.cpp
+++ b/test/basic/MatrixParitioned.cpp
@@ -1,103 +1,110 @@
 
 #include <catch2/catch.hpp>
 #include <xtensor/xrandom.hpp>
 #include <xtensor/xmath.hpp>
 #include <Eigen/Eigen>
 #include <GooseFEM/GooseFEM.h>
 
 #define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12));
 
 TEST_CASE("GooseFEM::MatrixPartitioned", "MatrixPartitioned.h")
 {
     SECTION("solve")
     {
         GooseFEM::Mesh::Quad4::Regular mesh(2, 2);
 
         size_t nne = mesh.nne();
         size_t ndim = mesh.ndim();
         size_t nelem = mesh.nelem();
         size_t nnode = mesh.nnode();
         auto dofs = mesh.dofs();
         size_t npp = xt::amax(dofs)();
         npp = (npp - npp % 2) / 2;
         xt::xtensor<size_t, 1> iip = xt::arange<size_t<(npp);
 
         xt::xtensor<double, 3> a = xt::empty<double>({nelem, nne * ndim, nne * ndim});
         xt::xtensor<double, 1> b = xt::random::rand<double>({nnode * ndim});
 
         for (size_t e = 0; e < nelem; ++e) {
             xt::xtensor<double, 2> ae = xt::random::rand<double>({nne * ndim, nne * ndim});
             ae = (ae + xt::transpose(ae)) / 2.0;
             xt::view(a, e, xt::all(), xt::all()) = ae;
         }
 
         GooseFEM::MatrixPartitioned A(mesh.conn(), dofs, iip);
         GooseFEM::MatrixPartitionedSolver<> Solver;
         A.assemble(a);
         xt::xtensor<double, 1> C = A.Dot(b);
         xt::xtensor<double, 1> B = Solver.Solve(A, C);
 
         REQUIRE(B.size() == b.size());
         REQUIRE(xt::allclose(B, b));
+
+        // check that allocating a different Solver instance still works
+        GooseFEM::MatrixPartitionedSolver<> NewSolver;
+        xt::xtensor<double, 1> NB = NewSolver.Solve(A, C);
+
+        REQUIRE(NB.size() == b.size());
+        REQUIRE(xt::allclose(NB, b));
     }
 
     SECTION("set/add/dot/solve - dofval")
     {
         xt::xtensor<double, 2> a = xt::random::rand<double>({10, 10});
         xt::xtensor<double, 1> x = xt::random::rand<double>({10});
         xt::xtensor<double, 1> b = xt::zeros<double>({10});
 
         xt::xtensor<double, 2> A = a + xt::transpose(a);
 
         for (size_t i = 0; i < A.shape(0); ++i) {
             for (size_t j = 0; j < A.shape(1); ++j) {
                 b(i) += A(i, j) * x(j);
             }
         }
 
         xt::xtensor<size_t, 2> conn = xt::zeros<size_t>({1, 5});
         xt::xtensor<size_t, 2> dofs = xt::arange<size_t>(10).reshape({5, 2});
         xt::xtensor<size_t, 1> iip = xt::arange<size_t<(5);
 
         GooseFEM::MatrixPartitioned K(conn, dofs, iip);
         GooseFEM::MatrixPartitionedSolver<> Solver;
         K.set(xt::arange<size_t>(10), xt::arange<size_t>(10), a);
         K.add(xt::arange<size_t>(10), xt::arange<size_t>(10), xt::transpose(a));
 
         REQUIRE(xt::allclose(A, K.Todense()));
         REQUIRE(xt::allclose(b, K.Dot(x)));
         REQUIRE(xt::allclose(x, Solver.Solve(K, b)));
     }
 
     SECTION("set/add/dot/solve - nodevec")
     {
         xt::xtensor<double, 2> a = xt::random::rand<double>({10, 10});
         xt::xtensor<double, 2> x = xt::random::rand<double>({5, 2});
         xt::xtensor<double, 2> b = xt::zeros<double>({5, 2});
 
         xt::xtensor<double, 2> A = a + xt::transpose(a);
 
         for (size_t m = 0; m < x.shape(0); ++m) {
             for (size_t n = 0; n < x.shape(0); ++n) {
                 for (size_t i = 0; i < x.shape(1); ++i) {
                     for (size_t j = 0; j < x.shape(1); ++j) {
                         b(m, i) += A(m * x.shape(1) + i, n * x.shape(1) + j) * x(n, j);
                     }
                 }
             }
         }
 
         xt::xtensor<size_t, 2> conn = xt::zeros<size_t>({1, 5});
         xt::xtensor<size_t, 2> dofs = xt::arange<size_t>(10).reshape({5, 2});
         xt::xtensor<size_t, 1> iip = xt::arange<size_t<(5);
 
         GooseFEM::MatrixPartitioned K(conn, dofs, iip);
         GooseFEM::MatrixPartitionedSolver<> Solver;
         K.set(xt::arange<size_t>(10), xt::arange<size_t>(10), a);
         K.add(xt::arange<size_t>(10), xt::arange<size_t>(10), xt::transpose(a));
 
         REQUIRE(xt::allclose(A, K.Todense()));
         REQUIRE(xt::allclose(b, K.Dot(x)));
         REQUIRE(xt::allclose(x, Solver.Solve(K, b)));
     }
 }