diff --git a/src/ASM/ASMbase.C b/src/ASM/ASMbase.C
index 052670b5..d44bbc6c 100644
--- a/src/ASM/ASMbase.C
+++ b/src/ASM/ASMbase.C
@@ -539,7 +539,7 @@ bool ASMbase::evalSolution (Matrix&, const Vector&,
 }
 
 
-bool ASMbase::evalSolution (Matrix&, const Integrand&, const int*) const
+bool ASMbase::evalSolution (Matrix&, const Integrand&, const int*, bool) const
 {
   std::cerr <<" *** ASMBase::evalSolution: Must be implemented in sub-class."
 	    << std::endl;
diff --git a/src/ASM/ASMbase.h b/src/ASM/ASMbase.h
index 26dfcc25..721e712b 100644
--- a/src/ASM/ASMbase.h
+++ b/src/ASM/ASMbase.h
@@ -270,14 +270,17 @@ public:
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods
   //! \param[in] npe Number of visualization nodes over each knot span
+  //! \param[in] project Flag indicating if the projected solution is wanted
   //!
   //! \details The secondary solution is derived from the primary solution,
   //! which is assumed to be stored within the \a integrand for current patch.
   //! If \a npe is NULL, the solution is evaluated at the Greville points and
   //! then projected onto the spline basis to obtain the control point values,
   //! which then are returned through \a sField.
+  //! If \a npe is not NULL and \a project is \e true, the solution is also
+  //! projected onto the spline basis, and then evaluated at the \a npe points.
   virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
-			    const int* npe = 0) const;
+			    const int* npe = 0, bool project = false) const;
 
   //! \brief Evaluates the secondary solution field at the given points.
   //! \param[out] sField Solution field
diff --git a/src/ASM/ASMs1D.C b/src/ASM/ASMs1D.C
index 260c6859..d858e233 100644
--- a/src/ASM/ASMs1D.C
+++ b/src/ASM/ASMs1D.C
@@ -681,21 +681,24 @@ bool ASMs1D::tesselate (ElementBlock& grid, const int* npe) const
 
   // Evaluate the spline curve at all points
   size_t nx = gpar.size();
-  RealArray XYZ(3*nx,0.0);
-
-  // Establish the block grid coordinates
-  size_t i, j;
-  for (i = j = 0; i < nx; i++, j += 3)
+  RealArray XYZ(curv->dimension()*nx);
+  //curv->gridEvaluator(XYZ,gpar);
+  //TODO: Replace this loop by the above line when available in GoTools
+  size_t m, jp;
+  Go::Point pt;
+  for (m = jp = 0; m < nx; m++)
   {
-    Go::Point pt;
-    curv->point(pt,gpar[i]);
+    curv->point(pt,gpar[m]);
     for (int k = 0; k < pt.size(); k++)
-      XYZ[j+k] = pt[k];
+      XYZ[jp++] = pt[k];
   }
 
+  // Establish the block grid coordinates
+  size_t i, j, l;
   grid.resize(nx);
-  for (i = j = 0; i < grid.getNoNodes(); i++, j += 3)
-    grid.setCoor(i,XYZ[j],XYZ[j+1],XYZ[j+2]);
+  for (i = l = 0; i < grid.getNoNodes(); i++, l += curv->dimension())
+    for (j = 0; j < nsd; j++)
+      grid.setCoor(i,j,XYZ[l+j]);
 
   // Establish the block grid topology
   int n[2], ip = 0;
@@ -703,8 +706,8 @@ bool ASMs1D::tesselate (ElementBlock& grid, const int* npe) const
   n[1] = n[0] + 1;
 
   for (i = 1; i < nx; i++)
-    for (j = 0; j < 2; j++)
-      grid.setNode(ip++,n[j]++);
+    for (l = 0; l < 2; l++)
+      grid.setNode(ip++,n[l]++);
 
   return true;
 }
