diff --git a/Apps/FiniteDefElasticity/LinearMaterial.C b/Apps/FiniteDefElasticity/LinearMaterial.C
index 5c5f053e..43a3ac61 100644
--- a/Apps/FiniteDefElasticity/LinearMaterial.C
+++ b/Apps/FiniteDefElasticity/LinearMaterial.C
@@ -43,7 +43,7 @@ bool LinearMaterial::evaluate (Matrix& C, SymmTensor& sigma, double& U,
   // Push-forward the constitutive matrix to current configuration
 
   size_t i, j, k, l;
-  size_t n = F.dim();
+  size_t n = eps.dim();
   size_t m = C.rows();
   Matrix T(m,m), Ctmp;
 
diff --git a/Apps/FiniteDefElasticity/LinearMaterial.h b/Apps/FiniteDefElasticity/LinearMaterial.h
index ffcb934f..7e5c8f93 100644
--- a/Apps/FiniteDefElasticity/LinearMaterial.h
+++ b/Apps/FiniteDefElasticity/LinearMaterial.h
@@ -18,13 +18,21 @@
 
 
 /*!
-  \brief Class representing a isotropic linear elastic material model.
+  \brief Class representing a general linear elastic material model.
+  \details This class is a wrapper for any linear-elastic material model,
+  when it is to be used in a nonlinear finite element formulation based on
+  the updated Lagrangian formulation. The \a evaluate method of this class
+  just invokes the corresponding method of the object pointed to by the
+  \a material member, and then transforms the resulting constitutive matrix
+  and stress tensor to the current updated reference frame, based on the
+  supplied deformation tensor, \a F.
 */
 
 class LinearMaterial : public Material
 {
 public:
   //! \brief The constructor initializes the material properties pointer.
+  //! \param[in] Pointer to material model for linear elastic analysis.
   LinearMaterial(const Material* mat) : material(mat) {}
   //! \brief The destructor deletes the material properties object.
   virtual ~LinearMaterial() { delete material; }
@@ -53,7 +61,7 @@ public:
 			char iop = 1, const TimeDomain* prm = 0) const;
 
 private:
-  const Material* material; //!< Pointer to actual material properties object
+  const Material* material; //!< Linear-elastic material properties object
 };
 
 #endif
diff --git a/Apps/FiniteDefElasticity/Material/cons2d.f b/Apps/FiniteDefElasticity/Material/cons2d.f
index 91bffcb6..7d1f47ca 100644
--- a/Apps/FiniteDefElasticity/Material/cons2d.f
+++ b/Apps/FiniteDefElasticity/Material/cons2d.f
@@ -1,29 +1,27 @@
-      subroutine cons2d (ipsw,iter,iwr,lfirst,mTYP,mVER,
-     &                   detF,F,Fp,pmat,HV,U,Sig,Cst,ierr)
+      subroutine cons2d (ipsw,iter,iwr,lfirst,mTYP,mVER,nHV,nTM,
+     &                   detF,nDF,F,pMAT,HV,U,Sig,Cst,ierr)
 C
 C     CONS2D: 2D wrapper on the 3D material of FENRIS.
 C
       implicit none
 C
-      integer  ipsw, iter, iwr, lfirst, mTYP, mVER, ierr
-      real*8   detF, F(4), Fp(4), pmat(*), HV, U, Sig(3), Cst(3,3)
-      real*8   F3D(3,3), F3P(3,3), S3D(6), C3D(6,6)
+      integer  ipsw, iter, iwr, lfirst, mTYP, mVER, nHV, nTM, nDF, ierr
+      real*8   detF, F(*), pMAT(*), HV(*), U, Sig(3), Cst(3,3)
+      real*8   F3D(3,3), S3D(6), C3D(6,6)
 C
-      F3D      = 0.0D0
-      F3D(1,1) = F(1)
-      F3D(2,1) = F(2)
-      F3D(1,2) = F(3)
-      F3D(2,2) = F(4)
-      F3D(3,3) = 1.0D0
-      F3P      = 0.0D0
-      F3P(1,1) = Fp(1)
-      F3P(2,1) = Fp(2)
-      F3P(1,2) = Fp(3)
-      F3P(2,2) = Fp(4)
-      F3P(3,3) = 1.0D0
-C
-      call cons3d (ipsw,iter,iwr,lfirst,mTYP,mVER,
-     &             detF,F3D,F3P,pmat,HV,U,S3D,C3D,ierr)
+      if (ndf .eq. 2) then
+         F3D      = 0.0D0
+         F3D(1,1) = F(1)
+         F3D(2,1) = F(2)
+         F3D(1,2) = F(3)
+         F3D(2,2) = F(4)
+         F3D(3,3) = 1.0D0
+         call cons3d (ipsw,iter,iwr,lfirst,mTYP,mVER,nHV,nTM,
+     &                detF,F3D,pMAT,HV,U,S3D,C3D,ierr)
+      else
+         call cons3d (ipsw,iter,iwr,lfirst,mTYP,mVER,nHV,nTM,
+     &                detF,F,pMAT,HV,U,S3D,C3D,ierr)
+      end if
 C
       Cst(1:2,1:2) = C3D(1:2,1:2)
       Cst(1:2,3)   = C3D(1:2,4)
diff --git a/Apps/FiniteDefElasticity/Material/plas2d.f b/Apps/FiniteDefElasticity/Material/plas2d.f
index 6dd69629..7ee8bd18 100644
--- a/Apps/FiniteDefElasticity/Material/plas2d.f
+++ b/Apps/FiniteDefElasticity/Material/plas2d.f
@@ -1,30 +1,34 @@
       subroutine plas2d (ipsw,iwr,iter,lfirst,ntm,pmat,detF,
-     &                   Fn1,Fn,be,Epp,Epl,Sig,Cst,ierr)
+     &                   ndf,Fn1,Fn,be,Epp,Epl,Sig,Cst,ierr)
 C
 C     PLAS2D: 2D wrapper on the 3D plastic material of FENRIS.
 C
       implicit none
 C
-      integer  ipsw, iwr, iter, lfirst, ntm, ierr
-      real*8   detF, Fn1(4), Fn(4), pmat(*), be(*)
+      integer  ipsw, iwr, iter, lfirst, ntm, ndf, ierr
+      real*8   detF, Fn1(*), Fn(*), pmat(*), be(*)
       real*8   Epp, Epl(*), Sig(3), Cst(3,3)
       real*8   F3D(3,3), F3P(3,3), S3D(6), C3D(6,6)
 C
-      F3D      = 0.0D0
-      F3D(1,1) = Fn(1)
-      F3D(2,1) = Fn(2)
-      F3D(1,2) = Fn(3)
-      F3D(2,2) = Fn(4)
-      F3D(3,3) = 1.0D0
-      F3P      = 0.0D0
-      F3P(1,1) = Fn1(1)
-      F3P(2,1) = Fn1(2)
-      F3P(1,2) = Fn1(3)
-      F3P(2,2) = Fn1(4)
-      F3P(3,3) = 1.0D0
-C
-      call plas3d (ipsw,iwr,iter,lfirst,ntm,pmat,detF,
-     &             F3P,F3D,be,Epp,Epl,S3D,C3D,ierr)
+      if (ndf .eq. 2) then
+         F3D      = 0.0D0
+         F3D(1,1) = Fn(1)
+         F3D(2,1) = Fn(2)
+         F3D(1,2) = Fn(3)
+         F3D(2,2) = Fn(4)
+         F3D(3,3) = 1.0D0
+         F3P      = 0.0D0
+         F3P(1,1) = Fn1(1)
+         F3P(2,1) = Fn1(2)
+         F3P(1,2) = Fn1(3)
+         F3P(2,2) = Fn1(4)
+         F3P(3,3) = 1.0D0
+         call plas3d (ipsw,iwr,iter,lfirst,ntm,pmat,detF,
+     &                F3P,F3D,be,Epp,Epl,S3D,C3D,ierr)
+      else
+         call plas3d (ipsw,iwr,iter,lfirst,ntm,pmat,detF,
+     &                Fn,Fn1,be,Epp,Epl,S3D,C3D,ierr)
+      end if
 C
       Cst(1:2,1:2) = C3D(1:2,1:2)
       Cst(1:2,3)   = C3D(1:2,4)
