diff --git a/include/GooseFEM/TyingsPeriodic.h b/include/GooseFEM/TyingsPeriodic.h
index a77d624..2e78926 100644
--- a/include/GooseFEM/TyingsPeriodic.h
+++ b/include/GooseFEM/TyingsPeriodic.h
@@ -1,258 +1,246 @@
 /**
 Tools to store and apply nodal/DOF tyings.
 
 \file TyingsPeriodic.h
 \copyright Copyright 2017. Tom de Geus. All rights reserved.
 \license This project is released under the GNU Public License (GPLv3).
 */
 
 #ifndef GOOSEFEM_TYINGSPERIODIC_H
 #define GOOSEFEM_TYINGSPERIODIC_H
 
 #include "config.h"
 
 #include <Eigen/Eigen>
 #include <Eigen/Sparse>
 
 namespace GooseFEM {
 
 /**
 Tools to store and apply nodal/DOF tyings.
 */
 namespace Tyings {
 
 /**
 Nodal tyings per periodic boundary conditions.
 The idea is that the displacement of all DOFs of a node are tied to another node
 and to the average displacement gradient.
 The latter is applied/measured using 'virtual' control nodes.
 
 Consider the DOF list \f$ u \f$ renumbered such that it is split up in
 independent and dependent DOFs as follows
 
 \f$ u = \begin{bmatrix} u_i \\ u_d \end{bmatrix}\f$
 
 whereby the independent DOFs are furthermore split up in unknown and prescribed nodes as follows
 
 \f$ u_i = \begin{bmatrix} u_u \\ u_p \end{bmatrix}\f$
 
 such that
 
 \f$ u = \begin{bmatrix} u_u \\ u_p \\ u_d \end{bmatrix}\f$
 
 \todo Document how the DOFs are tied to the control nodes, and what the has to do with the mean.
 */
 class Periodic {
 public:
 
     Periodic() = default;
 
     /**
     Constructor.
 
     \param coor Nodal coordinates [nnode, ndim].
     \param dofs DOF-numbers per node [nnode, ndim].
     \param control_dofs DOF-numbers per control node [ndim, ndim].
     \param nodal_tyings List of nodal tyings, see nodal_tyings(). [ntyings, 2].
     */
     Periodic(
         const xt::xtensor<double, 2>& coor,
         const xt::xtensor<size_t, 2>& dofs,
         const xt::xtensor<size_t, 2>& control_dofs,
         const xt::xtensor<size_t, 2>& nodal_tyings);
 
     /**
     Constructor.
 
     \param coor Nodal coordinates [nnode, ndim].
     \param dofs DOF-numbers per node [nnode, ndim].
     \param control_dofs DOF-numbers per control node [ndim, ndim].
     \param nodal_tyings List of nodal tyings, see nodal_tyings(). [ntyings, 2].
     \param iip List of prescribed DOF-numbers.
     */
     Periodic(
         const xt::xtensor<double, 2>& coor,
         const xt::xtensor<size_t, 2>& dofs,
         const xt::xtensor<size_t, 2>& control_dofs,
         const xt::xtensor<size_t, 2>& nodal_tyings,
         const xt::xtensor<size_t, 1>& iip);
 
     /**
-    Number of dependent DOFs.
-
-    \return unsigned int
+    \return Number of dependent DOFs.
     */
     size_t nnd() const;
 
     /**
-    Number of independent DOFs.
-
-    \return unsigned int
+    \return Number of independent DOFs.
     */
     size_t nni() const;
 
     /**
-    Number of independent unknown DOFs.
-
-    \return unsigned int
+    \return Number of independent unknown DOFs.
     */
     size_t nnu() const;
 
     /**
-    Number of independent prescribed DOFs.
-
-    \return unsigned int
+    \return Number of independent prescribed DOFs.
     */
     size_t nnp() const;
 
     /**
-    Return the DOF-numbers per node, as used internally (after renumbering).
-
-    \return [nnode, ndim].
+    \return DOF-numbers per node, as used internally (after renumbering), [nnode, ndim].
     */
     xt::xtensor<size_t, 2> dofs() const;
 
     /**
-    Return the DOF-numbers for each control node, as used internally (after renumbering).
-
-    \return [ndim, ndim].
+    \return DOF-numbers for each control node, as used internally (after renumbering), [ndim, ndim].
     */
     xt::xtensor<size_t, 2> control() const;
 
     /**
     Return the applied nodal tyings.
     Per tying (row) two node numbers are specified,
     according to the convention (independent, dependent).
 
     \return [ntyings, 2].
     */
     xt::xtensor<size_t, 2> nodal_tyings() const;
 
     /**
     Dependent DOFs.
 
     \return List of DOF numbers.
     */
     xt::xtensor<size_t, 1> iid() const;
 
     /**
     Independent DOFs.
 
     \return List of DOF numbers.
     */
     xt::xtensor<size_t, 1> iii() const;
 
     /**
     Independent unknown DOFs.
 
     \return List of DOF numbers.
     */
     xt::xtensor<size_t, 1> iiu() const;
 
     /**
     Independent prescribed DOFs.
 
     \return List of DOF numbers.
     */
     xt::xtensor<size_t, 1> iip() const;
 
     /**
     Return tying matrix such as to get the dependent DOFs \f$ u_d \f$ from
     the independent DOFs \f$ u_i \f$ as follows
 
     \f$ u_d = C_{di} u_i \f$
 
     Note that this can be further partitioned in
 
     \f$ u_d = C_{du} u_u + C_{dp} u_p \f$
 
     See Cdu() and Cdp().
 
     \return Sparse matrix.
     */
     Eigen::SparseMatrix<double> Cdi() const;
 
     /**
     Unknown part of the partitioned tying matrix, see Cdi().
 
     \return Sparse matrix.
     */
     Eigen::SparseMatrix<double> Cdu() const;
 
     /**
     Prescribed part of the partitioned tying matrix, see Cdi().
 
     \return Sparse matrix.
     */
     Eigen::SparseMatrix<double> Cdp() const;
 
 private:
     size_t m_nnu; ///< See nnu().
     size_t m_nnp; ///< See nnp().
     size_t m_nni; ///< See nni().
     size_t m_nnd; ///< See nnd().
     size_t m_ndim; ///< Number of dimensions.
     size_t m_nties; ///< Number of nodal ties.
     xt::xtensor<size_t, 2> m_dofs; ///< See dofs().
     xt::xtensor<size_t, 2> m_control; ///< See control().
     xt::xtensor<size_t, 2> m_tyings; ///< See nodal_tyings().
     xt::xtensor<double, 2> m_coor; ///< Nodal coordinates [nnode, ndim].
 };
 
 /**
 Add control nodes to an existing system.
 */
 class Control {
 public:
 
     Control() = default;
 
     /**
     Constructor.
 
     \param coor Nodal coordinates [nnode, ndim].
     \param dofs DOF-numbers per node [nnode, ndim].
     */
     Control(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& dofs);
 
     /**
     Nodal coordinates, for the system with control nodes added to it.
 
     \param [nnode + ndim, ndim], with nnode the number of nodes of the original system.
     */
     xt::xtensor<double, 2> coor() const;
 
     /**
     DOF-numbers per node, for the system with control nodes added to it.
 
     \param [nnode + ndim, ndim], with nnode the number of nodes of the original system.
     */
     xt::xtensor<size_t, 2> dofs() const;
 
     /**
     DOF-numbers of each control node.
 
     \param [ndim, ndim].
     */
     xt::xtensor<size_t, 2> controlDofs() const;
 
     /**
     Node-numbers of the control nodes.
 
     \param [ndim].
     */
     xt::xtensor<size_t, 1> controlNodes() const;
 
 private:
     xt::xtensor<double, 2> m_coor; ///< See coor().
     xt::xtensor<size_t, 2> m_dofs; ///< See dofs().
     xt::xtensor<size_t, 2> m_control_dofs;  ///< See controlDofs().
     xt::xtensor<size_t, 1> m_control_nodes; ///< See controlNodes().
 };
 
 } // namespace Tyings
 } // namespace GooseFEM
 
 #include "TyingsPeriodic.hpp"
 
 #endif