@@ -767,15 +770,38 @@ bool ASMs1D::evalSolution (Matrix& sField, const Vector& locSol,
 
 
 bool ASMs1D::evalSolution (Matrix& sField, const Integrand& integrand,
-			   const int* npe) const
+			   const int* npe, bool project) const
 {
   if (npe)
   {
     // Compute parameter values of the result sampling points
     RealArray gpar;
     if (this->getGridParameters(gpar,npe[0]-1))
-      // Evaluate the secondary solution at all sampling points
-      return this->evalSolution(sField,integrand,&gpar);
+      if (project)
+      {
+	// Project the secondary solution onto the spline basis
+	Go::SplineCurve* c = this->projectSolution(integrand);
+	if (c)
+	{
+	  // Evaluate the projected field at the result sampling points
+	  //const Vector& svec = sField; // using utl::matrix cast operator
+	  sField.resize(c->dimension(),gpar.size());
+	  //c->gridEvaluator(svec,gpar);
+	  //TODO: Replace this loop by the above line when available in GoTools
+	  Go::Point pt;
+	  for (size_t j = 0; j < sField.cols(); j++)
+	  {
+	    c->point(pt,gpar[j]);
+	    for (size_t i = 0; i < sField.rows(); i++)
+	      sField(i+1,j+1) = pt[i];
+	  }
+	  delete c;
+	  return true;
+	}
+      }
+      else
+	// Evaluate the secondary solution directly at all sampling points
+	return this->evalSolution(sField,integrand,&gpar);
   }
   else
   {
@@ -791,6 +817,7 @@ bool ASMs1D::evalSolution (Matrix& sField, const Integrand& integrand,
     }
   }
 
+  std::cerr <<" *** ASMs1D::evalSolution: Failure!"<< std::endl;
   return false;
 }
 
@@ -805,7 +832,7 @@ Go::SplineCurve* ASMs1D::projectSolution (const Integrand& integrand) const
 {
   //TODO: This requires the method Go::CurveInterpolator::regularInterpolation
   //      which is not yet present in GoTools.
-  std::cerr <<"ASMs1D::projectSolution: Not implemented yet!"<< std::endl;
+  std::cerr <<" *** ASMs1D::projectSolution: Not implemented yet!"<< std::endl;
   return 0;
 }
 
diff --git a/src/ASM/ASMs1D.h b/src/ASM/ASMs1D.h
index 59fa0bdf..e0f98c62 100644
--- a/src/ASM/ASMs1D.h
+++ b/src/ASM/ASMs1D.h
@@ -154,14 +154,17 @@ public:
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods
   //! \param[in] npe Number of visualization nodes over each knot span
+  //! \param[in] project Flag indicating if the projected solution is wanted
   //!
   //! \details The secondary solution is derived from the primary solution,
   //! which is assumed to be stored within the \a integrand for current patch.
   //! If \a npe is NULL, the solution is evaluated at the Greville points and
   //! then projected onto the spline basis to obtain the control point values,
   //! which then are returned through \a sField.
+  //! If \a npe is not NULL and \a project is \e true, the solution is also
+  //! projected onto the spline basis, and then evaluated at the \a npe points.
   virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
-			    const int* npe = 0) const;
+			    const int* npe = 0, bool project = false) const;
 
   //! \brief Projects the secondary solution field onto the primary basis.
   //! \param[in] integrand Object with problem-specific data and methods
diff --git a/src/ASM/ASMs1DLag.C b/src/ASM/ASMs1DLag.C
index 19f3100e..b6fcad9f 100644
--- a/src/ASM/ASMs1DLag.C
+++ b/src/ASM/ASMs1DLag.C
@@ -73,16 +73,16 @@ bool ASMs1DLag::generateFEMTopology ()
   coord.resize(nx);
 
   // Evaluate the nodal coordinates
+  Go::Point pt;
   for (size_t i = 0; i < nx; i++)
   {
-    Go::Point pt;
     curv->point(pt,gpar[i]);
     for (int k = 0; k < pt.size(); k++)
       coord[i][k] = pt[k];
     MLGN[i] = ++gNod;
   }
 
-  // Connectivety array: local --> global node relation
+  // Connectivity array: local --> global node relation
   MLGE.resize(nel);
   MNPC.resize(nel);
 
@@ -299,8 +299,16 @@ int ASMs1DLag::evalPoint (const double* xi, double* param, Vec3& X) const
 }
 
 
-bool ASMs1DLag::tesselate (ElementBlock& grid, const int*) const
+bool ASMs1DLag::tesselate (ElementBlock& grid, const int* npe) const
 {
+  const int p1 = curv->order();
+  if (p1 != npe[0])
+  {
+    std::cout <<"\nLagrange elements: The number of visualization points is "
+	      << p1 <<" by default\n"<< std::endl;
+    const_cast<int*>(npe)[0] = p1;
+  }
+
   size_t i, l;
 
   grid.resize(nx);
@@ -348,7 +356,7 @@ bool ASMs1DLag::evalSolution (Matrix&, const Vector&,
 
 
 bool ASMs1DLag::evalSolution (Matrix& sField, const Integrand& integrand,
-			      const int*) const
+			      const int*, bool) const
 {
   sField.resize(0,0);
 
diff --git a/src/ASM/ASMs1DLag.h b/src/ASM/ASMs1DLag.h
index 13df14ec..eec37291 100644
--- a/src/ASM/ASMs1DLag.h
+++ b/src/ASM/ASMs1DLag.h
@@ -39,7 +39,7 @@ public:
 
   //! \brief Generates the finite element topology data for the patch.
   //! \details The data generated are the element-to-node connectivity array,
-  //! and the global node and element numbers.
+  //! the nodal coordinate array, as well as global node and element numbers.
   virtual bool generateFEMTopology();
 
   //! \brief Clears the contents of the patch, making it empty.
@@ -86,8 +86,9 @@ public:
 
   //! \brief Creates a line element model of this patch for visualization.
   //! \param[out] grid The generated line grid
+  //! \param[in] npe Number of visualization nodes over each knot span
   //! \note The number of element nodes must be set in \a grid on input.
-  virtual bool tesselate(ElementBlock& grid, const int*) const;
+  virtual bool tesselate(ElementBlock& grid, const int* npe) const;
 
   //! \brief Evaluates the primary solution field at all visualization points.
   //! \details The number of visualization points is the same as the order of
@@ -110,7 +111,7 @@ public:
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods
   virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
-			    const int*) const;
+			    const int*, bool = false) const;
 
   //! \brief Evaluates the secondary solution field at the given points.
   //! \param[out] sField Solution field
