diff --git a/include/GooseFEM/MatrixPartitionedTyings.hpp b/include/GooseFEM/MatrixPartitionedTyings.hpp
index 32caec5..8affee1 100644
--- a/include/GooseFEM/MatrixPartitionedTyings.hpp
+++ b/include/GooseFEM/MatrixPartitionedTyings.hpp
@@ -1,393 +1,393 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
 #ifndef GOOSEFEM_MATRIXPARTITIONEDTYINGS_HPP
 #define GOOSEFEM_MATRIXPARTITIONEDTYINGS_HPP
 
 // -------------------------------------------------------------------------------------------------
 
 #include "MatrixPartitionedTyings.h"
 
 // =================================================================================================
 
 namespace GooseFEM {
 
 // -------------------------------------------------------------------------------------------------
 
 inline MatrixPartitionedTyings::MatrixPartitionedTyings(
-  const xt::xtensor<size_t,2> &conn,
-  const xt::xtensor<size_t,2> &dofs,
+  const xt::xtensor<size_t,2>& conn,
+  const xt::xtensor<size_t,2>& dofs,
   const Eigen::SparseMatrix<double>& Cdu,
   const Eigen::SparseMatrix<double>& Cdp) :
   m_conn(conn),
   m_dofs(dofs),
   m_Cdu(Cdu),
   m_Cdp(Cdp)
 {
   GOOSEFEM_ASSERT(Cdu.rows() == Cdp.rows());
 
   m_nnu   = static_cast<size_t>(m_Cdu.cols());
   m_nnp   = static_cast<size_t>(m_Cdp.cols());
   m_nnd   = static_cast<size_t>(m_Cdp.rows());
   m_nni   = m_nnu + m_nnp;
   m_ndof  = m_nni + m_nnd;
 
   m_iiu   = xt::arange<size_t>(m_nnu);
-  m_iip   = xt::arange<size_t>(m_nnp) + m_nnu;
-  m_iid   = xt::arange<size_t>(m_nnd) + m_nni;
+  m_iip   = xt::arange<size_t>(m_nnu, m_nnu + m_nnp);
+  m_iid   = xt::arange<size_t>(m_nni, m_nni + m_nnd);
 
   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_Cud   = m_Cdu.transpose();
   m_Cpd   = m_Cdp.transpose();
 
   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_Tud.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim);
   m_Tpd.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim);
   m_Tdu.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim);
   m_Tdp.reserve(m_nelem*m_nne*m_ndim*m_nne*m_ndim);
   m_Tdd.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);
   m_Aud.resize(m_nnu,m_nnd);
   m_Apd.resize(m_nnp,m_nnd);
   m_Adu.resize(m_nnd,m_nnu);
   m_Adp.resize(m_nnd,m_nnp);
   m_Add.resize(m_nnd,m_nnd);
 
   GOOSEFEM_ASSERT(xt::amax(m_conn)[0] + 1 == m_nnode);
   GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim);
   GOOSEFEM_ASSERT(m_ndof == xt::amax(m_dofs)[0] + 1);
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline size_t MatrixPartitionedTyings::nelem() const
 { return m_nelem; }
 
 inline size_t MatrixPartitionedTyings::nne() const
 { return m_nne; }
 
 inline size_t MatrixPartitionedTyings::nnode() const
 { return m_nnode; }
 
 inline size_t MatrixPartitionedTyings::ndim() const
 { return m_ndim; }
 
 inline size_t MatrixPartitionedTyings::ndof() const
 { return m_ndof; }
 
 inline size_t MatrixPartitionedTyings::nnu() const
 { return m_nnu; }
 
 inline size_t MatrixPartitionedTyings::nnp() const
 { return m_nnp; }
 
 inline xt::xtensor<size_t,2> MatrixPartitionedTyings::dofs() const
 { return m_dofs; }
 
 inline xt::xtensor<size_t,1> MatrixPartitionedTyings::iiu() const
 { return m_iiu; }
 
 inline xt::xtensor<size_t,1> MatrixPartitionedTyings::iip() const
 { return m_iip; }
 
 inline xt::xtensor<size_t,1> MatrixPartitionedTyings::iii() const
 { return xt::arange<size_t>(m_nni); }
 
 inline xt::xtensor<size_t,1> MatrixPartitionedTyings::iid() const
 { return m_iid; }
 
 // -------------------------------------------------------------------------------------------------
 
 inline void MatrixPartitionedTyings::factorize()
 {
   if (!m_factor) return;
 
   m_ACuu = m_Auu + m_Aud * m_Cdu + m_Cud * m_Adu + m_Cud * m_Add * m_Cdu;
   m_ACup = m_Aup + m_Aud * m_Cdp + m_Cud * m_Adp + m_Cud * m_Add * m_Cdp;
   // m_ACpu = m_Apu + m_Apd * m_Cdu + m_Cpd * m_Adu + m_Cpd * m_Add * m_Cdu;
   // m_ACpp = m_App + m_Apd * m_Cdp + m_Cpd * m_Adp + m_Cpd * m_Add * m_Cdp;
 
   m_solver.compute(m_ACuu);
 
   m_factor = false;
 }
 
 // -------------------------------------------------------------------------------------------------
 
