diff --git a/Apps/LinearElasticity/main_LinEl3D.C b/Apps/LinearElasticity/main_LinEl3D.C
index 021d7cf2..f6719fa0 100644
--- a/Apps/LinearElasticity/main_LinEl3D.C
+++ b/Apps/LinearElasticity/main_LinEl3D.C
@@ -233,9 +233,9 @@ int main (int argc, char** argv)
 
   model->setQuadratureRule(nGauss);
 
-  Matrix eNorm;
+  Matrix eNorm, ssol;
   Vector gNorm, load;
-  Vectors displ(1);
+  Vectors displ(1), projs;
   std::vector<Mode> modes;
   std::vector<Mode>::const_iterator it;
 
@@ -255,20 +255,40 @@ int main (int argc, char** argv)
     if (!model->solveSystem(displ.front(),1))
       return 3;
 
-    // Evaluate solution norms
     model->setMode(SIM::RECOVERY);
-    if (!model->solutionNorms(displ,eNorm,gNorm))
+    if (SIMbase::discretization == SIMbase::Spline)
+    {
+      // Project the FE stresses onto the splines basis
+      if (!model->project(ssol,displ.front(),SIMbase::GLOBAL))
+	return 4;
+      else
+	projs.push_back(ssol);
+    }
+
+    // Integrate solution norms and errors
+    if (!model->solutionNorms(displ,projs,eNorm,gNorm))
       return 4;
 
     if (linalg.myPid == 0)
     {
       std::cout <<"Energy norm |u^h| = a(u^h,u^h)^0.5   : "<< gNorm(1);
-      std::cout <<"\nExternal energy ((f,u^h)+(t,u^h)^0.5 : "<< gNorm(2);
-      if (gNorm.size() > 2)
-	std::cout <<"\nExact norm  |u|   = a(u,u)^0.5       : "<< gNorm(3);
-      if (gNorm.size() > 3)
-	std::cout <<"\nExact error a(e,e)^0.5, e=u-u^h      : "<< gNorm(4)
+      std::cout	<<"\nExternal energy ((f,u^h)+(t,u^h)^0.5 : "<< gNorm(2);
+      if (model->haveAnaSol() && gNorm.size() >= 4)
+	std::cout <<"\nExact norm  |u|   = a(u,u)^0.5       : "<< gNorm(3)
+		  <<"\nExact error a(e,e)^0.5, e=u-u^h      : "<< gNorm(4)
 		  <<"\nExact relative error (%) : "<< gNorm(4)/gNorm(3)*100.0;
+      size_t j = model->haveAnaSol() ? 5 : 3;
+      for (size_t i = 0; i < projs.size() && j < gNorm.size(); i++)
+      {
+        std::cout <<"\nEnergy norm |u^r| = a(u^r,u^r)^0.5   : "<< gNorm(j++);
+	std::cout <<"\nError norm a(e,e)^0.5, e=u^r-u^h     : "<< gNorm(j++);
+	std::cout <<"\n- relative error (% of |u^r|) : "
+		  << gNorm(j-1)/gNorm(j-2)*100.0;
+	if (model->haveAnaSol() && j++ <= gNorm.size())
+	  std::cout <<"\nExact error a(e,e)^0.5, e=u-u^r      : "<< gNorm(j-1)
+		    <<"\n- relative error (% of |u|)   : "
+		    << gNorm(j-1)/gNorm(3)*100.0;
+      }
       std::cout << std::endl;
     }
 
@@ -349,6 +369,10 @@ int main (int argc, char** argv)
     if (!model->writeGlvS(displ.front(),n,iStep,nBlock))
       return 10;
 
+    // Write projected solution fields to VTF-file
+    if (!model->writeGlvP(projs.front(),n,iStep,nBlock))
+      return 10;
+
     // Write eigenmodes
     for (it = modes.begin(); it != modes.end(); it++)
       if (!model->writeGlvM(*it, iop==3 || iop==4 || iop==6, n, nBlock))
@@ -359,7 +383,7 @@ int main (int argc, char** argv)
       if (!model->writeGlvN(eNorm,iStep,nBlock))
 	return 12;
 
-    model->writeGlvStep(1,0,1);
+    model->writeGlvStep(1);
     model->closeGlv();
   }
 
