diff --git a/include/GooseFEM/TyingsPeriodic.h b/include/GooseFEM/TyingsPeriodic.h
index 81cd680..d12ae4a 100644
--- a/include/GooseFEM/TyingsPeriodic.h
+++ b/include/GooseFEM/TyingsPeriodic.h
@@ -1,96 +1,226 @@
 /**
-Methods to store and apply DOF tyings.
+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:
-    // Constructors
+
     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); // (independent, dependent)
+        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, // (independent, dependent)
+        const xt::xtensor<size_t, 2>& nodal_tyings,
         const xt::xtensor<size_t, 1>& iip);
 
-    // Dimensions
-    size_t nnd() const; // dependent DOFs
-    size_t nni() const; // independent DOFs
-    size_t nnu() const; // independent, unknown DOFs
-    size_t nnp() const; // independent, prescribed DOFs
-
-    // DOF lists
-    xt::xtensor<size_t, 2> dofs() const;    // DOFs
-    xt::xtensor<size_t, 2> control() const; // control DOFs
-    xt::xtensor<size_t, 1> iid() const;     // dependent DOFs
-    xt::xtensor<size_t, 1> iii() const;     // independent DOFs
-    xt::xtensor<size_t, 1> iiu() const;     // independent, unknown DOFs
-    xt::xtensor<size_t, 1> iip() const;     // independent, prescribed DOFs
-
-    // Return the tying matrix
-    // u_d = C_di * u_i
-    // u_d = [C_du, C_dp]^T * [u_u, u_p] = C_du * u_u + C_dp * u_p
+    /**
+    Number of dependent DOFs.
+
+    \return unsigned int
+    */
+    size_t nnd() const;
+
+    /**
+    Number of independent DOFs.
+
+    \return unsigned int
+    */
+    size_t nni() const;
+
+    /**
+    Number of independent unknown DOFs.
+
+    \return unsigned int
+    */
+    size_t nnu() const;
+
+    /**
+    Number of independent prescribed DOFs.
+
+    \return unsigned int
+    */
+    size_t nnp() const;
+
+    /**
+    Return the DOF-numbers per node, as used internally (after renumbering).
+
+    \return [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].
+    */
+    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;
-    size_t m_nnp;
-    size_t m_nni;
-    size_t m_nnd;
-    size_t m_ndim;
-    size_t m_nties; // number of nodal ties
-    xt::xtensor<size_t, 2> m_dofs;
-    xt::xtensor<size_t, 2> m_control;
-    xt::xtensor<size_t, 2> m_tyings; // nodal ties: (independent, dependent)
-    xt::xtensor<double, 2> m_coor;
+    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].
 };
 
 class Control {
 public:
     // Constructors
     Control() = default;
     Control(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& dofs);
 
     // Extract new lists
     xt::xtensor<double, 2> coor() const;
     xt::xtensor<size_t, 2> dofs() const;
     xt::xtensor<size_t, 2> controlDofs() const;
     xt::xtensor<size_t, 1> controlNodes() const;
 
 private:
     xt::xtensor<double, 2> m_coor;
     xt::xtensor<size_t, 2> m_dofs;
     xt::xtensor<size_t, 2> m_control_dofs;
     xt::xtensor<size_t, 1> m_control_nodes;
 };
 
 } // namespace Tyings
 } // namespace GooseFEM
 
 #include "TyingsPeriodic.hpp"
 
 #endif
