diff --git a/include/GooseFEM/Matrix.h b/include/GooseFEM/Matrix.h
index 12831ee..9f6c035 100644
--- a/include/GooseFEM/Matrix.h
+++ b/include/GooseFEM/Matrix.h
@@ -1,239 +1,242 @@
 /**
 Sparse matrix.
 
 \file Matrix.h
 \copyright Copyright 2017. Tom de Geus. All rights reserved.
 \license This project is released under the GNU Public License (GPLv3).
 */
 
 #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;
 
 /**
 Sparse matrix.
 
 See Vector() for bookkeeping definitions.
 */
 class Matrix {
 public:
 
     Matrix() = default;
 
     /**
     Constructor.
 
     \param conn connectivity [#nelem, #nne].
     \param dofs DOFs per node [#nnode, #ndim].
     */
     Matrix(const xt::xtensor<size_t, 2>& conn, const xt::xtensor<size_t, 2>& dofs);
 
     /**
     \return Number of elements.
     */
     size_t nelem() const;
 
     /**
     \return Number of nodes per element.
     */
     size_t nne() const;
 
     /**
     \return Number of nodes.
     */
     size_t nnode() const;
 
     /**
     \return Number of dimensions.
     */
     size_t ndim() const;
 
     /**
     \return Number of DOFs.
     */
     size_t ndof() const;
 
     /**
     \return DOFs per node [#nnode, #ndim]
     */
     xt::xtensor<size_t, 2> dofs() const;
 
     /**
     Assemble from matrices stored per element.
 
     \param elemmat [#nelem, #nne * #ndim, #nne * #ndim].
     */
     virtual void assemble(const xt::xtensor<double, 3>& elemmat);
 
     /**
-    Overwrite matrix with dense matrix.
+    Overwrite matrix with (sparse) matrix.
 
     \param rows Row numbers in the matrix [n].
     \param cols Column numbers in the matrix [n].
     \param matrix Data entries on (rows, cols) [n].
     */
     virtual 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 matrix to the matrix.
+    Add a (sparse) matrix to the matrix.
 
     \param rows Row numbers in the matrix [n].
     \param cols Column numbers in the matrix [n].
     \param matrix Data entries on (rows, cols) [n].
     */
     virtual void add(
         const xt::xtensor<size_t, 1>& rows,
         const xt::xtensor<size_t, 1>& cols,
         const xt::xtensor<double, 2>& matrix);
 
     /**
     \return Matrix as dense matrix [#ndof, #ndof].
     */
     xt::xtensor<double, 2> Todense() const;
 
     /**
     Covert matrix to dense matrix.
 
     \param ret overwritten [#ndof, #ndof].
     */
     virtual void todense(xt::xtensor<double, 2>& ret) const;
 
     /**
     Dot-product \f$ b_i = A_{ij} x_j \f$.
 
     \param x nodevec [#nelem, #ndim].
     \return b nodevec overwritten [#nelem, #ndim].
     */
     xt::xtensor<double, 2> Dot(const xt::xtensor<double, 2>& x) const;
 
     /**
     Same as Dot(const xt::xtensor<double, 2>&, xt::xtensor<double, 2>& b)
     but writing to preallocated data.
 
     \param x nodevec [#nelem, #ndim].
     \param b nodevec overwritten [#nelem, #ndim].
     */
     virtual void dot(const xt::xtensor<double, 2>& x, xt::xtensor<double, 2>& b) const;
 
     /**
     Dot-product \f$ b_i = A_{ij} x_j \f$.
 
     \param x dofval [#ndof].
     \return b dofval overwritten [#ndof].
     */
     xt::xtensor<double, 1> Dot(const xt::xtensor<double, 1>& x) const;
 
     /**
     Same as Dot(const xt::xtensor<double, 1>&, xt::xtensor<double, 1>& b)
     but writing to preallocated data.
 
     \param x dofval [#ndof].
     \param b dofval overwritten [#ndof].
     */
     virtual void dot(const xt::xtensor<double, 1>& x, xt::xtensor<double, 1>& b) const;
 
 protected:
+
     xt::xtensor<size_t, 2> m_conn; ///< Connectivity [#nelem, #nne].
     xt::xtensor<size_t, 2> m_dofs; ///< DOF-numbers per node [#nnode, #ndim].
     size_t m_nelem; ///< See nelem().
     size_t m_nne; ///< See nne().
     size_t m_nnode; ///< See nnode().
     size_t m_ndim; ///< See ndim().
     size_t m_ndof; ///< See ndof().
     bool m_changed = true; ///< Signal changes to data.
 
 private:
+
     Eigen::SparseMatrix<double> m_A; ///< The matrix.
     std::vector<Eigen::Triplet<double>> m_T; ///< Matrix entries.
     template <class> friend class MatrixSolver; ///< Grant access to solver class
 
     /**
     Convert "nodevec" to "dofval" (overwrite entries that occur more than once).
 
     \param nodevec input [#nnode, #ndim]
     \return dofval output [#ndof]
     */
     Eigen::VectorXd AsDofs(const xt::xtensor<double, 2>& nodevec) const;
 
     /**
     Convert "dofval" to "nodevec" (overwrite entries that occur more than once).
 
     \param dofval input [#ndof]
     \param nodevec output [#nnode, #ndim]
     */
     void asNode(const Eigen::VectorXd& dofval, xt::xtensor<double, 2>& nodevec) const;
 };
 
 /**
 Solver for Matrix().
 The idea is that this solver class can be used to solve for multiple right-hand-sides
 using one factorisation.
 */
 template <class Solver = Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>>>
 class MatrixSolver {
 public:
 
     MatrixSolver() = default;
 
     /**
     Solve \f$ x = A^{-1} b \f$.
 
     \param A sparse matrix, see Matrix().
     \param b nodevec [nelem, ndim].
     \return x nodevec [nelem, ndim].
     */
     xt::xtensor<double, 2> Solve(Matrix& A, const xt::xtensor<double, 2>& b);
 
     /**
     Same as Solve(Matrix&, const xt::xtensor<double, 2>&)
     but writing to preallocated data.
 
     \param A sparse matrix, see Matrix().
     \param b nodevec [nelem, ndim].
     \param x nodevec overwritten [nelem, ndim].
     */
     void solve(Matrix& A, const xt::xtensor<double, 2>& b, xt::xtensor<double, 2>& x);
 
     /**
-    Solve \f$ x = A^{-1} b \f$.
+    Same as Solve(Matrix&, const xt::xtensor<double, 2>&)
+    but for "dofval" input and output.
 
     \param A sparse matrix, see Matrix().
     \param b dofval [ndof].
     \return x dofval [ndof].
     */
     xt::xtensor<double, 1> Solve(Matrix& A, const xt::xtensor<double, 1>& b);
 
     /**
     Same as Solve(Matrix&, const xt::xtensor<double, 1>&)
     but writing to preallocated data.
 
     \param A sparse matrix, see Matrix().
     \param b dofval [ndof].
     \param x dofval overwritten [ndof].
     */
     void solve(Matrix& A, const xt::xtensor<double, 1>& b, xt::xtensor<double, 1>& x);
 
 private:
     Solver m_solver; ///< Solver.
     bool m_factor = true; ///< Signal to force factorization.
     void factorize(Matrix& A); ///< Compute inverse (evaluated by "solve").
 };
 
 } // namespace GooseFEM
 
 #include "Matrix.hpp"
 
 #endif