@@ -135,7 +136,7 @@ protected:
   virtual void getNodalCoordinates(Matrix& X) const;
 
   //! \brief Returns the number of nodal points in the patch.
-  virtual int getSize() const { return nx; }
+  virtual int getSize(int = 0) const { return nx; }
 
 private:
   //! \brief The number of nodes in each direction for the patch
diff --git a/src/ASM/ASMs1DSpec.C b/src/ASM/ASMs1DSpec.C
index 3f7a61b1..656937ad 100644
--- a/src/ASM/ASMs1DSpec.C
+++ b/src/ASM/ASMs1DSpec.C
@@ -157,7 +157,7 @@ bool ASMs1DSpec::integrate (Integrand& integrand, int lIndex,
 
 
 bool ASMs1DSpec::evalSolution (Matrix& sField, const Integrand& integrand,
-			       const int*) const
+			       const int*, bool) const
 {
   sField.resize(0,0);
   if (!curv) return false;
diff --git a/src/ASM/ASMs1DSpec.h b/src/ASM/ASMs1DSpec.h
index 76e2ba7b..b34a401d 100644
--- a/src/ASM/ASMs1DSpec.h
+++ b/src/ASM/ASMs1DSpec.h
@@ -65,9 +65,8 @@ public:
   //! \brief Evaluates the secondary solution field at all visualization points.
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods
-  //! \param[in] npe Number of visualization nodes over each knot span
   virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
-			    const int* npe) const;
+			    const int*, bool = false) const;
 
 protected:
 
diff --git a/src/ASM/ASMs2D.C b/src/ASM/ASMs2D.C
index c1852ce7..d4c317b5 100644
--- a/src/ASM/ASMs2D.C
+++ b/src/ASM/ASMs2D.C
@@ -1178,14 +1178,10 @@ bool ASMs2D::tesselate (ElementBlock& grid, const int* npe) const
 
   // Establish the block grid coordinates
   size_t i, j, l;
-  double X[3] = { 0.0, 0.0, 0.0 };
   grid.resize(nx,ny);
-  for (i = j = 0; i < grid.getNoNodes(); i++, j += surf->dimension())
-  {
-    for (l = 0; l < nsd; l++)
-      X[l] = XYZ[j+l];
-    grid.setCoor(i,X[0],X[1],X[2]);
-  }
+  for (i = l = 0; i < grid.getNoNodes(); i++, l += surf->dimension())
+    for (j = 0; j < nsd; j++)
+      grid.setCoor(i,j,XYZ[l+j]);
 
   // Establish the block grid topology
   int n[4], ip = 0;
@@ -1288,20 +1284,31 @@ bool ASMs2D::evalSolution (Matrix& sField, const Vector& locSol,
 
 
 bool ASMs2D::evalSolution (Matrix& sField, const Integrand& integrand,
-			   const int* npe) const
+			   const int* npe, bool project) const
 {
-  bool retVal = true;
-
   if (npe)
   {
     // Compute parameter values of the result sampling points
     RealArray gpar[2];
-    for (int dir = 0; dir < 2 && retVal; dir++)
-      retVal = this->getGridParameters(gpar[dir],dir,npe[dir]-1);
-
-    // Evaluate the secondary solution at all sampling points
-    if (retVal)
-      retVal = this->evalSolution(sField,integrand,gpar);
+    if (this->getGridParameters(gpar[0],0,npe[0]-1) &&
+	this->getGridParameters(gpar[1],1,npe[1]-1))
+      if (project)
+      {
+	// Project the secondary solution onto the spline basis
+	Go::SplineSurface* s = this->projectSolution(integrand);
+	if (s)
+	{
+	  // Evaluate the projected field at the result sampling points
+	  const Vector& svec = sField; // using utl::matrix cast operator
+	  sField.resize(s->dimension(),gpar[0].size()*gpar[1].size());
+	  s->gridEvaluator(const_cast<Vector&>(svec),gpar[0],gpar[1]);
+	  delete s;
+	  return true;
+	}
+      }
+      else
+	// Evaluate the secondary solution directly at all sampling points
+	return this->evalSolution(sField,integrand,gpar);
   }
   else
   {
@@ -1313,12 +1320,12 @@ bool ASMs2D::evalSolution (Matrix& sField, const Integrand& integrand,
       sField.resize(s->dimension(),s->numCoefs_u()*s->numCoefs_v());
       sField.fill(&(*s->coefs_begin()));
       delete s;
+      return true;
     }
-    else
-      retVal = false;
   }
 
-  return retVal;
+  std::cerr <<" *** ASMs2D::evalSolution: Failure!"<< std::endl;
+  return false;
 }
 
 
