diff --git a/Apps/FiniteDefElasticity/main_NonLinEl.C b/Apps/FiniteDefElasticity/main_NonLinEl.C
index 80f399e9..e05d82a9 100644
--- a/Apps/FiniteDefElasticity/main_NonLinEl.C
+++ b/Apps/FiniteDefElasticity/main_NonLinEl.C
@@ -41,6 +41,7 @@
   \arg -nv \a nv : Number of visualization points per knot-span in v-direction
   \arg -nw \a nw : Number of visualization points per knot-span in w-direction
   \arg -skip2nd : Skip VTF-output of secondary solution fields
+  \arg -noEnergy : Skip integration of energy norms
   \arg -saveInc \a dtSave : Time increment between each result save to VTF
   \arg -dumpInc \a dtDump [raw] : Time increment between each solution dump
   \arg -ignore \a p1, \a p2, ... : Ignore these patches in the analysis
@@ -65,12 +66,14 @@ int main (int argc, char** argv)
   int form = SIM::TOTAL_LAGRANGE;
   int nGauss = 4;
   int format = -1;
+  int outPrec = 3;
   int n[3] = { 2, 2, 2 };
   std::vector<int> ignoredPatches;
   double dtSave = 0.0;
   double dtDump = 0.0;
   bool dumpWithID = true;
   bool skip2nd = false;
+  bool energy = true;
   bool checkRHS = false;
   bool fixDup = false;
   bool twoD = false;
@@ -108,8 +111,12 @@ int main (int argc, char** argv)
       n[2] = atoi(argv[++i]);
     else if (!strcmp(argv[i],"-skip2nd"))
       skip2nd = true;
+    else if (!strcmp(argv[i],"-noEnergy"))
+      energy = false;
     else if (!strcmp(argv[i],"-saveInc") && i < argc-1)
       dtSave = atof(argv[++i]);