diff --git a/include/GooseFEM/MatrixDiagonal.h b/include/GooseFEM/MatrixDiagonal.h
index 29329ce..1d6114b 100644
--- a/include/GooseFEM/MatrixDiagonal.h
+++ b/include/GooseFEM/MatrixDiagonal.h
@@ -1,85 +1,176 @@
 /**
 Diagonal matrix.
 
 \file MatrixDiagonal.h
 \copyright Copyright 2017. Tom de Geus. All rights reserved.
 \license This project is released under the GNU Public License (GPLv3).
 */
 
 #ifndef GOOSEFEM_MATRIXDIAGONAL_H
 #define GOOSEFEM_MATRIXDIAGONAL_H
 
 #include "config.h"
 
 namespace GooseFEM {
 
+/**
+Diagonal matrix.
+
+See Vector() for bookkeeping definitions.
+*/
 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
+    /**
+    Constructor.
 
-    // DOF lists
-    xt::xtensor<size_t, 2> dofs() const; // DOFs
+    \param conn connectivity [#nelem, #nne].
+    \param dofs DOFs per node [#nnode, #ndim].
+    */
+    MatrixDiagonal(const xt::xtensor<size_t, 2>& conn, const xt::xtensor<size_t, 2>& dofs);
 
-    // Set matrix components
+    /**
+    \return Number of elements.
+    */
+    size_t nelem() const;
+
+    /**
+    \return Number of nodes per element.
+    */
+    size_t nne() const;
+
+    /**
+    \return Number of nodes.
+    */
+    size_t nnode() const;
+
+    /**
+    \return Number of dimensions.
+    */
+    size_t ndim() const;
+
+    /**
+    \return Number of DOFs.
+    */
+    size_t ndof() const;
+
+    /**
+    \return DOFs per node [#nnode, #ndim]
+    */
+    xt::xtensor<size_t, 2> dofs() const;
+
+    /**
+    Assemble from matrices stored per element.
+    \warning Ignores any off-diagonal terms.
+
+    \param elemmat [#nelem, #nne * #ndim, #nne * #ndim].
+    */
+    virtual void assemble(const xt::xtensor<double, 3>& elemmat);
+
+    /**
+    Set all (diagonal) matrix components.
+
+    \param A The matrix [#ndof].
+    */
     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;
+    /**
+    Return matrix as diagonal matrix.
 
-    // 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);
+    \param [#ndof].
+    */
+    virtual xt::xtensor<double, 1> Todiagonal() const;
 
-    // Return matrix as diagonal matrix (column)
-    xt::xtensor<double, 1> Todiagonal() const;
+    /**
+    Dot-product \f$ b_i = A_{ij} x_j \f$.
 
-    // Auto-allocation of the functions above
+    \param x nodevec [#nelem, #ndim].
+    \return b nodevec overwritten [#nelem, #ndim].
+    */
     xt::xtensor<double, 2> Dot(const xt::xtensor<double, 2>& x) const;
+
+    /**
+    Same as Dot(const xt::xtensor<double, 2>&, xt::xtensor<double, 2>& b)
+    but writing to preallocated data.
+
+    \param x nodevec [#nelem, #ndim].
+    \param b nodevec overwritten [#nelem, #ndim].
+    */
+    virtual void dot(const xt::xtensor<double, 2>& x, xt::xtensor<double, 2>& b) const;
+
+    /**
+    Dot-product \f$ b_i = A_{ij} x_j \f$.
+
+    \param x dofval [#ndof].
+    \return b dofval overwritten [#ndof].
+    */
     xt::xtensor<double, 1> Dot(const xt::xtensor<double, 1>& x) const;
+
+    /**
+    Same as Dot(const xt::xtensor<double, 1>&, xt::xtensor<double, 1>& b)
+    but writing to preallocated data.
+
+    \param x dofval [#ndof].
+    \param b dofval overwritten [#ndof].
+    */
+    virtual void dot(const xt::xtensor<double, 1>& x, xt::xtensor<double, 1>& b) const;
+
+    /**
+    Solve \f$ x = A^{-1} b \f$.
+
+    \param b nodevec [nelem, ndim].
+    \return x nodevec [nelem, ndim].
+    */
     xt::xtensor<double, 2> Solve(const xt::xtensor<double, 2>& b);
+
+    /**
+    Same as Solve(const xt::xtensor<double, 2>&)
+    but writing to preallocated data.
+
+    \param b nodevec [nelem, ndim].
+    \param x nodevec overwritten [nelem, ndim].
+    */
+    virtual void solve(const xt::xtensor<double, 2>& b, xt::xtensor<double, 2>& x);
+
+    /**
+    Same as Solve(const xt::xtensor<double, 2>&)
+    but for "dofval" input and output.
+
+    \param b dofval [ndof].
+    \return x dofval [ndof].
+    */
     xt::xtensor<double, 1> Solve(const xt::xtensor<double, 1>& b);
 
+    /**
+    Same as Solve(const xt::xtensor<double, 1>&)
+    but writing to preallocated data.
+
+    \param b dofval [ndof].
+    \param x dofval overwritten [ndof].
+    */
+    virtual void solve(const xt::xtensor<double, 1>& b, xt::xtensor<double, 1>& x);
+
+protected:
+
+    xt::xtensor<size_t, 2> m_conn; ///< Connectivity [#nelem, #nne].
+    xt::xtensor<size_t, 2> m_dofs; ///< DOF-numbers per node [#nnode, #ndim].
+    size_t m_nelem; ///< See nelem().
+    size_t m_nne; ///< See nne().
+    size_t m_nnode; ///< See nnode().
+    size_t m_ndim; ///< See ndim().
+    size_t m_ndof; ///< See ndof().
+    bool m_changed = true; ///< Signal changes to data.
+
 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 = 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();
+    xt::xtensor<double, 1> m_A; ///< The matrix.
+    xt::xtensor<double, 1> m_inv; /// Inverse of the matrix.
+    void factorize(); ///< Compute inverse (automatically evaluated by "solve").
 };
 
 } // namespace GooseFEM
 
 #include "MatrixDiagonal.hpp"
 
 #endif