diff --git a/src/ASM/ASMs2D.h b/src/ASM/ASMs2D.h
index fd1f27f5..4f6e196e 100644
--- a/src/ASM/ASMs2D.h
+++ b/src/ASM/ASMs2D.h
@@ -235,14 +235,17 @@ public:
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods
   //! \param[in] npe Number of visualization nodes over each knot span
+  //! \param[in] project Flag indicating if the projected solution is wanted
   //!
   //! \details The secondary solution is derived from the primary solution,
   //! which is assumed to be stored within the \a integrand for current patch.
   //! If \a npe is NULL, the solution is evaluated at the Greville points and
   //! then projected onto the spline basis to obtain the control point values,
   //! which then are returned through \a sField.
+  //! If \a npe is not NULL and \a project is \e true, the solution is also
+  //! projected onto the spline basis, and then evaluated at the \a npe points
   virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
-			    const int* npe = 0) const;
+			    const int* npe = 0, bool project = false) const;
 
   //! \brief Projects the secondary solution field onto the primary basis.
   //! \param[in] integrand Object with problem-specific data and methods
diff --git a/src/ASM/ASMs2DLag.C b/src/ASM/ASMs2DLag.C
index 83a20823..cff1c666 100644
--- a/src/ASM/ASMs2DLag.C
+++ b/src/ASM/ASMs2DLag.C
@@ -93,7 +93,7 @@ bool ASMs2DLag::generateFEMTopology ()
   // Number of nodes per element
   const int nen = p1*p2;
 
-  // Connectivety array: local --> global node relation
+  // Connectivity array: local --> global node relation
   MLGE.resize(nel);
   MNPC.resize(nel);
 
@@ -453,7 +453,7 @@ bool ASMs2DLag::evalSolution (Matrix&, const Vector&,
 
 
 bool ASMs2DLag::evalSolution (Matrix& sField, const Integrand& integrand,
-			      const int*) const
+			      const int*, bool) const
 {
   sField.resize(0,0);
 
diff --git a/src/ASM/ASMs2DLag.h b/src/ASM/ASMs2DLag.h
index 6fa2c2a7..f17319db 100644
--- a/src/ASM/ASMs2DLag.h
+++ b/src/ASM/ASMs2DLag.h
@@ -112,7 +112,7 @@ public:
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods
   virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
-			    const int*) const;
+			    const int*, bool = false) const;
 
   //! \brief Evaluates the secondary solution field at the given points.
   //! \param[out] sField Solution field
diff --git a/src/ASM/ASMs2DSpec.C b/src/ASM/ASMs2DSpec.C
index 3b7c5556..5b0bfee3 100644
--- a/src/ASM/ASMs2DSpec.C
+++ b/src/ASM/ASMs2DSpec.C
@@ -264,7 +264,7 @@ bool ASMs2DSpec::integrate (Integrand& integrand, int lIndex,
 
 
 bool ASMs2DSpec::evalSolution (Matrix& sField, const Integrand& integrand,
-			       const int*) const
+			       const int*, bool) const
 {
   sField.resize(0,0);
 
diff --git a/src/ASM/ASMs2DSpec.h b/src/ASM/ASMs2DSpec.h
index 4ba53fa3..9ea877aa 100644
--- a/src/ASM/ASMs2DSpec.h
+++ b/src/ASM/ASMs2DSpec.h
@@ -65,9 +65,8 @@ public:
   //! \brief Evaluates the secondary solution field at all visualization points.
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods
-  //! \param[in] npe Number of visualization nodes over each knot span
   virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
-			    const int* npe) const;
+			    const int*, bool = false) const;
 
 protected:
 