diff --git a/src/ASM/ElmNorm.h b/src/ASM/ElmNorm.h
index ee39883c..bc0aaca6 100644
--- a/src/ASM/ElmNorm.h
+++ b/src/ASM/ElmNorm.h
@@ -28,7 +28,7 @@ class ElmNorm : public LocalIntegral
 {
 public:
   //! \brief The constructor assigns the internal pointer.
-  ElmNorm(double* p) : ptr(p) {}
+  ElmNorm(double* p, size_t n) : ptr(p), nnv(n) {}
   //! \brief Empty destructor.
   virtual ~ElmNorm() {}
 
@@ -37,8 +37,12 @@ public:
   //! \brief Indexing operator for referencing.
   const double& operator[](size_t i) const { return ptr[i]; }
 
+  //! \brief Returns the number of norm values.
+  size_t size() const { return nnv; }
+
 private:
   double* ptr; //!< Pointer to the actual norm values
+  size_t  nnv; //!< Number of norm values
 };
 
 #endif
diff --git a/src/ASM/GlbNorm.C b/src/ASM/GlbNorm.C
index 4a4ea34c..de2af75d 100644
--- a/src/ASM/GlbNorm.C
+++ b/src/ASM/GlbNorm.C
@@ -35,9 +35,10 @@ bool GlbNorm::assemble (const LocalIntegral* elmObj, int elmId)
   if (!ptr) return false;
 
   ElmNorm& elVals = *const_cast<ElmNorm*>(ptr);
-  for (size_t i = 0; i < myVals.size(); i++)
+  for (size_t i = 0; i < elVals.size(); i++)
   {
-    myVals[i] += elVals[i];
+    if (i < myVals.size())
+      myVals[i] += elVals[i];
     this->applyFinalOp(elVals[i]);
   }
 
diff --git a/src/ASM/IntegrandBase.C b/src/ASM/IntegrandBase.C
index f21233d3..25395b63 100644
--- a/src/ASM/IntegrandBase.C
+++ b/src/ASM/IntegrandBase.C
@@ -185,8 +185,8 @@ ElmNorm& NormBase::getElmNormBuffer (LocalIntegral*& elmInt, const size_t nn)
 
   static double* data = 0;
   if (!data) data = new double[nn];
-  static ElmNorm buf(data);
   memset(data,0,nn*sizeof(double));
+  static ElmNorm buf(data,nn);
   elmInt = &buf;
   return buf;
 }
diff --git a/src/ASM/IntegrandBase.h b/src/ASM/IntegrandBase.h
index 85401067..112d09d8 100644
--- a/src/ASM/IntegrandBase.h
+++ b/src/ASM/IntegrandBase.h
@@ -221,6 +221,9 @@ public:
   //! \brief Returns the number of field components.
   virtual size_t getNoFields() const { return 0; }
 
+  //! \brief Accesses a projected secondary solution vector of current patch.
+  virtual Vector& getProjection(size_t) { static Vector dummy; return dummy; }
+
 protected:
   //! \brief Returns the element norm object to use in the integration.
   //! \param elmInt The local integral object to receive norm contributions
diff --git a/src/Integrands/Elasticity.C b/src/Integrands/Elasticity.C
index 6f7f1db3..2c344181 100644
--- a/src/Integrands/Elasticity.C
+++ b/src/Integrands/Elasticity.C
@@ -550,9 +550,49 @@ NormBase* Elasticity::getNormIntegrand (AnaSol* asol) const
 }
 
 
+ElasticityNorm::ElasticityNorm (Elasticity& p, STensorFunc* a)
+  : NormBase(p), anasol(a)
+{
+  nrcmp = myProblem.getNoFields(2);
+}
+
+
 size_t ElasticityNorm::getNoFields () const
 {
-  return anasol ? 4 : 2;
+  size_t nf = anasol ? 4 : 2;
+  for (size_t i = 0; i < prjsol.size(); i++)
+    if (!prjsol.empty())
+       nf += anasol ? 3 : 2;
+
+  return nf;
+}
+
+
+Vector& ElasticityNorm::getProjection (size_t i)
+{
+  if (prjsol.size() < i)
+  {
+    prjsol.resize(i);
+    mySols.resize(i);
+  }
+
+  return prjsol[i-1];
+}
+
+
+bool ElasticityNorm::initElement (const std::vector<int>& MNPC)
+{
+  // Extract projected solution vectors for this element
+  int ierr = 0;
+  for (size_t i = 0; i < mySols.size() && ierr == 0; i++)
+    if (!prjsol[i].empty())
+      ierr = utl::gather(MNPC,nrcmp,prjsol[i],mySols[i]);
+
+  if (ierr == 0) return myProblem.initElement(MNPC);
+
+  std::cerr <<" *** ElasticityNorm::initElement: Detected "
+            << ierr <<" node numbers out of range."<< std::endl;
+  return false;
 }
 
 