diff --git a/include/GooseFEM/MatrixDiagonal.hpp b/include/GooseFEM/MatrixDiagonal.hpp
index 2cf9eb9..4f28b31 100644
--- a/include/GooseFEM/MatrixDiagonal.hpp
+++ b/include/GooseFEM/MatrixDiagonal.hpp
@@ -1,180 +1,180 @@
 /**
 Implementation of MatrixDiagonal.h
 
 \file MatrixDiagonal.hpp
 \copyright Copyright 2017. Tom de Geus. All rights reserved.
 \license This project is released under the GNU Public License (GPLv3).
 */
 
 #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(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) {
+    if (!m_changed) {
         return;
     }
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_ndof; ++d) {
         m_inv(d) = 1.0 / m_A(d);
     }
 
-    m_factor = false;
+    m_changed = 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;
+    m_changed = 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;
+    m_changed = 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::Todiagonal() 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.h b/include/GooseFEM/MatrixDiagonalPartitioned.h
index 257a8a1..9d64799 100644
--- a/include/GooseFEM/MatrixDiagonalPartitioned.h
+++ b/include/GooseFEM/MatrixDiagonalPartitioned.h
@@ -1,146 +1,226 @@
 /**
 Diagonal matrix that is partitioned in:
 -   unknown DOFs
 -   prescribed DOFs
 
 \file MatrixDiagonalPartitioned.h
 \copyright Copyright 2017. Tom de Geus. All rights reserved.
 \license This project is released under the GNU Public License (GPLv3).
 */
 
 #ifndef GOOSEFEM_MATRIXDIAGONALPARTITIONED_H
 #define GOOSEFEM_MATRIXDIAGONALPARTITIONED_H
 
 #include "config.h"