diff --git a/include/GooseFEM/Vector.h b/include/GooseFEM/Vector.h
index 1bb2d27..17afda4 100644
--- a/include/GooseFEM/Vector.h
+++ b/include/GooseFEM/Vector.h
@@ -1,392 +1,392 @@
 /**
 Methods to switch between storage types based on a mesh.
 
 \file Vector.h
 \copyright Copyright 2017. Tom de Geus. All rights reserved.
 \license This project is released under the GNU Public License (GPLv3).
 */
 
 #ifndef GOOSEFEM_VECTOR_H
 #define GOOSEFEM_VECTOR_H
 
 #include "config.h"
 
 namespace GooseFEM {
 
 /**
-Class to switch between storage types:
+Class to switch between storage types. In particular:
 
 -   "dofval": DOF values [#ndof].
 -   "nodevec": nodal vectors [#nnode, #ndim].
 -   "elemvec": nodal vectors stored per element [#nelem, #nne, #ndim].
 */
 class Vector {
 public:
 
     Vector() = default;
 
     /**
     Constructor.
 
     \param conn connectivity [#nelem, #nne].
     \param dofs DOFs per node [#nnode, #ndim].
     */
     Vector(const xt::xtensor<size_t, 2>& conn, const xt::xtensor<size_t, 2>& dofs);
 
     /**
     Number of elements.
 
     \return unsigned int
     */
     size_t nelem() const;
 
     /**
     Number of nodes per element.
 
     \return unsigned int
     */
     size_t nne() const;
 
     /**
     Number of nodes.
 
     \return unsigned int
     */
     size_t nnode() const;
 
     /**
     Number of dimensions.
 
     \return unsigned int
     */
     size_t ndim() const;
 
     /**
     Number of DOFs.
 
     \return unsigned int
     */
     size_t ndof() const;
 
     /**
     Connectivity (nodes per element).
 
     \return [#nelem, #nne]
     */
     xt::xtensor<size_t, 2> conn() const;
 
     /**
     DOFs per node.
 
     \return [#nnode, #ndim]
     */
     xt::xtensor<size_t, 2> dofs() const;
 
     /**
     Copy "nodevec" to another "nodevec".
 
     \param nodevec_src input [#nnode, #ndim]
     \param nodevec_dest input [#nnode, #ndim]
     \return nodevec output [#nnode, #ndim]
     */
     xt::xtensor<double, 2> Copy(
         const xt::xtensor<double, 2>& nodevec_src,
         const xt::xtensor<double, 2>& nodevec_dest) const;
 
     /**
     Copy "nodevec" to another "nodevec".
 
     \param nodevec_src input [#nnode, #ndim]
     \param nodevec_dest output [#nnode, #ndim]
     */
     void copy(const xt::xtensor<double, 2>& nodevec_src, xt::xtensor<double, 2>& nodevec_dest) const;
 
     /**
     Convert "nodevec" to "dofval" (overwrite entries that occur more than once).
 
     \param nodevec input [#nnode, #ndim]
     \return dofval output [#ndof]
     */
     xt::xtensor<double, 1> AsDofs(const xt::xtensor<double, 2>& nodevec) const;
 
     /**
     Convert "nodevec" to "dofval" (overwrite entries that occur more than once).
 
     \param nodevec input [#nnode, #ndim]
     \param dofval output [#ndof]
     */
     void asDofs(const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 1>& dofval) const;
 
     /**
     Convert "elemvec" to "dofval" (overwrite entries that occur more than once).
 
     \param elemvec input [#nelem, #nne, #ndim]
     \return dofval output [#ndof]
     */
     xt::xtensor<double, 1> AsDofs(const xt::xtensor<double, 3>& elemvec) const;
 
     /**
     Convert "elemvec" to "dofval" (overwrite entries that occur more than once).
 
     \param elemvec input [#nelem, #nne, #ndim]
     \param dofval output [#ndof]
     */
     void asDofs(const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 1>& dofval) const;
 
     /**
     Convert "dofval" to "nodevec" (overwrite entries that occur more than once).
 
     \param dofval input [#ndof]
     \return nodevec output [#nnode, #ndim]
     */
     xt::xtensor<double, 2> AsNode(const xt::xtensor<double, 1>& dofval) const;
 
     /**
     Convert "dofval" to "nodevec" (overwrite entries that occur more than once).
 
     \param dofval input [#ndof]
     \param nodevec input [#nnode, #ndim]
     */
     void asNode(const xt::xtensor<double, 1>& dofval, xt::xtensor<double, 2>& nodevec) const;
 
     /**
     Convert "elemvec" to "nodevec" (overwrite entries that occur more than once).
 
     \param elemvec input [#nelem, #nne, #ndim]
     \return nodevec output [#nnode, #ndim]
     */
     xt::xtensor<double, 2> AsNode(const xt::xtensor<double, 3>& elemvec) const;
 
     /**
     Convert "elemvec" to "nodevec" (overwrite entries that occur more than once).
 
     \param elemvec input [#nelem, #nne, #ndim]
     \param nodevec [#nnode, #ndim]
     */
     void asNode(const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 2>& nodevec) const;
 
     /**
     Convert "dofval" to "elemvec" (overwrite entries that occur more than once).
 
     \param dofval input [#ndof].
     \return elemvec output [#nelem, #nne, #ndim].
     */
     xt::xtensor<double, 3> AsElement(const xt::xtensor<double, 1>& dofval) const;
 
     /**
     Convert "dofval" to "elemvec" (overwrite entries that occur more than once).
 
     \param dofval input [#ndof].
     \param elemvec output [#nelem, #nne, #ndim].
     */
     void asElement(const xt::xtensor<double, 1>& dofval, xt::xtensor<double, 3>& elemvec) const;
 
     /**
     Convert "nodevec" to "elemvec" (overwrite entries that occur more than once).
 
     \param nodevec input [#nnode, #ndim].
     \return elemvec output [#nelem, #nne, #ndim].
     */
     xt::xtensor<double, 3> AsElement(const xt::xtensor<double, 2>& nodevec) const;
 
     /**
     Convert "nodevec" to "elemvec" (overwrite entries that occur more than once).
 
     \param nodevec input [#nnode, #ndim].
     \param elemvec output [#nelem, #nne, #ndim].
     */
     void asElement(const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 3>& elemvec) const;
 
     /**
     Assemble "nodevec" to "dofval" (adds entries that occur more that once).
 
     \param nodevec input [#nnode, #ndim]
     \return dofval output [#ndof]
     */
     xt::xtensor<double, 1> AssembleDofs(const xt::xtensor<double, 2>& nodevec) const;
 
     /**
     Assemble "nodevec" to "dofval" (adds entries that occur more that once).
 
     \param nodevec input [#nnode, #ndim]
     \param dofval output [#ndof]
     */
     void assembleDofs(const xt::xtensor<double, 2>& nodevec, xt::xtensor<double, 1>& dofval) const;
 
     /**
     Assemble "elemvec" to "dofval" (adds entries that occur more that once).
 
     \param elemvec input [#nelem, #nne, #ndim]
     \return dofval output [#ndof]
     */
     xt::xtensor<double, 1> AssembleDofs(const xt::xtensor<double, 3>& elemvec) const;
 
     /**
     Assemble "elemvec" to "dofval" (adds entries that occur more that once).
 
     \param elemvec input [#nelem, #nne, #ndim]
     \param dofval output [#ndof]
     */
     void assembleDofs(const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 1>& dofval) const;
 
     /**
     Assemble "elemvec" to "nodevec" (adds entries that occur more that once.
 
     \param elemvec input [#nelem, #nne, #ndim]
     \return nodevec output [#nnode, #ndim]
     */
     xt::xtensor<double, 2> AssembleNode(const xt::xtensor<double, 3>& elemvec) const;
 
     /**
     Assemble "elemvec" to "nodevec" (adds entries that occur more that once.
 
     \param elemvec input [#nelem, #nne, #ndim]
     \param nodevec output [#nnode, #ndim]
     */
     void assembleNode(const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 2>& nodevec) const;
 
     /**
     Shape of "dofval".
 
     \return [#ndof]
     */
     std::array<size_t, 1> shape_dofval() const;
 
     /**
     Shape of "nodevec".
 
     \return [#nnode, #ndim]
     */
     std::array<size_t, 2> shape_nodevec() const;
 
     /**
     Shape of "elemvec".
 
     \return [#nelem, #nne, #ndim]
     */
     std::array<size_t, 3> shape_elemvec() const;
 
     /**
     Shape of "elemmat".
 
     \return [#nelem, #nne * #ndim, #nne * #ndim]
     */
     std::array<size_t, 3> shape_elemmat() const;
 
     /**
     Allocated "dofval".
 
     \return [#ndof]
     */
     xt::xtensor<double, 1> allocate_dofval() const;
 
     /**
     Allocated and initialised "dofval".
 
     \param val value to which to initialise.
     \return [#ndof]
     */
     xt::xtensor<double, 1> allocate_dofval(double val) const;
 
     /**
     Allocated "nodevec".
 
     \return [#nnode, #ndim]
     */
     xt::xtensor<double, 2> allocate_nodevec() const;
 
     /**
     Allocated and initialised "nodevec".
 
     \param val value to which to initialise.
     \return [#nnode, #ndim]
     */
     xt::xtensor<double, 2> allocate_nodevec(double val) const;
 
     /**
     Allocated "elemvec".
 
     \return [#nelem, #nne, #ndim]
     */
     xt::xtensor<double, 3> allocate_elemvec() const;
 
     /**
     Allocated and initialised "elemvec".
 
     \param val value to which to initialise.
     \return [#nelem, #nne, #ndim]
     */
     xt::xtensor<double, 3> allocate_elemvec(double val) const;
 
     /**
     Allocated "elemmat".
 
     \return [#nelem, #nne * #ndim, #nne * #ndim]
     */
     xt::xtensor<double, 3> allocate_elemmat() const;
 
     /**
     Allocated and initialised "elemmat".
 
     \param val value to which to initialise.
     \return [#nelem, #nne * #ndim, #nne * #ndim]
     */
     xt::xtensor<double, 3> allocate_elemmat(double val) const;
 
     /**
     \cond
     */
     [[ deprecated ]]
     std::array<size_t, 1> ShapeDofval() const;
 
     [[ deprecated ]]
     std::array<size_t, 2> ShapeNodevec() const;
 
     [[ deprecated ]]
     std::array<size_t, 3> ShapeElemvec() const;
 
     [[ deprecated ]]
     std::array<size_t, 3> ShapeElemmat() const;
 
     [[ deprecated ]]
     xt::xtensor<double, 1> AllocateDofval() const;
 
     [[ deprecated ]]
     xt::xtensor<double, 1> AllocateDofval(double val) const;
 
     [[ deprecated ]]
     xt::xtensor<double, 2> AllocateNodevec() const;
 
     [[ deprecated ]]
     xt::xtensor<double, 2> AllocateNodevec(double val) const;
 
     [[ deprecated ]]
     xt::xtensor<double, 3> AllocateElemvec() const;
 
     [[ deprecated ]]
     xt::xtensor<double, 3> AllocateElemvec(double val) const;
 
     [[ deprecated ]]
     xt::xtensor<double, 3> AllocateElemmat() const;
 
     [[ deprecated ]]
     xt::xtensor<double, 3> AllocateElemmat(double val) const;
     /**
     \endcond
     */
 
 protected:
     xt::xtensor<size_t, 2> m_conn; ///< See conn()
     xt::xtensor<size_t, 2> m_dofs; ///< See dofs()
     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
 };
 
 } // namespace GooseFEM
 
 #include "Vector.hpp"
 
 #endif