@@ -560,18 +600,19 @@ bool ElasticityNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
 			      const Vec3& X) const
 {
   Elasticity& problem = static_cast<Elasticity&>(myProblem);
-  ElmNorm& pnorm = NormBase::getElmNormBuffer(elmInt);
+  ElmNorm& pnorm = NormBase::getElmNormBuffer(elmInt,this->getNoFields());
 
   // Evaluate the inverse constitutive matrix at this point
   Matrix Cinv;
   if (!problem.formCinverse(Cinv,X)) return false;
 
   // Evaluate the finite element stress field
-  Vector sigmah, sigma;
+  Vector sigmah, sigma, error;
   if (!problem.evalSol(sigmah,fe.N,fe.dNdX,X))
     return false;
-  else if (sigmah.size() == 4 && Cinv.rows() == 3)
-    sigmah.erase(sigmah.begin()+2); // Remove the sigma_zz if plane strain
+
+  bool planeStrain = sigmah.size() == 4 && Cinv.rows() == 3;
+  if (planeStrain) sigmah.erase(sigmah.begin()+2); // Remove the sigma_zz
 
   double detJW = fe.detJxW;
   if (problem.isAxiSymmetric())
@@ -590,6 +631,7 @@ bool ElasticityNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
     pnorm[1] += f*u*detJW;
   }
 
+  size_t ip = 2;
   if (anasol)
   {
     // Evaluate the analytical stress field
@@ -598,12 +640,36 @@ bool ElasticityNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
       sigma.erase(sigma.begin()+2); // Remove the sigma_zz if plane strain
 
     // Integrate the energy norm a(u,u)
-    pnorm[2] += sigma.dot(Cinv*sigma)*detJW;
+    pnorm[ip++] += sigma.dot(Cinv*sigma)*detJW;
     // Integrate the error in energy norm a(u-u^h,u-u^h)
-    sigma -= sigmah;
-    pnorm[3] += sigma.dot(Cinv*sigma)*detJW;
+    error = sigma - sigmah;
+    pnorm[ip++] += error.dot(Cinv*error)*detJW;
   }
 
+  size_t i, j, k;
+  for (i = 0; i < mySols.size(); i++)
+    if (!mySols[i].empty())
+    {
+      // Evaluate projected stress field
+      Vector sigmar(sigmah.size());
+      for (j = k = 0; j < nrcmp && k < sigmar.size(); j++)
+	if (!planeStrain || j != 2)
+	  sigmar[k++] = mySols[i].dot(fe.N,j,nrcmp);
+
+      // Integrate the energy norm a(u^r,u^r)
+      pnorm[ip++] += sigmar.dot(Cinv*sigmar)*detJW;
+      // Integrate the error in energy norm a(u^r-u^h,u^r-u^h)
+      error = sigmar - sigmah;
+      pnorm[ip++] += error.dot(Cinv*error)*detJW;
+
+      if (anasol)
+      {
+	// Integrate the error in the projected solution a(u-u^r,u-u^r)
+	error = sigma - sigmar;
+	pnorm[ip++] += error.dot(Cinv*error)*detJW;
+      }
+    }
+
   return true;
 }
 
diff --git a/src/Integrands/Elasticity.h b/src/Integrands/Elasticity.h
index e111a017..e4be9077 100644
--- a/src/Integrands/Elasticity.h
+++ b/src/Integrands/Elasticity.h
@@ -225,7 +225,7 @@ public:
   //! \brief The only constructor initializes its data members.
   //! \param[in] p The linear elasticity problem to evaluate norms for
   //! \param[in] a The analytical stress field (optional)
-  ElasticityNorm(Elasticity& p, STensorFunc* a = 0) : NormBase(p), anasol(a) {}
+  ElasticityNorm(Elasticity& p, STensorFunc* a = 0);
   //! \brief Empty destructor.
   virtual ~ElasticityNorm() {}
 
@@ -234,6 +234,14 @@ public:
   //! \brief Returns a 1-based index of the external energy norm.
   virtual size_t indExt() const { return 2; }
 
