diff --git a/include/GooseFEM/MatrixPartitioned.h b/include/GooseFEM/MatrixPartitioned.h
index 6136095..17d44d5 100644
--- a/include/GooseFEM/MatrixPartitioned.h
+++ b/include/GooseFEM/MatrixPartitioned.h
@@ -1,176 +1,182 @@
 /*
 
 (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);
 
+
+    // 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> 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
     void factorize(MatrixPartitioned& matrix); // compute inverse (evaluated by "solve")
 };
 
 } // namespace GooseFEM
 
 #include "MatrixPartitioned.hpp"
 
 #endif
diff --git a/include/GooseFEM/MatrixPartitioned.hpp b/include/GooseFEM/MatrixPartitioned.hpp
index 5859bba..4ffce36 100644
--- a/include/GooseFEM/MatrixPartitioned.hpp
+++ b/include/GooseFEM/MatrixPartitioned.hpp
@@ -1,475 +1,525 @@
 /*
 
 (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_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::set(
     const xt::xtensor<size_t, 1>& rows,
     const xt::xtensor<size_t, 1>& cols,
     const xt::xtensor<double, 2>& matrix)
 {
     GOOSEFEM_ASSERT(rows.size() == matrix.shape(0));
     GOOSEFEM_ASSERT(cols.size() == matrix.shape(1));
     GOOSEFEM_ASSERT(xt::amax(cols)() < m_ndof);
     GOOSEFEM_ASSERT(xt::amax(rows)() < m_ndof);
 
     std::vector<Eigen::Triplet<double>> Tuu;
     std::vector<Eigen::Triplet<double>> Tup;
     std::vector<Eigen::Triplet<double>> Tpu;
     std::vector<Eigen::Triplet<double>> Tpp;
 
     for (size_t i = 0; i < rows.size(); ++i) {
         for (size_t j = 0; j < cols.size(); ++j) {
             size_t di = rows(i);
             size_t dj = cols(j);
             double v = matrix(i, j);
 
             if (di < m_nnu && dj < m_nnu) {
                 Tuu.push_back(Eigen::Triplet<double>(di, dj, v));
             }
             else if (di < m_nnu) {
                 Tup.push_back(Eigen::Triplet<double>(di, dj - m_nnu, v));
             }
             else if (dj < m_nnu) {
                 Tpu.push_back(Eigen::Triplet<double>(di - m_nnu, dj, v));
             }
             else {
                 Tpp.push_back(Eigen::Triplet<double>(di - m_nnu, dj - m_nnu, v));
             }
         }
     }
 
     m_Auu.setFromTriplets(Tuu.begin(), Tuu.end());
     m_Aup.setFromTriplets(Tup.begin(), Tup.end());
     m_Apu.setFromTriplets(Tpu.begin(), Tpu.end());
     m_App.setFromTriplets(Tpp.begin(), Tpp.end());
     m_changed = true;
 }
 
 inline void MatrixPartitioned::add(
     const xt::xtensor<size_t, 1>& rows,
     const xt::xtensor<size_t, 1>& cols,
     const xt::xtensor<double, 2>& matrix)
 {
     GOOSEFEM_ASSERT(rows.size() == matrix.shape(0));
     GOOSEFEM_ASSERT(cols.size() == matrix.shape(1));
     GOOSEFEM_ASSERT(xt::amax(cols)() < m_ndof);
     GOOSEFEM_ASSERT(xt::amax(rows)() < m_ndof);
 
     std::vector<Eigen::Triplet<double>> Tuu;
     std::vector<Eigen::Triplet<double>> Tup;
     std::vector<Eigen::Triplet<double>> Tpu;
     std::vector<Eigen::Triplet<double>> Tpp;
 
-    Eigen::SparseMatrix<double> Auu;
-    Eigen::SparseMatrix<double> Aup;
-    Eigen::SparseMatrix<double> Apu;
-    Eigen::SparseMatrix<double> App;
+    Eigen::SparseMatrix<double> Auu(m_nnu, m_nnu);
+    Eigen::SparseMatrix<double> Aup(m_nnu, m_nnp);
+    Eigen::SparseMatrix<double> Apu(m_nnp, m_nnu);
+    Eigen::SparseMatrix<double> App(m_nnp, m_nnp);
 
     for (size_t i = 0; i < rows.size(); ++i) {
         for (size_t j = 0; j < cols.size(); ++j) {
             size_t di = rows(i);
             size_t dj = cols(j);
             double v = matrix(i, j);
 
             if (di < m_nnu && dj < m_nnu) {
                 Tuu.push_back(Eigen::Triplet<double>(di, dj, v));
             }
             else if (di < m_nnu) {
                 Tup.push_back(Eigen::Triplet<double>(di, dj - m_nnu, v));
             }
             else if (dj < m_nnu) {
                 Tpu.push_back(Eigen::Triplet<double>(di - m_nnu, dj, v));
             }
             else {
                 Tpp.push_back(Eigen::Triplet<double>(di - m_nnu, dj - m_nnu, v));
             }
         }
     }
 
     Auu.setFromTriplets(Tuu.begin(), Tuu.end());
     Aup.setFromTriplets(Tup.begin(), Tup.end());
     Apu.setFromTriplets(Tpu.begin(), Tpu.end());
     App.setFromTriplets(Tpp.begin(), Tpp.end());
     m_Auu += Auu;
     m_Aup += Aup;
     m_Apu += Apu;
     m_App += App;
     m_changed = true;
 }
+
+inline void MatrixPartitioned::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}));
+
+    Eigen::VectorXd X_u = this->AsDofs_u(x);
+    Eigen::VectorXd X_p = this->AsDofs_p(x);
+    Eigen::VectorXd B_u = m_Auu * X_u + m_Aup * X_p;
+    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_u(m_part(m, i));
+            }
+            else{
+                b(m, i) = B_p(m_part(m, i) - m_nnu);
+            }
+        }
+    }
+}
+
+inline void MatrixPartitioned::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::VectorXd X_u = this->AsDofs_u(x);
+    Eigen::VectorXd X_p = this->AsDofs_p(x);
+
+    Eigen::Map<Eigen::VectorXd>(b.data(), m_nnu).noalias() = m_Auu * X_u + m_Aup * X_p;
+    Eigen::Map<Eigen::VectorXd>(b.data() + m_nnu, m_ndof).noalias() = m_Auu * X_u + m_Aup * X_p;
+}
+
+inline xt::xtensor<double, 2> MatrixPartitioned::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> MatrixPartitioned::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 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> 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> 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 = 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 = 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> ret = x;
     this->solve(matrix, b, ret);
     return ret;
 }
 
 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> ret = x;
     this->solve(matrix, b, ret);
     return ret;
 }
 
 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/python/MatrixPartitioned.hpp b/python/MatrixPartitioned.hpp
index f0fadf0..5a1b75c 100644
--- a/python/MatrixPartitioned.hpp
+++ b/python/MatrixPartitioned.hpp
@@ -1,124 +1,152 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
 #include <Eigen/Eigen>
 #include <GooseFEM/GooseFEM.h>
 #include <pybind11/eigen.h>
 #include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
 
 namespace py = pybind11;
 
 void init_MatrixPartitioned(py::module& m)
 {
 
     py::class_<GooseFEM::MatrixPartitioned>(m, "MatrixPartitioned")
 
         .def(
             py::init<
                 const xt::xtensor<size_t, 2>&,
                 const xt::xtensor<size_t, 2>&,
                 const xt::xtensor<size_t, 1>&>(),
             "Sparse, partitioned, matrix",
             py::arg("conn"),
             py::arg("dofs"),
             py::arg("iip"))
 
         .def("nelem", &GooseFEM::MatrixPartitioned::nelem, "Number of element")
 
         .def("nne", &GooseFEM::MatrixPartitioned::nne, "Number of nodes per element")
 
         .def("nnode", &GooseFEM::MatrixPartitioned::nnode, "Number of nodes")
 
         .def("ndim", &GooseFEM::MatrixPartitioned::ndim, "Number of dimensions")
 
         .def("ndof", &GooseFEM::MatrixPartitioned::ndof, "Number of degrees-of-freedom")
 
         .def("nnu", &GooseFEM::MatrixPartitioned::nnu, "Number of unknown degrees-of-freedom")
 
         .def("nnp", &GooseFEM::MatrixPartitioned::nnp, "Number of prescribed degrees-of-freedom")
 
         .def(
             "assemble",
             &GooseFEM::MatrixPartitioned::assemble,
             "Assemble matrix from 'elemmat",
             py::arg("elemmat"))
 
+        .def(
+            "set",
+            &GooseFEM::MatrixPartitioned::set,
+            "Overwrite with a dense (sub-) matrix",
+            py::arg("rows"),
+            py::arg("cols"),
+            py::arg("matrix"))
+
+        .def(
+            "add",
+            &GooseFEM::MatrixPartitioned::add,
+            "Add a dense (sub-) matrix to the current matrix",
+            py::arg("rows"),
+            py::arg("cols"),
+            py::arg("matrix"))
+
         .def("dofs", &GooseFEM::MatrixPartitioned::dofs, "Return degrees-of-freedom")
 
         .def("iiu", &GooseFEM::MatrixPartitioned::iiu, "Return unknown degrees-of-freedom")
 
         .def("iip", &GooseFEM::MatrixPartitioned::iip, "Return prescribed degrees-of-freedom")
 
         .def(
             "Reaction",
             py::overload_cast<const xt::xtensor<double, 1>&, const xt::xtensor<double, 1>&>(
                 &GooseFEM::MatrixPartitioned::Reaction, py::const_),
             "Reaction",
             py::arg("x"),
             py::arg("b"))
 
         .def(
             "Reaction",
             py::overload_cast<const xt::xtensor<double, 2>&, const xt::xtensor<double, 2>&>(
                 &GooseFEM::MatrixPartitioned::Reaction, py::const_),
             "Reaction",
             py::arg("x"),
             py::arg("b"))
 
         .def(
             "Reaction_p",
             py::overload_cast<const xt::xtensor<double, 1>&, const xt::xtensor<double, 1>&>(
                 &GooseFEM::MatrixPartitioned::Reaction_p, py::const_),
             "Reaction_p",
             py::arg("x_u"),
             py::arg("x_p"))
 
+        .def(
+            "Dot",
+            py::overload_cast<const xt::xtensor<double, 1>&>(&GooseFEM::MatrixPartitioned::Dot, py::const_),
+            "Dot",
+            py::arg("x"))
+
+        .def(
+            "Dot",
+            py::overload_cast<const xt::xtensor<double, 2>&>(&GooseFEM::MatrixPartitioned::Dot, py::const_),
+            "Dot",
+            py::arg("x"))
+
         .def("__repr__", [](const GooseFEM::MatrixPartitioned&) {
             return "<GooseFEM.MatrixPartitioned>";
         });
 
     py::class_<GooseFEM::MatrixPartitionedSolver<>>(m, "MatrixPartitionedSolver")
 
         .def(py::init<>(), "Sparse, partitioned, matrix solver")
 
         .def(
             "Solve",
             py::overload_cast<
                 GooseFEM::MatrixPartitioned&,
                 const xt::xtensor<double, 1>&,
                 const xt::xtensor<double, 1>&>(&GooseFEM::MatrixPartitionedSolver<>::Solve),
             "Solve",
             py::arg("matrix"),
             py::arg("b"),
             py::arg("x"))
 
         .def(
             "Solve",
             py::overload_cast<
                 GooseFEM::MatrixPartitioned&,
                 const xt::xtensor<double, 2>&,
                 const xt::xtensor<double, 2>&>(&GooseFEM::MatrixPartitionedSolver<>::Solve),
             "Solve",
             py::arg("matrix"),
             py::arg("b"),
             py::arg("x"))
 
         .def(
             "Solve_u",
             py::overload_cast<
                 GooseFEM::MatrixPartitioned&,
                 const xt::xtensor<double, 1>&,
                 const xt::xtensor<double, 1>&>(&GooseFEM::MatrixPartitionedSolver<>::Solve_u),
             "Solve_u",
             py::arg("matrix"),
             py::arg("b_u"),
             py::arg("x_p"))
 
         .def("__repr__", [](const GooseFEM::MatrixPartitionedSolver<>&) {
             return "<GooseFEM.MatrixPartitionedSolver>";
         });
 }