diff --git a/Apps/FiniteDefElasticity/NeoHookeMaterial.C b/Apps/FiniteDefElasticity/NeoHookeMaterial.C
index be3401ab..0d262b78 100644
--- a/Apps/FiniteDefElasticity/NeoHookeMaterial.C
+++ b/Apps/FiniteDefElasticity/NeoHookeMaterial.C
@@ -20,14 +20,15 @@ extern "C" {
   void cons2d_(const int& ipsw, const int& iter, const int& iwr,
 	       const int& lfirst, const int& mTYP, const int& mVER,
 	       const int& nHV, const int& nTM, const double& detF,
-	       const double* F, const double* pMAT, double* HV, double& Engy,
-	       const double* Sig, double* Cst, int& ierr);
+	       const int& nDF,
+	       const double* F, const double* pMAT, double* HV,
+	       double& Engy, const double* Sig, double* Cst, int& ierr);
   //! \brief Interface to 3D nonlinear material routines (FORTRAN-77 code).
   void cons3d_(const int& ipsw, const int& iter, const int& iwr,
 	       const int& lfirst, const int& mTYP, const int& mVER,
 	       const int& nHV, const int& nTM, const double& detF,
-	       const double* F, const double* pMAT, double* HV, double& Engy,
-	       const double* Sig, double* Cst, int& ierr);
+	       const double* F, const double* pMAT, double* HV,
+	       double& Engy, const double* Sig, double* Cst, int& ierr);
 }
 #ifndef INT_DEBUG
 #define INT_DEBUG 0
@@ -83,14 +84,14 @@ bool NeoHookeMaterial::evaluate (Matrix& C, SymmTensor& sigma, double& U,
     return false;
   }
 
-  size_t ndim = F.dim();
+  size_t ndim = sigma.dim();
   size_t ncmp = ndim*(ndim+1)/2;
   C.resize(ncmp,ncmp);
   int ierr = -99;
 #ifdef USE_FTNMAT
   // Invoke the FORTRAN routine for Neo-Hookean hyperelastic material models.
   if (ndim == 2)