diff --git a/src/ASM/ASMs2DmxLag.C b/src/ASM/ASMs2DmxLag.C
index 90c83471..039035d7 100644
--- a/src/ASM/ASMs2DmxLag.C
+++ b/src/ASM/ASMs2DmxLag.C
@@ -126,7 +126,7 @@ bool ASMs2DmxLag::generateFEMTopology ()
   const int nelx = (nx2-1)/(p1-1);
   const int nely = (ny2-1)/(p2-1);
 
-  // Add connectivety for second basis: local --> global node relation
+  // Add connectivity for second basis: local --> global node relation
   int i, j, a, b, iel;
   for (j = iel = 0; j < nely; j++)
     for (i = 0; i < nelx; i++, iel++)
@@ -414,7 +414,7 @@ bool ASMs2DmxLag::evalSolution (Matrix& sField, const Vector& locSol,
 
 
 bool ASMs2DmxLag::evalSolution (Matrix& sField, const Integrand& integrand,
-				const int*) const
+				const int*, bool) const
 {
   sField.resize(0,0);
 
diff --git a/src/ASM/ASMs2DmxLag.h b/src/ASM/ASMs2DmxLag.h
index 51944b2c..83d8838d 100644
--- a/src/ASM/ASMs2DmxLag.h
+++ b/src/ASM/ASMs2DmxLag.h
@@ -118,7 +118,7 @@ public:
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods
   virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
-			    const int*) const;
+			    const int*, bool = false) const;
 
   //! \brief Extracts nodal results for this patch from the global vector.
   //! \param[in] globVec Global solution vector in DOF-order
diff --git a/src/ASM/ASMs3D.C b/src/ASM/ASMs3D.C
index c4f65fbd..13c95e13 100644
--- a/src/ASM/ASMs3D.C
+++ b/src/ASM/ASMs3D.C
@@ -1821,20 +1821,33 @@ bool ASMs3D::evalSolution (Matrix& sField, const Vector& locSol,
 
 
 bool ASMs3D::evalSolution (Matrix& sField, const Integrand& integrand,
-			   const int* npe) const
+			   const int* npe, bool project) const
 {
-  bool retVal = true;
-
   if (npe)
   {
     // Compute parameter values of the result sampling points
     RealArray gpar[3];
-    for (int dir = 0; dir < 3 && retVal; dir++)
-      retVal = this->getGridParameters(gpar[dir],dir,npe[dir]-1);
-
-    // Evaluate the secondary solution at all sampling points
-    if (retVal)
-      retVal = this->evalSolution(sField,integrand,gpar);
+    if (this->getGridParameters(gpar[0],0,npe[0]-1) &&
+	this->getGridParameters(gpar[1],1,npe[1]-1) &&
+	this->getGridParameters(gpar[2],2,npe[2]-1))
+      if (project)
+      {
+	// Project the secondary solution onto the spline basis
+	Go::SplineVolume* v = this->projectSolution(integrand);
+	if (v)
+	{
+	  // Evaluate the projected field at the result sampling points
+	  const Vector& svec = sField; // using utl::matrix cast operator
+	  sField.resize(v->dimension(),
+			gpar[0].size()*gpar[1].size()*gpar[2].size());
+	  v->gridEvaluator(const_cast<Vector&>(svec),gpar[0],gpar[1],gpar[2]);
+	  delete v;
+	  return true;
+	}
+      }
+      else
+	// Evaluate the secondary solution at all sampling points
+	return this->evalSolution(sField,integrand,gpar);
   }
   else
   {
@@ -1847,12 +1860,12 @@ bool ASMs3D::evalSolution (Matrix& sField, const Integrand& integrand,
 		    v->numCoefs(0)*v->numCoefs(1)*v->numCoefs(2));
       sField.fill(&(*v->coefs_begin()));
       delete v;
+      return true;
     }
-    else
-      retVal = false;
   }
 
-  return retVal;
+  std::cerr <<" *** ASMs3D::evalSolution: Failure!"<< std::endl;
+  return false;
 }
 
 
diff --git a/src/ASM/ASMs3D.h b/src/ASM/ASMs3D.h
index 0fd01ea4..bcbcab56 100644
--- a/src/ASM/ASMs3D.h
+++ b/src/ASM/ASMs3D.h
@@ -300,14 +300,17 @@ public:
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods
   //! \param[in] npe Number of visualization nodes over each knot span
+  //! \param[in] project Flag indicating if the projected solution is wanted
   //!
   //! \details The secondary solution is derived from the primary solution,
   //! which is assumed to be stored within the \a integrand for current patch.
   //! If \a npe is NULL, the solution is evaluated at the Greville points and
   //! then projected onto the spline basis to obtain the control point values,
   //! which then are returned through \a sField.
+  //! If \a npe is not NULL and \a project is \e true, the solution is also
+  //! projected onto the spline basis, and then evaluated at the \a npe points
   virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
