diff --git a/src/GooseFEM/DynamicsDiagonalPeriodic.h b/src/GooseFEM/DynamicsDiagonalPeriodic.h
index 720f8bb..89bed53 100644
--- a/src/GooseFEM/DynamicsDiagonalPeriodic.h
+++ b/src/GooseFEM/DynamicsDiagonalPeriodic.h
@@ -1,351 +1,353 @@
 /* ========================================== DESCRIPTION ==========================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
 #ifndef GOOSEFEM_DYNAMICS_DIAGONAL_PERIODIC_H
 #define GOOSEFEM_DYNAMICS_DIAGONAL_PERIODIC_H
 
 #include "Macros.h"
 #include <cppmat/cppmat.h>
 
 // -------------------------------------------------------------------------------------------------
 
 namespace GooseFEM {
 namespace Dynamics {
 namespace Diagonal {
 
 // =========================================== OVERVIEW ============================================
 
 template <class Element>
 class Periodic
 {
 public:
 
   // variables
   // ---------
 
   // element/quadrature/material definition (should match mesh)
   std::shared_ptr<Element> elem;
 
   // mesh : nodal quantities and connectivity
   MatS   dofs;      // DOF-numbers       of each node    (column of 'vectors')
   MatS   conn;      // node numbers      of each element (column of elements)
   MatD   x0;        // initial positions of each node    (column of vectors)
   MatD   u;         // displacements     of each node    (column of vectors)
   MatD   v;         // velocities        of each node    (column of vectors)
   MatD   a;         // accelerations     of each node    (column of vectors)
   size_t nnode;     // number of nodes
   size_t nelem;     // number of elements
   size_t nne;       // number of nodes-per-element
   size_t ndim;      // number of dimensions
   size_t ndof;      // number of DOFs (after eliminating periodic dependencies)
 
   // linear system
   ColD   Minv;      // inverse of mass matrix (constructed diagonal -> inverse == diagonal)
   ColD   Fu;        // force that depends on displacement (column)
   ColD   Fv;        // force that depends on velocity (and displacement) (column)
   ColD   V;         // velocity (column)
   ColD   V_n;       // velocity (column, last increment)
   ColD   A;         // acceleration (column)
   ColD   A_n;       // acceleration (column, last increment)
 
   // time integration
   double dt;        // time step
   double t = 0.0;   // current time
   double alpha;     // non-Galilean damping coefficient
 
   // constructor
   // -----------
 
+  Periodic(){};
+
   Periodic(std::shared_ptr<Element> elem, const MatD &x0, const MatS &conn, const MatS &dofs,
     double dt, double alpha=0.0 );
 
   // functions
   // ---------
 
   void velocityVerlet();  // one time step of time integrator
   void Verlet();          // one time step of time integrator (Fv == 0)
   void computeMinv();     // element loop: inverse mass matrix            <- "elem->computeM"
   void computeFu();       // element loop: displacement dependent forces  <- "elem->computeFu"
   void computeFv();       // element loop: velocity dependent forces      <- "elem->computeFv"
   void post();            // element loop: post-process                   <- "elem->post"
 
 };
 
 // ========================================== SOURCE CODE ==========================================
 
 template<class Element>
 Periodic<Element>::Periodic(
   std::shared_ptr<Element> _elem, const MatD &_x0, const MatS &_conn, const MatS &_dofs,
   double _dt, double _alpha
 )
 {
   // problem specific settings
   elem  = _elem;
 
   // mesh definition
   x0    = _x0;
   conn  = _conn;
   dofs  = _dofs;
 
   // time integration
   dt    = _dt;
   alpha = _alpha;
 
   // extract mesh size from definition
   nnode = static_cast<size_t>(x0  .rows());
   ndim  = static_cast<size_t>(x0  .cols());
   nelem = static_cast<size_t>(conn.rows());
   nne   = static_cast<size_t>(conn.cols());
   ndof  = dofs.maxCoeff()+1;
 
   // basic check (mostly the user is 'trusted')
   assert( dofs.size() == nnode * ndim );
   assert( ndof        <  nnode * ndim );
 
   // nodal quantities
   u.conservativeResize(nnode,ndim); u.setZero();
   v.conservativeResize(nnode,ndim); v.setZero();
   a.conservativeResize(nnode,ndim); a.setZero();
 
   // linear system (DOFs)
   Minv.conservativeResize(ndof); Minv.setZero();
   Fu  .conservativeResize(ndof); Fu  .setZero();
   Fv  .conservativeResize(ndof); Fv  .setZero();
   V   .conservativeResize(ndof); V   .setZero();
   V_n .conservativeResize(ndof); V_n .setZero();
   A   .conservativeResize(ndof); A   .setZero();
   A_n .conservativeResize(ndof); A_n .setZero();
 
   // compute inverse mass matrix : assumed constant in time
   computeMinv();
 }
 
 // =================================================================================================
 
 template<class Element>
 void Periodic<Element>::velocityVerlet()
 {
   // (1) new nodal positions (displacements)
   // - apply update (nodal) : x_{n+1} = x_n + dt * v_n + .5 * dt^2 * a_n"
   u += dt * v + ( .5 * std::pow(dt,2.) ) * a;
   // - compute forces that are dependent on the displacement, but not on the velocity
   computeFu();
 
   // (2a) estimate nodal velocities
   // - update velocities (DOFs)
   V.noalias() = V_n + dt * A;
   // - convert to nodal velocity (periodicity implies that several nodes depend on the same DOF)
   for ( size_t i=0; i<nnode*ndim; ++i ) v(i) = V(dofs(i));
   // - compute forces that are dependent on the velocity
   computeFv();
   // - solve for accelerations (DOFs)
   A.noalias() = Minv.cwiseProduct( - Fu - Fv - alpha*V );
   // - update velocities (DOFs)
   V.noalias() = V_n + ( .5 * dt ) * ( A_n + A );
   // - convert to nodal velocity (periodicity implies that several nodes depend on the same DOF)
   for ( size_t i=0; i<nnode*ndim; ++i ) v(i) = V(dofs(i));
   // - compute forces that are dependent on the velocity
   computeFv();
 
   // (2b) new nodal velocities
   // - solve for accelerations (DOFs)
   A.noalias() = Minv.cwiseProduct( - Fu - Fv - alpha*V );
   // - update velocities (DOFs)
   V.noalias() = V_n + ( .5 * dt ) * ( A_n + A );
   // - convert to nodal velocity (periodicity implies that several nodes depend on the same DOF)
   for ( size_t i=0; i<nnode*ndim; ++i ) v(i) = V(dofs(i));
   // - compute forces that are dependent on the velocity
   computeFv();
 
   // (3) new nodal accelerations
   // - solve for accelerations (DOFs)
   A.noalias() = Minv.cwiseProduct( - Fu - Fv - alpha*V );
   // - convert to nodal acceleration (periodicity implies that several nodes depend on the same DOF)
   for ( size_t i=0; i<nnode*ndim; ++i ) a(i) = A(dofs(i));
 
   // store history
   A_n = A;  // accelerations (DOFs)
   V_n = V;  // nodal velocities
   t  += dt; // current time instance
 
   // N.B. at this point:
   // "a" == "A" == "A_n"  ->  new nodal accelerations, their DOF equivalents, and a 'back-up'
   // "v" ==        "V_n"  ->  new nodal velocities,                           and a 'back-up'
   // "u"                  ->  new nodal displacements
   // The forces "Fu" and "Fv" correspond to this state of the system
 }
 
 // =================================================================================================
 
 template<class Element>
 void Periodic<Element>::Verlet()
 {
   // (1) new nodal positions (displacements)
   // - apply update (nodal) : x_{n+1} = x_n + dt * v_n + .5 * dt^2 * a_n"
   u += dt * v + ( .5 * std::pow(dt,2.) ) * a;
   // - compute forces that are dependent on the displacement, but not on the velocity
   computeFu();
   // - solve for accelerations (DOFs)
   A.noalias() = Minv.cwiseProduct( - Fu );
   // - convert to nodal acceleration (periodicity implies that several nodes depend on the same DOF)
   for ( size_t i=0; i<nnode*ndim; ++i ) a(i) = A(dofs(i));
 
   // (2) propagate velocities
   // - update velocities (DOFs)
   V.noalias() = V_n + ( .5 * dt ) * ( A_n + A );
   // - convert to nodal velocity (periodicity implies that several nodes depend on the same DOF)
   for ( size_t i=0; i<nnode*ndim; ++i ) v(i) = V(dofs(i));
 
   // store history
   A_n = A;  // accelerations (DOFs)
   V_n = V;  // nodal velocities
   t  += dt; // current time instance
 
   // N.B. at this point:
   // "a" == "A" == "A_n"  ->  new nodal accelerations, their DOF equivalents, and a 'back-up'
   // "v" ==        "V_n"  ->  new nodal velocities,                           and a 'back-up'
   // "u"                  ->  new nodal displacements
   // The forces "Fu" correspond to this state of the system
 }
 
 // =================================================================================================
 
 template<class Element>
 void Periodic<Element>::computeMinv()
 {
   // array with DOF numbers of each node
   cppmat::matrix2<size_t> dofe(nne,ndim);
 
   // zero-initialize mass matrix (only the inverse is stored)
   ColD M(ndof); M.setZero();
 
   // loop over all elements
   for ( size_t e = 0 ; e < nelem ; ++e )
   {
     // - nodal positions, displacement, and DOF numbers for element "e"
     for ( size_t m = 0 ; m < nne ; ++m ) {
       for ( size_t i = 0 ; i < ndim ; ++i ) {
         elem->xe(m,i) = x0  (conn(e,m),i);
         elem->ue(m,i) = u   (conn(e,m),i);
         dofe    (m,i) = dofs(conn(e,m),i);
       }
     }
 
     // - compute element mass matrix using problem specific "Element" class
     elem->computeM(e);
 
     // - assemble element mass matrix : take only the diagonal
     for ( size_t i = 0 ; i < nne*ndim ; ++i ) M(dofe[i]) += elem->M(i,i);
 
     // - check that the user provided a diagonal mass matrix
     #ifndef NDEBUG
-      for ( size_t i = 0 ; i < nne*ndim ; ++i ) {
-        for ( size_t j = 0 ; j < nne*ndim ; ++j ) {
-          if ( i != j ) assert( ! elem->M(i,j) );
-          else          assert(   elem->M(i,i) );
-        }
+    for ( size_t i = 0 ; i < nne*ndim ; ++i ) {
+      for ( size_t j = 0 ; j < nne*ndim ; ++j ) {
+        if ( i != j ) assert( ! elem->M(i,j) );
+        else          assert(   elem->M(i,i) );
       }
+    }
     #endif
   }
 
   // compute inverse of the mass matrix
   Minv = M.cwiseInverse();
 }
 
 // =================================================================================================
 
 template<class Element>
 void Periodic<Element>::computeFu()
 {
   // array with DOF numbers of each node
   cppmat::matrix2<size_t> dofe(nne,ndim);
 
   // zero-initialize displacement dependent force
   Fu.setZero();
 
   // loop over all elements
   for ( size_t e = 0 ; e < nelem ; ++e )
   {
     // - nodal positions, displacement, and DOF numbers for element "e"
     for ( size_t m = 0 ; m < nne ; ++m ) {
       for ( size_t i = 0 ; i < ndim ; ++i ) {
         elem->xe(m,i) = x0  (conn(e,m),i);
         elem->ue(m,i) = u   (conn(e,m),i);
         dofe    (m,i) = dofs(conn(e,m),i);
       }
     }
 
     // - compute element force using problem specific "Element" class
     elem->computeFu(e);
 
     // - assemble force to global system
     for ( size_t i = 0 ; i < nne*ndim ; ++i ) Fu(dofe[i]) += elem->fu(i);
   }
 }
 
 // =================================================================================================
 
 template<class Element>
 void Periodic<Element>::computeFv()
 {
   // array with DOF numbers of each node
   cppmat::matrix2<size_t> dofe(nne,ndim);
 
   // zero-initialize velocity dependent force
   Fv.setZero();
 
   // loop over all elements
   for ( size_t e = 0 ; e < nelem ; ++e )
   {
     // - nodal positions, displacement, velocity, and DOF numbers for element "e"
     for ( size_t m = 0 ; m < nne ; ++m ) {
       for ( size_t i = 0 ; i < ndim ; ++i ) {
         elem->xe(m,i) = x0  (conn(e,m),i);
         elem->ue(m,i) = u   (conn(e,m),i);
         elem->ve(m,i) = v   (conn(e,m),i);
         dofe    (m,i) = dofs(conn(e,m),i);
       }
     }
 
     // - compute element force using problem specific "Element" class
     elem->computeFv(e);
 
     // - assemble force to global system
     for ( size_t i = 0 ; i < nne*ndim ; ++i ) Fv(dofe[i]) += elem->fv(i);
   }
 }
 
 // =================================================================================================
 
 template<class Element>
 void Periodic<Element>::post()
 {
   // array with DOF numbers of each node
   cppmat::matrix2<size_t> dofe(nne,ndim);
 
   // loop over all elements
   for ( size_t e = 0 ; e < nelem ; ++e )
   {
     // - nodal positions, displacement, velocity, and DOF numbers for element "e"
     for ( size_t m = 0 ; m < nne ; ++m ) {
       for ( size_t i = 0 ; i < ndim ; ++i ) {
         elem->xe(m,i) = x0  (conn(e,m),i);
         elem->ue(m,i) = u   (conn(e,m),i);
         elem->ve(m,i) = v   (conn(e,m),i);
         dofe    (m,i) = dofs(conn(e,m),i);
       }
     }
 
     // - run post-process function of the problem specific "Element" class (output is handled there)
     elem->post(e);
   }
 }
 
 // =================================================================================================
 
 }}} // namespace GooseFEM::Dynamics::Diagonal
 
 #endif
diff --git a/src/GooseFEM/DynamicsDiagonalSemiPeriodic.h b/src/GooseFEM/DynamicsDiagonalSemiPeriodic.h
index 0708763..4fefbb8 100644
--- a/src/GooseFEM/DynamicsDiagonalSemiPeriodic.h
+++ b/src/GooseFEM/DynamicsDiagonalSemiPeriodic.h
@@ -1,377 +1,383 @@
 /* ========================================== DESCRIPTION ==========================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
 #ifndef GOOSEFEM_DYNAMICS_DIAGONAL_SEMIPERIODIC_H
 #define GOOSEFEM_DYNAMICS_DIAGONAL_SEMIPERIODIC_H
 
 #include "Macros.h"
 #include <cppmat/cppmat.h>
 
 // -------------------------------------------------------------------------------------------------
 
 namespace GooseFEM {
 namespace Dynamics {
 namespace Diagonal {
 
 // =========================================== OVERVIEW ============================================
 
 template <class Element>
 class SemiPeriodic
 {
 public:
 
   // variables
   // ---------
 
   // element/quadrature/material definition (should match mesh)
   std::shared_ptr<Element> elem;
 
   // mesh : nodal quantities and connectivity
   MatS   dofs;      // DOF-numbers       of each node    (column of 'vectors')
   MatS   conn;      // node numbers      of each element (column of elements)
   MatD   x0;        // initial positions of each node    (column of vectors)
   MatD   u;         // displacements     of each node    (column of vectors)
   MatD   v;         // velocities        of each node    (column of vectors)
   MatD   a;         // accelerations     of each node    (column of vectors)
   size_t nnode;     // number of nodes
   size_t nelem;     // number of elements
   size_t nne;       // number of nodes-per-element
   size_t ndim;      // number of dimensions
   size_t ndof;      // number of DOFs (after eliminating periodic dependencies)
 
   // fixed DOFs
   size_t nfixed;    // number of fixed DOFs
   ColS   fixedDofs; // DOF-number that are prescribed
   ColD   fixedV;    // prescribed velocity     for the fixed DOFs
   ColD   fixedA;    // prescribed acceleration for the fixed DOFs
 
   // linear system
   ColD   Minv;      // inverse of mass matrix (constructed diagonal -> inverse == diagonal)
   ColD   Fu;        // force that depends on displacement (column)
   ColD   Fv;        // force that depends on velocity (and displacement) (column)
   ColD   V;         // velocity (column)
   ColD   V_n;       // velocity (column, last increment)
   ColD   A;         // acceleration (column)
   ColD   A_n;       // acceleration (column, last increment)
 
   // time integration
   double dt;        // time step
   double t = 0.0;   // current time
   double alpha;     // non-Galilean damping coefficient
 
   // constructor
   // -----------
 
+  SemiPeriodic(){};
+
   SemiPeriodic(std::shared_ptr<Element> elem, const MatD &x0, const MatS &conn, const MatS &dofs,
     const ColS &fixedDofs, double dt, double alpha=0.0 );
 
   // functions
   // ---------
 
   void velocityVerlet();  // one time step of time integrator
   void Verlet();          // one time step of time integrator (Fv == 0)
   void computeMinv();     // element loop: inverse mass matrix            <- "elem->computeM"
   void computeFu();       // element loop: displacement dependent forces  <- "elem->computeFu"
   void computeFv();       // element loop: velocity dependent forces      <- "elem->computeFv"
   void post();            // element loop: post-process                   <- "elem->post"
 
 };
 
 // ========================================== SOURCE CODE ==========================================
 
 template<class Element>
 SemiPeriodic<Element>::SemiPeriodic(
   std::shared_ptr<Element> _elem, const MatD &_x0, const MatS &_conn, const MatS &_dofs,
   const ColS &_fixedDofs, double _dt, double _alpha
 )
 {
   // problem specific settings
   elem  = _elem;
 
   // mesh definition
   x0    = _x0;
   conn  = _conn;
   dofs  = _dofs;
 
   // time integration
   dt    = _dt;
   alpha = _alpha;
 
   // fixed DOFs
   fixedDofs = _fixedDofs;
   nfixed    = static_cast<size_t>(fixedDofs.size());
 
   // extract mesh size from definition
   nnode = static_cast<size_t>(x0  .rows());
   ndim  = static_cast<size_t>(x0  .cols());
   nelem = static_cast<size_t>(conn.rows());
   nne   = static_cast<size_t>(conn.cols());
   ndof  = dofs.maxCoeff()+1;
 
   // basic check (mostly the user is 'trusted')
   assert( dofs.size() == nnode * ndim );
   assert( ndof        <  nnode * ndim );
 
+  #ifndef NDEBUG
+  for ( size_t i = 0 ; i < nfixed ; ++i ) assert( fixedDofs(i) < ndof );
+  #endif
+
   // nodal quantities
   u.conservativeResize(nnode,ndim); u.setZero();
   v.conservativeResize(nnode,ndim); v.setZero();
   a.conservativeResize(nnode,ndim); a.setZero();
 
   // linear system (DOFs)
   Minv.conservativeResize(ndof); Minv.setZero();
   Fu  .conservativeResize(ndof); Fu  .setZero();
   Fv  .conservativeResize(ndof); Fv  .setZero();
   V   .conservativeResize(ndof); V   .setZero();
   V_n .conservativeResize(ndof); V_n .setZero();
   A   .conservativeResize(ndof); A   .setZero();
   A_n .conservativeResize(ndof); A_n .setZero();
 
   // fixed DOFs : default zero velocity and acceleration
   fixedV.conservativeResize(nfixed); fixedV.setZero();
   fixedA.conservativeResize(nfixed); fixedA.setZero();
 
   // compute inverse mass matrix : assumed constant in time
   computeMinv();
 }
 
 // =================================================================================================
 
 template<class Element>
 void SemiPeriodic<Element>::velocityVerlet()
 {
   // (1) new nodal positions (displacements)
   // - apply update (nodal) : x_{n+1} = x_n + dt * v_n + .5 * dt^2 * a_n"
   u += dt * v + ( .5 * std::pow(dt,2.) ) * a;
   // - compute forces that are dependent on the displacement, but not on the velocity
   computeFu();
 
   // (2a) estimate nodal velocities
   // - update velocities (DOFs)
   V.noalias() = V_n + dt * A;
+  // - set fixed velocity
+  for ( size_t i=0; i<nfixed; ++i ) V(fixedDofs(i)) = fixedV(i);
   // - convert to nodal velocity (periodicity implies that several nodes depend on the same DOF)
   for ( size_t i=0; i<nnode*ndim; ++i ) v(i) = V(dofs(i));
-  // - set fixed velocity
-  for ( size_t i=0; i<nfixed; ++i ) v(fixedDofs(i)) = fixedV(i);
   // - compute forces that are dependent on the velocity
   computeFv();
   // - solve for accelerations (DOFs)
   A.noalias() = Minv.cwiseProduct( - Fu - Fv - alpha*V );
   // - update velocities (DOFs)
   V.noalias() = V_n + ( .5 * dt ) * ( A_n + A );
+  // - set fixed velocity
+  for ( size_t i=0; i<nfixed; ++i ) V(fixedDofs(i)) = fixedV(i);
   // - convert to nodal velocity (periodicity implies that several nodes depend on the same DOF)
   for ( size_t i=0; i<nnode*ndim; ++i ) v(i) = V(dofs(i));
-  // - set fixed velocity
-  for ( size_t i=0; i<nfixed; ++i ) v(fixedDofs(i)) = fixedV(i);
   // - compute forces that are dependent on the velocity
   computeFv();
 
   // (2b) new nodal velocities
   // - solve for accelerations (DOFs)
   A.noalias() = Minv.cwiseProduct( - Fu - Fv - alpha*V );
   // - update velocities (DOFs)
   V.noalias() = V_n + ( .5 * dt ) * ( A_n + A );
+  // - set fixed velocity
+  for ( size_t i=0; i<nfixed; ++i ) V(fixedDofs(i)) = fixedV(i);
   // - convert to nodal velocity (periodicity implies that several nodes depend on the same DOF)
   for ( size_t i=0; i<nnode*ndim; ++i ) v(i) = V(dofs(i));
-  // - set fixed velocity
-  for ( size_t i=0; i<nfixed; ++i ) v(fixedDofs(i)) = fixedV(i);
   // - compute forces that are dependent on the velocity
   computeFv();
 
   // (3) new nodal accelerations
   // - solve for accelerations (DOFs)
   A.noalias() = Minv.cwiseProduct( - Fu - Fv - alpha*V );
+  // - set fixed acceleration
+  for ( size_t i=0; i<nfixed; ++i ) A(fixedDofs(i)) = fixedA(i);
   // - convert to nodal acceleration (periodicity implies that several nodes depend on the same DOF)
   for ( size_t i=0; i<nnode*ndim; ++i ) a(i) = A(dofs(i));
-  // - set fixed acceleration
-  for ( size_t i=0; i<nfixed; ++i ) a(fixedDofs(i)) = fixedA(i);
 
   // store history
   A_n = A;  // accelerations (DOFs)
   V_n = V;  // nodal velocities
   t  += dt; // current time instance
 
   // N.B. at this point:
   // "a" == "A" == "A_n"  ->  new nodal accelerations, their DOF equivalents, and a 'back-up'
   // "v" ==        "V_n"  ->  new nodal velocities,                           and a 'back-up'
   // "u"                  ->  new nodal displacements
   // The forces "Fu" and "Fv" correspond to this state of the system
 }
 
 // =================================================================================================
 
 template<class Element>
 void SemiPeriodic<Element>::Verlet()
 {
   // (1) new nodal positions (displacements)
   // - apply update (nodal) : x_{n+1} = x_n + dt * v_n + .5 * dt^2 * a_n"
   u += dt * v + ( .5 * std::pow(dt,2.) ) * a;
   // - compute forces that are dependent on the displacement, but not on the velocity
   computeFu();
   // - solve for accelerations (DOFs)
   A.noalias() = Minv.cwiseProduct( - Fu );
+  // - set fixed acceleration
+  for ( size_t i=0; i<nfixed; ++i ) A(fixedDofs(i)) = fixedA(i);
   // - convert to nodal acceleration (periodicity implies that several nodes depend on the same DOF)
   for ( size_t i=0; i<nnode*ndim; ++i ) a(i) = A(dofs(i));
-  // - set fixed acceleration
-  for ( size_t i=0; i<nfixed; ++i ) a(fixedDofs(i)) = fixedA(i);
 
   // (2) propagate velocities
   // - update velocities (DOFs)
   V.noalias() = V_n + ( .5 * dt ) * ( A_n + A );
+  // - set fixed velocity
+  for ( size_t i=0; i<nfixed; ++i ) V(fixedDofs(i)) = fixedV(i);
   // - convert to nodal velocity (periodicity implies that several nodes depend on the same DOF)
   for ( size_t i=0; i<nnode*ndim; ++i ) v(i) = V(dofs(i));
-  // - set fixed velocity
-  for ( size_t i=0; i<nfixed; ++i ) v(fixedDofs(i)) = fixedV(i);
 
   // store history
   A_n = A;  // accelerations (DOFs)
   V_n = V;  // nodal velocities
   t  += dt; // current time instance
 
   // N.B. at this point:
   // "a" == "A" == "A_n"  ->  new nodal accelerations, their DOF equivalents, and a 'back-up'
   // "v" ==        "V_n"  ->  new nodal velocities,                           and a 'back-up'
   // "u"                  ->  new nodal displacements
   // The forces "Fu" correspond to this state of the system
 }
 
 // =================================================================================================
 
 template<class Element>
 void SemiPeriodic<Element>::computeMinv()
 {
   // array with DOF numbers of each node
   cppmat::matrix2<size_t> dofe(nne,ndim);
 
   // zero-initialize mass matrix (only the inverse is stored)
   ColD M(ndof); M.setZero();
 
   // loop over all elements
   for ( size_t e = 0 ; e < nelem ; ++e )
   {
     // - nodal positions, displacement, and DOF numbers for element "e"
     for ( size_t m = 0 ; m < nne ; ++m ) {
       for ( size_t i = 0 ; i < ndim ; ++i ) {
         elem->xe(m,i) = x0  (conn(e,m),i);
         elem->ue(m,i) = u   (conn(e,m),i);
         dofe    (m,i) = dofs(conn(e,m),i);
       }
     }
 
     // - compute element mass matrix using problem specific "Element" class
     elem->computeM(e);
 
     // - assemble element mass matrix : take only the diagonal
     for ( size_t i = 0 ; i < nne*ndim ; ++i ) M(dofe[i]) += elem->M(i,i);
 
     // - check that the user provided a diagonal mass matrix
     #ifndef NDEBUG
-      for ( size_t i = 0 ; i < nne*ndim ; ++i ) {
-        for ( size_t j = 0 ; j < nne*ndim ; ++j ) {
-          if ( i != j ) assert( ! elem->M(i,j) );
-          else          assert(   elem->M(i,i) );
-        }
+    for ( size_t i = 0 ; i < nne*ndim ; ++i ) {
+      for ( size_t j = 0 ; j < nne*ndim ; ++j ) {
+        if ( i != j ) assert( ! elem->M(i,j) );
+        else          assert(   elem->M(i,i) );
       }
+    }
     #endif
   }
 
   // compute inverse of the mass matrix
   Minv = M.cwiseInverse();
 }
 
 // =================================================================================================
 
 template<class Element>
 void SemiPeriodic<Element>::computeFu()
 {
   // array with DOF numbers of each node
   cppmat::matrix2<size_t> dofe(nne,ndim);
 
   // zero-initialize displacement dependent force
   Fu.setZero();
 
   // loop over all elements
   for ( size_t e = 0 ; e < nelem ; ++e )
   {
     // - nodal positions, displacement, and DOF numbers for element "e"
     for ( size_t m = 0 ; m < nne ; ++m ) {
       for ( size_t i = 0 ; i < ndim ; ++i ) {
         elem->xe(m,i) = x0  (conn(e,m),i);
         elem->ue(m,i) = u   (conn(e,m),i);
         dofe    (m,i) = dofs(conn(e,m),i);
       }
     }
 
     // - compute element force using problem specific "Element" class
     elem->computeFu(e);
 
     // - assemble force to global system
     for ( size_t i = 0 ; i < nne*ndim ; ++i ) Fu(dofe[i]) += elem->fu(i);
   }
 }
 
 // =================================================================================================
 
 template<class Element>
 void SemiPeriodic<Element>::computeFv()
 {
   // array with DOF numbers of each node
   cppmat::matrix2<size_t> dofe(nne,ndim);
 
   // zero-initialize velocity dependent force
   Fv.setZero();
 
   // loop over all elements
   for ( size_t e = 0 ; e < nelem ; ++e )
   {
     // - nodal positions, displacement, velocity, and DOF numbers for element "e"
     for ( size_t m = 0 ; m < nne ; ++m ) {
       for ( size_t i = 0 ; i < ndim ; ++i ) {
         elem->xe(m,i) = x0  (conn(e,m),i);
         elem->ue(m,i) = u   (conn(e,m),i);
         elem->ve(m,i) = v   (conn(e,m),i);
         dofe    (m,i) = dofs(conn(e,m),i);
       }
     }
 
     // - compute element force using problem specific "Element" class
     elem->computeFv(e);
 
     // - assemble force to global system
     for ( size_t i = 0 ; i < nne*ndim ; ++i ) Fv(dofe[i]) += elem->fv(i);
   }
 }
 
 // =================================================================================================
 
 template<class Element>
 void SemiPeriodic<Element>::post()
 {
   // array with DOF numbers of each node
   cppmat::matrix2<size_t> dofe(nne,ndim);
 
   // loop over all elements
   for ( size_t e = 0 ; e < nelem ; ++e )
   {
     // - nodal positions, displacement, velocity, and DOF numbers for element "e"
     for ( size_t m = 0 ; m < nne ; ++m ) {
       for ( size_t i = 0 ; i < ndim ; ++i ) {
         elem->xe(m,i) = x0  (conn(e,m),i);
         elem->ue(m,i) = u   (conn(e,m),i);
         elem->ve(m,i) = v   (conn(e,m),i);
         dofe    (m,i) = dofs(conn(e,m),i);
       }
     }
 
     // - run post-process function of the problem specific "Element" class (output is handled there)
     elem->post(e);
   }
 }
 
 // =================================================================================================
 
 }}} // namespace GooseFEM::Dynamics::Diagonal
 
 #endif
diff --git a/src/GooseFEM/Macros.h b/src/GooseFEM/Macros.h
index 0f9abbe..598e4cf 100644
--- a/src/GooseFEM/Macros.h
+++ b/src/GooseFEM/Macros.h
@@ -1,57 +1,57 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
 #ifndef GOOSEFEM_MACROS_H
 #define GOOSEFEM_MACROS_H
 
 #define _USE_MATH_DEFINES // to use "M_PI" from "math.h"
 
 #include <assert.h>
 #include <cstdlib>
 #include <vector>
 #include <string>
 #include <memory>
 #include <iostream>
 #include <iomanip>
 #include <math.h>
 
 #include <Eigen/Dense>
 
 #include <cppmat/tensor.h>
 #include <cppmat/tensor2.h>
 #include <cppmat/tensor3.h>
 
 // =================================================================================================
 
 #define GOOSEFEM_WORLD_VERSION 0
 #define GOOSEFEM_MAJOR_VERSION 0
-#define GOOSEFEM_MINOR_VERSION 2
+#define GOOSEFEM_MINOR_VERSION 3
 
 #define GOOSEFEM_VERSION_AT_LEAST(x,y,z) \
   (GOOSEFEM_WORLD_VERSION>x || (GOOSEFEM_WORLD_VERSION>=x && \
   (GOOSEFEM_MAJOR_VERSION>y || (GOOSEFEM_MAJOR_VERSION>=y && \
                                 GOOSEFEM_MINOR_VERSION>=z))))
 
 #define GOOSEFEM_VERSION(x,y,z) \
   (GOOSEFEM_WORLD_VERSION==x && \
    GOOSEFEM_MAJOR_VERSION==y && \
    GOOSEFEM_MINOR_VERSION==z)
 
 // =================================================================================================
 
 namespace GooseFEM {
 
 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> MatD;
 typedef Eigen::Matrix<size_t, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> MatS;
 typedef Eigen::Matrix<double, Eigen::Dynamic,              1, Eigen::ColMajor> ColD;
 typedef Eigen::Matrix<size_t, Eigen::Dynamic,              1, Eigen::ColMajor> ColS;
 typedef Eigen::Matrix<int   , Eigen::Dynamic,              1, Eigen::ColMajor> ColI;
 
 // -------------------------------------------------------------------------------------------------
 
 }
 
 #endif