-inline void MatrixPartitionedTyings::assemble(const xt::xtensor<double,3> &elemmat)
+inline void MatrixPartitionedTyings::assemble(const xt::xtensor<double,3>& elemmat)
 {
   GOOSEFEM_ASSERT(elemmat.shape() ==\
     std::decay_t<decltype(elemmat)>::shape_type({m_nelem, m_nne*m_ndim, m_nne*m_ndim}));
 
   m_Tuu.clear();
   m_Tup.clear();
   m_Tpu.clear();
   m_Tpp.clear();
   m_Tud.clear();
   m_Tpd.clear();
   m_Tdu.clear();
   m_Tdp.clear();
   m_Tdd.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_dofs(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_dofs(m_conn(e,n),j);
 
             if      (di < m_nnu and 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 and dj < m_nni)
               m_Tup.push_back(Eigen::Triplet<double>(di      ,dj-m_nnu,elemmat(e,m*m_ndim+i,n*m_ndim+j)));
             else if (di < m_nnu)
               m_Tud.push_back(Eigen::Triplet<double>(di      ,dj-m_nni,elemmat(e,m*m_ndim+i,n*m_ndim+j)));
             else if (di < m_nni and 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 if (di < m_nni and dj < m_nni)
               m_Tpp.push_back(Eigen::Triplet<double>(di-m_nnu,dj-m_nnu,elemmat(e,m*m_ndim+i,n*m_ndim+j)));
             else if (di < m_nni)
               m_Tpd.push_back(Eigen::Triplet<double>(di-m_nnu,dj-m_nni,elemmat(e,m*m_ndim+i,n*m_ndim+j)));
             else if (dj < m_nnu)
               m_Tdu.push_back(Eigen::Triplet<double>(di-m_nni,dj      ,elemmat(e,m*m_ndim+i,n*m_ndim+j)));
             else if (dj < m_nni)
               m_Tdp.push_back(Eigen::Triplet<double>(di-m_nni,dj-m_nnu,elemmat(e,m*m_ndim+i,n*m_ndim+j)));
             else
               m_Tdd.push_back(Eigen::Triplet<double>(di-m_nni,dj-m_nni,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_Aud.setFromTriplets(m_Tud.begin(), m_Tud.end());
   m_Apd.setFromTriplets(m_Tpd.begin(), m_Tpd.end());
   m_Adu.setFromTriplets(m_Tdu.begin(), m_Tdu.end());
   m_Adp.setFromTriplets(m_Tdp.begin(), m_Tdp.end());
   m_Add.setFromTriplets(m_Tdd.begin(), m_Tdd.end());
 
   m_factor = true;
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline void MatrixPartitionedTyings::solve(
-  const xt::xtensor<double,2> &b,
-        xt::xtensor<double,2> &x)
+  const xt::xtensor<double,2>& b,
+        xt::xtensor<double,2>& x)
 {
   GOOSEFEM_ASSERT(b.shape() ==\
     std::decay_t<decltype(b)>::shape_type({m_nnode, m_ndim}));
   GOOSEFEM_ASSERT(x.shape() ==\
     std::decay_t<decltype(x)>::shape_type({m_nnode, m_ndim}));
 
   this->factorize();
 
   Eigen::VectorXd B_u = this->asDofs_u(b);
   Eigen::VectorXd B_d = this->asDofs_d(b);
   Eigen::VectorXd X_p = this->asDofs_p(x);
 
   B_u += m_Cud * B_d;
 
   Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - m_ACup * X_p));
   Eigen::VectorXd X_d = m_Cdu * X_u + m_Cdp * 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_dofs(m,i) < m_nnu)
         x(m,i) = X_u(m_dofs(m,i));
       else if (m_dofs(m,i) >= m_nni)
         x(m,i) = X_d(m_dofs(m,i)-m_nni);
     }
   }
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline void MatrixPartitionedTyings::solve(
-  const xt::xtensor<double,1> &b,
-        xt::xtensor<double,1> &x)
+  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();
 
   Eigen::VectorXd B_u = this->asDofs_u(b);
   Eigen::VectorXd B_d = this->asDofs_d(b);
   Eigen::VectorXd X_p = this->asDofs_p(x);
 
   Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - m_ACup * X_p));
   Eigen::VectorXd X_d = m_Cdu * X_u + m_Cdp * X_p;
 
   #pragma omp parallel for
   for (size_t d = 0 ; d < m_nnu ; ++d)
     x(m_iiu(d)) = X_u(d);
 
   #pragma omp parallel for
   for (size_t d = 0 ; d < m_nnd ; ++d)
     x(m_iid(d)) = X_d(d);
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline void MatrixPartitionedTyings::solve_u(
-  const xt::xtensor<double,1> &b_u,
-  const xt::xtensor<double,1> &b_d,
-  const xt::xtensor<double,1> &x_p,
-        xt::xtensor<double,1> &x_u)
+  const xt::xtensor<double,1>& b_u,
+  const xt::xtensor<double,1>& b_d,
+  const xt::xtensor<double,1>& x_p,
+        xt::xtensor<double,1>& x_u)
 {
   GOOSEFEM_ASSERT(b_u.size() == m_nnu);
   GOOSEFEM_ASSERT(b_d.size() == m_nnd);
   GOOSEFEM_ASSERT(x_p.size() == m_nnp);
   GOOSEFEM_ASSERT(x_u.size() == m_nnu);
 
   this->factorize();
 
   Eigen::VectorXd B_u(m_nnu,1);
   Eigen::VectorXd B_d(m_nnd,1);
   Eigen::VectorXd X_p(m_nnp,1);
 
   std::copy(b_u.begin(), b_u.end(), B_u.data());
   std::copy(b_d.begin(), b_d.end(), B_d.data());
   std::copy(x_p.begin(), x_p.end(), X_p.data());
 
   Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - m_ACup * X_p));
 
   std::copy(X_u.data(), X_u.data()+m_nnu, x_u.begin());
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_u(
-  const xt::xtensor<double,1> &dofval) const
+  const xt::xtensor<double,1>& dofval) const
 {
   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 MatrixPartitionedTyings::asDofs_u(
-  const xt::xtensor<double,2> &nodevec) const
+  const xt::xtensor<double,2>& nodevec) const
 {
   assert(nodevec.shape() ==\
     std::decay_t<decltype(nodevec)>::shape_type({m_nnode, m_ndim}));
 
   Eigen::VectorXd dofval_u(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_dofs(m,i) < m_nnu)
         dofval_u(m_dofs(m,i)) = nodevec(m,i);
 
   return dofval_u;
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_p(
-  const xt::xtensor<double,1> &dofval) const
+  const xt::xtensor<double,1>& dofval) const
 {
   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 MatrixPartitionedTyings::asDofs_p(
-  const xt::xtensor<double,2> &nodevec) const
+  const xt::xtensor<double,2>& nodevec) const
 {
   assert(nodevec.shape() ==\
     std::decay_t<decltype(nodevec)>::shape_type({m_nnode, m_ndim}));
 
   Eigen::VectorXd dofval_p(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_dofs(m,i) >= m_nnu and m_dofs(m,i) < m_nni)
         dofval_p(m_dofs(m,i)-m_nnu) = nodevec(m,i);
 
   return dofval_p;
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_d(
-  const xt::xtensor<double,1> &dofval) const
+  const xt::xtensor<double,1>& dofval) const
 {
   assert(dofval.size() == m_ndof);
 
   Eigen::VectorXd dofval_d(m_nnd,1);
 
   #pragma omp parallel for
   for (size_t d = 0 ; d < m_nnd ; ++d)
     dofval_d(d) = dofval(m_iip(d));
 
   return dofval_d;
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline Eigen::VectorXd MatrixPartitionedTyings::asDofs_d(
-  const xt::xtensor<double,2> &nodevec) const
+  const xt::xtensor<double,2>& nodevec) const
 {
   assert(nodevec.shape() ==\
     std::decay_t<decltype(nodevec)>::shape_type({m_nnode, m_ndim}));
 
   Eigen::VectorXd dofval_d(m_nnd,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_dofs(m,i) >= m_nni)
         dofval_d(m_dofs(m,i)-m_nni) = nodevec(m,i);
 
   return dofval_d;
 }
 
 // -------------------------------------------------------------------------------------------------
 
 } // namespace ...
 
 // =================================================================================================
 
 #endif