-			    const int* npe = 0) const;
+			    const int* npe = 0, bool project = false) const;
 
   //! \brief Projects the secondary solution field onto the primary basis.
   //! \param[in] integrand Object with problem-specific data and methods
diff --git a/src/ASM/ASMs3DLag.C b/src/ASM/ASMs3DLag.C
index a44e9b89..8af891c2 100644
--- a/src/ASM/ASMs3DLag.C
+++ b/src/ASM/ASMs3DLag.C
@@ -98,7 +98,7 @@ bool ASMs3DLag::generateFEMTopology()
   // Number of nodes in a xy-surface of an element
   const int ct  = p1*p2;
 
-  // Connectivety array: local --> global node relation
+  // Connectivity array: local --> global node relation
   MLGE.resize(nel);
   MNPC.resize(nel);
 
@@ -614,7 +614,7 @@ bool ASMs3DLag::evalSolution (Matrix&, const Vector&,
 
 
 bool ASMs3DLag::evalSolution (Matrix& sField, const Integrand& integrand,
-			      const int*) const
+			      const int*, bool) const
 {
   sField.resize(0,0);
 
diff --git a/src/ASM/ASMs3DLag.h b/src/ASM/ASMs3DLag.h
index 29cacb85..a299ef34 100644
--- a/src/ASM/ASMs3DLag.h
+++ b/src/ASM/ASMs3DLag.h
@@ -120,7 +120,7 @@ public:
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods
   virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
-			    const int*) const;
+			    const int*, bool = false) const;
 
   //! \brief Evaluates the secondary solution field at the given points.
   //! \param[out] sField Solution field
diff --git a/src/ASM/ASMs3DSpec.C b/src/ASM/ASMs3DSpec.C
index 992841f4..a8eb3b5e 100644
--- a/src/ASM/ASMs3DSpec.C
+++ b/src/ASM/ASMs3DSpec.C
@@ -408,7 +408,7 @@ bool ASMs3DSpec::integrateEdge (Integrand& integrand, int lEdge,
 
 
 bool ASMs3DSpec::evalSolution (Matrix& sField, const Integrand& integrand,
-			       const int*) const
+			       const int*, bool) const
 {
   sField.resize(0,0);
 
diff --git a/src/ASM/ASMs3DSpec.h b/src/ASM/ASMs3DSpec.h
index 5d104492..6f087d21 100644
--- a/src/ASM/ASMs3DSpec.h
+++ b/src/ASM/ASMs3DSpec.h
@@ -73,9 +73,8 @@ public:
   //! \brief Evaluates the secondary solution field at all visualization points.
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods
-  //! \param[in] npe Number of visualization nodes over each knot span
   virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
-			    const int* npe) const;
+			    const int*, bool = false) const;
 
 protected:
 
diff --git a/src/ASM/ASMs3DmxLag.C b/src/ASM/ASMs3DmxLag.C
index 5f1faacb..a5e3ed40 100644
--- a/src/ASM/ASMs3DmxLag.C
+++ b/src/ASM/ASMs3DmxLag.C
@@ -131,7 +131,7 @@ bool ASMs3DmxLag::generateFEMTopology ()
   const int nely = (ny2-1)/(p2-1);
   const int nelz = (nz2-1)/(p3-1);
 
-  // Add connectivety for second basis: local --> global node relation
+  // Add connectivity for second basis: local --> global node relation
   int i, j, k, a, b, c, iel;
   for (k = iel = 0; k < nelz; k++)
     for (j = 0; j < nely; j++)
@@ -451,7 +451,7 @@ bool ASMs3DmxLag::evalSolution (Matrix& sField, const Vector& locSol,
 
 
 bool ASMs3DmxLag::evalSolution (Matrix& sField, const Integrand& integrand,
-				const int*) const
+				const int*, bool) const
 {
   sField.resize(0,0);
 
diff --git a/src/ASM/ASMs3DmxLag.h b/src/ASM/ASMs3DmxLag.h
index 05d51625..fa9c69bf 100644
--- a/src/ASM/ASMs3DmxLag.h
+++ b/src/ASM/ASMs3DmxLag.h
@@ -118,7 +118,7 @@ public:
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods
   virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
-			    const int*) const;
+			    const int*, bool = false) const;
 
   //! \brief Extracts nodal results for this patch from the global vector.
   //! \param[in] globVec Global solution vector in DOF-order
diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C
index 95fd4b32..cdbcca0b 100644
--- a/src/SIM/SIMbase.C
+++ b/src/SIM/SIMbase.C
@@ -926,6 +926,7 @@ bool SIMbase::writeGlvBC (const int* nViz, int& nBlock) const
 
     if (msgLevel > 1)
       std::cout <<"Writing boundary conditions for patch "<< i+1 << std::endl;