+  //! \brief Initializes current element for numerical integration.
+  //! \param[in] MNPC Matrix of nodal point correspondance for current element
+  virtual bool initElement(const std::vector<int>& MNPC);
+  //! \brief Initializes current element for numerical integration.
+  //! \param[in] MNPC Matrix of nodal point correspondance for current element
+  virtual bool initElement(const std::vector<int>& MNPC, const Vec3&, size_t)
+  { return this->initElement(MNPC); }
+
   //! \brief Evaluates the integrand at an interior point.
   //! \param elmInt The local integral object to receive the contributions
   //! \param[in] fe Finite element data of current integration point
@@ -252,8 +260,16 @@ public:
   //! \brief Returns the number of norm quantities.
   virtual size_t getNoFields() const;
 
+  //! \brief Accesses a projected secondary solution vector of current patch.
+  virtual Vector& getProjection(size_t i);
+
 private:
   STensorFunc* anasol; //!< Analytical stress field
+
+  Vectors prjsol; //!< Projected secondary solution vectors for current patch
+  Vectors mySols; //!< Local element vectors
+
+  size_t nrcmp; //!< Number of projected solution components
 };
 
 #endif
diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C
index a3645a9d..d83c3357 100644
--- a/src/SIM/SIMbase.C
+++ b/src/SIM/SIMbase.C
@@ -697,7 +697,8 @@ double SIMbase::solutionNorms (const Vector& x, double* inf,
 }
 
 