diff --git a/include/GooseFEM/TyingsPeriodic.h b/include/GooseFEM/TyingsPeriodic.h
index 6eb26b6..07e5401 100644
--- a/include/GooseFEM/TyingsPeriodic.h
+++ b/include/GooseFEM/TyingsPeriodic.h
@@ -1,124 +1,124 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
 #ifndef GOOSEFEM_TYINGSPERIODIC_H
 #define GOOSEFEM_TYINGSPERIODIC_H
 
 // -------------------------------------------------------------------------------------------------
 
 #include "config.h"
 
 #include <Eigen/Eigen>
 #include <Eigen/Sparse>
 
 // =================================================================================================
 
 namespace GooseFEM {
 namespace Tyings {
 
 // -------------------------------------------------------------------------------------------------
 
 class Periodic
 {
 public:
 
   // Constructors
 
   Periodic() = default;
 
   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,2>& nodal_tyings); // (independent, dependent)
 
   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,2>& nodal_tyings, // (independent, dependent)
     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
 
   Eigen::SparseMatrix<double> Cdi() const;
   Eigen::SparseMatrix<double> Cdu() const;
   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;
+  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;
+  xt::xtensor<size_t,2> m_tyings;  // nodal ties: (independent, dependent)
   xt::xtensor<double,2> m_coor;
 
 };
 
 // -------------------------------------------------------------------------------------------------
 
 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 ...
 
 // -------------------------------------------------------------------------------------------------
 
 #include "TyingsPeriodic.hpp"
 
 // =================================================================================================
 
 #endif
diff --git a/include/GooseFEM/TyingsPeriodic.hpp b/include/GooseFEM/TyingsPeriodic.hpp
index adeccc7..8259e66 100644
--- a/include/GooseFEM/TyingsPeriodic.hpp
+++ b/include/GooseFEM/TyingsPeriodic.hpp
@@ -1,236 +1,285 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
 #ifndef GOOSEFEM_TYINGSPERIODIC_HPP
 #define GOOSEFEM_TYINGSPERIODIC_HPP
 
 // -------------------------------------------------------------------------------------------------
 
 #include "TyingsPeriodic.h"
 #include "Mesh.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];
 
   xt::xtensor<size_t,1> dependent = xt::view(m_tyings, xt::all(), 1);
-
-  xt::xtensor<size_t,1> iid = xt::flatten(xt::view(dofs, xt::keep(dependent), xt::all()));
+  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; }
+{
+  return m_dofs;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline xt::xtensor<size_t,2> Periodic::control() const
-{ return m_control; }
+{
+  return m_control;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline size_t Periodic::nnu() const
-{ return m_nnu; }
+{
+  return m_nnu;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline size_t Periodic::nnp() const
-{ return m_nnp; }
+{
+  return m_nnp;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline size_t Periodic::nni() const
-{ return m_nni; }
+{
+  return m_nni;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline size_t Periodic::nnd() const
-{ return m_nnd; }
+{
+  return m_nnd;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline xt::xtensor<size_t,1> Periodic::iiu() const
-{ return xt::arange<size_t>(m_nnu); }
+{
+  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; }
+{
+  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); }
+{
+  return xt::arange<size_t>(m_nni);
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline xt::xtensor<size_t,1> Periodic::iid() const
-{ return xt::arange<size_t>(m_nnd) + m_nni; }
+{
+  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)
-    {
+  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.));
 
       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)
-    {
+  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 )
+      if (m_dofs(ni,j) < m_nnu)
         data.push_back(Eigen::Triplet<double>(i*m_ndim+j, m_dofs(ni,j), +1.));
 
       for (size_t k = 0; k < m_ndim; ++k)
-        if ( m_control(j,k) < m_nnu )
+        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)
-    {
+  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 )
+      if (m_dofs(ni,j) >= m_nnu)
         data.push_back(Eigen::Triplet<double>(i*m_ndim+j, m_dofs(ni,j)-m_nnu, +1.));
 
       for (size_t k = 0; k < m_ndim; ++k)
-        if ( m_control(j,k) >= m_nnu )
+        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)[0] + 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; }
+{
+  return m_coor;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline xt::xtensor<size_t,2> Control::dofs() const
-{ return m_dofs; }
+{
+  return m_dofs;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline xt::xtensor<size_t,2> Control::controlDofs() const
-{ return m_control_dofs; }
+{
+  return m_control_dofs;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline xt::xtensor<size_t,1> Control::controlNodes() const
-{ return m_control_nodes; }
+{
+  return m_control_nodes;
+}
 
 // -------------------------------------------------------------------------------------------------
 
 }} // namespace ...
 
 // =================================================================================================
 
 #endif
diff --git a/include/GooseFEM/VectorPartitionedTyings.hpp b/include/GooseFEM/VectorPartitionedTyings.hpp
index 6b4c6b6..f2a436d 100644
--- a/include/GooseFEM/VectorPartitionedTyings.hpp
+++ b/include/GooseFEM/VectorPartitionedTyings.hpp
@@ -1,275 +1,329 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
 #ifndef GOOSEFEM_VECTORPARTITIONEDTYINGS_HPP
 #define GOOSEFEM_VECTORPARTITIONEDTYINGS_HPP
 
 // -------------------------------------------------------------------------------------------------
 
 #include "VectorPartitionedTyings.h"
 
 // =================================================================================================
 
 namespace GooseFEM {
 
 // -------------------------------------------------------------------------------------------------
 
 inline VectorPartitionedTyings::VectorPartitionedTyings(
-  const xt::xtensor<size_t,2> &conn,
-  const xt::xtensor<size_t,2> &dofs,
+  const xt::xtensor<size_t,2>& conn,
+  const xt::xtensor<size_t,2>& dofs,
   const Eigen::SparseMatrix<double>& Cdu,
   const Eigen::SparseMatrix<double>& Cdp,
   const Eigen::SparseMatrix<double>& Cdi) :
   m_conn(conn),
   m_dofs(dofs),
   m_Cdu(Cdu),
   m_Cdp(Cdp),
   m_Cdi(Cdi)
 {
   GOOSEFEM_ASSERT(Cdu.rows() == Cdp.rows());
   GOOSEFEM_ASSERT(Cdi.rows() == Cdp.rows());
 
   m_nnu   = static_cast<size_t>(m_Cdu.cols());
   m_nnp   = static_cast<size_t>(m_Cdp.cols());
   m_nnd   = static_cast<size_t>(m_Cdp.rows());
   m_nni   = m_nnu + m_nnp;
   m_ndof  = m_nni + m_nnd;
 
   m_iiu   = xt::arange<size_t>(m_nnu);
-  m_iip   = xt::arange<size_t>(m_nnp) + m_nnu;
-  m_iid   = xt::arange<size_t>(m_nnd) + m_nni;
+  m_iip   = xt::arange<size_t>(m_nnu, m_nnu + m_nnp);
+  m_iid   = xt::arange<size_t>(m_nni, m_nni + m_nnd);
 
   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_Cud   = m_Cdu.transpose();
   m_Cpd   = m_Cdp.transpose();
   m_Cid   = m_Cdi.transpose();
 
   GOOSEFEM_ASSERT(static_cast<size_t>(m_Cdi.cols()) == m_nni);
   GOOSEFEM_ASSERT(xt::amax(m_conn)[0] + 1 == m_nnode);
   GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim);
   GOOSEFEM_ASSERT(m_ndof == xt::amax(m_dofs)[0] + 1);
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline size_t VectorPartitionedTyings::nelem() const
-{ return m_nelem; }
+{
+  return m_nelem;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline size_t VectorPartitionedTyings::nne() const
-{ return m_nne; }
+{
+  return m_nne;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline size_t VectorPartitionedTyings::nnode() const
-{ return m_nnode; }
+{
+  return m_nnode;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline size_t VectorPartitionedTyings::ndim() const
-{ return m_ndim; }
+{
+  return m_ndim;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline size_t VectorPartitionedTyings::ndof() const
-{ return m_ndof; }
+{
+  return m_ndof;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline size_t VectorPartitionedTyings::nnu() const
-{ return m_nnu; }
+{
+  return m_nnu;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline size_t VectorPartitionedTyings::nnp() const
-{ return m_nnp; }
+{
+  return m_nnp;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline size_t VectorPartitionedTyings::nni() const
-{ return m_nni; }
+{
+  return m_nni;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline size_t VectorPartitionedTyings::nnd() const
-{ return m_nnd; }
+{
+  return m_nnd;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline xt::xtensor<size_t,2> VectorPartitionedTyings::dofs() const
-{ return m_dofs; }
+{
+  return m_dofs;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline xt::xtensor<size_t,1> VectorPartitionedTyings::iiu() const
-{ return m_iiu; }
+{
+  return m_iiu;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline xt::xtensor<size_t,1> VectorPartitionedTyings::iip() const
-{ return m_iip; }
+{
+  return m_iip;
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline xt::xtensor<size_t,1> VectorPartitionedTyings::iii() const
-{ return xt::arange<size_t>(m_nni); }
+{
+  return xt::arange<size_t>(m_nni);
+}
+
+// -------------------------------------------------------------------------------------------------
 
 inline xt::xtensor<size_t,1> VectorPartitionedTyings::iid() const
-{ return m_iid; }
+{
+  return m_iid;
+}
 
 // -------------------------------------------------------------------------------------------------
 
 inline void VectorPartitionedTyings::copy_p(
   const xt::xtensor<double,1>& dofval_src,
         xt::xtensor<double,1>& dofval_dest) const
 {
   GOOSEFEM_ASSERT(dofval_src.size() == m_ndof || dofval_src.size() == m_nni);
   GOOSEFEM_ASSERT(dofval_dest.size() == m_ndof || dofval_dest.size() == m_nni);
 
   #pragma omp parallel for
   for (size_t i = m_nnu; i < m_nni; ++i)
     dofval_dest(i) = dofval_src(i);
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline void VectorPartitionedTyings::asDofs_i(
   const xt::xtensor<double,2>& nodevec,
         xt::xtensor<double,1>& dofval_i,
         bool apply_tyings) const
 {
   GOOSEFEM_ASSERT(nodevec.shape() ==\
     std::decay_t<decltype(nodevec)>::shape_type({m_nnode, m_ndim}));
   GOOSEFEM_ASSERT(dofval_i.size() == m_nni);
 
   #pragma omp parallel for
   for (size_t m = 0 ; m < m_nnode ; ++m)
     for (size_t i = 0 ; i < m_ndim ; ++i)
       if (m_dofs(m,i) < m_nni)
         dofval_i(m_dofs(m,i)) = nodevec(m,i);
 
   if (!apply_tyings)
     return;
 
   Eigen::VectorXd Dofval_d = this->Eigen_asDofs_d(nodevec);
   Eigen::VectorXd Dofval_i = m_Cid * Dofval_d;
 
   #pragma omp parallel for
   for (size_t i = 0 ; i < m_nni ; ++i)
     dofval_i(i) += Dofval_i(i);
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline void VectorPartitionedTyings::asNode(
   const xt::xtensor<double,1>& dofval,
         xt::xtensor<double,2>& nodevec) const
 {
   GOOSEFEM_ASSERT(dofval.size() == m_ndof);
   GOOSEFEM_ASSERT(nodevec.shape() ==\
     std::decay_t<decltype(nodevec)>::shape_type({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)
       nodevec(m,i) = dofval(m_dofs(m,i));
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline void VectorPartitionedTyings::asElement(
   const xt::xtensor<double,2>& nodevec,
         xt::xtensor<double,3>& elemvec) const
 {
   GOOSEFEM_ASSERT(nodevec.shape() ==\
     std::decay_t<decltype(nodevec)>::shape_type({m_nnode, m_ndim}));
   GOOSEFEM_ASSERT(elemvec.shape() ==\
     std::decay_t<decltype(elemvec)>::shape_type({m_nelem, m_nne, m_ndim}));
 
   #pragma omp parallel for
   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)
         elemvec(e,m,i) = nodevec(m_conn(e,m),i);
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline void VectorPartitionedTyings::assembleDofs(
   const xt::xtensor<double,3>& elemvec,
         xt::xtensor<double,1>& dofval) const
 {
   GOOSEFEM_ASSERT(elemvec.shape() ==\
     std::decay_t<decltype(elemvec)>::shape_type({m_nelem, m_nne, m_ndim}));
   GOOSEFEM_ASSERT(dofval.size() == m_ndof);
 
   dofval.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)
           dofval(m_dofs(m_conn(e,m),i)) += elemvec(e,m,i);
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline void VectorPartitionedTyings::assembleNode(
   const xt::xtensor<double,3>& elemvec,
         xt::xtensor<double,2>& nodevec) const
 {
   GOOSEFEM_ASSERT(elemvec.shape() ==\
     std::decay_t<decltype(elemvec)>::shape_type({m_nelem, m_nne, m_ndim}));
   GOOSEFEM_ASSERT(nodevec.shape() ==\
     std::decay_t<decltype(nodevec)>::shape_type({m_nnode, m_ndim}));
 
   xt::xtensor<double,1> dofval = this->AssembleDofs(elemvec);
 
   this->asNode(dofval, nodevec);
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline xt::xtensor<double,3> VectorPartitionedTyings::AsElement(
   const xt::xtensor<double,2>& nodevec) const
 {
   xt::xtensor<double,3> elemvec = xt::empty<double>({m_nelem, m_nne, m_ndim});
 
   this->asElement(nodevec, elemvec);
 
   return elemvec;
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline xt::xtensor<double,1> VectorPartitionedTyings::AssembleDofs(
   const xt::xtensor<double,3>& elemvec) const
 {
   xt::xtensor<double,1> dofval = xt::empty<double>({m_ndof});
 
   this->assembleDofs(elemvec, dofval);
 
   return dofval;
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline xt::xtensor<double,2> VectorPartitionedTyings::AssembleNode(
   const xt::xtensor<double,3>& elemvec) const
 {
   xt::xtensor<double,2> nodevec = xt::empty<double>({m_nnode, m_ndim});
 
   this->assembleNode(elemvec, nodevec);
 
   return nodevec;
 }
 
 // -------------------------------------------------------------------------------------------------
 
 inline Eigen::VectorXd VectorPartitionedTyings::Eigen_asDofs_d(
   const xt::xtensor<double,2>& nodevec) const
 {
   GOOSEFEM_ASSERT(nodevec.shape() ==\
     std::decay_t<decltype(nodevec)>::shape_type({m_nnode, m_ndim}));
 
   Eigen::VectorXd dofval_d(m_nnd,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_dofs(m,i) >= m_nni)
         dofval_d(m_dofs(m,i)-m_nni) = nodevec(m,i);
 
   return dofval_d;
 }
 
 // -------------------------------------------------------------------------------------------------
 
 } // namespace ...
 
 // =================================================================================================
 
 #endif