+    else if (!strcmp(argv[i],"-outPrec") && i < argc-1)
+      outPrec = atoi(argv[++i]);
     else if (!strcmp(argv[i],"-dumpInc") && i < argc-1)
     {
       dtDump = atof(argv[++i]);
@@ -291,11 +298,11 @@ int main (int argc, char** argv)
   while (simulator.advanceStep(params))
   {
     // Solve the nonlinear FE problem at this load step
-    if (!simulator.solveStep(params,SIM::STATIC))
+    if (!simulator.solveStep(params,SIM::STATIC,"displacement",energy))
       return 5;
 
     // Print solution components at the user-defined points
-    simulator.dumpResults(params.time.t,std::cout);
+    simulator.dumpResults(params.time.t,std::cout,outPrec);
 
     if (params.time.t + epsT*params.time.dt > nextDump)
     {
diff --git a/src/ASM/ASMbase.C b/src/ASM/ASMbase.C
index 18b71e20..052670b5 100644
--- a/src/ASM/ASMbase.C
+++ b/src/ASM/ASMbase.C
@@ -501,6 +501,27 @@ bool ASMbase::tesselate (ElementBlock&, const int*) const
 }
 
 
+bool ASMbase::getSolution (Matrix& sField, const Vector& locSol,
+			   const IntVec& nodes) const
+{
+  sField.resize(nf,nodes.size());
+  for (size_t i = 0; i < nodes.size(); i++)
+    if (nodes[i] < 1 || (size_t)nodes[i] > MLGN.size())
+    {
+      std::cerr <<" *** ASMbase::getSolution: Node #"<< nodes[i]
+		<<" is out of range [1,"<< MLGN.size() <<"]."<< std::endl;
+      return false;
+    }
+    else
+    {
+      for (unsigned char j = 0; j < nf; j++)
+	sField(j+1,i+1) = locSol[nf*(nodes[i]-1)+j];
+    }
+
+  return true;
+}
+
+
 bool ASMbase::evalSolution (Matrix&, const Vector&, const int*) const
 {
   std::cerr <<" *** ASMBase::evalSolution: Must be implemented in sub-class."
diff --git a/src/ASM/ASMbase.h b/src/ASM/ASMbase.h
index 8eb35f92..0e11e0cb 100644
--- a/src/ASM/ASMbase.h
+++ b/src/ASM/ASMbase.h
@@ -229,7 +229,8 @@ public:
   //! \param[in] xi Dimensionless parameters in range [0.0,1.0] of the point
   //! \param[out] param The parameters of the point in the knot-span domain
   //! \param[out] X The Cartesian coordinates of the point
-  virtual bool evalPoint(const double* xi, double* param, Vec3& X) const = 0;
+  //! \return Local node number within the patch that matches the point
+  virtual int evalPoint(const double* xi, double* param, Vec3& X) const = 0;
 
   //! \brief Creates a standard FE model of this patch for visualization.
   //! \param[out] grid The generated finite element grid
@@ -237,6 +238,13 @@ public:
   //! \note The number of element nodes must be set in \a grid on input.
   virtual bool tesselate(ElementBlock& grid, const int* npe) const;
 
+  //! \brief Extract the primary solution field at the specified nodes.
+  //! \param[out] sField Solution field
+  //! \param[in] locSol Solution vector local to current patch
+  //! \param[in] nodes 1-based local node numbers to extract solution for
+  bool getSolution(Matrix& sField, const Vector& locSol,
+		   const IntVec& nodes) const;
+
   //! \brief Evaluates the primary solution field at all visualization points.
   //! \param[out] sField Solution field
   //! \param[in] locSol Solution vector local to current patch
diff --git a/src/ASM/ASMs1D.C b/src/ASM/ASMs1D.C
index 4b4e2ed2..b1285b59 100644
--- a/src/ASM/ASMs1D.C
+++ b/src/ASM/ASMs1D.C
@@ -603,18 +603,35 @@ bool ASMs1D::integrate (Integrand& integrand, int lIndex,
 }
 
 
-bool ASMs1D::evalPoint (const double* xi, double* param, Vec3& X) const
+int ASMs1D::evalPoint (const double* xi, double* param, Vec3& X) const
 {
-  if (!curv) return false;
+  if (!curv) return -1;
 
   param[0] = (1.0-xi[0])*curv->startparam() + xi[0]*curv->endparam();
 
   Go::Point X0;
   curv->point(X0,param[0]);
-  for (unsigned char i = 0; i < nsd; i++)
-    X[i] = X0[i];
+  for (unsigned char d = 0; d < nsd; d++)
+    X[d] = X0[d];
 
-  return true;
+  // Check if this point matches any of the control points (nodes)
+  Vec3 Xnod;
+  size_t inod = 1;
+  RealArray::const_iterator cit = curv->coefs_begin();
+  for (int i = 0; cit != curv->coefs_end(); cit++, i++)
+  {
+    if (i < nsd) Xnod[i] = *cit;
+    if (i+1 == curv->dimension())
+      if (X.equal(Xnod,0.001))
+	return inod;
+      else
+      {
+	inod++;
+	i = -1;
+      }
+  }
+
+  return 0;
 }
 
 
diff --git a/src/ASM/ASMs1D.h b/src/ASM/ASMs1D.h
index 3625e6e6..9b4663dc 100644
--- a/src/ASM/ASMs1D.h
+++ b/src/ASM/ASMs1D.h
@@ -126,7 +126,9 @@ public:
   //! \param[in] xi Dimensionless parameter in range [0.0,1.0] of the point
   //! \param[out] param The parameter of the point in the knot-span domain
   //! \param[out] X The Cartesian coordinates of the point
-  virtual bool evalPoint(const double* xi, double* param, Vec3& X) const;
+  //! \return Local node number within the patch that matches the point, if any
+  //! \return 0 if no node (control point) matches this point
+  virtual int evalPoint(const double* xi, double* param, Vec3& X) const;
 
   //! \brief Creates a line element model of this patch for visualization.
   //! \param[out] grid The generated line grid
diff --git a/src/ASM/ASMs1DLag.C b/src/ASM/ASMs1DLag.C
index ac28374a..19f3100e 100644
--- a/src/ASM/ASMs1DLag.C
+++ b/src/ASM/ASMs1DLag.C
@@ -281,6 +281,24 @@ bool ASMs1DLag::integrate (Integrand& integrand, int lIndex,
 }
 
 
+int ASMs1DLag::evalPoint (const double* xi, double* param, Vec3& X) const
+{
+  if (!curv) return -1;
+
+  param[0] = (1.0-xi[0])*curv->startparam() + xi[0]*curv->endparam();
+
+  // Evaluate the parametric values of the nodes
+  RealArray u;
+  if (!this->getGridParameters(u,curv->order()-1)) return -1;
+
+  // Search for the closest node
+  size_t n = utl::find_closest(u,param[0]);
+  X = coord[n];
+
+  return 1+n;
+}
+
+
 bool ASMs1DLag::tesselate (ElementBlock& grid, const int*) const
 {
   size_t i, l;
diff --git a/src/ASM/ASMs1DLag.h b/src/ASM/ASMs1DLag.h
index 51015a46..f11dc2d1 100644
--- a/src/ASM/ASMs1DLag.h
+++ b/src/ASM/ASMs1DLag.h
@@ -79,6 +79,13 @@ public:
   // Post-processing methods
   // =======================
 
+  //! \brief Evaluates the geometry at a specified point.
+  //! \param[in] xi Dimensionless parameters in range [0.0,1.0] of the point
+  //! \param[out] param The parameters of the point in the knot-span domain
+  //! \param[out] X The Cartesian coordinates of the point
+  //! \return Local node number within the patch that is closest to the point
+  virtual int evalPoint(const double* xi, double* param, Vec3& X) const;
+
   //! \brief Creates a line element model of this patch for visualization.
   //! \param[out] grid The generated line grid
   //! \note The number of element nodes must be set in \a grid on input.
diff --git a/src/ASM/ASMs2D.C b/src/ASM/ASMs2D.C
index f600cf27..abe2ded7 100644
--- a/src/ASM/ASMs2D.C
+++ b/src/ASM/ASMs2D.C
@@ -1078,19 +1078,36 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex,
 }
 
 
-bool ASMs2D::evalPoint (const double* xi, double* param, Vec3& X) const
+int ASMs2D::evalPoint (const double* xi, double* param, Vec3& X) const
 {
-  if (!surf) return false;
+  if (!surf) return -2;
 
   param[0] = (1.0-xi[0])*surf->startparam_u() + xi[0]*surf->endparam_u();
   param[1] = (1.0-xi[1])*surf->startparam_v() + xi[1]*surf->endparam_v();
 
   Go::Point X0;
   surf->point(X0,param[0],param[1]);
-  for (unsigned char i = 0; i < nsd; i++)
-    X[i] = X0[i];
+  for (unsigned char d = 0; d < nsd; d++)
+    X[d] = X0[d];
 
-  return true;
+  // Check if this point matches any of the control points (nodes)
+  Vec3 Xnod;
+  size_t inod = 1;
+  RealArray::const_iterator cit = surf->coefs_begin();
+  for (int i = 0; cit != surf->coefs_end(); cit++, i++)
+  {
+    if (i < nsd) Xnod[i] = *cit;
+    if (i+1 == surf->dimension())
+      if (X.equal(Xnod,0.001))
+	return inod;
+      else
+      {
+	inod++;
+	i = -1;
+      }
+  }
+
+  return 0;
 }
 
 
diff --git a/src/ASM/ASMs2D.h b/src/ASM/ASMs2D.h
index 9ee14ad2..d20920e1 100644
--- a/src/ASM/ASMs2D.h
+++ b/src/ASM/ASMs2D.h
@@ -199,7 +199,9 @@ public:
   //! \param[in] xi Dimensionless parameters in range [0.0,1.0] of the point
   //! \param[out] param The (u,v) parameters of the point in knot-span domain
   //! \param[out] X The Cartesian coordinates of the point
-  virtual bool evalPoint(const double* xi, double* param, Vec3& X) const;
+  //! \return Local node number within the patch that matches the point, if any
+  //! \return 0 if no node (control point) matches this point
+  virtual int evalPoint(const double* xi, double* param, Vec3& X) const;
 
   //! \brief Creates a quad element model of this patch for visualization.
   //! \param[out] grid The generated quadrilateral grid
diff --git a/src/ASM/ASMs2DLag.C b/src/ASM/ASMs2DLag.C
index 314d53fe..83a20823 100644
--- a/src/ASM/ASMs2DLag.C
+++ b/src/ASM/ASMs2DLag.C
@@ -386,6 +386,28 @@ bool ASMs2DLag::integrate (Integrand& integrand, int lIndex,
 }
 
 
+int ASMs2DLag::evalPoint (const double* xi, double* param, Vec3& X) const
+{
+  if (!surf) return -2;
+
+  param[0] = (1.0-xi[0])*surf->startparam_u() + xi[0]*surf->endparam_u();
+  param[1] = (1.0-xi[1])*surf->startparam_v() + xi[1]*surf->endparam_v();
+
+  // Evaluate the parametric values of the nodes
+  RealArray u, v;
+  if (!this->getGridParameters(u,0,surf->order_u()-1)) return -2;
+  if (!this->getGridParameters(v,1,surf->order_v()-1)) return -2;
+
+  // Search for the closest node
+  size_t i = utl::find_closest(u,param[0]);
+  size_t j = utl::find_closest(v,param[1]);
+  size_t n = u.size()*j + i;
+  X = coord[n];
+
+  return 1+n;
+}
+
+
 bool ASMs2DLag::tesselate (ElementBlock& grid, const int* npe) const
 {
   const int p1 = surf->order_u();
diff --git a/src/ASM/ASMs2DLag.h b/src/ASM/ASMs2DLag.h
index ccd5fbc3..d14984b7 100644
--- a/src/ASM/ASMs2DLag.h
+++ b/src/ASM/ASMs2DLag.h
@@ -79,6 +79,13 @@ public:
   // Post-processing methods
   // =======================
 
+  //! \brief Evaluates the geometry at a specified point.
+  //! \param[in] xi Dimensionless parameters in range [0.0,1.0] of the point
+  //! \param[out] param The parameters of the point in the knot-span domain
+  //! \param[out] X The Cartesian coordinates of the point
+  //! \return Local node number within the patch that is closest to the point
+  virtual int evalPoint(const double* xi, double* param, Vec3& X) const;
+
   //! \brief Creates a quad element model of this patch for visualization.
   //! \param[out] grid The generated quadrilateral grid
   //! \param[in] npe Number of visualization nodes over each knot span
diff --git a/src/ASM/ASMs2DmxLag.C b/src/ASM/ASMs2DmxLag.C
index 262ec07d..8cc35311 100644
--- a/src/ASM/ASMs2DmxLag.C
+++ b/src/ASM/ASMs2DmxLag.C
@@ -247,7 +247,7 @@ bool ASMs2DmxLag::integrate (Integrand& integrand,
 
 	  // Evaluate the integrand and accumulate element contributions
 	  fe.detJxW *= wg[i]*wg[j];
-	  if (!integrand.evalInt(elmInt,fe,time,X))
+	  if (!integrand.evalIntMx(elmInt,fe,time,X))
 	    return false;
 	}
 
@@ -358,7 +358,7 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, int lIndex,
 
 	// Evaluate the integrand and accumulate element contributions
 	fe.detJxW *= wg[i];
-	if (!integrand.evalBou(elmInt,fe,time,X,normal))
+	if (!integrand.evalBouMx(elmInt,fe,time,X,normal))
 	  return false;
       }
 
diff --git a/src/ASM/ASMs3D.C b/src/ASM/ASMs3D.C
index 5ba48e6b..15551107 100644
--- a/src/ASM/ASMs3D.C
+++ b/src/ASM/ASMs3D.C
@@ -1593,20 +1593,37 @@ bool ASMs3D::integrateEdge (Integrand& integrand, int lEdge,
 }
 
 
-bool ASMs3D::evalPoint (const double* xi, double* param, Vec3& X) const
+int ASMs3D::evalPoint (const double* xi, double* param, Vec3& X) const
 {
-  if (!svol) return false;
+  if (!svol) return -3;
 
-  for (int i = 0; i < 3; i++)
+  int i;
+  for (i = 0; i < 3; i++)
     param[i] = (1.0-xi[i])*svol->startparam(i) + xi[i]*svol->endparam(i);
 
   Go::Point X0;
   svol->point(X0,param[0],param[1],param[2]);
-  X.x = X0[0];
-  X.y = X0[1];
-  X.z = X0[2];
+  for (i = 0; i < 3 && i < svol->dimension(); i++)
+    X[i] = X0[i];
 
-  return true;
+  // Check if this point matches any of the control points (nodes)
+  Vec3 Xnod;
+  size_t inod = 1;
+  RealArray::const_iterator cit = svol->coefs_begin();
+  for (i = 0; cit != svol->coefs_end(); cit++, i++)
+  {
+    if (i < 3) Xnod[i] = *cit;
+    if (i+1 == svol->dimension())
+      if (X.equal(Xnod,0.001))
+	return inod;
+      else
+      {
+	inod++;
+	i = -1;
+      }
+  }
+
+  return 0;
 }
 
 
diff --git a/src/ASM/ASMs3D.h b/src/ASM/ASMs3D.h
index 65b54485..ca3ca002 100644
--- a/src/ASM/ASMs3D.h
+++ b/src/ASM/ASMs3D.h
@@ -264,7 +264,9 @@ public:
   //! \param[in] xi Dimensionless parameters in range [0.0,1.0] of the point
   //! \param[out] param The (u,v,w) parameters of the point in knot-span domain
   //! \param[out] X The Cartesian coordinates of the point
-  virtual bool evalPoint(const double* xi, double* param, Vec3& X) const;
+  //! \return Local node number within the patch that matches the point, if any
+  //! \return 0 if no node (control point) matches this point
+  virtual int evalPoint(const double* xi, double* param, Vec3& X) const;
 
   //! \brief Creates a hexahedron element model of this patch for visualization.
   //! \param[out] grid The generated hexahedron grid
diff --git a/src/ASM/ASMs3DLag.C b/src/ASM/ASMs3DLag.C
index 49d92f3c..a44e9b89 100644
--- a/src/ASM/ASMs3DLag.C
+++ b/src/ASM/ASMs3DLag.C
@@ -544,6 +544,29 @@ bool ASMs3DLag::integrateEdge (Integrand& integrand, int lEdge,
 }
 
 
+int ASMs3DLag::evalPoint (const double* xi, double* param, Vec3& X) const
+{
+  if (!svol) return -3;
+
+  // Evaluate the parametric values of the point and nodes
+  RealArray u[3];
+  for (int d = 0; d < 3; d++)
+  {
+    param[d] = (1.0-xi[d])*svol->startparam(d) + xi[d]*svol->endparam(d);
+    if (!this->getGridParameters(u[d],d,svol->order(d)-1)) return -3;
+  }
+
+  // Search for the closest node
+  size_t i = utl::find_closest(u[0],param[0]);
+  size_t j = utl::find_closest(u[1],param[1]);
+  size_t k = utl::find_closest(u[2],param[2]);
+  size_t n = u[0].size()*(u[1].size()*k + j) + i;
+  X = coord[n];
+
+  return 1+n;
+}
+
+
 bool ASMs3DLag::tesselate (ElementBlock& grid, const int* npe) const
 {
   const int p1 = svol->order(0);
diff --git a/src/ASM/ASMs3DLag.h b/src/ASM/ASMs3DLag.h
index b4c81e4f..f74cf3f3 100644
--- a/src/ASM/ASMs3DLag.h
+++ b/src/ASM/ASMs3DLag.h
@@ -87,6 +87,13 @@ public:
   // Post-processing methods
   // =======================
 
+  //! \brief Evaluates the geometry at a specified point.
+  //! \param[in] xi Dimensionless parameters in range [0.0,1.0] of the point
+  //! \param[out] param The parameters of the point in the knot-span domain
+  //! \param[out] X The Cartesian coordinates of the point
+  //! \return Local node number within the patch that is closest to the point
+  virtual int evalPoint(const double* xi, double* param, Vec3& X) const;
+
   //! \brief Creates a hexahedron element model of this patch for visualization.
   //! \param[out] grid The generated hexahedron grid
   //! \param[in] npe Number of visualization nodes over each knot span
diff --git a/src/LinAlg/LinAlgInit.C b/src/LinAlg/LinAlgInit.C
index 5388eef5..3351c7ad 100644
--- a/src/LinAlg/LinAlgInit.C
+++ b/src/LinAlg/LinAlgInit.C
@@ -19,11 +19,12 @@
 #include "slepceps.h"
 #endif
 
-LinAlgInit* LinAlgInit::instance;
-int LinAlgInit::refs;
+
+LinAlgInit* LinAlgInit::instance = 0;
+int         LinAlgInit::refs = 0;
 
 
-LinAlgInit& LinAlgInit::Init(int argc, char** argv)
+LinAlgInit& LinAlgInit::Init (int argc, char** argv)
 {
   if (!instance)
     instance = new LinAlgInit(argc,argv);
@@ -31,18 +32,6 @@ LinAlgInit& LinAlgInit::Init(int argc, char** argv)
   return *instance;
 }
 
-void LinAlgInit::increfs()
-{
-  refs++;
-}
-
-void LinAlgInit::decrefs()
-{
-  refs--;
-  if (!refs)
-    delete instance;
-}
-
 
 LinAlgInit::LinAlgInit (int argc, char** argv)
 {
diff --git a/src/LinAlg/LinAlgInit.h b/src/LinAlg/LinAlgInit.h
index 09b4f9be..8fbd5c1c 100644
--- a/src/LinAlg/LinAlgInit.h
+++ b/src/LinAlg/LinAlgInit.h
@@ -1,4 +1,4 @@
-// $Id: LinAlgInit.h,v 1.2 2011-01-02 15:50:35 kmo Exp $
+// $Id$
 //==============================================================================
 //!
 //! \file LinAlgInit.h
@@ -17,31 +17,37 @@
 
 /*!
   \brief Initializer for linear algebra packages.
+
+  \details This class is a reference-counted singleton. We need control over the
+  destruction order since the destructor calls \a PetscFinalize() which shuts
+  down MPI. Thus, the LinAlgInit object must be destroyed after all other
+  objects that do MPI calls in their destructor.
+  This cannot be guaranteed without the reference-counting since this class is
+  instanciated as a local variable in the \a main() scope of the applications.
 */
 
 class LinAlgInit
 {
+  //! \brief The constructor uses the command-line arguments to set up things.
+  LinAlgInit(int argc, char** argv);
 public:
-  //! \brief This class is a ref-counted singleton. We need control over 
-  //         the destruction order since the destructor calls PetscFinalize() 
-  //         which shuts down MPI. Thus this object must be destroyed after 
-  //         objects that do MPI calls in their destructor.
-  //         This we cannot guarantee without the ref-counting since this
-  //         class is instanced as a local variable in the main() scope of 
-  //         applications.
-  static LinAlgInit& Init(int argc, char** argv);
   //! \brief The destructor finalizes the linear algebra packages.
   ~LinAlgInit();
 
-  int myPid; //!< Processor ID in parallel simulations
-  static void increfs();
-  static void decrefs();
-private:
-  static int refs;
-  //! \brief The constructor uses the command-line arguments to set up things.
-  LinAlgInit(int argc, char** argv);
+  //! \brief Increments the reference counter.
+  static void increfs() { ++refs; }
+  //! \brief Decrements the reference counter.
+  static void decrefs() { if (--refs == 0) delete instance; }
 
-  static LinAlgInit* instance;
+  //! \brief Instanciates the singleton LinAlgInit object.
+  static LinAlgInit& Init(int argc, char** argv);
+
+  int myPid; //!< Processor ID in parallel simulations
+
+private:
+  static int refs; //!< Reference counter
+
+  static LinAlgInit* instance; //!< Pointer to the singleton object
 };
 
 #endif
diff --git a/src/LinAlg/LinSolParams.C b/src/LinAlg/LinSolParams.C
index cb96b564..7e8c98f1 100644
--- a/src/LinAlg/LinSolParams.C
+++ b/src/LinAlg/LinSolParams.C
@@ -1,4 +1,4 @@
-// $Id: LinSolParams.C,v 1.6 2011-02-08 12:45:08 rho Exp $
+// $Id$
 //==============================================================================
 //!
 //! \file LinSolParams.C
@@ -25,7 +25,7 @@ void LinSolParams::setDefault ()
   // Use GMRES with ILU preconditioner as default
   method  = KSPGMRES;
   prec    = PCILU;
-  package = MAT_SOLVER_PETSC; 
+  package = MAT_SOLVER_PETSC;
   levels  = 0;
   overlap = 0;
   nullspc = NONE;
@@ -126,8 +126,7 @@ bool LinSolParams::read (std::istream& is, int nparam)
 	nullspc = NONE;
     }
 
-    else
-    {
+    else {
       std::cerr <<" *** LinSolParams::read: Unknown keyword: "
 		<< cline << std::endl;
       return false;
diff --git a/src/LinAlg/SAM.C b/src/LinAlg/SAM.C
index c0f01cbf..c5814048 100644
--- a/src/LinAlg/SAM.C
+++ b/src/LinAlg/SAM.C
@@ -1,4 +1,4 @@
-// $Id: SAM.C,v 1.27 2011-02-05 18:07:50 kmo Exp $
+// $Id$
 //==============================================================================
 //!
 //! \file SAM.C
@@ -729,3 +729,21 @@ real SAM::normReact (const Vector& u, const Vector& rf) const
 
   return retVal;
 }
+
+
+bool SAM::getNodalReactions (int inod, const Vector& rf, Vector& nrf) const
+{
+  if (inod < 1 || inod > nnod) return false;
+
+  bool haveRF = false;
+  int ip = madof[inod-1]-1;
+  nrf.resize(madof[inod]-1-ip,true);
+  for (size_t i = 0; i < nrf.size(); i++, ip++)
+    if (msc[ip] < 0 && -msc[ip] <= (int)rf.size())
+    {
+      haveRF = true;
+      nrf[i] = rf(-msc[ip]);
+    }
+
+  return haveRF;
+}
diff --git a/src/LinAlg/SAM.h b/src/LinAlg/SAM.h
index 39cb9ebf..5116982e 100644
--- a/src/LinAlg/SAM.h
+++ b/src/LinAlg/SAM.h
@@ -1,4 +1,4 @@
-// $Id: SAM.h,v 1.25 2011-02-08 12:45:46 rho Exp $
+// $Id$
 //==============================================================================
 //!
 //! \file SAM.h
@@ -101,7 +101,8 @@ public:
   //! \param sysRHS The system right-hand-side load vector to be initialized
   //! \param reactionForces Pointer to vector of nodal reaction forces
   //! \return \e false if no free DOFs in the system, otherwise \e true
-  virtual bool initForAssembly(SystemVector& sysRHS, Vector* reactionForces = 0) const;
+  virtual bool initForAssembly(SystemVector& sysRHS,
+			       Vector* reactionForces = 0) const;
 
   //! \brief Adds an element stiffness matrix into the system stiffness matrix.
   //! \param sysK    The left-hand-side system stiffness matrix
@@ -204,6 +205,12 @@ public:
   //! \param[in] u The (incremental) nodal displacement vector
   //! \param[in] rf Compressed reaction force vector
   virtual real normReact(const Vector& u, const Vector& rf) const;
+  //! \brief Returns a vector of reaction forces for a given node.
+  //! \param[in] inod Identifier for the node to get the reaction forces for
+  //! \param[in] rf Compressed reaction force vector for the entire model
+  //! \param[out] nrf Nodal reaction forces
+  //! \return \e true of the specified node has reaction forces
+  bool getNodalReactions(int inod, const Vector& rf, Vector& nrf) const;
 
 protected:
   //! \brief Initializes the DOF-to-equation connectivity array \a MEQN.
diff --git a/src/SIM/NonLinSIM.C b/src/SIM/NonLinSIM.C
index 01fb9ced..dbd11a2f 100644
--- a/src/SIM/NonLinSIM.C
+++ b/src/SIM/NonLinSIM.C
@@ -125,7 +125,8 @@ bool NonLinSIM::advanceStep (SolvePrm& param)
 }
 
 
-bool NonLinSIM::solveStep (SolvePrm& param, SIM::SolutionMode mode)
+bool NonLinSIM::solveStep (SolvePrm& param, SIM::SolutionMode mode,
+			   const char* compName, bool energyNorm)
 {
   PROFILE1("NonLinSIM::solveStep");
 
@@ -160,7 +161,7 @@ bool NonLinSIM::solveStep (SolvePrm& param, SIM::SolutionMode mode)
 	  return false;
 
 	model->setMode(SIM::RECOVERY);
-	if (!this->solutionNorms(param.time))
+	if (!this->solutionNorms(param.time,compName,energyNorm))
 	  return false;
 
 	param.time.first = false;
@@ -324,7 +325,8 @@ bool NonLinSIM::updateConfiguration (SolvePrm& param)
 }
 
 
-bool NonLinSIM::solutionNorms (const TimeDomain& time, const char* compName)
+bool NonLinSIM::solutionNorms (const TimeDomain& time, const char* compName,
+			       bool energyNorm)
 {
   if (msgLevel < 0 || solution.empty()) return true;
 
@@ -334,8 +336,9 @@ bool NonLinSIM::solutionNorms (const TimeDomain& time, const char* compName)
   double dMax[nsd];
   double normL2 = model->solutionNorms(solution.front(),dMax,iMax);
 
-  if (!model->solutionNorms(time,solution,eNorm,gNorm))
-    gNorm.clear();
+  if (energyNorm)
+    if (!model->solutionNorms(time,solution,eNorm,gNorm))
+      gNorm.clear();
 
   if (myPid == 0)
   {
@@ -426,9 +429,9 @@ void NonLinSIM::dumpStep (int iStep, double time, std::ostream& os,
 }
 
 
-bool NonLinSIM::dumpResults (double time, std::ostream& os) const
+bool NonLinSIM::dumpResults (double time, std::ostream& os, int precision) const
 {
-  return model->dumpResults(solution.front(),time,os);
+  return model->dumpResults(solution.front(),time,os,precision);
 }
 
 
diff --git a/src/SIM/NonLinSIM.h b/src/SIM/NonLinSIM.h
index 84ffb088..6e8d145e 100644
--- a/src/SIM/NonLinSIM.h
+++ b/src/SIM/NonLinSIM.h
@@ -74,7 +74,11 @@ public:
   //! \brief Solves the nonlinear equations by Newton-Raphson iterations.
   //! \param param Solution algorithm parameters
   //! \param[in] mode Solution mode to use for this step
-  bool solveStep(SolvePrm& param, SIM::SolutionMode mode = SIM::STATIC);
+  //! \param[in] compName Solution name to be used in the norm output
+  //! \param[in] energyNorm If \e true, integrate energy norm of the solution
+  bool solveStep(SolvePrm& param, SIM::SolutionMode mode = SIM::STATIC,
+		 const char* compName = "displacement",
+		 bool energyNorm = false);
 
   //! \brief Enum describing the norm used for convergence checks.
   enum CNORM { L2, ENERGY };
@@ -86,8 +90,9 @@ public:
   //! \brief Computes and prints some solution norm quantities.
   //! \param[in] time Parameters for nonlinear/time-dependent simulations
   //! \param[in] compName Solution name to be used in the norm output
+  //! \param[in] energyNorm If \e true, integrate energy norm of the solution
   bool solutionNorms(const TimeDomain& time,
-		     const char* compName = "displacement");
+		     const char* compName, bool energyNorm);
 
   //! \brief Saves the converged results to VTF file of a given load/time step.
   //! \param[in] iStep Load/time step identifier
@@ -107,7 +112,8 @@ public:
   //! \brief Dumps solution variables at user-defined points.
   //! \param[in] time Current time/load parameter
   //! \param[in] os The output stream to write the solution to
-  bool dumpResults(double time, std::ostream& os) const;
+  //! \param[in] precision Number of digits after the decimal point
+  bool dumpResults(double time, std::ostream& os, int precision = 3) const;
 
   //! \brief Returns a const reference to current solution vector.
   const Vector& getSolution() const { return solution.front(); }
diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C
index eb8120c6..95fd4b32 100644
--- a/src/SIM/SIMbase.C
+++ b/src/SIM/SIMbase.C
@@ -142,9 +142,9 @@ bool SIMbase::parse (char* keyWord, std::istream& is)
 	std::cout <<"\n\tPoint "<< i+1 <<": P"<< myPoints[i].patch <<" xi =";
       for (int j = 0; j < 3 && (cline = strtok(NULL," ")); j++)
       {
-        myPoints[i].param[j] = atof(cline);
+        myPoints[i].par[j] = atof(cline);
 	if (myPid == 0)
-	  std::cout <<' '<< myPoints[i].param[j];
+	  std::cout <<' '<< myPoints[i].par[j];
       }
     }
     if (myPid == 0 && nres > 0)
@@ -306,11 +306,11 @@ bool SIMbase::preprocess (const std::vector<int>& ignoredPatches, bool fixDup)
     int pid = this->getLocalPatchIndex(p->patch);
     if (pid < 1 || myModel[pid-1]->empty())
       p = myPoints.erase(p);
+    else if ((p->inod = myModel[pid-1]->evalPoint(p->par,p->par,p->X)) < 0)
+      p = myPoints.erase(p);
     else
     {
       p->npar = myModel[pid-1]->getNoParamDim();
-      myModel[pid-1]->evalPoint(p->param,p->param,p->X);
-
       std::cout <<"\nResult point #"<< 1+(int)(p-myPoints.begin())
 		<<": patch #"<< p->patch;
       switch (p->npar) {
@@ -318,11 +318,12 @@ bool SIMbase::preprocess (const std::vector<int>& ignoredPatches, bool fixDup)
       case 2: std::cout <<" (u,v)=("; break;
       case 3: std::cout <<" (u,v,w)=("; break;
       }
-      std::cout << p->param[0];
+      std::cout << p->par[0];
       for (unsigned char c = 1; c < p->npar; c++)
-	std::cout <<','<< p->param[c];
+	std::cout <<','<< p->par[c];
       if (p->npar > 1) std::cout <<')';
-	std::cout <<", X = "<< p->X;
+      if (p->inod > 0) std::cout <<", node #"<< p->inod;
+      std::cout <<", X = "<< p->X;
       p++;
     }
   }
@@ -1355,35 +1356,47 @@ bool SIMbase::dumpResults (const Vector& psol, double time, std::ostream& os,
 
   size_t i, j, k;
   Matrix sol1, sol2;
+  Vector reactionFS;
+  const Vector* reactionForces = myEqSys->getReactions();
 
   for (i = 0; i < myModel.size(); i++)
   {
     if (myModel[i]->empty()) continue; // skip empty patches
 
     ResPointVec::const_iterator p;
-    std::vector<size_t> points;
+    std::vector<int> points;
     RealArray params[3];
 
     // Find all evaluation points within this patch, if any
     for (j = 0, p = myPoints.begin(); p != myPoints.end(); j++, p++)
       if (this->getLocalPatchIndex(p->patch) == (int)(i+1))
-      {
-	points.push_back(j+1);
-	for (k = 0; k < myModel[i]->getNoParamDim(); k++)
-	  params[k].push_back(p->param[k]);
-      }
+	if (discretization == Spline)
+	{
+	  points.push_back(p->inod > 0 ? p->inod : -(j+1));
+	  for (k = 0; k < myModel[i]->getNoParamDim(); k++)
+	    params[k].push_back(p->par[k]);
+	}
+	else if (p->inod > 0)
+	  points.push_back(p->inod);
 
     if (points.empty()) continue; // no points in this patch
 
-    // Evaluate the primary solution variables
     myModel[i]->extractNodeVec(psol,myProblem->getSolution());
-    if (!myModel[i]->evalSolution(sol1,myProblem->getSolution(),params,false))
-      return false;
+    if (discretization == Spline)
+    {
+      // Evaluate the primary solution variables
+      if (!myModel[i]->evalSolution(sol1,myProblem->getSolution(),params,false))
+	return false;
 
-    // Evaluate the secondary solution variables
-    LocalSystem::patch = i;
-    if (!myModel[i]->evalSolution(sol2,*myProblem,params,false))
-      return false;
+      // Evaluate the secondary solution variables
+      LocalSystem::patch = i;
+      if (!myModel[i]->evalSolution(sol2,*myProblem,params,false))
+	return false;
+    }
+    else
+      // Extract nodal primary solution variables
+      if (!myModel[i]->getSolution(sol1,myProblem->getSolution(),points))
+	return false;
 
     // Formatted output, use scientific notation with fixed field width
     std::streamsize flWidth = 8 + outputPrecision;
@@ -1391,13 +1404,32 @@ bool SIMbase::dumpResults (const Vector& psol, double time, std::ostream& os,
     std::ios::fmtflags oldF = os.flags(std::ios::scientific | std::ios::right);
     for (j = 0; j < points.size(); j++)
     {
-      os <<"  Result point #"<< points[j];
-      os <<":  sol1 =";
+      if (points[j] > 0)
+      {
+	points[j] = myModel[i]->getNodeID(points[j]);
+	os <<"  Node #"<< points[j] <<":\tsol1 =";
+      }
+      else
+	os <<"  Point #"<< -points[j] <<":\tsol1 =";
       for (k = 1; k <= sol1.rows(); k++)
 	os << std::setw(flWidth) << sol1(k,j+1);
-      os <<"   sol2 =";
-      for (k = 1; k <= sol2.rows(); k++)
-	os << std::setw(flWidth) << sol2(k,j+1);
+
+      if (discretization == Spline)
+      {
+	os <<"   sol2 =";
+	for (k = 1; k <= sol2.rows(); k++)
+	  os << std::setw(flWidth) << sol2(k,j+1);
+      }
+
+      if (reactionForces && points[j] > 0)
+	// Print nodal reaction forces for nodes with prescribed DOFs
+	if (mySam->getNodalReactions(points[j],*reactionForces,reactionFS))
+	{
+	  os <<"\n\t\treac =";
+	  for (k = 0; k < reactionFS.size(); k++)
+	    os << std::setw(flWidth) << reactionFS[k];
+	}
+
       os << std::endl;
     }
     os.precision(oldPrec);
diff --git a/src/SIM/SIMbase.h b/src/SIM/SIMbase.h
index 28f0c03d..404ef5dd 100644
--- a/src/SIM/SIMbase.h
+++ b/src/SIM/SIMbase.h
@@ -51,12 +51,13 @@ struct Mode
 
 struct ResultPoint
 {
-  unsigned char npar;     //!< Number of parameters
-  size_t        patch;    //!< Patch index [0,nPatch>
-  double        param[3]; //!< Parameters of the point (u,v,w)
-  Vec3          X;        //!< Spatial coordinates of the point
+  unsigned char npar;   //!< Number of parameters
+  size_t        patch;  //!< Patch index [0,nPatch>
+  int           inod;   //!< Local node number of the closest node
+  double        par[3]; //!< Parameters of the point (u,v,w)
+  Vec3          X;      //!< Spatial coordinates of the point
   // \brief Default constructor.
-  ResultPoint() : npar(0), patch(0) { param[0] = param[1] = param[2] = 0.0; }
+  ResultPoint() : npar(0), patch(0), inod(0) { par[0] = par[1] = par[2] = 0.0; }
 };
 
 typedef std::vector<ResultPoint> ResPointVec; //!< Result point container
@@ -154,11 +155,6 @@ public:
   bool assembleSystem(const TimeDomain& time, const Vectors& pSol = Vectors(),
 		      bool newLHSmatrix = true);
 
-  //! \brief Finalizes system assembly.
-  //TODO: This HAS to be called from subclasses if assembleSystem is overridden.
-  //      This need to be fixed, probably through a wrapper function.
-  bool finalizeAssembly(bool newLHSmatrix);
-
   //! \brief Administers assembly of the linear equation system.
   //! \param[in] pSol Previous primary solution vectors in DOF-order
   //!
@@ -238,6 +234,7 @@ public:
   //! spline basis to obtain the control point values of the secondary solution.
   bool project(Vector& psol);
 
+
   // Post-processing methods
   // =======================
 
@@ -363,15 +360,18 @@ protected:
   //! \brief Initializes for integration of Neumann terms for a given property.
   virtual bool initNeumann(size_t) { return true; }
 
-  //! \brief Extract local solution vector(s) for a given patch.
+  //! \brief Reads a LinSolParams object from the given stream.
+  //! \details This method helps with encapsulating PETSc in libIFEM.
+  void readLinSolParams(std::istream& is, int npar);
+
+  //! \brief Finalizes the system assembly.
+  bool finalizeAssembly(bool newLHSmatrix);
+
+  //! \brief Extracts local solution vector(s) for the given patch.
   //! \param[in] patch Pointer to the patch to extract solution vectors for
   //! \param[in] sol Global primary solution vectors in DOF-order
   void extractPatchSolution(const ASMbase* patch, const Vectors& sol);
 
-  //! \brief Read a LinSolParams object from the given stream.
-  //! \details This method helps with encapsulating PETSc in libIFEM.
-  void readLinSolParams(std::istream& is, int npar);
-
 public:
   //! \brief Enum defining the available discretization methods.
   enum Discretization { Spline, Lagrange, Spectral };
diff --git a/src/Utility/Tensor.C b/src/Utility/Tensor.C
index 8282e710..f9de50ca 100644
--- a/src/Utility/Tensor.C
+++ b/src/Utility/Tensor.C
@@ -71,7 +71,7 @@ Tensor& Tensor::operator= (const Tensor& T)
     else
       ndim = n;
     for (t_ind i = 1; i <= ndim; i++)
-      for (t_ind j = (this->symmetric() ? i : 1); j <= n; j++)
+      for (t_ind j = (this->symmetric() ? i : 1); j <= ndim; j++)
 	v[this->index(i,j)] = T(i,j);
   }
 
diff --git a/src/Utility/Utilities.C b/src/Utility/Utilities.C
index c195ceb3..f6941503 100644
--- a/src/Utility/Utilities.C
+++ b/src/Utility/Utilities.C
@@ -189,3 +189,17 @@ void utl::Pascal (int p, real x, real y, real z, std::vector<real>& phi)
 	  phi.push_back(a);
 	}
 }
+
+
+size_t utl::find_closest (const std::vector<real>& a, real v)
+{
+  // The lower_bound function uses binary search to find the index of the value.
+  // Thus, this works only when the vector a is sorted in increasing order.
+  size_t i = std::lower_bound(a.begin(),a.end(),v) - a.begin();
+  if (i > 0 && i < a.size())
+    return a[i]-v < v-a[i-1] ? i : i-1;
+  else if (i == 0)
+    return 0;
+  else
+    return i-1;
+}
diff --git a/src/Utility/Utilities.h b/src/Utility/Utilities.h
index 911f3cf5..7a60220a 100644
--- a/src/Utility/Utilities.h
+++ b/src/Utility/Utilities.h
@@ -89,6 +89,14 @@ namespace utl
   void Pascal(int p, real x, real y, std::vector<real>& phi);
   //! \brief Evaluates the monomials of Pascal's triangle in 3D for order \a p.
   void Pascal(int p, real x, real y, real z, std::vector<real>& phi);
+
+  //! \brief Searches for a real value in an ordered array of reals.
+  //! \param[in] a The array of real values
+  //! \param[in] v The value to search for
+  //! \return 0-based index of the closest array value.
+  //! \note It is assumed that the array \a a is sorted in encreasing order on
+  //! input. If this is not the case the method will deliver incorrect result.
+  size_t find_closest(const std::vector<real>& a, real v);
 }
 
 #endif