diff --git a/include/GooseFEM/TyingsPeriodic.hpp b/include/GooseFEM/TyingsPeriodic.hpp
index 7b6619a..7acb022 100644
--- a/include/GooseFEM/TyingsPeriodic.hpp
+++ b/include/GooseFEM/TyingsPeriodic.hpp
@@ -1,241 +1,255 @@
 /**
 Implementation of TyingsPeriodic.h
 
 \file TyingsPeriodic.hpp
 \copyright Copyright 2017. Tom de Geus. All rights reserved.
 \license This project is released under the GNU Public License (GPLv3).
 */
 
 #ifndef GOOSEFEM_TYINGSPERIODIC_HPP
 #define GOOSEFEM_TYINGSPERIODIC_HPP
 
 #include "Mesh.h"
 #include "TyingsPeriodic.h"
 
 #include <Eigen/Eigen>
 #include <Eigen/Sparse>
 
 namespace GooseFEM {
 namespace Tyings {
 
 inline Periodic::Periodic(
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& dofs,
     const xt::xtensor<size_t, 2>& control,
     const xt::xtensor<size_t, 2>& nodal_tyings)
     : Periodic(coor, dofs, control, nodal_tyings, xt::empty<size_t>({0}))
 {
 }
 
 inline Periodic::Periodic(
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& dofs,
     const xt::xtensor<size_t, 2>& control,
     const xt::xtensor<size_t, 2>& nodal_tyings,
     const xt::xtensor<size_t, 1>& iip)
     : m_tyings(nodal_tyings), m_coor(coor)
 {
     m_ndim = m_coor.shape(1);
     m_nties = m_tyings.shape(0);
 
+    GOOSEFEM_ASSERT(xt::has_shape(m_tyings, {m_nties, size_t(2)}));
+    GOOSEFEM_ASSERT(xt::has_shape(control, {m_ndim, m_ndim}));
+    GOOSEFEM_ASSERT(xt::has_shape(dofs, m_coor.shape()));
+    GOOSEFEM_ASSERT(xt::amax(control)() <= xt::amax(dofs)());
+    GOOSEFEM_ASSERT(xt::amin(control)() >= xt::amin(dofs)());
+    GOOSEFEM_ASSERT(xt::amax(iip)() <= xt::amax(dofs)());
+    GOOSEFEM_ASSERT(xt::amin(iip)() >= xt::amin(dofs)());
+    GOOSEFEM_ASSERT(xt::amax(nodal_tyings)() < m_coor.shape(0));
+
     xt::xtensor<size_t, 1> dependent = xt::view(m_tyings, xt::all(), 1);
     xt::xtensor<size_t, 2> dependent_dofs = xt::view(dofs, xt::keep(dependent), xt::all());
     xt::xtensor<size_t, 1> iid = xt::flatten(dependent_dofs);
     xt::xtensor<size_t, 1> iii = xt::setdiff1d(dofs, iid);
     xt::xtensor<size_t, 1> iiu = xt::setdiff1d(iii, iip);
 
     m_nnu = iiu.size();
     m_nnp = iip.size();
     m_nni = iii.size();
     m_nnd = iid.size();
 
     GooseFEM::Mesh::Reorder reorder({iiu, iip, iid});
 
     m_dofs = reorder.apply(dofs);
     m_control = reorder.apply(control);
 }
 
 inline xt::xtensor<size_t, 2> Periodic::dofs() const
 {
     return m_dofs;
 }
 
 inline xt::xtensor<size_t, 2> Periodic::control() const
 {
     return m_control;
 }
 
+inline xt::xtensor<size_t, 2> Periodic::nodal_tyings() const
+{
+    return m_tyings;
+}
+
 inline size_t Periodic::nnu() const
 {
     return m_nnu;
 }
 
 inline size_t Periodic::nnp() const
 {
     return m_nnp;
 }
 
 inline size_t Periodic::nni() const
 {
     return m_nni;
 }
 
 inline size_t Periodic::nnd() const
 {
     return m_nnd;
 }
 
 inline xt::xtensor<size_t, 1> Periodic::iiu() const
 {
     return xt::arange<size_t>(m_nnu);
 }
 
 inline xt::xtensor<size_t, 1> Periodic::iip() const
 {
     return xt::arange<size_t>(m_nnp) + m_nnu;
 }
 
 inline xt::xtensor<size_t, 1> Periodic::iii() const
 {
     return xt::arange<size_t>(m_nni);
 }
 
 inline xt::xtensor<size_t, 1> Periodic::iid() const
 {
     return xt::arange<size_t>(m_nni, m_nni + m_nnd);
 }
 
 inline Eigen::SparseMatrix<double> Periodic::Cdi() const
 {
     std::vector<Eigen::Triplet<double>> data;
 
     data.reserve(m_nties * m_ndim * (m_ndim + 1));
 
     for (size_t i = 0; i < m_nties; ++i) {
         for (size_t j = 0; j < m_ndim; ++j) {
 
             size_t ni = m_tyings(i, 0);
             size_t nd = m_tyings(i, 1);
 
             data.push_back(Eigen::Triplet<double>(i * m_ndim + j, m_dofs(ni, j), +1.0));
 
             for (size_t k = 0; k < m_ndim; ++k) {
                 data.push_back(Eigen::Triplet<double>(
                     i * m_ndim + j, m_control(j, k), m_coor(nd, k) - m_coor(ni, k)));
             }
         }
     }
 
     Eigen::SparseMatrix<double> Cdi;
     Cdi.resize(m_nnd, m_nni);
     Cdi.setFromTriplets(data.begin(), data.end());
 
     return Cdi;
 }
 
 inline Eigen::SparseMatrix<double> Periodic::Cdu() const
 {
     std::vector<Eigen::Triplet<double>> data;
 
     data.reserve(m_nties * m_ndim * (m_ndim + 1));
 
     for (size_t i = 0; i < m_nties; ++i) {
         for (size_t j = 0; j < m_ndim; ++j) {
 
             size_t ni = m_tyings(i, 0);
             size_t nd = m_tyings(i, 1);
 
             if (m_dofs(ni, j) < m_nnu) {
                 data.push_back(Eigen::Triplet<double>(i * m_ndim + j, m_dofs(ni, j), +1.0));
             }
 
             for (size_t k = 0; k < m_ndim; ++k) {
                 if (m_control(j, k) < m_nnu) {
                     data.push_back(Eigen::Triplet<double>(
                         i * m_ndim + j, m_control(j, k), m_coor(nd, k) - m_coor(ni, k)));
                 }
             }
         }
     }
 
     Eigen::SparseMatrix<double> Cdu;
     Cdu.resize(m_nnd, m_nnu);
     Cdu.setFromTriplets(data.begin(), data.end());
 
     return Cdu;
 }
 
 inline Eigen::SparseMatrix<double> Periodic::Cdp() const
 {
     std::vector<Eigen::Triplet<double>> data;
 
     data.reserve(m_nties * m_ndim * (m_ndim + 1));
 
     for (size_t i = 0; i < m_nties; ++i) {
         for (size_t j = 0; j < m_ndim; ++j) {
 
             size_t ni = m_tyings(i, 0);
             size_t nd = m_tyings(i, 1);
 
             if (m_dofs(ni, j) >= m_nnu) {
                 data.push_back(Eigen::Triplet<double>(i * m_ndim + j, m_dofs(ni, j) - m_nnu, +1.0));
             }
 
             for (size_t k = 0; k < m_ndim; ++k) {
                 if (m_control(j, k) >= m_nnu) {
                     data.push_back(Eigen::Triplet<double>(
                         i * m_ndim + j, m_control(j, k) - m_nnu, m_coor(nd, k) - m_coor(ni, k)));
                 }
             }
         }
     }
 
     Eigen::SparseMatrix<double> Cdp;
     Cdp.resize(m_nnd, m_nnp);
     Cdp.setFromTriplets(data.begin(), data.end());
 
     return Cdp;
 }
 
 inline Control::Control(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& dofs)
     : m_coor(coor), m_dofs(dofs)
 {
     GOOSEFEM_ASSERT(coor.shape().size() == 2);
     GOOSEFEM_ASSERT(coor.shape() == dofs.shape());
 
     size_t nnode = coor.shape(0);
     size_t ndim = coor.shape(1);
 
     m_control_dofs = xt::arange<size_t>(ndim * ndim).reshape({ndim, ndim});
     m_control_dofs += xt::amax(dofs)() + 1;
 
     m_control_nodes = nnode + xt::arange<size_t>(ndim);
 
     m_coor = xt::concatenate(xt::xtuple(coor, xt::zeros<double>({ndim, ndim})));
     m_dofs = xt::concatenate(xt::xtuple(dofs, m_control_dofs));
 }
 
 inline xt::xtensor<double, 2> Control::coor() const
 {
     return m_coor;
 }
 
 inline xt::xtensor<size_t, 2> Control::dofs() const
 {
     return m_dofs;
 }
 
 inline xt::xtensor<size_t, 2> Control::controlDofs() const
 {
     return m_control_dofs;
 }
 
 inline xt::xtensor<size_t, 1> Control::controlNodes() const
 {
     return m_control_nodes;
 }
 
 } // namespace Tyings
 } // namespace GooseFEM
 
 #endif