-bool SIMbase::solutionNorms (const TimeDomain& time, const Vectors& psol,
+bool SIMbase::solutionNorms (const TimeDomain& time,
+			     const Vectors& psol, const Vectors& ssol,
 			     Vector& gNorm, Matrix* eNorm)
 {
   NormBase* norm = myProblem->getNormIntegrand(mySol);
@@ -714,6 +715,12 @@ bool SIMbase::solutionNorms (const TimeDomain& time, const Vectors& psol,
   myProblem->initIntegration(time);
   const Vector& primsol = psol.front();
 
+  size_t i, j, k;
+  for (i = 0; i < ssol.size(); i++)
+    if (!ssol[i].empty())
+      norm->getProjection(i+1);
+
+  size_t nCmp = ssol.empty() ? 0 : ssol.front().size() / mySam->getNoNodes();
   size_t nNorms = norm->getNoFields();
   gNorm.resize(nNorms,true);
 
@@ -724,21 +731,24 @@ bool SIMbase::solutionNorms (const TimeDomain& time, const Vectors& psol,
   {
     eNorm->resize(nNorms,mySam->getNoElms(),true);
     elementNorms.reserve(eNorm->cols());
-    for (size_t i = 0; i < eNorm->cols(); i++)
-      elementNorms.push_back(new ElmNorm(eNorm->ptr(i)));
+    for (i = 0; i < eNorm->cols(); i++)
+      elementNorms.push_back(new ElmNorm(eNorm->ptr(i),nNorms));
   }
 
   // Loop over the different material regions, integrating solution norm terms
   // for the patch domain associated with each material
   bool ok = true;
-  size_t i, j = 0, lp = 0;
-  for (i = 0; i < myProps.size() && ok; i++)
+  size_t lp = 0;
+  for (i = j = 0; i < myProps.size() && ok; i++)
     if (myProps[i].pcode == Property::MATERIAL)
       if ((j = myProps[i].patch) < 1 || j > myModel.size())
 	ok = false;
       else if (this->initMaterial(myProps[i].pindx))
       {
 	myModel[j-1]->extractNodeVec(primsol,myProblem->getSolution());
+	for (k = 0; k < ssol.size(); k++)
+	  if (!ssol[k].empty())
+	    myModel[j-1]->extractNodeVec(ssol[k],norm->getProjection(k+1),nCmp);
 	ok = myModel[j-1]->integrate(*norm,globalNorm,time,elementNorms);
 	lp = j;
       }
@@ -751,8 +761,11 @@ bool SIMbase::solutionNorms (const TimeDomain& time, const Vectors& psol,
     for (i = 0; i < myModel.size() && ok; i++)
     {
       myModel[i]->extractNodeVec(primsol,myProblem->getSolution());
+      for (k = 0; k < ssol.size(); k++)
+	if (!ssol[k].empty())
+	  myModel[i]->extractNodeVec(ssol[k],norm->getProjection(k+1),nCmp);
       ok = myModel[i]->integrate(*norm,globalNorm,time,elementNorms);
-      lp = j;
+      lp = i+1;
     }
 
   // Integrate norm contributions due to Neumann boundary conditions, if any.
@@ -900,13 +913,15 @@ bool SIMbase::writeGlv (const char* inpFile, const int* nViz, int format)
 {
   if (myVtf) return false;
 
+#if HAS_VTFAPI == 2
+  const char* ext = ".vtfx";
+#else
+  const char* ext = ".vtf";
+#endif
+
   // Open a new VTF-file
   char* vtfName = new char[strlen(inpFile)+10];
   strtok(strcpy(vtfName,inpFile),".");
-  const char* ext = ".vtf";
-#if HAS_VTFAPI == 2
-  ext = ".vtfx";
-#endif
   if (nProc > 1)
     sprintf(vtfName+strlen(vtfName),"_p%04d%s",myPid,ext);
   else
@@ -1220,6 +1235,52 @@ bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
 }
 
 
+bool SIMbase::writeGlvP (const Vector& ssol, const int* nViz,
+			 int iStep, int& nBlock, double time,
+			 int idBlock, const char* prefix)
+{
+  if (ssol.empty())
+    return true;
+  else if (!myVtf)
+    return false;
+
+  Matrix field;
+  Vector lovec;
+  const size_t nf = myProblem->getNoFields(2);
+  std::vector<int> sID[nf];
+
+  size_t i, j;
+  int geomID = 0;
+  for (i = 0; i < myModel.size(); i++)
+  {
+    if (myModel[i]->empty()) continue; // skip empty patches
+
+    if (msgLevel > 1)
+      std::cout <<"Writing projected solution for patch "<< i+1 << std::endl;
+
+    // Evaluate the solution variables at the visualization points
+    myModel[i]->extractNodeVec(ssol,lovec,nf);
+    if (!myModel[i]->evalSolution(field,lovec,nViz))
+      return false;
+
+    // Write out to VTF-file as scalar fields
+    geomID++;
+    for (j = 0; j < field.rows(); j++)
+      if (!myVtf->writeNres(field.getRow(1+j),++nBlock,geomID))
+	return false;
+      else
+	sID[j].push_back(nBlock);
+  }
+
+  // Write result block identifications
+  for (j = 0; j < nf && !sID[j].empty(); j++)
+    if (!myVtf->writeSblk(sID[j],myProblem->getField2Name(j,prefix),
+			  ++idBlock,iStep)) return false;
+
+  return true;
+}
+
+
 bool SIMbase::writeGlvStep (int iStep, double value, int itype)
 {
   if (itype == 0)
@@ -1280,7 +1341,7 @@ bool SIMbase::writeGlvN (const Matrix& norms, int iStep, int& nBlock)
 
   Matrix field;
   int geomID = 0;
-  std::vector<int> sID[3];
+  std::vector<int> sID[10];
 
   size_t i, j, k;
   for (i = 0; i < myModel.size(); i++)
@@ -1293,25 +1354,39 @@ bool SIMbase::writeGlvN (const Matrix& norms, int iStep, int& nBlock)
     geomID++;
     myModel[i]->extractElmRes(norms,field);
 
-    for (j = k = 0; j < field.rows() && j < 4; j++)
-      if (field.rows()%2 || j != 1) // Skip the external norms (always zero)
+    for (j = k = 0; j < field.rows() && j < 10; j++)
+      if (j != 1) // Skip the external norms (always zero)
 	if (!myVtf->writeEres(field.getRow(1+j),++nBlock,geomID))
 	  return false;
 	else
 	  sID[k++].push_back(nBlock);
   }
 
-  const char* label[3] = {
+  const char* label[9] = {
     "a(u^h,u^h)^0.5",
     "a(u,u)^0.5",
-    "a(e,e)^0.5, e=u-u^h"
+    "a(e,e)^0.5, e=u-u^h",
+    "a(u^r,u^r)^0.5",
+    "a(e,e)^0.5, e=u^r-u^h",
+    "a(e,e)^0.5, e=u-u^r",
+    "a(u^rr,u^rr)^0.5",
+    "a(e,e)^0.5, e=u^rr-u^h",
+    "a(e,e)^0.5, e=u-u^rr"
   };
 
-  int idBlock = 100;
-  for (j = 0; j < 3; j++)
-    if (!sID[j].empty())
-      if (!myVtf->writeSblk(sID[j],label[j],++idBlock,iStep))
-	return false;
+  int idBlock = 200;
+  for (j = k = 0; !sID[j].empty(); j++, k++)
+  {
+    if (!mySol)
+    {
+      if (k == 1)
+	k = 3;
+      else if (k%3 == 2)
+        k ++;
+    }
+    if (!myVtf->writeSblk(sID[j],label[k],++idBlock,iStep))
+      return false;
+  }
 
   return true;
 }
@@ -1375,9 +1450,9 @@ bool SIMbase::dumpGeometry (std::ostream& os) const
 
 bool SIMbase::dumpBasis (std::ostream& os, int basis, size_t patch) const
 {
-  size_t start = patch?patch-1:0;
-  size_t end = patch?start+1:myModel.size();
-  for (size_t i = start; i < end; i++)
+  size_t start = patch ? patch-1 : 0;
+  size_t end = patch ? start : myModel.size();
+  for (size_t i = start; i < end && i < myModel.size(); i++)
     if (!myModel[i]->empty())
       if (!myModel[i]->write(os,basis))
 	return false;
@@ -1438,7 +1513,7 @@ bool SIMbase::dumpSolution (const Vector& psol, std::ostream& os) const
 
 
 bool SIMbase::dumpResults (const Vector& psol, double time, std::ostream& os,
-			   bool formatted, std::streamsize outputPrecision) const
+			   bool formatted, std::streamsize precision) const
 {
   if (psol.empty() || myPoints.empty())
     return true;
@@ -1491,8 +1566,8 @@ bool SIMbase::dumpResults (const Vector& psol, double time, std::ostream& os,
 	return false;
 
     // Formatted output, use scientific notation with fixed field width
-    std::streamsize flWidth = 8 + outputPrecision;
-    std::streamsize oldPrec = os.precision(outputPrecision);
+    std::streamsize flWidth = 8 + precision;
+    std::streamsize oldPrec = os.precision(precision);
     std::ios::fmtflags oldF = os.flags(std::ios::scientific | std::ios::right);
     for (j = 0; j < points.size(); j++)
     {
@@ -1537,21 +1612,81 @@ bool SIMbase::dumpResults (const Vector& psol, double time, std::ostream& os,
 }
 
 
-bool SIMbase::project (Vector& psol)
+bool SIMbase::project (Matrix& ssol, const Vector& psol,
+		       ProjectionMethod pMethod) const
 {
-  Matrix values;
-  Vector ssol;
-  ssol.resize(psol.size()*myProblem->getNoFields());
-  for (size_t i = 0; i < myModel.size(); i++)
-  {
-    myModel[i]->extractNodeVec(psol,myProblem->getSolution());
-    if (!myModel[i]->evalSolution(values,*myProblem))
-      return false;
-    if (!myModel[i]->injectNodeVec(values,ssol,values.rows()))
-      return false;
-  }
-  psol = ssol;
+  PROFILE1("Solution projection");
 
+  if (msgLevel > 1)
+    std::cout <<"\nProjecting secondary solution ...\n"<< std::endl;
+
+  ssol.resize(0,0);
+
+  size_t i, j, n;
+  size_t ngNodes = mySam->getNoNodes();
+
+  Matrix values;
+  Vector count(myModel.size() > 1 ? ngNodes : 0);
+
+  for (i = 0; i < myModel.size(); i++)
+  {
+    if (myModel[i]->empty()) continue; // skip empty patches
+
+    // Extract the primary solution control point values for this patch
+    myModel[i]->extractNodeVec(psol,myProblem->getSolution());
+
+    // Project the secondary solution and retrieve control point values
+    switch (pMethod) {
+    case GLOBAL:
+      if (!myModel[i]->evalSolution(values,*myProblem))
+	return false;
+      break;
+
+    case LOCAL:
+      // Annette, add your local projection stuff here...
+
+    default:
+      std::cerr <<" *** SIMbase::project: Projection method "<< pMethod
+		<<" not implemented."<< std::endl;
+      return false;
+    }
+
+    size_t nComps = values.rows();
+    size_t nNodes = myModel[i]->getNoNodes();
+    if (ssol.empty())
+      ssol.resize(nComps,ngNodes);
+
+    // Nodal averaging for nodes referred to by two or more patches
+    // (these are typically the interface nodes)
+    for (n = 1; n <= nNodes; n++)
+      if (count.empty())
+	ssol.fillColumn(myModel[i]->getNodeID(n),values.getColumn(n));
+      else
+      {
+	int inod = myModel[i]->getNodeID(n);
+	for (j = 1; j <= nComps; j++)
+	  ssol(j,inod) += values(j,n);
+	count(inod) ++;
+      }
+  }
+
+  // Divide through by count(n) to get the nodal average at the interface nodes
+  for (n = 1; n <= count.size(); n++)
+    if (count(n) > 1.0)
+      for (j = 1; j <= ssol.rows(); j++)
+	ssol(j,n) /= count(n);
+
+  return true;
+}
+
+
+bool SIMbase::project (Vector& sol) const
+{
+  Matrix secsol;
+  if (!this->project(secsol,sol))
+    return false;
+
+  sol = secsol;
   return true;
 }
 
diff --git a/src/SIM/SIMbase.h b/src/SIM/SIMbase.h
index 54221f5b..4e9e1363 100644
--- a/src/SIM/SIMbase.h
+++ b/src/SIM/SIMbase.h
@@ -90,16 +90,16 @@ public:
   //! \param is The file stream to read from
   virtual bool parse(char* keyWord, std::istream& is);
 
+  //! \brief Reads patches from a given stream.
+  //! \param[in] isp The stream to read from
+  virtual void readPatches(std::istream& isp) {}
+
   //! \brief Performs some pre-processing tasks on the FE model.
   //! \param[in] ignoredPatches Indices of patches to ignore in the analysis
   //! \param[in] fixDup Merge duplicated FE nodes on patch interfaces?
   bool preprocess(const std::vector<int>& ignoredPatches = std::vector<int>(),
 		  bool fixDup = false);
 
-  //! \brief Read patches from given stream
-  //! \param[in] isp The stream to read from
-  virtual void readPatches(std::istream& isp) {};
-
   //! \brief Allocates the system matrices of the FE problem to be solved.
   //! \param[in] mType The matrix format to use
   //! \param[in] nMats Number of system matrices
@@ -215,24 +215,47 @@ public:
 
   //! \brief Integrates some solution norm quantities.
   //! \details If an analytical solution is provided, norms of the exact
-  //! error in the solution are computed as well.
+  //! error in the solution are computed as well. If projected secondary
+  //! solutions are provided (i.e., \a ssol is not empty), norms of the
+  //! difference between these solutions and the directly evaluated secondary
+  //! solution are computed as well.
+  //! \param[in] time Parameters for nonlinear/time-dependent simulations.
+  //! \param[in] psol Primary solution vectors
+  //! \param[in] ssol Secondary solution vectors
+  //! \param[out] gNorm Global norm quantities
+  //! \param[out] eNorm Element-wise norm quantities
+  bool solutionNorms(const TimeDomain& time,
+		     const Vectors& psol, const Vectors& ssol,
+		     Vector& gNorm, Matrix* eNorm = 0);
+  //! \brief Integrates some solution norm quantities.
   //! \param[in] time Parameters for nonlinear/time-dependent simulations.
   //! \param[in] psol Primary solution vectors
   //! \param[out] gNorm Global norm quantities
   //! \param[out] eNorm Element-wise norm quantities
-  virtual bool solutionNorms(const TimeDomain& time, const Vectors& psol,
-			     Vector& gNorm, Matrix* eNorm = 0);
-
+  //!
+  //! \details Use this version if no projected solutions are needed/available.
+  bool solutionNorms(const TimeDomain& time, const Vectors& psol,
+		     Vector& gNorm, Matrix* eNorm = 0)
+  { return this->solutionNorms(time,psol,Vectors(),gNorm,eNorm); }
   //! \brief Integrates some solution norm quantities.
-  //! \details If an analytical solution is provided, norms of the exact
-  //! error in the solution are computed as well.
   //! \param[in] psol Primary solution vectors
+  //! \param[in] ssol Secondary solution vectors
   //! \param[out] eNorm Element-wise norm quantities
   //! \param[out] gNorm Global norm quantities
   //!
   //! \details Use this version for linear/stationary problems only.
-  virtual bool solutionNorms(const Vectors& psol, Matrix& eNorm, Vector& gNorm)
-  { return this->solutionNorms(TimeDomain(),psol,gNorm,&eNorm); }
+  bool solutionNorms(const Vectors& psol, const Vectors& ssol,
+		     Matrix& eNorm, Vector& gNorm)
+  { return this->solutionNorms(TimeDomain(),psol,ssol,gNorm,&eNorm); }
+  //! \brief Integrates some solution norm quantities.
+  //! \param[in] psol Primary solution vectors
+  //! \param[out] eNorm Element-wise norm quantities
+  //! \param[out] gNorm Global norm quantities
+  //!
+  //! \details Use this version for linear/stationary problems,
+  //! and when no projected solutions are needed/available.
+  bool solutionNorms(const Vectors& psol, Matrix& eNorm, Vector& gNorm)
+  { return this->solutionNorms(TimeDomain(),psol,Vectors(),gNorm,&eNorm); }
 
   //! \brief Computes the total reaction forces in the model.
   //! \param[out] RF Reaction force in each spatial direction + energy
@@ -251,13 +274,21 @@ public:
 		   int nev, int ncv, int iop, double shift,
 		   size_t iA = 0, size_t iB = 1);
 
-  //! \brief Projects the secondary solution associated with \a psol.
-  //! \param psol Control point values of primary(in)/secondary(out) solution
+  //! \brief Enum defining the available projection methods.
+  enum ProjectionMethod { GLOBAL, LOCAL };
+  //! \brief Projects the secondary solution associated with a primary solution.
+  //! \param[out] ssol Control point values of the secondary solution
+  //! \param[in] psol Control point values of the primary solution
+  //! \param[in] pMethod Projection method to use
   //!
   //! \details The secondary solution, defined through the Integrand object,
   //! corresponding to the primary solution \a psol is projected onto the
   //! spline basis to obtain the control point values of the secondary solution.
-  bool project(Vector& psol);
+  bool project(Matrix& ssol, const Vector& psol,
+	       ProjectionMethod pMethod = GLOBAL) const;
+  //! \brief Projects the secondary solution associated with a primary solution.
+  //! \param sol Control point values of primary(in)/secondary(out) solution
+  bool project(Vector& sol) const;
 
   //! \brief Evaluates the secondary solution field for specified patch.
   //! \param[out] field Control point values of the secondary solution field
@@ -324,6 +355,18 @@ public:
 			 const char* pvecName = 0, int idBlock = 10,
 			 int psolComps = 0);
 
+  //! \brief Writes projected solutions for a given time step to the VTF-file.
+  //! \param[in] ssol Secondary solution vector (control point values)
+  //! \param[in] nViz Number of visualization points over each knot-span
+  //! \param[in] iStep Load/time step identifier
+  //! \param nBlock Running result block counter
+  //! \param[in] time Load/time step parameter
+  //! \param[in] idBlock Starting value of result block numbering
+  //! \param[in] prefix Common prefix for the field components
+  bool writeGlvP(const Vector& ssol, const int* nViz, int iStep,
+		 int& nBlock, double time = 0.0, int idBlock = 100,
+		 const char* pvecName = "Global projected");
+
   //! \brief Writes a mode shape and associated eigenvalue to the VTF-file.
   //! \details The eigenvalue is used as a label on the step state info
   //! that is associated with the eigenvector.
@@ -344,7 +387,7 @@ public:
   //! \brief Writes time/load step info to the VTF-file.
   //! \param[in] iStep Load/time step identifier
   //! \param[in] value Time or load parameter of the step
-  //! \param[in] itype Type indentifier of the step
+  //! \param[in] itype Type identifier of the step
   bool writeGlvStep(int iStep, double value = 0.0, int itype = 0);
 
   //! \brief Closes the current VTF-file.
@@ -357,7 +400,7 @@ public:
   //! \param os Output stream to write the spline data to
   //! \param[in] basis Which basis to dump for mixed methods (0 = geometry)
   //! \param[in] patch Which patch to dump for (0 = all)
-  bool dumpBasis(std::ostream& os, int basis = 0, size_t patch=0) const;
+  bool dumpBasis(std::ostream& os, int basis = 0, size_t patch = 0) const;
   //! \brief Dumps the entire solution in ASCII format.
   //! \param[in] psol Primary solution vector to derive other quantities from
   //! \param os Output stream to write the solution data to
@@ -378,6 +421,9 @@ public:
   void dumpPrimSol(const Vector& psol, std::ostream& os,
 		   bool withID = true) const;
 
+  //! \brief Returns whether an analytical solution is available or not.
+  bool haveAnaSol() const { return mySol ? true : false; }
+
 protected:
   //! \brief Defines the type of a property set.
   //! \param[in] code The property code to be associated with the property type