diff --git a/src/GooseFEM/DynamicsDiagonalPeriodic.h b/src/GooseFEM/DynamicsDiagonalPeriodic.h
index 04117a9..351621b 100644
--- a/src/GooseFEM/DynamicsDiagonalPeriodic.h
+++ b/src/GooseFEM/DynamicsDiagonalPeriodic.h
@@ -1,343 +1,345 @@
 /* ========================================== 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(std::shared_ptr<Element> elem, const MatD &x0,const MatS &conn,const MatS &dofs,
     double dt, double alpha=0.0 );
 
   // functions
   // ---------
 
   void velocityVerlet();  // 1 time step of time integrator
   void Verlet();          // 1 time step of time integrator (Fv == 0)
-  void computeMinv();     // element loop <- "elem->computeM"
-  void computeFu();       // element loop <- "elem->computeFu"
-  void computeFv();
-  void post();
+  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;
 
   // check
   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
+  // 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.noalias() = A;  // accelerations (DOFs)
   V_n.noalias() = 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.noalias() = A;  // accelerations (DOFs)
   V_n.noalias() = V;  // nodal velocities
   t            += dt; // current time instance
 }
 
 // =================================================================================================
 
 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) );
+      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/DynamicsDiagonalQuad4.h b/src/GooseFEM/DynamicsDiagonalQuad4.h
index 40b42ae..3e0e22d 100644
--- a/src/GooseFEM/DynamicsDiagonalQuad4.h
+++ b/src/GooseFEM/DynamicsDiagonalQuad4.h
@@ -1,326 +1,349 @@
-/* ========================================== DESCRIPTION ==========================================
+/* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
 #ifndef GOOSEFEM_DYNAMICS_DIAGONAL_QUAD4_H
 #define GOOSEFEM_DYNAMICS_DIAGONAL_QUAD4_H
 
 #include "Macros.h"
 #include <cppmat/cppmat.h>
 
 // -------------------------------------------------------------------------------------------------
 
 namespace GooseFEM {
 namespace Dynamics {
 namespace Diagonal {
 
 // -------------------------------------------------------------------------------------------------
 
 using vec = cppmat::cartesian2d::vector  <double>;
 using T2  = cppmat::cartesian2d::tensor2 <double>;
 using T2s = cppmat::cartesian2d::tensor2s<double>;
 using T2d = cppmat::cartesian2d::tensor2d<double>;
 
-// =================================================================================================
+// ========================== N.B. most loops are unrolled for efficiency ==========================
 
 template <class QuadraturePoint>
 class Quad4
 {
 public:
 
   // variables
   // ---------
 
   // arrays / matrices
   cppmat::tiny::matrix2<double,4,2> xe, ue, ve, xi, xi_n, dNdxi, dNdx;
   cppmat::tiny::vector <double,4>   w, w_n;
   cppmat::tiny::matrix2<double,8,8> M;
   cppmat::tiny::vector <double,8>   fu, fv;
   // tensors
   cppmat::cartesian2d::tensor2<double> J, Jinv, gradu, gradv;
   // scalars
   double Jdet, V;
   // sizes (nodes per element, dimensions, quadrature points)
   size_t nne=4, ndim=2, nk=4;
   // quadrature point
   std::shared_ptr<QuadraturePoint> quad;
 
   // constructor
   // -----------
 
   Quad4(std::shared_ptr<QuadraturePoint> quad);
 
   // functions
   // ---------
 
   void computeM (size_t elem); // mass matrix                     <- quad->density
   void computeFu(size_t elem); // displacement dependent forces   <- quad->stressStrain
   void computeFv(size_t elem); // displacement dependent forces   <- quad->stressStrainRate
   void post     (size_t elem); // post-process                    <- quad->stressStrain(Rate)
 
 };
 
 // =================================================================================================
 
 template <class QuadraturePoint>
 Quad4<QuadraturePoint>::Quad4(std::shared_ptr<QuadraturePoint> _quad)
 {
   // quadrature point routines
   quad = _quad;
 
   // integration point coordinates/weights: normal Gauss integration
   xi(0,0) = -1./std::sqrt(3.); xi(0,1) = -1./std::sqrt(3.); w(0) = 1.;
   xi(1,0) = +1./std::sqrt(3.); xi(1,1) = -1./std::sqrt(3.); w(1) = 1.;
   xi(2,0) = +1./std::sqrt(3.); xi(2,1) = +1./std::sqrt(3.); w(2) = 1.;
   xi(3,0) = -1./std::sqrt(3.); xi(3,1) = +1./std::sqrt(3.); w(3) = 1.;
 
   // integration point coordinates/weights: integration at the nodes
   xi_n(0,0) = -1.; xi_n(0,1) = -1.; w_n(0) = 1.;
   xi_n(1,0) = +1.; xi_n(1,1) = -1.; w_n(1) = 1.;
   xi_n(2,0) = +1.; xi_n(2,1) = +1.; w_n(2) = 1.;
   xi_n(3,0) = -1.; xi_n(3,1) = +1.; w_n(3) = 1.;
 }
 
 // =================================================================================================
 
 template <class QuadraturePoint>
 void Quad4<QuadraturePoint>::computeM(size_t elem)
 {
   // zero-initialize element mass matrix
   M.zeros();
 
   // loop over integration points (coincide with the nodes to get a diagonal mass matrix)
   for ( size_t k = 0 ; k < nne ; ++k )
   {
     // - shape function gradients (local coordinates)
     dNdxi(0,0) = -.25*(1.-xi_n(k,1)); dNdxi(0,1) = -.25*(1.-xi_n(k,0));
     dNdxi(1,0) = +.25*(1.-xi_n(k,1)); dNdxi(1,1) = -.25*(1.+xi_n(k,0));
     dNdxi(2,0) = +.25*(1.+xi_n(k,1)); dNdxi(2,1) = +.25*(1.+xi_n(k,0));
     dNdxi(3,0) = -.25*(1.+xi_n(k,1)); dNdxi(3,1) = +.25*(1.-xi_n(k,0));
 
     // - Jacobian
+    //   J(i,j) += dNdxi(m,i) * xe(m,j)
     J.zeros();
-    for ( size_t m = 0 ; m < nne ; ++m )
-      for ( size_t i = 0 ; i < ndim ; ++i )
-        for ( size_t j = 0 ; j < ndim ; ++j )
-          J(i,j) += dNdxi(m,i) * xe(m,j);
+    for ( size_t m = 0 ; m < nne ; ++m ) {
+      J(0,0) += dNdxi(m,0) * xe(m,0);
+      J(0,1) += dNdxi(m,0) * xe(m,1);
+      J(1,0) += dNdxi(m,1) * xe(m,0);
+      J(1,1) += dNdxi(m,1) * xe(m,1);
+    }
 
     // - determinant and inverse of the Jacobian
     Jdet = J.det();
     Jinv = J.inv();
 
     // - integration point volume (== volume associated with the node, in this element)
     V = w_n(k) * Jdet;
 
     // - assemble to element mass matrix (use the delta properties of the shape functions)
-    for ( size_t i = 0 ; i < ndim ; ++i )
-      M(i,i) += quad->density(elem,k,V) * V;
+    //   M(m+i,n+i) = N(m) * rho * V * N(n);
+    M(k*2  ,k*2  ) = quad->density(elem,k,V) * V;
+    M(k*2+1,k*2+1) = quad->density(elem,k,V) * V;
   }
 }
 
 // =================================================================================================
 
 template <class QuadraturePoint>
 void Quad4<QuadraturePoint>::computeFu(size_t elem)
 {
   // zero-initialize element force
   fu.zeros();
 
   // loop over integration points
   for ( size_t k = 0 ; k < nk ; ++k )
   {
     // - shape function gradients (local coordinates)
     dNdxi(0,0) = -.25*(1.-xi(k,1)); dNdxi(0,1) = -.25*(1.-xi(k,0));
     dNdxi(1,0) = +.25*(1.-xi(k,1)); dNdxi(1,1) = -.25*(1.+xi(k,0));
     dNdxi(2,0) = +.25*(1.+xi(k,1)); dNdxi(2,1) = +.25*(1.+xi(k,0));
     dNdxi(3,0) = -.25*(1.+xi(k,1)); dNdxi(3,1) = +.25*(1.-xi(k,0));
 
     // - Jacobian
+    //   J(i,j) += dNdxi(m,i) * xe(m,j)
     J.zeros();
-    for ( size_t m = 0 ; m < nne ; ++m )
-      for ( size_t i = 0 ; i < ndim ; ++i )
-        for ( size_t j = 0 ; j < ndim ; ++j )
-          J(i,j) += dNdxi(m,i) * xe(m,j);
+    for ( size_t m = 0 ; m < nne ; ++m ) {
+      J(0,0) += dNdxi(m,0) * xe(m,0);
+      J(0,1) += dNdxi(m,0) * xe(m,1);
+      J(1,0) += dNdxi(m,1) * xe(m,0);
+      J(1,1) += dNdxi(m,1) * xe(m,1);
+    }
 
     // - determinant and inverse of the Jacobian
     Jdet = J.det();
     Jinv = J.inv();
 
     // - integration point volume
     V = w(k) * Jdet;
 
     // - shape function gradients (global coordinates)
-    dNdx.zeros();
-    for ( size_t m = 0 ; m < nne ; ++m )
-      for ( size_t i = 0 ; i < ndim ; ++i )
-        for ( size_t j = 0 ; j < ndim ; ++j )
-          dNdx(m,i) += Jinv(i,j) * dNdxi(m,j);
+    //   dNdx(m,i) += Jinv(i,j) * dNdxi(m,j)
+    for ( size_t m = 0 ; m < nne ; ++m ) {
+      dNdx(m,0) = Jinv(0,0) * dNdxi(m,0) + Jinv(0,1) * dNdxi(m,1);
+      dNdx(m,1) = Jinv(1,0) * dNdxi(m,0) + Jinv(1,1) * dNdxi(m,1);
+    }
 
     // - displacement gradient
+    //   gradu(i,j) += dNdx(m,i) * ue(m,j)
     gradu.zeros();
-    for ( size_t m = 0 ; m < nne ; ++m )
-      for ( size_t i = 0 ; i < ndim ; ++i )
-        for ( size_t j = 0 ; j < ndim ; ++j )
-          gradu(i,j) += dNdx(m,i) * ue(m,j);
+    for ( size_t m = 0 ; m < nne ; ++m ) {
+      gradu(0,0) += dNdx(m,0) * ue(m,0);
+      gradu(0,1) += dNdx(m,0) * ue(m,1);
+      gradu(1,0) += dNdx(m,1) * ue(m,0);
+      gradu(1,1) += dNdx(m,1) * ue(m,1);
+    }
 
     // - strain tensor (symmetric part of "gradu")
-    for ( size_t i = 0 ; i < ndim ; ++i ) {
-      for ( size_t j = i ; j < ndim ; ++j ) {
-        quad->eps(i,j) = .5 * ( gradu(i,j) + gradu(j,i) );
-        quad->eps(j,i) = quad->eps(i,j);
-      }
-    }
+    //   quad->eps(i,j) = .5 * ( gradu(i,j) + gradu(j,i) )
+    quad->eps(0,0) = gradu(0,0);
+    quad->eps(0,1) = .5 * ( gradu(0,1) + gradu(1,0) );
+    quad->eps(1,0) = quad->eps(0,1);
+    quad->eps(1,1) = gradu(1,1);
 
     // - constitutive response
     quad->stressStrain(elem,k,V);
 
     // - assemble to element force
-    for ( size_t m = 0 ; m < nne ; ++m )
-      for ( size_t i = 0 ; i < ndim ; ++i )
-        for ( size_t j = 0 ; j < ndim ; ++j )
-          fu(m*ndim+j) += dNdx(m,i) * quad->sig(i,j) * V;
+    //   fu(m*ndim+j) += dNdx(m,i) * quad->sig(i,j) * V;
+    for ( size_t m = 0 ; m < nne ; ++m ) {
+      fu(m*ndim+0) += dNdx(m,0) * quad->sig(0,0) * V + dNdx(m,1) * quad->sig(1,0) * V;
+      fu(m*ndim+1) += dNdx(m,0) * quad->sig(0,1) * V + dNdx(m,1) * quad->sig(1,1) * V;
+    }
   }
 }
 
 // =================================================================================================
 
 template <class QuadraturePoint>
 void Quad4<QuadraturePoint>::computeFv(size_t elem)
 {
   // zero-initialize element force
   fv.zeros();
 
   // loop over integration points
   for ( size_t k = 0 ; k < nk ; ++k )
   {
     // - shape function gradients (local coordinates)
     dNdxi(0,0) = -.25*(1.-xi(k,1)); dNdxi(0,1) = -.25*(1.-xi(k,0));
     dNdxi(1,0) = +.25*(1.-xi(k,1)); dNdxi(1,1) = -.25*(1.+xi(k,0));
     dNdxi(2,0) = +.25*(1.+xi(k,1)); dNdxi(2,1) = +.25*(1.+xi(k,0));
     dNdxi(3,0) = -.25*(1.+xi(k,1)); dNdxi(3,1) = +.25*(1.-xi(k,0));
 
     // - Jacobian
+    //   J(i,j) += dNdxi(m,i) * xe(m,j)
     J.zeros();
-    for ( size_t m = 0 ; m < nne ; ++m )
-      for ( size_t i = 0 ; i < ndim ; ++i )
-        for ( size_t j = 0 ; j < ndim ; ++j )
-          J(i,j) += dNdxi(m,i) * xe(m,j);
+    for ( size_t m = 0 ; m < nne ; ++m ) {
+      J(0,0) += dNdxi(m,0) * xe(m,0);
+      J(0,1) += dNdxi(m,0) * xe(m,1);
+      J(1,0) += dNdxi(m,1) * xe(m,0);
+      J(1,1) += dNdxi(m,1) * xe(m,1);
+    }
 
     // - determinant and inverse of the Jacobian
     Jdet = J.det();
     Jinv = J.inv();
 
     // - integration point volume
     V = w(k) * Jdet;
 
     // - shape function gradients (global coordinates)
-    dNdx.zeros();
-    for ( size_t m = 0 ; m < nne ; ++m )
-      for ( size_t i = 0 ; i < ndim ; ++i )
-        for ( size_t j = 0 ; j < ndim ; ++j )
-          dNdx(m,i) += Jinv(i,j) * dNdxi(m,j);
+		//   dNdx(m,i) += Jinv(i,j) * dNdxi(m,j)
+    for ( size_t m = 0 ; m < nne ; ++m ) {
+      dNdx(m,0) = Jinv(0,0) * dNdxi(m,0) + Jinv(0,1) * dNdxi(m,1);
+      dNdx(m,1) = Jinv(1,0) * dNdxi(m,0) + Jinv(1,1) * dNdxi(m,1);
+    }
 
     // - velocity gradient
+		//   gradv(i,j) += dNdx(m,i) * ve(m,j)
     gradv.zeros();
-    for ( size_t m = 0 ; m < nne ; ++m )
-      for ( size_t i = 0 ; i < ndim ; ++i )
-        for ( size_t j = 0 ; j < ndim ; ++j )
-          gradv(i,j) += dNdx(m,i) * ve(m,j);
-
-    // - strain-rate tensor (symmetric part of "gradu")
-    for ( size_t i = 0 ; i < ndim ; ++i ) {
-      for ( size_t j = i ; j < ndim ; ++j ) {
-        quad->epsdot(i,j) = .5 * ( gradu(i,j) + gradu(j,i) );
-        quad->epsdot(j,i) = quad->epsdot(i,j);
-      }
+    for ( size_t m = 0 ; m < nne ; ++m ) {
+      gradv(0,0) += dNdx(m,0) * ve(m,0);
+      gradv(0,1) += dNdx(m,0) * ve(m,1);
+      gradv(1,0) += dNdx(m,1) * ve(m,0);
+      gradv(1,1) += dNdx(m,1) * ve(m,1);
     }
 
+    // - strain-rate tensor (symmetric part of "gradv")
+		//   quad->epsdot(i,j) = .5 * ( gradv(i,j) + gradv(j,i) )
+    quad->epsdot(0,0) = gradv(0,0);
+    quad->epsdot(0,1) = .5 * ( gradv(0,1) + gradv(1,0) );
+    quad->epsdot(1,0) = quad->epsdot(0,1);
+    quad->epsdot(1,1) = gradv(1,1);
+
     // - constitutive response
     quad->stressStrainRate(elem,k,V);
 
     // - assemble to element force
-    for ( size_t m = 0 ; m < nne ; ++m )
-      for ( size_t i = 0 ; i < ndim ; ++i )
-        for ( size_t j = 0 ; j < ndim ; ++j )
-          fv(m*ndim+j) += dNdx(m,i) * quad->sig(i,j) * V;
+		//   fv(m*ndim+j) += dNdx(m,i) * quad->sig(i,j) * V;
+    for ( size_t m = 0 ; m < nne ; ++m ) {
+      fv(m*2+0) += dNdx(m,0) * quad->sig(0,0) * V + dNdx(m,1) * quad->sig(1,0) * V;
+      fv(m*2+1) += dNdx(m,0) * quad->sig(0,1) * V + dNdx(m,1) * quad->sig(1,1) * V;
+    }
   }
 }
 
 // =================================================================================================
 
 template <class QuadraturePoint>
 void Quad4<QuadraturePoint>::post(size_t elem)
 {
   // loop over integration points
   for ( size_t k = 0 ; k < nk ; ++k )
   {
     // - shape function gradients (local coordinates)
     dNdxi(0,0) = -.25*(1.-xi(k,1)); dNdxi(0,1) = -.25*(1.-xi(k,0));
     dNdxi(1,0) = +.25*(1.-xi(k,1)); dNdxi(1,1) = -.25*(1.+xi(k,0));
     dNdxi(2,0) = +.25*(1.+xi(k,1)); dNdxi(2,1) = +.25*(1.+xi(k,0));
     dNdxi(3,0) = -.25*(1.+xi(k,1)); dNdxi(3,1) = +.25*(1.-xi(k,0));
 
     // - Jacobian
+    //   J(i,j) += dNdxi(m,i) * xe(m,j)
     J.zeros();
-    for ( size_t m = 0 ; m < nne ; ++m )
-      for ( size_t i = 0 ; i < ndim ; ++i )
-        for ( size_t j = 0 ; j < ndim ; ++j )
-          J(i,j) += dNdxi(m,i) * xe(m,j);
+    for ( size_t m = 0 ; m < nne ; ++m ) {
+      J(0,0) += dNdxi(m,0) * xe(m,0);
+      J(0,1) += dNdxi(m,0) * xe(m,1);
+      J(1,0) += dNdxi(m,1) * xe(m,0);
+      J(1,1) += dNdxi(m,1) * xe(m,1);
+    }
 
     // - determinant and inverse of the Jacobian
     Jdet = J.det();
     Jinv = J.inv();
 
     // - integration point volume
     V = w(k) * Jdet;
 
     // - shape function gradients (global coordinates)
-    dNdx.zeros();
-    for ( size_t m = 0 ; m < nne ; ++m )
-      for ( size_t i = 0 ; i < ndim ; ++i )
-        for ( size_t j = 0 ; j < ndim ; ++j )
-          dNdx(m,i) += Jinv(i,j) * dNdxi(m,j);
+    //   dNdx(m,i) += Jinv(i,j) * dNdxi(m,j)
+    for ( size_t m = 0 ; m < nne ; ++m ) {
+      dNdx(m,0) = Jinv(0,0) * dNdxi(m,0) + Jinv(0,1) * dNdxi(m,1);
+      dNdx(m,1) = Jinv(1,0) * dNdxi(m,0) + Jinv(1,1) * dNdxi(m,1);
+    }
 
     // - displacement gradient
+    //   gradu(i,j) += dNdx(m,i) * ue(m,j)
     gradu.zeros();
-    for ( size_t m = 0 ; m < nne ; ++m )
-      for ( size_t i = 0 ; i < ndim ; ++i )
-        for ( size_t j = 0 ; j < ndim ; ++j )
-          gradu(i,j) += dNdx(m,i) * ue(m,j);
+    for ( size_t m = 0 ; m < nne ; ++m ) {
+      gradu(0,0) += dNdx(m,0) * ue(m,0);
+      gradu(0,1) += dNdx(m,0) * ue(m,1);
+      gradu(1,0) += dNdx(m,1) * ue(m,0);
+      gradu(1,1) += dNdx(m,1) * ue(m,1);
+    }
 
     // - strain tensor (symmetric part of "gradu")
-    for ( size_t i = 0 ; i < ndim ; ++i ) {
-      for ( size_t j = i ; j < ndim ; ++j ) {
-        quad->eps(i,j) = .5 * ( gradu(i,j) + gradu(j,i) );
-        quad->eps(j,i) = quad->eps(i,j);
-      }
-    }
+    //   quad->eps(i,j) = .5 * ( gradu(i,j) + gradu(j,i) )
+    quad->eps(0,0) = gradu(0,0);
+    quad->eps(0,1) = .5 * ( gradu(0,1) + gradu(1,0) );
+    quad->eps(1,0) = quad->eps(0,1);
+    quad->eps(1,1) = gradu(1,1);
 
     // - velocity gradient
+    //   gradv(i,j) += dNdx(m,i) * ve(m,j)
     gradv.zeros();
-    for ( size_t m = 0 ; m < nne ; ++m )
-      for ( size_t i = 0 ; i < ndim ; ++i )
-        for ( size_t j = 0 ; j < ndim ; ++j )
-          gradv(i,j) += dNdx(m,i) * ve(m,j);
-
-    // - strain-rate tensor (symmetric part of "gradu")
-    for ( size_t i = 0 ; i < ndim ; ++i ) {
-      for ( size_t j = i ; j < ndim ; ++j ) {
-        quad->epsdot(i,j) = .5 * ( gradu(i,j) + gradu(j,i) );
-        quad->epsdot(j,i) = quad->epsdot(i,j);
-      }
+    for ( size_t m = 0 ; m < nne ; ++m ) {
+      gradv(0,0) += dNdx(m,0) * ve(m,0);
+      gradv(0,1) += dNdx(m,0) * ve(m,1);
+      gradv(1,0) += dNdx(m,1) * ve(m,0);
+      gradv(1,1) += dNdx(m,1) * ve(m,1);
     }
 
+    // - strain-rate tensor (symmetric part of "gradv")
+    //   quad->epsdot(i,j) = .5 * ( gradv(i,j) + gradv(j,i) )
+    quad->epsdot(0,0) = gradv(0,0);
+    quad->epsdot(0,1) = .5 * ( gradv(0,1) + gradv(1,0) );
+    quad->epsdot(1,0) = quad->epsdot(0,1);
+    quad->epsdot(1,1) = gradv(1,1);
+
     // - constitutive response
     quad->stressStrainPost    (elem,k,V);
     quad->stressStrainRatePost(elem,k,V);
   }
 }
 
 // =================================================================================================
 
 }}} // namespace GooseFEM::Dynamics::Diagonal
 
 #endif