+
     if (!myModel[i]->evalSolution(field,bc,nViz))
       return false;
 
@@ -980,7 +981,7 @@ bool SIMbase::writeGlvV (const Vector& vec, const char* fieldName,
 
   Matrix field;
   Vector lovec;
-  int geomID = 0, idBlock = 2;
+  int geomID = 0;
   std::vector<int> vID;
 
   for (size_t i = 0; i < myModel.size(); i++)
@@ -989,6 +990,7 @@ bool SIMbase::writeGlvV (const Vector& vec, const char* fieldName,
 
     if (msgLevel > 1)
       std::cout <<"Writing vector field for patch "<< i+1 << std::endl;
+
     myModel[i]->extractNodeVec(vec,lovec);
     if (!myModel[i]->evalSolution(field,lovec,nViz))
       return false;
@@ -999,6 +1001,7 @@ bool SIMbase::writeGlvV (const Vector& vec, const char* fieldName,
       vID.push_back(nBlock);
   }
 
+  int idBlock = 2;
   return myVtf->writeVblk(vID,fieldName,idBlock,iStep);
 }
 
@@ -1016,8 +1019,9 @@ bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
 
   Matrix field;
   size_t i, j;
-  int geomID = 0, idBlock = 10;
-  std::vector<int> vID, dID[14], sID[14];
+  int geomID = 0;
+  const size_t sMAX = 21;
+  std::vector<int> vID, dID[3], sID[sMAX];
   bool haveAsol = false;
   bool scalarEq = myModel.front()->getNoFields() == 1;
   if (mySol)
@@ -1032,6 +1036,9 @@ bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
 
     if (msgLevel > 1)
       std::cout <<"Writing FE solution for patch "<< i+1 << std::endl;
+
+    // 1. Evaluate primary solution variables
+
     myModel[i]->extractNodeVec(psol,myProblem->getSolution());
     if (!myModel[i]->evalSolution(field,myProblem->getSolution(),nViz))
       return false;
@@ -1059,6 +1066,8 @@ bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
 
     if (psolOnly) continue; // skip secondary solution
 
+    // 2. Direct evaluation of secondary solution variables
+
     LocalSystem::patch = i;
     if (!myModel[i]->evalSolution(field,*myProblem,nViz))
       return false;
@@ -1070,16 +1079,33 @@ bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
 	vID.push_back(nBlock);
 
     size_t k = 0;
-    for (j = 1; j <= field.rows() && k < 14; j++)
+    for (j = 1; j <= field.rows() && k < sMAX; j++)
       if (!myVtf->writeNres(field.getRow(j),++nBlock,geomID))
 	return false;
       else
 	sID[k++].push_back(nBlock);
 
+    if (discretization == Spline)
+    {
+      // 3. Projection of secondary solution variables (splines only)
+
+      if (!myModel[i]->evalSolution(field,*myProblem,nViz,true))
+	return false;
+
+      for (j = 1; j <= field.rows() && k < sMAX; j++)
+	if (!myVtf->writeNres(field.getRow(j),++nBlock,geomID))
+	  return false;
+	else
+	  sID[k++].push_back(nBlock);
+    }
+
     if (haveAsol)
     {
+      // 4. Evaluate analytical solution variables
+
       if (msgLevel > 1)
 	std::cout <<"Writing exact solution for patch "<< i+1 << std::endl;
+
       const ElementBlock* grid = myVtf->getBlock(geomID);
       std::vector<Vec3>::const_iterator cit = grid->begin_XYZ();
       Vector solPt; field.fill(0.0);
@@ -1095,7 +1121,8 @@ bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
 	if (ok)
 	  field.fillColumn(j,solPt);
       }
-      for (j = 1; j <= field.rows() && k < 14 && ok; j++)
+
+      for (j = 1; j <= field.rows() && k < sMAX && ok; j++)
 	if (!myVtf->writeNres(field.getRow(j),++nBlock,geomID))
 	  return false;
 	else
@@ -1103,15 +1130,17 @@ bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
     }
   }
 