+#include "MatrixDiagonal.h"
 
 namespace GooseFEM {
 
-class MatrixDiagonalPartitioned {
+/**
+Diagonal and partitioned matrix.
+
+See Vector() for bookkeeping definitions.
+*/
+class MatrixDiagonalPartitioned : public MatrixDiagonal {
 public:
-    // Constructors
+
     MatrixDiagonalPartitioned() = default;
 
+    /**
+    Constructor.
+
+    \param conn connectivity [#nelem, #nne].
+    \param dofs DOFs per node [#nnode, #ndim].
+    \param iip prescribed DOFs [#nnp].
+    */
     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;
+    /**
+    Number of unknown DOFs.
+    */
+    size_t nnu() const;
+
+    /**
+    Number of prescribed DOFs.
+    */
+    size_t nnp() const;
+
+    /**
+    Unknown DOFs [#nnu].
+    */
+    xt::xtensor<size_t, 1> iiu() const;
+
+    /**
+    Prescribed DOFs [#nnp].
+    */
+    xt::xtensor<size_t, 1> iip() const;
+
+    void assemble(const xt::xtensor<double, 3>& elemmat) override;
 
+    xt::xtensor<double, 1> Todiagonal() const override;
+
+    void dot(const xt::xtensor<double, 2>& x, xt::xtensor<double, 2>& b) const override;
+    void dot(const xt::xtensor<double, 1>& x, xt::xtensor<double, 1>& b) const override;
+
+    /**
+    \param x_u dofval [#nnu].
+    \param x_p dofval [#nnp].
+    \return b_u dofval [#nnu].
+    */
+    xt::xtensor<double, 1> Dot_u(
+        const xt::xtensor<double, 1>& x_u,
+        const xt::xtensor<double, 1>& x_p) const;
+
+    /**
+    \param x_u dofval [#nnu].
+    \param x_p dofval [#nnp].
+    \param b_u (overwritten) dofval [#nnu].
+    */
     void dot_u(
         const xt::xtensor<double, 1>& x_u,
         const xt::xtensor<double, 1>& x_p,
-        xt::xtensor<double, 1>& b_u) const;
+              xt::xtensor<double, 1>& b_u) const;
 
+    /**
+    \param x_u dofval [#nnu].
+    \param x_p dofval [#nnp].
+    \return b_p dofval [#nnp].
+    */
+    xt::xtensor<double, 1> Dot_p(
+        const xt::xtensor<double, 1>& x_u,
+        const xt::xtensor<double, 1>& x_p) const;
+
+    /**
+    \param x_u dofval [#nnu].
+    \param x_p dofval [#nnp].
+    \param b_p (overwritten) dofval [#nnp].
+    */
     void dot_p(
         const xt::xtensor<double, 1>& x_u,
         const xt::xtensor<double, 1>& x_p,
-        xt::xtensor<double, 1>& b_p) const;
+              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, 2>& b, xt::xtensor<double, 2>& x) override; // modified with "x_u"
+    void solve(const xt::xtensor<double, 1>& b, xt::xtensor<double, 1>& x) override; // modified with "x_u"
 
-    void solve(const xt::xtensor<double, 1>& b, xt::xtensor<double, 1>& x); // modified with "x_u"
+    /**
+    \param b_u dofval [#nnu].
+    \param x_p dofval [#nnp].
+    \return x_u dofval [#nnu].
+    */
+    xt::xtensor<double, 1> Solve_u(
+        const xt::xtensor<double, 1>& b_u,
+        const xt::xtensor<double, 1>& x_p);
 
+    /**
+    \param b_u dofval [#nnu].
+    \param x_p dofval [#nnp].
+    \param x_u (overwritten) dofval [#nnu].
+    */
     void solve_u(
         const xt::xtensor<double, 1>& b_u,
         const xt::xtensor<double, 1>& x_p,
-        xt::xtensor<double, 1>& x_u);
+              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"
+    /**
+    Get right-hand-size for corresponding to the prescribed DOFs.
 
-    void reaction(
-        const xt::xtensor<double, 1>& x, xt::xtensor<double, 1>& b) const; // modified with "b_p"
+    \f$ b_p = A_{pu} * x_u + A_{pp} * x_p = A_{pp} * x_p \equiv A_{pp} * x_p \f$
 
-    void reaction_p(
-        const xt::xtensor<double, 1>& x_u,
-        const xt::xtensor<double, 1>& x_p,
-        xt::xtensor<double, 1>& b_p) const;
+    and assemble them to the appropriate places in "nodevec".
 
-    // Return matrix as diagonal matrix (column)
-    xt::xtensor<double, 1> Todiagonal() const;
+    \param x "nodevec" [#nnode, #ndim].
+    \param b "nodevec" [#nnode, #ndim].
+    \return Copy of `b` with \f$ b_p \f$ overwritten.
+    */
+    xt::xtensor<double, 2> Reaction(
+        const xt::xtensor<double, 2>& x,
+        const xt::xtensor<double, 2>& b) 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;
+    /**
+    Same as Reaction(const xt::xtensor<double, 2>&, const xt::xtensor<double, 2>&),
+    but inserting in-place.
 
-    xt::xtensor<double, 1> Dot_u(
-        const xt::xtensor<double, 1>& x_u, const xt::xtensor<double, 1>& x_p) const;
+    \param x "nodevec" [#nnode, #ndim].
+    \param b "nodevec" [#nnode, #ndim], \f$ b_p \f$ overwritten.
+    */
+    void reaction(
+        const xt::xtensor<double, 2>& x,
+              xt::xtensor<double, 2>& b) const;
 
-    xt::xtensor<double, 1> Dot_p(
-        const xt::xtensor<double, 1>& x_u, const xt::xtensor<double, 1>& x_p) const;
+    /**
+    Same as Reaction(const xt::xtensor<double, 2>&, const xt::xtensor<double, 2>&),
+    but of "dofval" input and output.
 
-    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);
+    \param x "dofval" [#ndof].
+    \param b "dofval" [#ndof].
+    \return Copy of `b` with \f$ b_p \f$ overwritten.
+    */
+    xt::xtensor<double, 1> Reaction(
+        const xt::xtensor<double, 1>& x,
+        const xt::xtensor<double, 1>& b) const;
 
-    xt::xtensor<double, 1> Solve_u(
-        const xt::xtensor<double, 1>& b_u, const xt::xtensor<double, 1>& x_p);
+    /**
+    Same as Reaction(const xt::xtensor<double, 1>&, const xt::xtensor<double, 1>&),
+    but inserting in-place.
 
-    xt::xtensor<double, 2> Reaction(
-        const xt::xtensor<double, 2>& x, const xt::xtensor<double, 2>& b) const;
+    \param x "dofval" [#ndof].
+    \param b "dofval" [#ndof], \f$ b_p \f$ overwritten.
+    */
+    void reaction(
+        const xt::xtensor<double, 1>& x,
+              xt::xtensor<double, 1>& b) const; // modified with "b_p"
 
-    xt::xtensor<double, 1> Reaction(
-        const xt::xtensor<double, 1>& x, const xt::xtensor<double, 1>& b) const;
+    /**
+    Same as Reaction(const xt::xtensor<double, 1>&, const xt::xtensor<double, 1>&),
+    but with partitioned input and output.
 
+    \param x_u unknown "dofval" [#nnu].
+    \param x_p prescribed "dofval" [#nnp].
+    \return b_p prescribed "dofval" [#nnp].
+    */
     xt::xtensor<double, 1> Reaction_p(
-        const xt::xtensor<double, 1>& x_u, const xt::xtensor<double, 1>& x_p) const;
+        const xt::xtensor<double, 1>& x_u,
+        const xt::xtensor<double, 1>& x_p) const;
+
+    /**
+    Same as Reaction_p(const xt::xtensor<double, 1>&, const xt::xtensor<double, 1>&),
+    but writing to preallocated output.
+
+    \param x_u unknown "dofval" [#nnu].
+    \param x_p prescribed "dofval" [#nnp].
+    \param b_p (overwritten) prescribed "dofval" [#nnp].
+    */
+    void reaction_p(
+        const xt::xtensor<double, 1>& x_u,
+        const xt::xtensor<double, 1>& x_p,
+              xt::xtensor<double, 1>& b_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 = 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/MatrixDiagonalPartitioned.hpp b/include/GooseFEM/MatrixDiagonalPartitioned.hpp
index 8b01c24..786829e 100644
--- a/include/GooseFEM/MatrixDiagonalPartitioned.hpp
+++ b/include/GooseFEM/MatrixDiagonalPartitioned.hpp
@@ -1,398 +1,332 @@
 /**
 Implementation of MatrixDiagonalPartitioned.h
 
 \file MatrixDiagonalPartitioned.hpp
 \copyright Copyright 2017. Tom de Geus. All rights reserved.
 \license This project is released under the GNU Public License (GPLv3).
 */
 
 #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)
+    : MatrixDiagonal(conn, 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_iip = iip;
     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}).apply(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_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) {
+    if (!m_changed) {
         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;
+    m_changed = 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;
+    m_changed = 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::Todiagonal() const
 {
     xt::xtensor<double, 1> ret = xt::zeros<double>({m_ndof});
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnu; ++d) {
         ret(m_iiu(d)) = m_Auu(d);
     }
 
     #pragma omp parallel for
     for (size_t d = 0; d < m_nnp; ++d) {
         ret(m_iip(d)) = m_App(d);
     }
 
     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> 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> 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> 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> 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