-    cons2d_(INT_DEBUG,0,6,0,mTYP,mVER,0,3,J,F.ptr(),pmat,
+    cons2d_(INT_DEBUG,0,6,0,mTYP,mVER,0,3,J,F.dim(),F.ptr(),pmat,
 	    &U,U,sigma.ptr(),C.ptr(),ierr);
   else
     cons3d_(INT_DEBUG,0,6,0,mTYP,mVER,0,6,J,F.ptr(),pmat,
diff --git a/Apps/FiniteDefElasticity/NeoHookeMaterial.h b/Apps/FiniteDefElasticity/NeoHookeMaterial.h
index ce1a1a2b..ee76e387 100644
--- a/Apps/FiniteDefElasticity/NeoHookeMaterial.h
+++ b/Apps/FiniteDefElasticity/NeoHookeMaterial.h
@@ -19,6 +19,11 @@
 
 /*!
   \brief Class representing a Neo-Hookean hyperelastic material model.
+  \details This class is a wrapper for the FORTRAN material routines of FENRIS,
+  implemented in by the subroutine CONS3D and subroutines called from CONS3D.
+  Which material routine to use is governed by the \a mVER parameter.
+  Only isotropic material properties are supported by this class.
+  \note In 2D, only plane strain is supported by this class.
 */
 
 class NeoHookeMaterial : public LinIsotropic
diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityUL.C b/Apps/FiniteDefElasticity/NonlinearElasticityUL.C
index 15923a4b..7e6b10a4 100644
--- a/Apps/FiniteDefElasticity/NonlinearElasticityUL.C
+++ b/Apps/FiniteDefElasticity/NonlinearElasticityUL.C
@@ -15,6 +15,7 @@
 #include "MaterialBase.h"
 #include "ElmMats.h"
 #include "ElmNorm.h"
+#include "TimeDomain.h"
 #include "Tensor.h"
 #include "Vec3Oper.h"
 
@@ -85,7 +86,22 @@ void NonlinearElasticityUL::setMode (SIM::SolutionMode mode)
 }
 
 
-bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt, double detJW,
+void NonlinearElasticityUL::initIntegration (const TimeDomain& prm)
+{
+  if (material)
+    material->initIntegration(prm);
+}
+
+
+void NonlinearElasticityUL::initResultPoints ()
+{
+  if (material)
+    material->initResultPoints();
+}
+
+
+bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt,
+				     const TimeDomain& prm, double detJW,
 				     const Vector& N, const Matrix& dNdX,
 				     const Vec3& X) const
 {
@@ -128,7 +144,7 @@ bool NonlinearElasticityUL::evalInt (LocalIntegral*& elmInt, double detJW,
   if (eKm || eKg || iS)
   {
     double U = 0.0;
-    if (!material->evaluate(Cmat,sigma,U,X,F,E, (eKg || iS) && lHaveStrains))
+    if (!material->evaluate(Cmat,sigma,U,X,F,E,(eKg || iS),&prm))
       return false;
   }
 
@@ -226,7 +242,7 @@ bool NonlinearElasticityUL::evalBou (LocalIntegral*& elmInt, double detJW,
 bool NonlinearElasticityUL::formDefGradient (const Matrix& dNdX,
 					     Tensor& F) const
 {
-  SymmTensor dummy(0);
+  static SymmTensor dummy(0);
   return this->kinematics(dNdX,F,dummy);
 }
 
@@ -250,6 +266,7 @@ bool NonlinearElasticityUL::kinematics (const Matrix& dNdX,
     return false;
   }
 
+
   // Compute the deformation gradient, [F] = [I] + [dudX] = [I] + [dNdX]*[u],
   // and the Green-Lagrange strains, E_ij = 0.5(F_ij+F_ji+F_ki*F_kj).
 
@@ -257,12 +274,17 @@ bool NonlinearElasticityUL::kinematics (const Matrix& dNdX,
   // element displacement vector, *eV, as a matrix whose number of columns
   // equals the number of rows in the matrix dNdX.
   Matrix dUdX;
-  if (dUdX.multiplyMat(*eV,dNdX)) // dUdX = Grad{u} = eV*dNdX
-    F = dUdX;
-  else
+  if (!dUdX.multiplyMat(*eV,dNdX)) // dUdX = Grad{u} = eV*dNdX
     return false;
 
   unsigned short int i, j, k;
+  if (dUdX.rows() < F.dim())
+    F.zero();
+
+  // Cannot use operator= here, in case F is of higher dimension than dUdX
+  for (i = 1; i <= dUdX.rows(); i++)
+    for (j = 1; j <= dUdX.cols(); j++)
+      F(i,j) = dUdX(i,j);
 
   // Now form the Green-Lagrange strain tensor.
   // Note that for the shear terms (i/=j) we actually compute 2*E_ij
@@ -294,7 +316,14 @@ NormBase* NonlinearElasticityUL::getNormIntegrand (AnaSol*) const
 }
 
 
-bool ElasticityNormUL::evalInt (LocalIntegral*& elmInt, double detJW,
+void ElasticityNormUL::initIntegration (const TimeDomain& prm)
+{
+  problem.initIntegration(prm);
+}
+
+
+bool ElasticityNormUL::evalInt (LocalIntegral*& elmInt,
+				const TimeDomain& prm, double detJW,
 				const Vector&, const Matrix& dNdX,
 				const Vec3& X) const
 {
@@ -312,18 +341,12 @@ bool ElasticityNormUL::evalInt (LocalIntegral*& elmInt, double detJW,
   Matrix C;
   SymmTensor S(E.dim());
   double U = 0.0;
-  if (!ulp->material->evaluate(C,S,U,X,F,E,2))
+  if (!ulp->material->evaluate(C,S,U,X,F,E,2,&prm))
     return false;
 
-  if (U == 0.0)
-  {
-    // Integrate the energy norm a(u^h,u^h) = Int_Omega0 (S:E) dV0
-    const RealArray& sigma = S;
-    const RealArray& epsil = E;
-    for (size_t i = 0; i < sigma.size(); i++)
-      U += sigma[i]*epsil[i];
-  }
-
+  // Integrate the energy norm a(u^h,u^h) = Int_Omega0 (S:E) dV0
+  if (U == 0.0) U = S.innerProd(E);
   pnorm[0] += U*detJW;
+
   return true;
 }
diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityUL.h b/Apps/FiniteDefElasticity/NonlinearElasticityUL.h
index 1e2ce8f9..6b959d07 100644
--- a/Apps/FiniteDefElasticity/NonlinearElasticityUL.h
+++ b/Apps/FiniteDefElasticity/NonlinearElasticityUL.h
@@ -43,13 +43,21 @@ public:
   //! \param[in] mode The solution mode to use
   virtual void setMode(SIM::SolutionMode mode);
 
+  //! \brief Initializes the integrand for a new integration loop.
+  //! \param[in] prm Nonlinear solution algorithm parameters
+  virtual void initIntegration(const TimeDomain& prm);
+  //! \brief Initializes the integrand for a result point loop.
+  virtual void initResultPoints();
+
   //! \brief Evaluates the integrand at an interior point.
   //! \param elmInt The local integral object to receive the contributions
+  //! \param[in] prm Nonlinear solution algorithm parameters
   //! \param[in] detJW Jacobian determinant times integration point weight
   //! \param[in] N Basis function values
   //! \param[in] dNdX Basis function gradients
   //! \param[in] X Cartesian coordinates of current integration point
-  virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
+  virtual bool evalInt(LocalIntegral*& elmInt,
+		       const TimeDomain& prm, double detJW,
 		       const Vector& N, const Matrix& dNdX,
 		       const Vec3& X) const;
 
@@ -115,6 +123,10 @@ public:
   //! \brief Empty destructor.
   virtual ~ElasticityNormUL() {}
 
+  //! \brief Initializes the integrand for a new integration loop.
+  //! \param[in] prm Nonlinear solution algorithm parameters
+  virtual void initIntegration(const TimeDomain& prm);
+
   //! \brief Initializes current element for numerical integration (mixed).
   //! \param[in] MNPC1 Nodal point correspondance for the displacement field
   virtual bool initElement(const std::vector<int>& MNPC1,
@@ -132,25 +144,29 @@ public:
 
   //! \brief Evaluates the integrand at an interior point.
   //! \param elmInt The local integral object to receive the contributions
+  //! \param[in] prm Nonlinear solution algorithm parameters
   //! \param[in] detJW Jacobian determinant times integration point weight
   //! \param[in] N Basis function values
   //! \param[in] dNdX Basis function gradients
   //! \param[in] X Cartesian coordinates of current integration point
-  virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
+  virtual bool evalInt(LocalIntegral*& elmInt,
+		       const TimeDomain& prm, double detJW,
 		       const Vector& N, const Matrix& dNdX,
 		       const Vec3& X) const;
   //! \brief Evaluates the integrand at an interior point (mixed).
   //! \param elmInt The local integral object to receive the contributions
+  //! \param[in] prm Nonlinear solution algorithm parameters
   //! \param[in] detJW Jacobian determinant times integration point weight
   //! \param[in] N1 Basis function values for the displacement field
   //! \param[in] dN1dX Basis function gradients for the displacement field
   //! \param[in] X Cartesian coordinates of current integration point
-  virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
+  virtual bool evalInt(LocalIntegral*& elmInt,
+		       const TimeDomain& prm, double detJW,
                        const Vector& N1, const Vector&,
                        const Matrix& dN1dX, const Matrix&,
                        const Vec3& X) const
   {
-    return this->evalInt(elmInt,detJW,N1,dN1dX,X);
+    return this->evalInt(elmInt,prm,detJW,N1,dN1dX,X);
   }
 
   //! \brief Evaluates the integrand at a boundary point (mixed).
diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityULMX.C b/Apps/FiniteDefElasticity/NonlinearElasticityULMX.C
index 90eb1f22..b323f0e9 100644
--- a/Apps/FiniteDefElasticity/NonlinearElasticityULMX.C
+++ b/Apps/FiniteDefElasticity/NonlinearElasticityULMX.C
@@ -15,6 +15,7 @@
 #include "MaterialBase.h"
 #include "ElmMats.h"
 #include "ElmNorm.h"
+#include "TimeDomain.h"
 #include "Utilities.h"
 #include "Vec3Oper.h"
 
@@ -133,19 +134,13 @@ bool NonlinearElasticityULMX::initElement (const std::vector<int>& MNPC,
   size_t nPm = utl::Pascal(p,nsd);
   if (Hh) Hh->resize(nPm,nPm,true);
 
-  // Extract the element level displacements of previous increment
-  int ierr = 0;
-  if (!primsol[1].empty())
-    if ((ierr = utl::gather(MNPC,nsd,primsol[1],myMats->b[2])))
-      std::cerr <<" *** NonlinearElasticityULMX::initElement: Detected "
-		<< ierr <<" node numbers out of range."<< std::endl;
-
   // The other element matrices are initialized by the parent class method
-  return this->NonlinearElasticityUL::initElement(MNPC) && ierr == 0;
+  return this->NonlinearElasticityUL::initElement(MNPC);
 }
 
 
-bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt, double detJW,
+bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt,
+				       const TimeDomain&, double detJW,
 				       const Vector& N, const Matrix& dNdX,
 				       const Vec3& X) const
 {
@@ -172,7 +167,9 @@ bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt, double detJW,
 
   // Invert the deformation gradient ==> Fi
   Matrix Fi(nsd,nsd);
-  Fi.fill(ptData.F.ptr());
+  for (unsigned short int i = 1; i <= nsd; i++)
+    for (unsigned short int j = 1; j <= nsd; j++)
+      Fi(i,j) = ptData.F(i,j);
   double J = Fi.inverse();
   if (J == 0.0) return false;
 
@@ -207,7 +204,8 @@ bool NonlinearElasticityULMX::evalInt (LocalIntegral*& elmInt, double detJW,
 }
 
 
-bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt)
+bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt,
+					       const TimeDomain& prm)
 {
   if (!iS && !eKm && !eKg) return true;
   if (!Hh) return false;
@@ -365,21 +363,20 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt)
     // Evaluate the constitutive relation
     double U = 0.0;
     Sig.push_back(SymmTensor(3));
-    if (!material->evaluate(Cmat,Sig.back(),U,pt.X,pt.F,Sig.back()))
+    if (!material->evaluate(Cmat,Sig.back(),U,pt.X,pt.F,Sig.back(),true,&prm))
       return false;
 
 #ifdef USE_FTNMAT
     // Project the constitutive matrix for the mixed model
     D[iP].resize(7,7);
-    int ipsw = INT_DEBUG > 1 ? 9 : 0;
-    pcst3d_(Cmat.ptr(),D[iP].ptr(),ipsw,6);
+    pcst3d_(Cmat.ptr(),D[iP].ptr(),INT_DEBUG,6);
 #else
     std::cerr <<" *** NonlinearElasticityULMX::finalizeElement: "
 	      <<" Does not work when compiled without USE_FTNMAT"<< std::endl;
     return false;
 #endif
 
-    if (U > 0.0)
+    if (!Sig.back().isZero())
     {
       // Integrate mean pressure
       lHaveStress = true;
@@ -417,8 +414,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt)
     if (eKm)
     {
 #ifdef USE_FTNMAT
-      int ipsw = INT_DEBUG > 1 ? 9 : 0;
-      mdma3d_(Bpres,Press,Sig[iP].ptr(),D[iP].ptr(),ipsw,6);
+      mdma3d_(Bpres,Press,Sig[iP].ptr(),D[iP].ptr(),INT_DEBUG,6);
 
       // Integrate the material stiffness matrix
       D[iP] *= dFmix*pt.detJW;
@@ -473,16 +469,18 @@ bool ElasticityNormULMX::initElement (const std::vector<int>& MNPC,
 }
 
 
-bool ElasticityNormULMX::evalInt (LocalIntegral*& elmInt, double detJW,
+bool ElasticityNormULMX::evalInt (LocalIntegral*& elmInt,
+				  const TimeDomain& prm, double detJW,
 				  const Vector& N, const Matrix& dNdX,
 				  const Vec3& X) const
 {
   NonlinearElasticityULMX& elp = static_cast<NonlinearElasticityULMX&>(problem);
-  return elp.evalInt(elmInt,detJW,N,dNdX,X);
+  return elp.evalInt(elmInt,prm,detJW,N,dNdX,X);
 }
 
 
-bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt)
+bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt,
+					  const TimeDomain& prm)
 {
   NonlinearElasticityULMX& elp = static_cast<NonlinearElasticityULMX&>(problem);
   if (!elp.Hh) return false;
@@ -574,7 +572,7 @@ bool ElasticityNormULMX::finalizeElement (LocalIntegral*& elmInt)
 
     // Compute the strain energy density
     double U = 0.0;
-    if (!elp.material->evaluate(C,Sig,U,pt.X,pt.F,Sig,false))
+    if (!elp.material->evaluate(C,Sig,U,pt.X,pt.F,Sig,false,&prm))
       return false;
 
     // Integrate strain energy
diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityULMX.h b/Apps/FiniteDefElasticity/NonlinearElasticityULMX.h
index 6fe268ca..852729bd 100644
--- a/Apps/FiniteDefElasticity/NonlinearElasticityULMX.h
+++ b/Apps/FiniteDefElasticity/NonlinearElasticityULMX.h
@@ -63,6 +63,7 @@ public:
 
   //! \brief Evaluates the integrand at an interior point.
   //! \param elmInt The local integral object to receive the contributions
+  //! \param[in] prm Nonlinear solution algorithm parameters
   //! \param[in] detJW Jacobian determinant times integration point weight
   //! \param[in] N Basis function values
   //! \param[in] dNdX Basis function gradients
@@ -72,7 +73,8 @@ public:
   //! depending on the basis function values in the internal member \a myData.
   //! The actual numerical integration of the tangent stiffness and internal
   //! forces is then performed by the \a finalizeElement method.
-  virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
+  virtual bool evalInt(LocalIntegral*& elmInt,
+		       const TimeDomain& prm, double detJW,
 		       const Vector& N, const Matrix& dNdX,
 		       const Vec3& X) const;
 
@@ -82,7 +84,8 @@ public:
   //! modes. All data needed during the integration has been stored internally
   //! in the \a myData member by the \a evalInt method.
   //! \param elmInt The local integral object to receive the contributions
-  virtual bool finalizeElement(LocalIntegral*&);
+  //! \param[in] prm Nonlinear solution algorithm parameters
+  virtual bool finalizeElement(LocalIntegral*& elmInt, const TimeDomain& prm);
 
   //! \brief Returns a pointer to an Integrand for solution norm evaluation.
   //! \note The Integrand object is allocated dynamically and has to be deleted
@@ -128,6 +131,7 @@ public:
 
   //! \brief Evaluates the integrand at an interior point.
   //! \param elmInt The local integral object to receive the contributions
+  //! \param[in] prm Nonlinear solution algorithm parameters
   //! \param[in] detJW Jacobian determinant times integration point weight
   //! \param[in] N Basis function values
   //! \param[in] dNdX Basis function gradients
@@ -137,7 +141,8 @@ public:
   //! corresponding NonlinearElasticityULMX object, for which this class is
   //! declared as friend such that it can access the data members. The actual
   //! norm integration us then performed by the \a finalizeElement method.
-  virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
+  virtual bool evalInt(LocalIntegral*& elmInt,
+		       const TimeDomain& prm, double detJW,
 		       const Vector& N, const Matrix& dNdX,
 		       const Vec3& X) const;
 
@@ -146,7 +151,8 @@ public:
   //! loops within an element required by the mixed formulation. It performs
   //! a subset of the tasks done by NonlinearElasticityULMX::finalizeElement.
   //! \param elmInt The local integral object to receive the contributions
-  virtual bool finalizeElement(LocalIntegral*& elmInt);
+  //! \param[in] prm Nonlinear solution algorithm parameters
+  virtual bool finalizeElement(LocalIntegral*& elmInt, const TimeDomain& prm);
 
   //! \brief Returns which integrand to be used.
   virtual int getIntegrandType() const { return 4; }
diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C
index 7df1355b..e6c8d9a4 100644
--- a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C
+++ b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C
@@ -190,7 +190,7 @@ const Vector& NonlinearElasticityULMixed::MixedElmMats::getRHSVector () const
 
 
 NonlinearElasticityULMixed::NonlinearElasticityULMixed (unsigned short int n)
-  : NonlinearElasticityUL(n), Fbar(n), Dmat(7,7)
+  : NonlinearElasticityUL(n), Fbar(3), Dmat(7,7)
 {
   if (myMats) delete myMats;
   myMats = new MixedElmMats();
@@ -290,7 +290,8 @@ bool NonlinearElasticityULMixed::initElementBou (const std::vector<int>& MNPC1,
 }
 
 
-bool NonlinearElasticityULMixed::evalInt (LocalIntegral*& elmInt, double detJW,
+bool NonlinearElasticityULMixed::evalInt (LocalIntegral*& elmInt,
+					  const TimeDomain& prm, double detJW,
 					  const Vector& N1,
 					  const Vector& N2,
 					  const Matrix& dNdX1,
@@ -298,8 +299,9 @@ bool NonlinearElasticityULMixed::evalInt (LocalIntegral*& elmInt, double detJW,
 					  const Vec3& X) const
 {
   // Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
+  Tensor F(nsd);
   SymmTensor E(nsd);
-  if (!this->kinematics(dNdX1,Fbar,E))
+  if (!this->kinematics(dNdX1,F,E))
     return false;
 
   bool lHaveStrains = !E.isZero();
@@ -317,10 +319,8 @@ bool NonlinearElasticityULMixed::evalInt (LocalIntegral*& elmInt, double detJW,
   double dVol = Theta*detJW;
 
   // Compute the mixed model deformation gradient, F_bar
-  double r1 = pow(fabs(Theta/Fbar.det()),1.0/3.0);
-
-  Matrix Fi(nsd,nsd);
-  Fi.fill(Fbar.ptr());
+  Fbar = F; // notice that F_bar always has dimension 3
+  double r1 = pow(fabs(Theta/F.det()),1.0/3.0);
 
   Fbar *= r1;
   if (nsd == 2) // In 2D we always assume plane strain so set F(3,3)=1
@@ -329,6 +329,8 @@ bool NonlinearElasticityULMixed::evalInt (LocalIntegral*& elmInt, double detJW,
     return false;
 
   // Invert the deformation gradient ==> Fi
+  Matrix Fi(nsd,nsd);
+  Fi.fill(F.ptr());
   double J = Fi.inverse();
   if (J == 0.0) return false;
 
@@ -351,7 +353,7 @@ bool NonlinearElasticityULMixed::evalInt (LocalIntegral*& elmInt, double detJW,
   // Evaluate the constitutive relation
   SymmTensor Eps(3), Sig(3), Sigma(nsd);
   double U, Bpres = 0.0, Mpres = 0.0;
-  if (!material->evaluate(Cmat,Sig,U,X,Fbar,Eps=E,lHaveStrains))
+  if (!material->evaluate(Cmat,Sig,U,X,Fbar,Eps=E,lHaveStrains,&prm))
     return false;
 
 #ifdef USE_FTNMAT
diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h
index ecb5287e..4ea67f39 100644
--- a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h
+++ b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h
@@ -70,13 +70,15 @@ public:
 
   //! \brief Evaluates the mixed field problem integrand at an interior point.
   //! \param elmInt The local integral object to receive the contributions
+  //! \param[in] prm Nonlinear solution algorithm parameters
   //! \param[in] detJW Jacobian determinant times integration point weight
   //! \param[in] N1 Basis function values, field 1 (displacements)
   //! \param[in] N2 Basis function values, field 2 (pressure and volume change)
   //! \param[in] dN1dX Basis function gradients, field 1
   //! \param[in] dN2dX Basis function gradients, field 2
   //! \param[in] X Cartesian coordinates of current integration point
-  virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
+  virtual bool evalInt(LocalIntegral*& elmInt,
+		       const TimeDomain& prm, double detJW,
 		       const Vector& N1, const Vector& N2,
 		       const Matrix& dN1dX, const Matrix& dN2dX,
 		       const Vec3& X) const;
diff --git a/Apps/FiniteDefElasticity/PlasticMaterial.C b/Apps/FiniteDefElasticity/PlasticMaterial.C
index 9667e027..941150c0 100644
--- a/Apps/FiniteDefElasticity/PlasticMaterial.C
+++ b/Apps/FiniteDefElasticity/PlasticMaterial.C
@@ -19,7 +19,8 @@ extern "C" {
   //! \brief Interface to 2D elasto-plastic material routines (FORTRAN-77 code).
   void plas2d_(const int& ipsw, const int& iwr, const int& iter,
 	       const int& lfirst, const int& ntm, const double* pMAT,
-	       const double& detF, const double* Fn1, const double* Fn,
+	       const double& detF, const int& ndF,
+	       const double* Fn1, const double* Fn,
 	       const double* be, const double& Epp, const double* Epl,
 	       const double* Sig, double* Cst, int& ierr);
   //! \brief Interface to 3D elasto-plastic material routines (FORTRAN-77 code).
@@ -36,10 +37,10 @@ extern "C" {
 
 
 /*!
-  See the Fortran routine Material/plas3d.f for interpretation of \a pMAT.
+  See the Fortran the routine Material/plas3d.f for interpretation of \a pMAT.
 */
 
-PlasticPrm::PlasticPrm (const RealArray& p) : pMAT(p)
+PlasticMaterial::PlasticMaterial (const RealArray& p) : pMAT(p), iP1(0), iP2(0)
 {
   if (pMAT.size() < 11) pMAT.resize(11,0.0);
 
@@ -57,13 +58,13 @@ PlasticPrm::PlasticPrm (const RealArray& p) : pMAT(p)
   else if (nu < 0.5)
   {
     // Calculate the Bulk and Shear moduli from E and nu
-    Smod = Emod / (2.0 + 2.0*nu);
     Bmod = Emod / (3.0 - 6.0*nu);
+    Smod = Emod / (2.0 + 2.0*nu);
   }
   else
   {
-    Smod = Emod / 3.0;
     Bmod = 0.0;
+    Smod = Emod / 3.0;
   }
 
   pMAT[0] = Emod;
@@ -73,7 +74,17 @@ PlasticPrm::PlasticPrm (const RealArray& p) : pMAT(p)
 }
 
 
-void PlasticPrm::print (std::ostream& os) const
+PlasticMaterial::~PlasticMaterial ()
+{
+  size_t i;
+  for (i = 0; i < itgPoints.size(); i++)
+    delete itgPoints[i];
+  for (i = 0; i < resPoints.size(); i++)
+    delete resPoints[i];
+}
+
+
+void PlasticMaterial::print (std::ostream& os) const
 {
   std::cout <<"PlasticMaterial: pMAT =";
   for (size_t i = 0; i < pMAT.size(); i++)
@@ -82,17 +93,103 @@ void PlasticPrm::print (std::ostream& os) const
 }
 
 
-PlasticMaterial::PlasticMaterial (const PlasticPrm* prm, unsigned short int n)
-  : pMAT(prm->pMAT), Fp(n)
+void PlasticMaterial::initIntegration (const TimeDomain& prm)
 {
-  memset(HVc,0,sizeof(HVc));
-  memset(HVp,0,sizeof(HVc));
+#if INT_DEBUG > 0
+  std::cout <<"PlasticMaterial::initIntegration: "<< iP1 << std::endl;
+#endif
+
+  iP1 = 0;
+  iAmIntegrating = true;
+
+  if (prm.it > 0 || prm.first) return;
+
+  for (size_t i = 0; i < itgPoints.size(); i++)
+    itgPoints[i]->updateHistoryVars();
+
+#if INT_DEBUG > 0
+  std::cout <<"PlasticMaterial::initIntegration: History updated "
+	    << itgPoints.size() << std::endl;
+#endif
+}
+
+
+void PlasticMaterial::initResultPoints ()
+{
+#if INT_DEBUG > 0
+  std::cout <<"PlasticMaterial::initResultPoints: "<< iP1 << std::endl;
+#endif
+
+  iP2 = 0;
+  iAmIntegrating = false;
+
+  for (size_t i = 0; i < resPoints.size(); i++)
+    resPoints[i]->updateHistoryVars();
+
+#if INT_DEBUG > 0
+  std::cout <<"PlasticMaterial::initResultPoints: History updated "
+	    << resPoints.size() << std::endl;
+#endif
 }
 
 
 bool PlasticMaterial::evaluate (Matrix& C, SymmTensor& sigma, double&,
 				const Vec3&, const Tensor& F, const SymmTensor&,
 				char, const TimeDomain* prm) const
+{
+  if (iAmIntegrating)
+  {
+    if (!prm)
+    {
+      std::cerr <<" *** PlasticMaterial::evaluate: No time domain."<< std::endl;
+      return false;
+    }
+
+    while (itgPoints.size() <= iP1)
+      itgPoints.push_back(new PlasticPoint(this,F.dim()));
+
+    if (prm->it == 0 && !prm->first)
+      itgPoints[iP1]->Fp = F;
+
+    return itgPoints[iP1++]->evaluate(C,sigma,F,*prm);
+  }
+  else // Result evaluation
+  {
+    while (resPoints.size() <= iP2)
+      resPoints.push_back(new PlasticPoint(this,F.dim()));
+
+    if (!resPoints[iP2]->evaluate(C,sigma,F,TimeDomain(0,false)))
+      return false;
+
+    // Assume only one evaluation per increment; always update Fp
+    resPoints[iP2++]->Fp = F;
+    return true;
+  }
+}
+
+
+PlasticMaterial::PlasticPoint::PlasticPoint (const PlasticMaterial* prm,
+					     unsigned short int n)
+  : pMAT(prm->pMAT), Fp(n)
+{
+  // Initialize the history variables
+  memset(HVc,0,sizeof(HVc));
+  memset(HVp,0,sizeof(HVp));
+  Fp = HVp[0] = HVp[1] = HVp[2] = 1.0;
+}
+
+
+void PlasticMaterial::PlasticPoint::updateHistoryVars ()
+{
+  // Update history variables with values of the new converged solution
+  memcpy(HVp,HVc,sizeof(HVc));
+}
+
+
+bool PlasticMaterial::PlasticPoint::evaluate (Matrix& C,
+					      SymmTensor& sigma,
+					      const Tensor& F,
+					      const TimeDomain& prm) const
 {
   double J = F.det();
   if (J == 0.0)
@@ -101,25 +198,23 @@ bool PlasticMaterial::evaluate (Matrix& C, SymmTensor& sigma, double&,
 	      <<" Singular/zero deformation gradient\n"<< F;
     return false;
   }
-  else if (!prm)
-  {
-    std::cerr <<" *** PlasticMaterial::evaluate: No time domain."<< std::endl;
-    return false;
-  }
 
-  size_t ndim = F.dim();
+  // Restore history variables from the previous, converged configuration
+  memcpy(const_cast<double*>(HVc),HVp,sizeof(HVc));
+
+  size_t ndim = sigma.dim();
   size_t ncmp = ndim*(ndim+1)/2;
   C.resize(ncmp,ncmp);
+
   int ierr = -99;
 #ifdef USE_FTNMAT
-  // Invoke the FORTRAN routine for plasticity material models.
-  memcpy(const_cast<double*>(HVc),HVp,sizeof(HVc));
+  // Invoke the FORTRAN routine for plasticity material models
   if (ndim == 2)
-    plas2d_(INT_DEBUG,6,prm->it,prm->first,4,&pMAT.front(),J,F.ptr(),Fp.ptr(),
-	    HVc,HVc[6],HVc+7,sigma.ptr(),C.ptr(),ierr);
+    plas2d_(INT_DEBUG,6,prm.it,prm.first,4,&pMAT.front(),J,F.dim(),
+	    F.ptr(),Fp.ptr(),HVc,HVc[6],HVc+7,sigma.ptr(),C.ptr(),ierr);
   else
-    plas3d_(INT_DEBUG,6,prm->it,prm->first,6,&pMAT.front(),J,F.ptr(),Fp.ptr(),
-	    HVc,HVc[6],HVc+7,sigma.ptr(),C.ptr(),ierr);
+    plas3d_(INT_DEBUG,6,prm.it,prm.first,6,&pMAT.front(),J,
+	    F.ptr(),Fp.ptr(),HVc,HVc[6],HVc+7,sigma.ptr(),C.ptr(),ierr);
 #if INT_DEBUG > 0
   std::cout <<"PlasticMaterial::sigma =\n"<< sigma;
   std::cout <<"PlasticMaterial::C ="<< C;
@@ -130,9 +225,3 @@ bool PlasticMaterial::evaluate (Matrix& C, SymmTensor& sigma, double&,
 
   return ierr == 0;
 }
-
-
-void PlasticMaterial::updateHistoryVars ()
-{
-  memcpy(HVp,HVc,sizeof(HVc));
-}
diff --git a/Apps/FiniteDefElasticity/PlasticMaterial.h b/Apps/FiniteDefElasticity/PlasticMaterial.h
index 61d17dfc..f3dcf532 100644
--- a/Apps/FiniteDefElasticity/PlasticMaterial.h
+++ b/Apps/FiniteDefElasticity/PlasticMaterial.h
@@ -17,55 +17,84 @@
 #include "MaterialBase.h"
 #include "Tensor.h"
 
-
-/*!
-  \brief Class containing parameters of an elasto-plastic material model.
-*/
-
-class PlasticPrm : public Material
-{
-public:
-  //! \brief Constructor initializing the material parameters.
-  PlasticPrm(const RealArray&);
-  //! \brief Empty destructor.
-  virtual ~PlasticPrm() {}
-
-  //! \brief Prints out material parameters to the given output stream.
-  virtual void print(std::ostream&) const;
-
-  //! \brief Evaluates the mass density at current point.
-  virtual double getMassDensity(const Vec3&) const { return pMAT[3]; }
-
-  //! \brief Dummy method (should not be invoked).
-  virtual bool evaluate(Matrix&, SymmTensor&, double&,
-                        const Vec3&, const Tensor&, const SymmTensor&,
-                        char, const TimeDomain*) const { return false; }
-
-private:
-  RealArray pMAT; //!< Material property parameters
-
-  friend class PlasticMaterial;
-};
+class PlasticMaterial;
 
 
 /*!
-  \brief Class representing an elasto-plastic material point.
+  \brief Class representing an history-dependent elasto-plastic material model.
+  \details This class is a wrapper for the plasticity routines of FENRIS,
+  implemented in by the FORTRAN subroutine PLAS3D.
+
+  The plasticity models have history variables at the integration points that
+  need to be "remembered" from one iteration/increment to the next. This is
+  maintained inside this class as a vector of integration point objects, and
+  a global counter keeping track of which integration point we are calculating.
+  This counter is initialized in the beginning of each iteration, and then
+  incremented by each invokation of the \a evaluate method. Therefore, it is
+  paramount that this method is invoked only once per integration point per
+  iteration, and in the same order in every iteration. Otherwise it will
+  output incorrect results.
+
+  A separate "integration point" vector is dedicated for results points.
+  This is needed because the results points (for visualization, etc.) are not
+  the same as the integration points used in the tangent evaluation.
 */
 
 class PlasticMaterial : public Material
 {
+  //! \brief Class representing an elasto-plastic material point.
+  class PlasticPoint
+  {
+  public:
+    //! \brief Constructor initializing the material parameters.
+    //! \param[in] prm Pointer to actual material model object.
+    //! \param[in] n Number of space dimensions
+    PlasticPoint(const PlasticMaterial* prm, unsigned short int n);
+
+    //! \brief Updates the internal history variables after convergence.
+    void updateHistoryVars();
+
+    //! \brief Evaluates the constitutive relation at this point.
+    //! \param[out] C Constitutive matrix at current point
+    //! \param[out] sigma Stress tensor at current point
+    //! \param[in] F Deformation gradient at current point
+    //! \param[in] prm Nonlinear solution algorithm parameters
+    bool evaluate(Matrix& C, SymmTensor& sigma,
+		  const Tensor& F, const TimeDomain& prm) const;
+
+  private:
+    const RealArray& pMAT; //!< Material property parameters
+
+    double HVc[10]; //!< History variables, current configuration
+    double HVp[10]; //!< History variables, previous configuration
+
+  public:
+    Tensor Fp; //!< Deformation gradient, previous configuration
+  };
+
 public:
   //! \brief Constructor initializing the material parameters.
-  //! \param[in] prm Pointer to actual material parameters object.
-  //! \param[in] n Number of space dimensions
-  PlasticMaterial(const PlasticPrm* prm, unsigned short int n);
-  //! \brief Empty destructor.
-  virtual ~PlasticMaterial() {}
+  PlasticMaterial(const RealArray&);
+  //! \brief The destructor frees the dynamically allocated data.
+  virtual ~PlasticMaterial();
+
+  //! \brief Prints out material parameters to the given output stream.
+  virtual void print(std::ostream&) const;
+
+  //! \brief Initializes the material model for a new integration loop.
+  //! \param[in] prm Nonlinear solution algorithm parameters
+  //!
+  //! \details This method is invoked once before starting the numerical
+  //! integration over the entire spatial domain, and is used to reset the
+  //! global integration point counter, and to update the history variables.
+  virtual void initIntegration(const TimeDomain& prm);
+  //! \brief Initializes the material model for a new result point loop.
+  virtual void initResultPoints();
 
   //! \brief Evaluates the mass density at current point.
   virtual double getMassDensity(const Vec3&) const { return pMAT[3]; }
 
-  //! \brief Evaluates the constitutive relation at an integration point.
+  //! \brief Evaluates the constitutive relation at current integration point.
   //! \param[out] C Constitutive matrix at current point
   //! \param[out] sigma Stress tensor at current point
   //! \param[in] F Deformation gradient at current point
@@ -74,18 +103,17 @@ public:
                         const Vec3&, const Tensor& F, const SymmTensor&,
                         char, const TimeDomain* prm) const;
 
-  //! \brief Returns a reference to the previous deformation gradient.
-  Tensor& defGrad() { return Fp; }
-
-  //! \brief Updates the internal history variables after convergence.
-  void updateHistoryVars();
-
 private:
-  const RealArray& pMAT; //!< Material property parameters
+  friend class PlasticPoint;
 
-  double HVc[10]; //!< History variables, current configuration
-  double HVp[10]; //!< History variables, previous configuration
-  Tensor Fp;      //!< Deformation gradient, previous configuration
+  RealArray     pMAT; //!< Material property parameters
+  mutable size_t iP1; //!< Global integration point counter
+  mutable size_t iP2; //!< Global result point counter
+
+  bool iAmIntegrating; //!< Flag indicating integration or result evaluation
+
+  mutable std::vector<PlasticPoint*> itgPoints; //!< Integration point data
+  mutable std::vector<PlasticPoint*> resPoints; //!< Result point data
 };
 
 #endif
diff --git a/Apps/FiniteDefElasticity/PlasticityUL.C b/Apps/FiniteDefElasticity/PlasticityUL.C
deleted file mode 100644
index 2d164729..00000000
--- a/Apps/FiniteDefElasticity/PlasticityUL.C
+++ /dev/null
@@ -1,144 +0,0 @@
-// $Id$
-//==============================================================================
-//!
-//! \file PlasticityUL.C
-//!
-//! \date Mar 16 2011
-//!
-//! \author Knut Morten Okstad / SINTEF
-//!
-//! \brief Integrand implementations for finite deformation plasticity problems.
-//!
-//==============================================================================
-
-#include "PlasticityUL.h"
-#include "PlasticMaterial.h"
-#include "ElmMats.h"
-#include "TimeDomain.h"
-
-
-PlasticityUL::PlasticityUL (unsigned short int n, char lop)
-  : NonlinearElasticityUL(n,lop)
-{
-  iP = 0;
-
-  // Both the current and previous solutions are needed
-  primsol.resize(2);
-}
-
-
-PlasticityUL::~PlasticityUL ()
-{
-  for (size_t i = 0; i < pBuf.size(); i++)
-    delete pBuf[i];
-}
-
-
-void PlasticityUL::print (std::ostream& os) const
-{
-  std::cout <<"PlasticityUL: Finite deformation plasticity"<< std::endl;
-
-  this->NonlinearElasticityUL::print(os);
-}
-
-
-void PlasticityUL::initIntegration (const TimeDomain& prm)
-{
-  iP = 0;
-
-  if (prm.it == 0 && !prm.first)
-    for (size_t i = 0; i < pBuf.size(); i++)
-      pBuf[i]->updateHistoryVars();
-}
-
-
-bool PlasticityUL::evalInt (LocalIntegral*& elmInt,
-			    const TimeDomain& prm, double detJW,
-			    const Vector& N, const Matrix& dNdX,
-			    const Vec3& X) const
-{
-  while (pBuf.size() <= iP)
-    pBuf.push_back(new PlasticMaterial(static_cast<PlasticPrm*>(material),nsd));
-  PlasticMaterial* ptData = pBuf[iP++];
-
-  if (prm.it == 0)
-  {
-    // Evaluate the deformation gradient, Fp, at previous configuration
-    const_cast<PlasticityUL*>(this)->eV = &myMats->b[2];
-    if (!this->formDefGradient(dNdX,ptData->defGrad()))
-      return false;
-  }
-
-  // Evaluate the deformation gradient, F, and the Green-Lagrange strains, E
-  Tensor F(nsd);
-  SymmTensor E(nsd);
-  const_cast<PlasticityUL*>(this)->eV = &myMats->b[1];
-  if (!this->kinematics(dNdX,F,E))
-    return false;
-
-  bool lHaveStrains = !E.isZero();
-  if (lHaveStrains)
-  {
-    // Invert the deformation gradient ==> Fi
-    Matrix Fi(nsd,nsd);
-    Fi.fill(F.ptr());
-    double J = Fi.inverse();
-    if (J == 0.0) return false;
-
-    // Scale with J=|F| since we are integrating over current configuration
-    detJW *= J;
-
-    if (eKm || iS)
-    {
-      // Push-forward the basis function gradients to current configuration
-      dNdx.multiply(dNdX,Fi); // dNdx = dNdX * F^-1
-      // Compute the small-deformation strain-displacement matrix B from dNdx
-      if (!this->formBmatrix(dNdx)) return false;
-#if INT_DEBUG > 0
-      std::cout <<"PlasticityUL::dNdx ="<< dNdx;
-      std::cout <<"PlasticityUL::B ="<< Bmat;
-#endif
-    }
-  }
-  else if (eKm || iS)
-    // Initial state, no deformation yet
-    if (!this->formBmatrix(dNdX)) return false;
-
-  // Evaluate the constitutive relation
-  SymmTensor sigma(nsd);
-  if (eKm || eKg || iS)
-  {
-    double U = 0.0;
-    if (!ptData->evaluate(Cmat,sigma,U,X,F,E,true,&prm))
-      return false;
-  }
-
-  if (eKm)
-  {
-    // Integrate the material stiffness matrix
-    CB.multiply(Cmat,Bmat).multiply(detJW); // CB = C*B*|J|*w
-    eKm->multiply(Bmat,CB,true,false,true); // EK += B^T * CB
-  }
-
-  if (eKg && lHaveStrains)
-    // Integrate the geometric stiffness matrix
-    this->formKG(*eKg,dNdx,sigma,detJW);
-
-  if (eM)
-    // Integrate the mass matrix
-    this->formMassMatrix(*eM,N,X,detJW);
-
-  if (iS && lHaveStrains)
-  {
-    // Integrate the internal forces
-    sigma *= -detJW;
-    if (!Bmat.multiply(sigma,*iS,true,true)) // ES -= B^T*sigma
-      return false;
-  }
-
-  if (eS)
-    // Integrate the load vector due to gravitation and other body forces
-    this->formBodyForce(*eS,N,X,detJW);
-
-  return this->getIntegralResult(elmInt);
-}
diff --git a/Apps/FiniteDefElasticity/PlasticityUL.h b/Apps/FiniteDefElasticity/PlasticityUL.h
deleted file mode 100644
index 72f831e6..00000000
--- a/Apps/FiniteDefElasticity/PlasticityUL.h
+++ /dev/null
@@ -1,72 +0,0 @@
-// $Id$
-//==============================================================================
-//!
-//! \file PlasticityUL.h
-//!
-//! \date Sep 21 2010
-//!
-//! \author Knut Morten Okstad / SINTEF
-//!
-//! \brief Integrand implementations for finite deformation plasticity problems.
-//!
-//==============================================================================
-
-#ifndef _PLASTICITY_UL_H
-#define _PLASTICITY_UL_H
-
-#include "NonlinearElasticityUL.h"
-
-class PlasticMaterial;
-struct TimeDomain;
-
-
-/*!
-  \brief Class representing the integrand of the nonlinear plasticity problem.
-  \details This class implements an Updated Lagrangian formulation. It inherits
-  most of the NonlinearElasticityUL methods, but reimplements \a evalInt using
-  an internal buffer to store history variables in every integration point
-  between the iterations.
-*/
-
-class PlasticityUL : public NonlinearElasticityUL
-{
-public:
-  //! \brief The default constructor invokes the parent class constructor.
-  //! \param[in] n Number of spatial dimensions
-  //! \param[in] lop Load option (0=on initial length, 1=on updated length)
-  PlasticityUL(unsigned short int n = 3, char lop = 0);
-  //! \brief The destructor frees the dynamically allocated integration buffers.
-  virtual ~PlasticityUL();
-
-  //! \brief Prints out problem definition to the given output stream.
-  virtual void print(std::ostream& os) const;
-
-  //! \brief Initializes the integrand for a new integration loop.
-  //! \param[in] prm Nonlinear solution algorithm parameters
-  //!
-  //! \details This method is invoked once before starting the numerical
-  //! integration over the entire spatial domain, and is used to reset the
-  //! internal integration point counter, and to update the history variables.
-  virtual void initIntegration(const TimeDomain& prm);
-
-  //! \brief Evaluates the integrand at an interior point.
-  //! \param elmInt The local integral object to receive the contributions
-  //! \param[in] prm Nonlinear solution algorithm parameters
-  //! \param[in] detJW Jacobian determinant times integration point weight
-  //! \param[in] N Basis function values
-  //! \param[in] dNdX Basis function gradients
-  //! \param[in] X Cartesian coordinates of current integration point
-  virtual bool evalInt(LocalIntegral*& elmInt,
-		       const TimeDomain& prm, double detJW,
-		       const Vector& N, const Matrix& dNdX,
-		       const Vec3& X) const;
-
-  //! \brief Returns a null pointer for solution norm evaluation.
-  virtual NormBase* getNormIntegrand(AnaSol* = 0) const { return 0; }
-
-private:
-  mutable size_t                        iP;   //!< Integration point counter
-  mutable std::vector<PlasticMaterial*> pBuf; //!< Integration point buffers
-};
-
-#endif
diff --git a/Apps/FiniteDefElasticity/SIMFiniteDefEl.C b/Apps/FiniteDefElasticity/SIMFiniteDefEl.C
index 1eabb441..c4803222 100644
--- a/Apps/FiniteDefElasticity/SIMFiniteDefEl.C
+++ b/Apps/FiniteDefElasticity/SIMFiniteDefEl.C
@@ -16,7 +16,6 @@
 #include "LinearMaterial.h"
 #include "NeoHookeMaterial.h"
 #include "PlasticMaterial.h"
-#include "PlasticityUL.h"
 #include "NonlinearElasticityULMixed.h"
 #include "NonlinearElasticityULMX.h"
 #include "NeoHookeElasticity.h"
@@ -32,9 +31,6 @@ SIMFiniteDefEl2D::SIMFiniteDefEl2D (const std::vector<int>& options)
 
   switch (form)
     {
-    case SIM::PLASTICITY:
-      myProblem = new PlasticityUL(2);
-      break;
     case SIM::MIXED_QnQn1:
       nf[1] = 2; // continuous volumetric change and pressure fields
       myProblem = new NonlinearElasticityULMixed(2);
@@ -109,7 +105,7 @@ bool SIMFiniteDefEl2D::parse (char* keyWord, std::istream& is)
       RealArray pMAT;
       while ((cline = strtok(NULL," ")))
 	pMAT.push_back(atof(cline));
-      mVec.push_back(new PlasticPrm(pMAT));
+      mVec.push_back(new PlasticMaterial(pMAT));
 
       if (myPid == 0)
       {
@@ -150,9 +146,6 @@ SIMFiniteDefEl3D::SIMFiniteDefEl3D (bool checkRHS,
 
   switch (form)
     {
-    case SIM::PLASTICITY:
-      myProblem = new PlasticityUL();
-      break;
     case SIM::MIXED_QnQn1:
       nf[1] = 2; // continuous volumetric change and pressure fields
       myProblem = new NonlinearElasticityULMixed();
@@ -233,7 +226,7 @@ bool SIMFiniteDefEl3D::parse (char* keyWord, std::istream& is)
       RealArray pMAT;
       while ((cline = strtok(NULL," ")))
 	pMAT.push_back(atof(cline));
-      mVec.push_back(new PlasticPrm(pMAT));
+      mVec.push_back(new PlasticMaterial(pMAT));
 
       if (myPid == 0)
       {
diff --git a/Apps/FiniteDefElasticity/SIMFiniteDefEl.h b/Apps/FiniteDefElasticity/SIMFiniteDefEl.h
index a1e97dcd..08f48369 100644
--- a/Apps/FiniteDefElasticity/SIMFiniteDefEl.h
+++ b/Apps/FiniteDefElasticity/SIMFiniteDefEl.h
@@ -28,8 +28,7 @@ namespace SIM
     NEOHOOKE_IV      = 5,
     UPDATED_LAGRANGE = 6,
     MIXED_QnPn1      = 7,
-    MIXED_QnQn1      = 8,
-    PLASTICITY       = 9
+    MIXED_QnQn1      = 8
   };
 };
 
diff --git a/Apps/FiniteDefElasticity/main_NonLinEl.C b/Apps/FiniteDefElasticity/main_NonLinEl.C
index d101a940..07dd76e8 100644
--- a/Apps/FiniteDefElasticity/main_NonLinEl.C
+++ b/Apps/FiniteDefElasticity/main_NonLinEl.C
@@ -51,7 +51,6 @@
   \arg -2Dpstrain : Use two-parametric simulation driver (plane strain)
   \arg -Tensor : Use tensorial total Lagrangian formulation (slow...)
   \arg -UL : Use updated Lagrangian formulation with nonlinear material
-  \arg -PL : Use placticity solution driver
   \arg -MX \a pord : Mixed formulation with internal discontinuous pressure
   \arg -mixed : Mixed formulation with continuous pressure and volumetric change
   \arg -lag : Use Lagrangian basis functions instead of splines/NURBS
@@ -146,8 +145,6 @@ int main (int argc, char** argv)
       twoD = SIMLinEl2D::planeStrain = true;
     else if (!strncmp(argv[i],"-2D",3))
       twoD = true;
-    else if (!strcmp(argv[i],"-PL"))
-      form = SIM::PLASTICITY;
     else if (!strcmp(argv[i],"-UL"))
     {
       if (form < SIM::UPDATED_LAGRANGE)
@@ -183,8 +180,8 @@ int main (int argc, char** argv)
 	      <<" [-lag] [-spec] [-2D[pstrain]] [-UL] [-MX [<p>]|-mixed]"
 	      <<" [-nGauss <n>]\n       [-vtf <format>] [-nviz <nviz>]"
 	      <<" [-nu <nu>] [-nv <nv>] [-nw <nw>]\n      "
-	      <<" [-saveInc <dtSave>] [-dumpInc <dtDump> [raw]]\n      "
-	      <<" [-ignore <p1> <p2> ...] [-fixDup] [-checkRHS] [-check]\n";
+	      <<" [-saveInc <dtSave>] [-skip2nd] [-dumpInc <dtDump> [raw]]\n   "
+	      <<"    [-ignore <p1> <p2> ...] [-fixDup] [-checkRHS] [-check]\n";
     return 0;
   }