+  // Write result block identifications
+
+  int idBlock = 10;
   if (scalarEq)
   {
+    if (!psolOnly)
+      if (!myVtf->writeVblk(vID,"Flux",idBlock,iStep))
+	return false;
+
     if (!myVtf->writeSblk(dID[0],"u",++idBlock,iStep))
       return false;
-
-    if (psolOnly) return true;
-
-    if (!myVtf->writeVblk(vID,"Flux",idBlock,iStep))
-      return false;
   }
   else
   {
@@ -1129,13 +1158,18 @@ bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
   if (psolOnly) return true;
 
   size_t nf = myProblem->getNoFields();
-  for (j = 0; j < nf && !sID[j].empty(); j++)
-    if (!myVtf->writeSblk(sID[j],myProblem->getFieldLabel(j,haveAsol?"FE":0),
+  for (i = j = 0; i < nf && !sID[j].empty(); i++, j++)
+    if (!myVtf->writeSblk(sID[j],myProblem->getFieldLabel(i,haveAsol?"FE":0),
 			  ++idBlock,iStep)) return false;
 
+  if (discretization == Spline)
+    for (i = 0; i < nf && !sID[j].empty(); i++, j++)
+      if (!myVtf->writeSblk(sID[j],myProblem->getFieldLabel(i,"Projected"),
+			    ++idBlock,iStep)) return false;
+
   if (haveAsol)
-    for (j = 0; j < nf && !sID[nf+j].empty(); j++)
-      if (!myVtf->writeSblk(sID[nf+j],myProblem->getFieldLabel(j,"Exact"),
+    for (i = 0; i < nf && !sID[j].empty(); i++, j++)
+      if (!myVtf->writeSblk(sID[j],myProblem->getFieldLabel(i,"Exact"),
 			    ++idBlock,iStep)) return false;
 
   return true;
@@ -1164,7 +1198,7 @@ bool SIMbase::writeGlvM (const Mode& mode, bool freq,
 
   Vector displ;
   Matrix field;
-  int geomID = 0, idBlock = 10;
+  int geomID = 0;
   std::vector<int> vID;
 
   for (size_t i = 0; i < myModel.size(); i++)
@@ -1184,6 +1218,7 @@ bool SIMbase::writeGlvM (const Mode& mode, bool freq,
   }
   if (myModel.size() > 1 && msgLevel > 1) std::cout << std::endl;
 
+  int idBlock = 10;
   if (!myVtf->writeDblk(vID,"Mode Shape",idBlock,mode.eigNo))
     return false;
 
@@ -1200,7 +1235,7 @@ bool SIMbase::writeGlvN (const Matrix& norms, int iStep, int& nBlock)
     return false;
 
   Matrix field;
-  int idBlock = 100, geomID = 0;
+  int geomID = 0;
   std::vector<int> sID[3];
 
   size_t i, j, k;
@@ -1208,9 +1243,10 @@ bool SIMbase::writeGlvN (const Matrix& norms, int iStep, int& nBlock)
   {
     if (myModel[i]->empty()) continue; // skip empty patches
 
-    geomID++;
     if (msgLevel > 1)
       std::cout <<"Writing element norms for patch "<< i+1 << std::endl;
+
+    geomID++;
     myModel[i]->extractElmRes(norms,field);
 
     for (j = k = 0; j < field.rows() && j < 4; j++)
@@ -1227,6 +1263,7 @@ bool SIMbase::writeGlvN (const Matrix& norms, int iStep, int& nBlock)
     "a(e,e)^0.5, e=u-u^h"
   };
 
+  int idBlock = 100;
   for (j = 0; j < 3; j++)
     if (!sID[j].empty())
       if (!myVtf->writeSblk(sID[j],label[j],++idBlock,iStep))
diff --git a/src/Utility/ElementBlock.C b/src/Utility/ElementBlock.C
index 9c7d93b2..675b4044 100644
--- a/src/Utility/ElementBlock.C
+++ b/src/Utility/ElementBlock.C
@@ -1,4 +1,4 @@
-// $Id: ElementBlock.C,v 1.4 2010-09-05 11:29:18 kmo Exp $
+// $Id$
 //==============================================================================
 //!
 //! \file ElementBlock.C
@@ -53,6 +53,15 @@ bool ElementBlock::setCoor (size_t i, real x, real y, real z)
 }
 
 
+bool ElementBlock::setCoor (size_t i, size_t j, real x)
+{
+  if (i >= coord.size() || j >= 3) return false;
+
+  coord[i][j] = x;
+  return true;
+}
+
+
 bool ElementBlock::setNode (size_t i, int nodeNumb)
 {
   if (i >= MMNPC.size()) return false;
diff --git a/src/Utility/ElementBlock.h b/src/Utility/ElementBlock.h
index 8805b2a2..82aed01e 100644
--- a/src/Utility/ElementBlock.h
+++ b/src/Utility/ElementBlock.h
@@ -1,4 +1,4 @@
-// $Id: ElementBlock.h,v 1.2 2010-03-10 13:59:02 kmo Exp $
+// $Id$
 //==============================================================================
 //!
 //! \file ElementBlock.h
@@ -35,6 +35,8 @@ public:
 
   //! \brief Defines the coordinates of node \a i
   bool setCoor(size_t i, real x, real y, real z);
+  //! \brief Defines the \a j'th coordinate of node \a i
+  bool setCoor(size_t i, size_t j, real x);
 
   //! \brief Defines the global number of element node \a i
   bool setNode(size_t i, int nodeNumb);