diff --git a/src/ASM/ASMbase.C b/src/ASM/ASMbase.C
index db5a1fbc..18b71e20 100644
--- a/src/ASM/ASMbase.C
+++ b/src/ASM/ASMbase.C
@@ -1,4 +1,4 @@
-// $Id: ASMbase.C,v 1.13 2010-12-29 18:02:10 kmo Exp $
+// $Id$
 //==============================================================================
 //!
 //! \file ASMbase.C
@@ -509,9 +509,27 @@ bool ASMbase::evalSolution (Matrix&, const Vector&, const int*) const
 }
 
 
+bool ASMbase::evalSolution (Matrix&, const Vector&,
+			    const RealArray*, bool) const
+{
+  std::cerr <<" *** ASMBase::evalSolution: Must be implemented in sub-class."
+	    << std::endl;
+  return false;
+}
+
+
 bool ASMbase::evalSolution (Matrix&, const Integrand&, const int*) const
 {
   std::cerr <<" *** ASMBase::evalSolution: Must be implemented in sub-class."
 	    << std::endl;
   return false;
 }
+
+
+bool ASMbase::evalSolution (Matrix&, const Integrand&,
+			    const RealArray*, bool) const
+{
+  std::cerr <<" *** ASMBase::evalSolution: Must be implemented in sub-class."
+	    << std::endl;
+  return false;
+}
diff --git a/src/ASM/ASMbase.h b/src/ASM/ASMbase.h
index aef8575c..915359d2 100644
--- a/src/ASM/ASMbase.h
+++ b/src/ASM/ASMbase.h
@@ -1,4 +1,4 @@
-// $Id: ASMbase.h,v 1.19 2011-01-28 13:34:28 kmo Exp $
+// $Id$
 //==============================================================================
 //!
 //! \file ASMbase.h
@@ -238,6 +238,20 @@ public:
   virtual bool evalSolution(Matrix& sField, const Vector& locSol,
 			    const int* npe) const;
 
+  //! \brief Evaluates the primary solution field at the given points.
+  //! \param[out] sField Solution field
+  //! \param[in] locSol Solution vector local to current patch
+  //! \param[in] gpar Parameter values of the result sampling points
+  //! \param[in] regular Flag indicating how the sampling points are defined
+  //!
+  //! \details When \a regular is \e true, it is assumed that the parameter
+  //! value array \a gpar forms a regular tensor-product point grid of dimension
+  //! \a gpar[0].size() \a X \a gpar[1].size() \a X \a gpar[2].size().
+  //! Otherwise, we assume that it contains the \a u, \a v and \a w parameters
+  //! directly for each sampling point.
+  virtual bool evalSolution(Matrix& sField, const Vector& locSol,
+			    const RealArray* gpar, bool regular = true) const;
+
   //! \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
@@ -251,6 +265,20 @@ public:
   virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
 			    const int* npe = 0) const;
 
+  //! \brief Evaluates the secondary solution field at the given points.
+  //! \param[out] sField Solution field
+  //! \param[in] integrand Object with problem-specific data and methods
+  //! \param[in] gpar Parameter values of the result sampling points
+  //! \param[in] regular Flag indicating how the sampling points are defined
+  //!
+  //! \details When \a regular is \e true, it is assumed that the parameter
+  //! value array \a gpar forms a regular tensor-product point grid of dimension
+  //! \a gpar[0].size() \a X \a gpar[1].size() \a X \a gpar[2].size().
+  //! Otherwise, we assume that it contains the \a u, \a v and \a w parameters
+  //! directly for each sampling point.
+  virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
+			    const RealArray* gpar, bool regular = true) const;
+
 
   // Methods for result extraction
   // =============================
diff --git a/src/ASM/ASMs1D.C b/src/ASM/ASMs1D.C
index 1c5e1e20..8e3d6324 100644
--- a/src/ASM/ASMs1D.C
+++ b/src/ASM/ASMs1D.C
@@ -27,9 +27,6 @@
 #include <ctype.h>
 #include <fstream>
 
-typedef std::vector<double>       DoubleVec;  //!< 1D double array
-typedef DoubleVec::const_iterator DoubleIter; //!< Iterator over DoubleVec
-
 
 ASMs1D::ASMs1D (const char* fileName, unsigned char n_s, unsigned char n_f)
   : ASMstruct(1,n_s,n_f)
@@ -137,8 +134,8 @@ bool ASMs1D::refine (const RealArray& xi)
   if (!curv || xi.empty()) return false;
   if (xi.front() < 0.0 || xi.back() > 1.0) return false;
 
-  DoubleVec  extraKnots;
-  DoubleIter uit = curv->basis().begin();
+  RealArray extraKnots;
+  RealArray::const_iterator uit = curv->basis().begin();
   double ucurr, uprev = *(uit++);
   while (uit != curv->basis().end())
   {
@@ -162,8 +159,8 @@ bool ASMs1D::uniformRefine (int nInsert)
 {
   if (!curv || nInsert < 1) return false;
 
-  DoubleVec  extraKnots;
-  DoubleIter uit = curv->basis().begin();
+  RealArray extraKnots;
+  RealArray::const_iterator uit = curv->basis().begin();
   double ucurr, uprev = *(uit++);
   while (uit != curv->basis().end())
   {
@@ -360,7 +357,7 @@ double ASMs1D::getKnotSpan (int i) const
   if (i < 0 || i >= curv->numCoefs() + curv->order() - 1)
     return 0.0;
 
-  DoubleIter uit = curv->basis().begin() + i;
+  RealArray::const_iterator uit = curv->basis().begin() + i;
   return *(uit+1) - *uit;
 }
 
@@ -379,7 +376,7 @@ bool ASMs1D::getElementCoordinates (Matrix& X, int iel) const
   const IntVec& mnpc = MNPC[iel-1];
   X.resize(nsd,mnpc.size());
 
-  DoubleIter cit = curv->coefs_begin();
+  RealArray::const_iterator cit = curv->coefs_begin();
   for (size_t n = 0; n < mnpc.size(); n++)
   {
     int ip = mnpc[n]*curv->dimension();
@@ -402,7 +399,7 @@ void ASMs1D::getNodalCoordinates (Matrix& X) const
 
   X.resize(nsd,n1);
 
-  DoubleIter cit = curv->coefs_begin();
+  RealArray::const_iterator cit = curv->coefs_begin();
   for (int inod = 0; inod < n1; inod++)
   {
     int ip = inod*curv->dimension();
@@ -418,7 +415,7 @@ Vec3 ASMs1D::getCoord (size_t inod) const
   int ip = (inod-1)*curv->dimension();
   if (ip < 0) return X;
 
-  DoubleIter cit = curv->coefs_begin() + ip;
+  RealArray::const_iterator cit = curv->coefs_begin() + ip;
   for (size_t i = 0; i < nsd; i++, cit++)
     X[i] = *cit;
 
@@ -438,7 +435,7 @@ void ASMs1D::extractBasis (double u, Vector& N, Matrix& dNdu) const
 {
   int p1 = curv->order();
 
-  DoubleVec bas(p1*2);
+  RealArray bas(p1*2);
   curv->basis().computeBasisValues(u,&bas.front(),1);
 
   N.resize(p1);
@@ -471,7 +468,7 @@ bool ASMs1D::integrate (Integrand& integrand,
   int pm1 = p1 - 1;
   int nCol = n1 - pm1;
   Matrix gpar(nGauss,nCol);
-  DoubleIter uit = curv->basis().begin() + pm1;
+  RealArray::const_iterator uit = curv->basis().begin() + pm1;
   double ucurr, uprev = *(uit++);
   for (int j = 1; j <= nCol; uit++, j++)
   {
@@ -551,18 +548,18 @@ bool ASMs1D::integrate (Integrand& integrand, int lIndex,
 
   // Integration of boundary point
 
-  int    iel;
-  double param;
+  FiniteElement fe;
+  int iel;
   switch (lIndex)
     {
     case 1:
+      fe.u = curv->startparam();
       iel = 1;
-      param = curv->startparam();
       break;
 
     case 2:
+      fe.u = curv->endparam();
       iel = this->getNoElms();
-      param = curv->endparam();
       break;
 
     default:
@@ -585,9 +582,8 @@ bool ASMs1D::integrate (Integrand& integrand, int lIndex,
   LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[MLGE[iel-1]-1];
 
   // Evaluate basis functions and corresponding derivatives
-  FiniteElement fe;
   Matrix dNdu;
-  this->extractBasis(param,fe.N,dNdu);
+  this->extractBasis(fe.u,fe.N,dNdu);
 
   // Cartesian coordinates of current integration point
   Vec4 X(Xnod*fe.N,time.t);
@@ -623,7 +619,7 @@ bool ASMs1D::getGridParameters (RealArray& prm, int nSegPerSpan) const
     return false;
   }
 
-  DoubleIter uit = curv->basis().begin();
+  RealArray::const_iterator uit = curv->basis().begin();
   double ucurr = 0.0, uprev = *(uit++);
   while (uit != curv->basis().end())
   {
@@ -648,13 +644,13 @@ bool ASMs1D::getGridParameters (RealArray& prm, int nSegPerSpan) const
 bool ASMs1D::tesselate (ElementBlock& grid, const int* npe) const
 {
   // Compute parameter values of the nodal points
-  DoubleVec gpar;
+  RealArray gpar;
   if (!this->getGridParameters(gpar,npe[0]-1))
     return false;
 
   // Evaluate the spline curve at all points
   size_t nx = gpar.size();
-  DoubleVec XYZ(3*nx,0.0);
+  RealArray XYZ(3*nx,0.0);
 
   // Establish the block grid coordinates
   size_t i, j;
@@ -699,10 +695,18 @@ bool ASMs1D::evalSolution (Matrix& sField, const Vector& locSol,
 			   const int* npe) const
 {
   // Compute parameter values of the result sampling points
-  DoubleVec gpar;
+  RealArray gpar;
   if (!this->getGridParameters(gpar,npe[0]-1))
     return false;
 
+  // Evaluate the primary solution at all sampling points
+  return this->evalSolution(sField,locSol,&gpar);
+}
+
+
+bool ASMs1D::evalSolution (Matrix& sField, const Vector& locSol,
+			   const RealArray* gpar, bool) const
+{
   const int p1 = curv->order();
   size_t nComp = locSol.size() / this->getNoNodes();
   if (nComp*this->getNoNodes() != locSol.size())
@@ -712,11 +716,12 @@ bool ASMs1D::evalSolution (Matrix& sField, const Vector& locSol,
   Vector Ytmp, basis(p1);
 
   // Evaluate the primary solution field at each point
-  size_t nPoints = gpar.size();
+  const RealArray& upar = *gpar;
+  size_t nPoints = upar.size();
   sField.resize(nComp,nPoints);
   for (size_t i = 0; i < nPoints; i++)
   {
-    curv->basis().computeBasisValues(gpar[i],&basis.front());
+    curv->basis().computeBasisValues(upar[i],&basis.front());
 
     IntVec ip;
     scatterInd(p1,curv->basis().lastKnotInterval(),ip);
@@ -739,13 +744,15 @@ bool ASMs1D::evalSolution (Matrix& sField, const Integrand& integrand,
     RealArray gpar;
     if (this->getGridParameters(gpar,npe[0]-1))
       // Evaluate the secondary solution at all sampling points
-      return this->evalSolution(sField,integrand,gpar);
+      return this->evalSolution(sField,integrand,&gpar);
   }
   else
   {
+    // Project the secondary solution onto the spline basis
     Go::SplineCurve* c = this->projectSolution(integrand);
     if (c)
     {
+      // Extract control point values from the spline object
       sField.resize(c->dimension(),c->numCoefs());
       sField.fill(&(*c->coefs_begin()));
       delete c;
@@ -771,7 +778,7 @@ Go::SplineCurve* ASMs1D::projectSolution (const Integrand& integrand) const
 
 
 bool ASMs1D::evalSolution (Matrix& sField, const Integrand& integrand,
-			   const RealArray& gpar) const
+			   const RealArray* gpar, bool) const
 {
   sField.resize(0,0);
 
@@ -785,11 +792,12 @@ bool ASMs1D::evalSolution (Matrix& sField, const Integrand& integrand,
   Matrix dNdu, dNdX, Jac;
 
   // Evaluate the secondary solution field at each point
-  size_t nPoints = gpar.size();
+  const RealArray& upar = *gpar;
+  size_t nPoints = upar.size();
   for (size_t i = 0; i < nPoints; i++)
   {
     // Fetch basis function derivatives at current integration point
-    this->extractBasis(gpar[i],N,dNdu);
+    this->extractBasis(upar[i],N,dNdu);
 
     // Fetch indices of the non-zero basis functions at this point
     IntVec ip;
diff --git a/src/ASM/ASMs1D.h b/src/ASM/ASMs1D.h
index 6b50f2fd..031138db 100644
--- a/src/ASM/ASMs1D.h
+++ b/src/ASM/ASMs1D.h
@@ -1,4 +1,4 @@
-// $Id: ASMs1D.h,v 1.4 2011-01-05 12:49:41 kmo Exp $
+// $Id$
 //==============================================================================
 //!
 //! \file ASMs1D.h
@@ -135,6 +135,13 @@ public:
   virtual bool evalSolution(Matrix& sField, const Vector& locSol,
 			    const int* npe) const;
 
+  //! \brief Evaluates the primary solution field at the given points.
+  //! \param[out] sField Solution field
+  //! \param[in] locSol Solution vector local to current patch
+  //! \param[in] gpar Parameter values of the result sampling points
+  virtual bool evalSolution(Matrix& sField, const Vector& locSol,
+			    const RealArray* gpar, bool = true) const;
+
   //! \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
@@ -150,19 +157,20 @@ public:
   //! \param[in] integrand Object with problem-specific data and methods
   Go::SplineCurve* projectSolution(const Integrand& integrand) const;
   //! \brief Projects the secondary solution field onto the primary basis.
+  //! \param[in] integrand Object with problem-specific data and methods
   virtual Go::GeomObject* evalSolution(const Integrand& integrand) const;
 
-protected:
-
-  // Internal utility methods
-  // ========================
-
   //! \brief Evaluates the secondary solution field at the given points.
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods
   //! \param[in] gpar Parameter values of the result sampling points
-  bool evalSolution(Matrix& sField, const Integrand& integrand,
-		    const RealArray& gpar) const;
+  virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
+			    const RealArray* gpar, bool = true) const;
+
+protected:
+
+  // Internal utility methods
+  // ========================
 
   //! \brief Calculates parameter values for the visualization nodal points.
   //! \param[out] prm Parameter values for all points
diff --git a/src/ASM/ASMs1DLag.C b/src/ASM/ASMs1DLag.C
index 84f10792..ac28374a 100644
--- a/src/ASM/ASMs1DLag.C
+++ b/src/ASM/ASMs1DLag.C
@@ -220,28 +220,26 @@ bool ASMs1DLag::integrate (Integrand& integrand, int lIndex,
 {
   if (!curv) return true; // silently ignore empty patches
 
-  // Number of elements
-  const int p1  = curv->order();
-  const int nel = (nx-1)/(p1-1);
-
   // Integration of boundary point
 
-  int iel;
+  FiniteElement fe;
   double x;
+  int iel;
   switch (lIndex)
     {
     case 1:
-      iel = 1;
+      fe.u = curv->startparam();
       x = -1.0;
+      iel = 1;
       break;
 
     case 2:
-      iel = nel;
+      fe.u = curv->endparam();
       x = 1.0;
+      iel = this->getNoElms();
 
     default:
-      iel = 0;
-      x = 0.0;
+      return false;
     }
 
   // Set up control point coordinates for current element
@@ -258,9 +256,8 @@ bool ASMs1DLag::integrate (Integrand& integrand, int lIndex,
   LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[MLGE[iel-1]-1];
 
   // Evaluate basis functions and corresponding derivatives
-  FiniteElement fe;
   Matrix dNdu, Jac;
-  if (!Lagrange::computeBasis(fe.N,dNdu,p1,x)) return false;
+  if (!Lagrange::computeBasis(fe.N,dNdu,curv->order(),x)) return false;
 
   // Cartesian coordinates of current integration point
   Vec4 X(Xnod*fe.N,time.t);
@@ -323,6 +320,15 @@ bool ASMs1DLag::evalSolution (Matrix& sField, const Vector& locSol,
 }
 
 
+bool ASMs1DLag::evalSolution (Matrix&, const Vector&,
+			      const RealArray*, bool) const
+{
+  std::cerr <<" *** ASMs1DLag::evalSolution(Matrix&,const Vector&,"
+	    <<"const RealArray*,bool): Not implemented."<< std::endl;
+  return false;
+}
+
+
 bool ASMs1DLag::evalSolution (Matrix& sField, const Integrand& integrand,
 			      const int*) const
 {
@@ -373,3 +379,12 @@ bool ASMs1DLag::evalSolution (Matrix& sField, const Integrand& integrand,
 
   return true;
 }
+
+
+bool ASMs1DLag::evalSolution (Matrix&, const Integrand&,
+			      const RealArray*, bool) const
+{
+  std::cerr <<" *** ASMs1DLag::evalSolution(Matrix&,const Integrand&,"
+	    <<"const RealArray*,bool): Not implemented."<< std::endl;
+  return false;
+}
diff --git a/src/ASM/ASMs1DLag.h b/src/ASM/ASMs1DLag.h
index a090b138..51015a46 100644
--- a/src/ASM/ASMs1DLag.h
+++ b/src/ASM/ASMs1DLag.h
@@ -1,4 +1,4 @@
-// $Id: ASMs1DLag.h,v 1.2 2010-09-07 07:51:46 kmo Exp $
+// $Id$
 //==============================================================================
 //!
 //! \file ASMs1DLag.h
@@ -79,29 +79,40 @@ public:
   // Post-processing methods
   // =======================
 
-  //! \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
+  //! \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.
-  virtual bool tesselate(ElementBlock& grid, const int* npe) const;
+  virtual bool tesselate(ElementBlock& grid, const int*) const;
 
   //! \brief Evaluates the primary solution field at all visualization points.
   //! \details The number of visualization points is the same as the order of
   //! the Lagrange elements by default.
   //! \param[out] sField Solution field
   //! \param[in] locSol Solution vector in DOF-order
-  //! \param[in] npe Number of visualization nodes over each knot span
   virtual bool evalSolution(Matrix& sField, const Vector& locSol,
-			    const int* npe) const;
+			    const int*) const;
+
+  //! \brief Evaluates the primary solution field at the given points.
+  //! \param[out] sField Solution field
+  //! \param[in] locSol Solution vector local to current patch
+  //! \param[in] gpar Parameter values of the result sampling points
+  virtual bool evalSolution(Matrix& sField, const Vector& locSol,
+			    const RealArray* gpar, bool = true) const;
 
   //! \brief Evaluates the secondary solution field at all visualization points.
   //! \details The number of visualization points is the same as the order of
   //! the Lagrange elements by default.
   //! \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*) const;
+
+  //! \brief Evaluates the secondary solution field at the given points.
+  //! \param[out] sField Solution field
+  //! \param[in] integrand Object with problem-specific data and methods
+  //! \param[in] gpar Parameter values of the result sampling points
+  virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
+			    const RealArray* gpar, bool = true) const;
 
 protected:
 
diff --git a/src/ASM/ASMs2D.C b/src/ASM/ASMs2D.C
index 8b122612..89c598ad 100644
--- a/src/ASM/ASMs2D.C
+++ b/src/ASM/ASMs2D.C
@@ -29,9 +29,6 @@
 #include <ctype.h>
 #include <fstream>
 
-typedef Go::SplineSurface::Dvector DoubleVec; //!< 1D double array
-typedef DoubleVec::const_iterator DoubleIter; //!< Iterator over DoubleVec
-
 
 ASMs2D::ASMs2D (const char* fileName, unsigned char n_s, unsigned char n_f)
   : ASMstruct(2,n_s,n_f)
@@ -139,8 +136,8 @@ bool ASMs2D::refine (int dir, const RealArray& xi)
   if (!surf || dir < 0 || dir > 1 || xi.empty()) return false;
   if (xi.front() < 0.0 || xi.back() > 1.0) return false;
 
-  DoubleVec  extraKnots;
-  DoubleIter uit = surf->basis(dir).begin();
+  RealArray extraKnots;
+  RealArray::const_iterator uit = surf->basis(dir).begin();
   double ucurr, uprev = *(uit++);
   while (uit != surf->basis(dir).end())
   {
@@ -168,8 +165,8 @@ bool ASMs2D::uniformRefine (int dir, int nInsert)
 {
   if (!surf || dir < 0 || dir > 1 || nInsert < 1) return false;
 
-  DoubleVec  extraKnots;
-  DoubleIter uit = surf->basis(dir).begin();
+  RealArray extraKnots;
+  RealArray::const_iterator uit = surf->basis(dir).begin();
   double ucurr, uprev = *(uit++);
   while (uit != surf->basis(dir).end())
   {
@@ -620,7 +617,7 @@ bool ASMs2D::getElementCoordinates (Matrix& X, int iel) const
   const IntVec& mnpc = MNPC[iel-1];
   X.resize(nsd,mnpc.size());
 
-  DoubleIter cit = surf->coefs_begin();
+  RealArray::const_iterator cit = surf->coefs_begin();
   for (size_t n = 0; n < mnpc.size(); n++)
   {
     int ip = this->coeffInd(mnpc[n])*surf->dimension();
@@ -643,7 +640,7 @@ void ASMs2D::getNodalCoordinates (Matrix& X) const
   const int n2 = surf->numCoefs_v();
   X.resize(nsd,n1*n2);
 
-  DoubleIter cit = surf->coefs_begin();
+  RealArray::const_iterator cit = surf->coefs_begin();
   size_t inod = 1;
   for (int i2 = 0; i2 < n2; i2++)
     for (int i1 = 0; i1 < n1; i1++, inod++)
@@ -661,7 +658,7 @@ Vec3 ASMs2D::getCoord (size_t inod) const
   int ip = this->coeffInd(inod-1)*surf->dimension();
   if (ip < 0) return X;
 
-  DoubleIter cit = surf->coefs_begin() + ip;
+  RealArray::const_iterator cit = surf->coefs_begin() + ip;
   for (size_t i = 0; i < nsd; i++, cit++)
     X[i] = *cit;
 
@@ -788,7 +785,7 @@ bool ASMs2D::integrate (Integrand& integrand,
   for (dir = 0; dir < 2; dir++)
   {
     int pm1 = (dir == 0 ? surf->order_u() : surf->order_v()) - 1;
-    DoubleIter uit = surf->basis(dir).begin() + pm1;
+    RealArray::const_iterator uit = surf->basis(dir).begin() + pm1;
     double ucurr, uprev = *(uit++);
     int nCol = (dir == 0 ? surf->numCoefs_u() : surf->numCoefs_v()) - pm1;
     gpar[dir].resize(nGauss,nCol);
@@ -982,7 +979,7 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex,
     else
     {
       int pm1 = (d == 0 ? surf->order_u() : surf->order_v()) - 1;
-      DoubleIter uit = surf->basis(d).begin() + pm1;
+      RealArray::const_iterator uit = surf->basis(d).begin() + pm1;
       double ucurr, uprev = *(uit++);
       int nCol = (d == 0 ? surf->numCoefs_u() : surf->numCoefs_v()) - pm1;
       gpar[d].resize(nGauss,nCol);
@@ -1005,6 +1002,9 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex,
   const int n2 = surf->numCoefs_v();
 
   FiniteElement fe(p1*p2);
+  fe.u = gpar[0](1,1);
+  fe.v = gpar[1](1,1);
+
   Matrix dNdu, Xnod, Jac;
   Vec4   X;
   Vec3   normal;
@@ -1030,7 +1030,7 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex,
       if (skipMe) continue;
 
       // Get element edge length in the parameter space
-      double dS = this->getParametricLength(iel,abs(edgeDir));
+      double dS = this->getParametricLength(iel,t1);
       if (dS < 0.0) return false; // topology error (probably logic error)
 
       // Set up control point coordinates for current element
@@ -1048,9 +1048,13 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex,
 
       // --- Integration loop over all Gauss points along the edge -------------
 
-      int ip = (abs(edgeDir) == 1 ? i2-p2 : i1-p1)*nGauss;
+      int ip = (t1 == 1 ? i2-p2 : i1-p1)*nGauss;
       for (int i = 0; i < nGauss; i++, ip++)
       {
+	// Parameter values of current integration point
+	if (gpar[0].size() > 1) fe.u = gpar[0](i+1,i1-p1+1);
+	if (gpar[1].size() > 1) fe.v = gpar[1](i+1,i2-p2+1);
+
 	// Fetch basis function derivatives at current integration point
 	extractBasis(spline[ip],fe.N,dNdu);
 
@@ -1090,7 +1094,7 @@ bool ASMs2D::getGridParameters (RealArray& prm, int dir, int nSegPerSpan) const
     return false;
   }
 
-  DoubleIter uit = surf->basis(dir).begin();
+  RealArray::const_iterator uit = surf->basis(dir).begin();
   double ucurr = 0.0, uprev = *(uit++);
   while (uit != surf->basis(dir).end())
   {
@@ -1129,7 +1133,7 @@ bool ASMs2D::getGrevilleParameters (RealArray& prm, int dir) const
 bool ASMs2D::tesselate (ElementBlock& grid, const int* npe) const
 {
   // Compute parameter values of the nodal points
-  DoubleVec gpar[2];
+  RealArray gpar[2];
   for (int dir = 0; dir < 2; dir++)
     if (!this->getGridParameters(gpar[dir],dir,npe[dir]-1))
       return false;
@@ -1137,7 +1141,7 @@ bool ASMs2D::tesselate (ElementBlock& grid, const int* npe) const
   // Evaluate the spline surface at all points
   size_t nx = gpar[0].size();
   size_t ny = gpar[1].size();
-  DoubleVec XYZ(surf->dimension()*nx*ny);
+  RealArray XYZ(surf->dimension()*nx*ny);
   surf->gridEvaluator(XYZ,gpar[0],gpar[1]);
 
   // Establish the block grid coordinates
@@ -1183,14 +1187,45 @@ bool ASMs2D::evalSolution (Matrix& sField, const Vector& locSol,
 			   const int* npe) const
 {
   // Compute parameter values of the result sampling points
-  DoubleVec gpar[2];
+  RealArray gpar[2];
   for (int dir = 0; dir < 2; dir++)
     if (!this->getGridParameters(gpar[dir],dir,npe[dir]-1))
       return false;
 
+  // Evaluate the primary solution at all sampling points
+  return this->evalSolution(sField,locSol,gpar);
+}
+
+
+bool ASMs2D::evalSolution (Matrix& sField, const Vector& locSol,
+			   const RealArray* gpar, bool regular) const
+{
   // Evaluate the basis functions at all points
   std::vector<Go::BasisPtsSf> spline;
-  surf->computeBasisGrid(gpar[0],gpar[1],spline);
+  if (regular)
+    surf->computeBasisGrid(gpar[0],gpar[1],spline);
+  else if (gpar[0].size() == gpar[1].size())
+  {
+    spline.resize(gpar[0].size());
+    std::vector<Go::BasisPtsSf> tmpSpline(1);
+    for (size_t i = 0; i < spline.size(); i++)
+    {
+      surf->computeBasisGrid(RealArray(1,gpar[0][i]),
+			     RealArray(1,gpar[1][i]),
+			     tmpSpline);
+      spline[i] = tmpSpline.front();
+    }
+    // TODO: Request a GoTools method replacing the above:
+    // void SplineSurface::computeBasisGrid(double param_u, double param_v,
+    //                                      BasisPtsSf& result) const
+    /*
+    spline.resize(gpar[0].size());
+    for (size_t i = 0; i < spline.size(); i++)
+      surf->computeBasis(gpar[0][i],gpar[1][i],spline[i]);
+    */
+  }
+  else
+    return false;
 
   const int p1 = surf->order_u();
   const int p2 = surf->order_v();
@@ -1238,9 +1273,11 @@ bool ASMs2D::evalSolution (Matrix& sField, const Integrand& integrand,
   }
   else
   {
+    // Project the secondary solution onto the spline basis
     Go::SplineSurface* s = this->projectSolution(integrand);
     if (s)
     {
+      // Extract control point values from the spline object
       sField.resize(s->dimension(),s->numCoefs_u()*s->numCoefs_v());
       sField.fill(&(*s->coefs_begin()));
       delete s;
@@ -1262,7 +1299,7 @@ Go::GeomObject* ASMs2D::evalSolution (const Integrand& integrand) const
 Go::SplineSurface* ASMs2D::projectSolution (const Integrand& integrand) const
 {
   // Compute parameter values of the result sampling points (Greville points)
-  DoubleVec gpar[2];
+  RealArray gpar[2];
   for (int dir = 0; dir < 2; dir++)
     if (!this->getGrevilleParameters(gpar[dir],dir))
       return 0;
@@ -1279,7 +1316,7 @@ Go::SplineSurface* ASMs2D::projectSolution (const Integrand& integrand) const
   // the result array. Think that is always the case, but beware if trying
   // other projection schemes later.
 
-  DoubleVec weights;
+  RealArray weights;
   if (surf->rational())
     surf->getWeights(weights);
 
@@ -1295,13 +1332,36 @@ Go::SplineSurface* ASMs2D::projectSolution (const Integrand& integrand) const
 
 
 bool ASMs2D::evalSolution (Matrix& sField, const Integrand& integrand,
-			   const RealArray* gpar) const
+			   const RealArray* gpar, bool regular) const
 {
   sField.resize(0,0);
 
   // Evaluate the basis functions and their derivatives at all points
   std::vector<Go::BasisDerivsSf> spline;
-  surf->computeBasisGrid(gpar[0],gpar[1],spline);
+  if (regular)
+    surf->computeBasisGrid(gpar[0],gpar[1],spline);
+  else if (gpar[0].size() == gpar[1].size())
+  {
+    spline.resize(gpar[0].size());
+    std::vector<Go::BasisDerivsSf> tmpSpline(1);
+    for (size_t i = 0; i < spline.size(); i++)
+    {
+      surf->computeBasisGrid(RealArray(1,gpar[0][i]),
+			     RealArray(1,gpar[1][i]),
+			     tmpSpline);
+      spline[i] = tmpSpline.front();
+    }
+    // TODO: Request a GoTools method replacing the above:
+    // void SplineSurface::computeBasisGrid(double param_u, double param_v,
+    //                                      BasisDerivsSf& result) const
+    /*
+    spline.resize(gpar[0].size());
+    for (size_t i = 0; i < spline.size(); i++)
+      surf->computeBasis(gpar[0][i],gpar[1][i],spline[i]);
+    */
+  }
+  else
+    return false;
 
   const int p1 = surf->order_u();
   const int p2 = surf->order_v();
diff --git a/src/ASM/ASMs2D.h b/src/ASM/ASMs2D.h
index 9056310a..4aa89db4 100644
--- a/src/ASM/ASMs2D.h
+++ b/src/ASM/ASMs2D.h
@@ -1,4 +1,4 @@
-// $Id: ASMs2D.h,v 1.11 2011-01-28 13:34:28 kmo Exp $
+// $Id$
 //==============================================================================
 //!
 //! \file ASMs2D.h
@@ -208,6 +208,20 @@ public:
   virtual bool evalSolution(Matrix& sField, const Vector& locSol,
 			    const int* npe) const;
 
+  //! \brief Evaluates the primary solution field at the given points.
+  //! \param[out] sField Solution field
+  //! \param[in] locSol Solution vector local to current patch
+  //! \param[in] gpar Parameter values of the result sampling points
+  //! \param[in] regular Flag indicating how the sampling points are defined
+  //!
+  //! \details When \a regular is \e true, it is assumed that the parameter
+  //! value array \a gpar forms a regular tensor-product point grid of dimension
+  //! \a gpar[0].size() \a X \a gpar[1].size().
+  //! Otherwise, we assume that it contains the \a u and \a v parameters
+  //! directly for each sampling point.
+  virtual bool evalSolution(Matrix& sField, const Vector& locSol,
+			    const RealArray* gpar, bool regular = true) const;
+
   //! \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
@@ -223,19 +237,27 @@ public:
   //! \param[in] integrand Object with problem-specific data and methods
   Go::SplineSurface* projectSolution(const Integrand& integrand) const;
   //! \brief Projects the secondary solution field onto the primary basis.
+  //! \param[in] integrand Object with problem-specific data and methods
   virtual Go::GeomObject* evalSolution(const Integrand& integrand) const;
 
-protected:
-
-  // Internal utility methods
-  // ========================
-
   //! \brief Evaluates the secondary solution field at the given points.
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods
   //! \param[in] gpar Parameter values of the result sampling points
-  bool evalSolution(Matrix& sField, const Integrand& integrand,
-		    const RealArray* gpar) const;
+  //! \param[in] regular Flag indicating how the sampling points are defined
+  //!
+  //! \details When \a regular is \e true, it is assumed that the parameter
+  //! value array \a gpar forms a regular tensor-product point grid of dimension
+  //! \a gpar[0].size() \a X \a gpar[1].size().
+  //! Otherwise, we assume that it contains the \a u and \a v parameters
+  //! directly for each sampling point.
+  virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
+			    const RealArray* gpar, bool regular = true) const;
+
+protected:
+
+  // Internal utility methods
+  // ========================
 
   //! \brief Calculates parameter values for visualization nodal points.
   //! \param[out] prm Parameter values in given direction for all points
diff --git a/src/ASM/ASMs2DLag.C b/src/ASM/ASMs2DLag.C
index 6ed962df..314d53fe 100644
--- a/src/ASM/ASMs2DLag.C
+++ b/src/ASM/ASMs2DLag.C
@@ -292,7 +292,20 @@ bool ASMs2DLag::integrate (Integrand& integrand, int lIndex,
   const int nelx = (nx-1)/(p1-1);
   const int nely = (ny-1)/(p2-1);
 
+  // Get parametric coordinates of the elements
   FiniteElement fe(p1*p2);
+  RealArray upar, vpar;
+  if (t1 == 1)
+  {
+    fe.u = edgeDir < 0 ? surf->startparam_u() : surf->endparam_u();
+    this->getGridParameters(vpar,1,1);
+  }
+  else if (t1 == 2)
+  {
+    this->getGridParameters(upar,0,1);
+    fe.v = edgeDir < 0 ? surf->startparam_v() : surf->endparam_v();
+  }
+
   Matrix dNdu, Xnod, Jac;
   Vec4   X;
   Vec3   normal;
@@ -337,6 +350,12 @@ bool ASMs2DLag::integrate (Integrand& integrand, int lIndex,
 	xi[t1-1] = edgeDir < 0 ? -1.0 : 1.0;
 	xi[t2-1] = xg[i];
 
+	// Parameter values of current integration point
+	if (upar.size() > 1)
+	  fe.u = 0.5*(upar[i1]*(1.0-xg[i]) + upar[i1+1]*(1.0+xg[i]));
+	if (vpar.size() > 1)
+	  fe.v = 0.5*(vpar[i2]*(1.0-xg[i]) + vpar[i2+1]*(1.0+xg[i]));
+
 	// Compute the basis functions and their derivatives, using
 	// tensor product of one-dimensional Lagrange polynomials
 	if (!Lagrange::computeBasis(fe.N,dNdu,p1,xi[0],p2,xi[1]))
@@ -402,6 +421,15 @@ bool ASMs2DLag::evalSolution (Matrix& sField, const Vector& locSol,
 }
 
 
+bool ASMs2DLag::evalSolution (Matrix&, const Vector&,
+			      const RealArray*, bool) const
+{
+  std::cerr <<" *** ASMs2DLag::evalSolution(Matrix&,const Vector&,"
+	    <<"const RealArray*,bool): Not implemented."<< std::endl;
+  return false;
+}
+
+
 bool ASMs2DLag::evalSolution (Matrix& sField, const Integrand& integrand,
 			      const int*) const
 {
@@ -457,3 +485,12 @@ bool ASMs2DLag::evalSolution (Matrix& sField, const Integrand& integrand,
 
   return true;
 }
+
+
+bool ASMs2DLag::evalSolution (Matrix&, const Integrand&,
+			      const RealArray*, bool) const
+{
+  std::cerr <<" *** ASMs2DLag::evalSolution(Matrix&,const Integrand&,"
+	    <<"const RealArray*,bool): Not implemented."<< std::endl;
+  return false;
+}
diff --git a/src/ASM/ASMs2DLag.h b/src/ASM/ASMs2DLag.h
index 0a94e449..ccd5fbc3 100644
--- a/src/ASM/ASMs2DLag.h
+++ b/src/ASM/ASMs2DLag.h
@@ -1,4 +1,4 @@
-// $Id: ASMs2DLag.h,v 1.6 2010-12-29 18:41:38 kmo Exp $
+// $Id$
 //==============================================================================
 //!
 //! \file ASMs2DLag.h
@@ -90,18 +90,30 @@ public:
   //! the Lagrange elements by default.
   //! \param[out] sField Solution field
   //! \param[in] locSol Solution vector in DOF-order
-  //! \param[in] npe Number of visualization nodes over each knot span
   virtual bool evalSolution(Matrix& sField, const Vector& locSol,
-			    const int* npe) const;
+			    const int*) const;
+
+  //! \brief Evaluates the primary solution field at the given points.
+  //! \param[out] sField Solution field
+  //! \param[in] locSol Solution vector local to current patch
+  //! \param[in] gpar Parameter values of the result sampling points
+  virtual bool evalSolution(Matrix& sField, const Vector& locSol,
+			    const RealArray* gpar, bool = true) const;
 
   //! \brief Evaluates the secondary solution field at all visualization points.
   //! \details The number of visualization points is the same as the order of
   //! the Lagrange elements by default.
   //! \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*) const;
+
+  //! \brief Evaluates the secondary solution field at the given points.
+  //! \param[out] sField Solution field
+  //! \param[in] integrand Object with problem-specific data and methods
+  //! \param[in] gpar Parameter values of the result sampling points
+  virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
+			    const RealArray* gpar, bool = true) const;
 
 protected:
 
diff --git a/src/ASM/ASMs2Dmx.C b/src/ASM/ASMs2Dmx.C
index 0642f75c..78eb321e 100644
--- a/src/ASM/ASMs2Dmx.C
+++ b/src/ASM/ASMs2Dmx.C
@@ -25,9 +25,6 @@
 #include "Vec3Oper.h"
 #include "Vec3.h"
 
-typedef Go::SplineSurface::Dvector DoubleVec; //!< 1D double array
-typedef DoubleVec::const_iterator DoubleIter; //!< Iterator over DoubleVec
-
 
 ASMs2Dmx::ASMs2Dmx (const char* fileName, unsigned char n_s,
 		    char n_f1, unsigned char n_f2)
@@ -289,7 +286,7 @@ bool ASMs2Dmx::getElementCoordinates (Matrix& X, int iel) const
   X.resize(nsd,nenod);
   const IntVec& mnpc = MNPC[iel-1];
 
-  DoubleIter cit = surf->coefs_begin();
+  RealArray::const_iterator cit = surf->coefs_begin();
   for (size_t n = 0; n < nenod; n++)
   {
     int iI = nodeInd[mnpc[lnod0+n]].I;
@@ -317,7 +314,7 @@ Vec3 ASMs2Dmx::getCoord (size_t inod) const
   }
 #endif
 
-  DoubleIter cit;
+  RealArray::const_iterator cit;
   const int I = nodeInd[inod-1].I;
   const int J = nodeInd[inod-1].J;
   if (inod <= nb1)
@@ -377,7 +374,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
   for (dir = 0; dir < 2; dir++)
   {
     int pm1 = (dir == 0 ? surf->order_u() : surf->order_v()) - 1;
-    DoubleIter uit = surf->basis(dir).begin() + pm1;
+    RealArray::const_iterator uit = surf->basis(dir).begin() + pm1;
     double ucurr, uprev = *(uit++);
     int nCol = (dir == 0 ? surf->numCoefs_u() : surf->numCoefs_v()) - pm1;
     gpar[dir].resize(nGauss,nCol);
@@ -519,7 +516,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
     else
     {
       int pm1 = (d == 0 ? surf->order_u() : surf->order_v()) - 1;
-      DoubleIter uit = surf->basis(d).begin() + pm1;
+      RealArray::const_iterator uit = surf->basis(d).begin() + pm1;
       double ucurr, uprev = *(uit++);
       int nCol = (d == 0 ? surf->numCoefs_u() : surf->numCoefs_v()) - pm1;
       gpar[d].resize(nGauss,nCol);
@@ -544,6 +541,9 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
 
   MxFiniteElement fe(basis1->order_u()*basis1->order_v(),
 		     basis2->order_u()*basis2->order_v());
+  fe.u = gpar[0](1,1);
+  fe.v = gpar[1](1,1);
+
   Matrix dN1du, dN2du, Xnod, Jac;
   Vec4   X;
   Vec3   normal;
@@ -569,7 +569,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
       if (skipMe) continue;
 
       // Get element edge length in the parameter space
-      double dS = this->getParametricLength(iel,abs(edgeDir));
+      double dS = this->getParametricLength(iel,t1);
       if (dS < 0.0) return false; // topology error (probably logic error)
 
       // Set up control point coordinates for current element
@@ -590,9 +590,13 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
 
       // --- Integration loop over all Gauss points along the edge -------------
 
-      int ip = (abs(edgeDir) == 1 ? i2-p2 : i1-p1)*nGauss;
+      int ip = (t1 == 1 ? i2-p2 : i1-p1)*nGauss;
       for (int i = 0; i < nGauss; i++, ip++)
       {
+	// Parameter values of current integration point
+	if (gpar[0].size() > 1) fe.u = gpar[0](i+1,i1-p1+1);
+	if (gpar[1].size() > 1) fe.v = gpar[1](i+1,i2-p2+1);
+
 	// Fetch basis function derivatives at current integration point
 	extractBasis(spline1[ip],fe.N1,dN1du);
 	extractBasis(spline2[ip],fe.N2,dN2du);
@@ -633,20 +637,35 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
 
 
 bool ASMs2Dmx::evalSolution (Matrix& sField, const Vector& locSol,
-			     const int* npe) const
+			     const RealArray* gpar, bool regular) const
 {
-  // Compute parameter values of the result sampling points
-  DoubleVec gpar[2];
-  for (int dir = 0; dir < 2; dir++)
-    if (!this->getGridParameters(gpar[dir],dir,npe[dir]-1))
-      return false;
-
   if (!basis1 || !basis2) return false;
 
   // Evaluate the basis functions at all points
   std::vector<Go::BasisPtsSf> spline1, spline2;
-  basis1->computeBasisGrid(gpar[0],gpar[1],spline1);
-  basis2->computeBasisGrid(gpar[0],gpar[1],spline2);
+  if (regular)
+  {
+    basis1->computeBasisGrid(gpar[0],gpar[1],spline1);
+    basis2->computeBasisGrid(gpar[0],gpar[1],spline2);
+  }
+  else if (gpar[0].size() == gpar[1].size())
+  {
+    std::vector<Go::BasisPtsSf> tmpSpline(1);
+    spline1.resize(gpar[0].size());
+    spline2.resize(gpar[0].size());
+    for (size_t i = 0; i < spline1.size(); i++)
+    {
+      basis1->computeBasisGrid(RealArray(1,gpar[0][i]),
+			       RealArray(1,gpar[1][i]),
+			       tmpSpline);
+      spline1[i] = tmpSpline.front();
+      basis2->computeBasisGrid(RealArray(1,gpar[0][i]),
+			       RealArray(1,gpar[1][i]),
+			       tmpSpline);
+      spline2[i] = tmpSpline.front();
+    }
+    // TODO: Request a GoTools method replacing the above (see ASMs2D)
+  }
 
   const int p1 = basis1->order_u();
   const int p2 = basis1->order_v();
@@ -699,22 +718,37 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const Vector& locSol,
 
 
 bool ASMs2Dmx::evalSolution (Matrix& sField, const Integrand& integrand,
-			     const int* npe) const
+			     const RealArray* gpar, bool regular) const
 {
   sField.resize(0,0);
 
-  // Compute parameter values of the result sampling points
-  DoubleVec gpar[2];
-  for (int dir = 0; dir < 2; dir++)
-    if (!this->getGridParameters(gpar[dir],dir,npe[dir]-1))
-      return false;
-
   if (!basis1 || !basis2) return false;
 
   // Evaluate the basis functions and their derivatives at all points
   std::vector<Go::BasisDerivsSf> spline1, spline2;
-  basis1->computeBasisGrid(gpar[0],gpar[1],spline1);
-  basis2->computeBasisGrid(gpar[0],gpar[1],spline2);
+  if (regular)
+  {
+    basis1->computeBasisGrid(gpar[0],gpar[1],spline1);
+    basis2->computeBasisGrid(gpar[0],gpar[1],spline2);
+  }
+  else if (gpar[0].size() == gpar[1].size())
+  {
+    std::vector<Go::BasisDerivsSf> tmpSpline(1);
+    spline1.resize(gpar[0].size());
+    spline2.resize(gpar[0].size());
+    for (size_t i = 0; i < spline1.size(); i++)
+    {
+      basis1->computeBasisGrid(RealArray(1,gpar[0][i]),
+			       RealArray(1,gpar[1][i]),
+			       tmpSpline);
+      spline1[i] = tmpSpline.front();
+      basis2->computeBasisGrid(RealArray(1,gpar[0][i]),
+			       RealArray(1,gpar[1][i]),
+			       tmpSpline);
+      spline2[i] = tmpSpline.front();
+    }
+    // TODO: Request a GoTools method replacing the above (see ASMs2D)
+  }
 
   const int p1 = basis1->order_u();
   const int p2 = basis1->order_v();
diff --git a/src/ASM/ASMs2Dmx.h b/src/ASM/ASMs2Dmx.h
index 06334ae9..5dc77b91 100644
--- a/src/ASM/ASMs2Dmx.h
+++ b/src/ASM/ASMs2Dmx.h
@@ -1,4 +1,4 @@
-// $Id: ASMs2Dmx.h,v 1.2 2010-12-30 15:02:02 kmo Exp $
+// $Id$
 //==============================================================================
 //!
 //! \file ASMs2Dmx.h
@@ -95,19 +95,33 @@ public:
   // Post-processing methods
   // =======================
 
-  //! \brief Evaluates the primary solution field at all visualization points.
+  //! \brief Evaluates the primary solution field at the given points.
   //! \param[out] sField Solution field
-  //! \param[in] locSol Solution vector in DOF-order
-  //! \param[in] npe Number of visualization nodes over each knot span
+  //! \param[in] locSol Solution vector local to current patch
+  //! \param[in] gpar Parameter values of the result sampling points
+  //! \param[in] regular Flag indicating how the sampling points are defined
+  //!
+  //! \details When \a regular is \e true, it is assumed that the parameter
+  //! value array \a gpar forms a regular tensor-product point grid of dimension
+  //! \a gpar[0].size() \a X \a gpar[1].size().
+  //! Otherwise, we assume that it contains the \a u and \a v parameters
+  //! directly for each sampling point.
   virtual bool evalSolution(Matrix& sField, const Vector& locSol,
-			    const int* npe) const;
+			    const RealArray* gpar, bool regular = true) const;
 
-  //! \brief Evaluates the secondary solution field at all visualization points.
+  //! \brief Evaluates the secondary solution field at the given 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
+  //! \param[in] gpar Parameter values of the result sampling points
+  //! \param[in] regular Flag indicating how the sampling points are defined
+  //!
+  //! \details When \a regular is \e true, it is assumed that the parameter
+  //! value array \a gpar forms a regular tensor-product point grid of dimension
+  //! \a gpar[0].size() \a X \a gpar[1].size().
+  //! Otherwise, we assume that it contains the \a u and \a v parameters
+  //! directly for each sampling point.
   virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
-			    const int* npe) const;
+			    const RealArray* gpar, bool regular = true) 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 734588f7..3690ce2b 100644
--- a/src/ASM/ASMs3D.C
+++ b/src/ASM/ASMs3D.C
@@ -29,9 +29,6 @@
 #include <ctype.h>
 #include <fstream>
 
-typedef Go::SplineVolume::Dvector DoubleVec;  //!< 1D double array
-typedef DoubleVec::const_iterator DoubleIter; //!< Iterator over DoubleVec
-
 
 ASMs3D::ASMs3D (const char* fileName, bool checkRHS, unsigned char n_f)
   : ASMstruct(3,3,n_f)
@@ -127,10 +124,10 @@ bool ASMs3D::checkRightHandSystem ()
   if (!svol) return false;
 
   // Evaluate the spline volume at its center
-  DoubleVec u(1,0.5*(svol->startparam(0) + svol->endparam(0)));
-  DoubleVec v(1,0.5*(svol->startparam(1) + svol->endparam(1)));
-  DoubleVec w(1,0.5*(svol->startparam(2) + svol->endparam(2)));
-  DoubleVec X(3), dXdu(3), dXdv(3), dXdw(3);
+  RealArray u(1,0.5*(svol->startparam(0) + svol->endparam(0)));
+  RealArray v(1,0.5*(svol->startparam(1) + svol->endparam(1)));
+  RealArray w(1,0.5*(svol->startparam(2) + svol->endparam(2)));
+  RealArray X(3), dXdu(3), dXdv(3), dXdw(3);
   svol->gridEvaluator(u,v,w,X,dXdu,dXdv,dXdw);
 
   // Check that |J| = (dXdu x dXdv) * dXdw > 0.0
@@ -149,8 +146,8 @@ bool ASMs3D::refine (int dir, const RealArray& xi)
   if (!svol || dir < 0 || dir > 2 || xi.empty()) return false;
   if (xi.front() < 0.0 || xi.back() > 1.0) return false;
 
-  DoubleVec  extraKnots;
-  DoubleIter uit = svol->basis(dir).begin();
+  RealArray extraKnots;
+  RealArray::const_iterator uit = svol->basis(dir).begin();
   double ucurr, uprev = *(uit++);
   while (uit != svol->basis(dir).end())
   {
@@ -174,8 +171,8 @@ bool ASMs3D::uniformRefine (int dir, int nInsert)
 {
   if (!svol || dir < 0 || dir > 2 || nInsert < 1) return false;
 
-  DoubleVec  extraKnots;
-  DoubleIter uit = svol->basis(dir).begin();
+  RealArray extraKnots;
+  RealArray::const_iterator uit = svol->basis(dir).begin();
   double ucurr, uprev = *(uit++);
   while (uit != svol->basis(dir).end())
   {
@@ -922,7 +919,7 @@ bool ASMs3D::getElementCoordinates (Matrix& X, int iel) const
   const IntVec& mnpc = MNPC[iel-1];
   X.resize(3,mnpc.size());
 
-  DoubleIter cit = svol->coefs_begin();
+  RealArray::const_iterator cit = svol->coefs_begin();
   for (size_t n = 0; n < mnpc.size(); n++)
   {
     int ip = this->coeffInd(mnpc[n])*svol->dimension();
@@ -946,7 +943,7 @@ void ASMs3D::getNodalCoordinates (Matrix& X) const
   const int n3 = svol->numCoefs(2);
   X.resize(3,n1*n2*n3);
 
-  DoubleIter cit = svol->coefs_begin();
+  RealArray::const_iterator cit = svol->coefs_begin();
   size_t inod = 1;
   for (int i3 = 0; i3 < n3; i3++)
     for (int i2 = 0; i2 < n2; i2++)
@@ -965,7 +962,7 @@ Vec3 ASMs3D::getCoord (size_t inod) const
   int ip = this->coeffInd(inod-1)*svol->dimension();
   if (ip < 0) return Vec3();
 
-  DoubleIter cit = svol->coefs_begin() + ip;
+  RealArray::const_iterator cit = svol->coefs_begin() + ip;
   return Vec3(*cit,*(cit+1),*(cit+2));
 }
 
@@ -1098,7 +1095,7 @@ bool ASMs3D::integrate (Integrand& integrand,
   for (dir = 0; dir < 3; dir++)
   {
     int pm1 = svol->order(dir) - 1;
-    DoubleIter uit = svol->basis(dir).begin() + pm1;
+    RealArray::const_iterator uit = svol->basis(dir).begin() + pm1;
     double ucurr, uprev = *(uit++);
     int nCol = svol->numCoefs(dir) - pm1;
     gpar[dir].resize(nGauss,nCol);
@@ -1302,7 +1299,7 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex,
     else
     {
       int pm1 = svol->order(d) - 1;
-      DoubleIter uit = svol->basis(d).begin() + pm1;
+      RealArray::const_iterator uit = svol->basis(d).begin() + pm1;
       double ucurr, uprev = *(uit++);
       int nCol = svol->numCoefs(d) - pm1;
       gpar[d].resize(nGauss,nCol);
@@ -1334,6 +1331,10 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex,
   const int nel2 = n2 - p2 + 1;
 
   FiniteElement fe(p1*p2*p3);
+  fe.u = gpar[0](1,1);
+  fe.v = gpar[1](1,1);
+  fe.w = gpar[2](1,1);
+
   Matrix dNdu, Xnod, Jac;
   Vec4   X;
   Vec3   normal;
@@ -1390,10 +1391,22 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex,
 
 	// --- Integration loop over all Gauss points in each direction --------
 
+	int k1, k2, k3;
 	int ip = (j2*nGauss*nf1 + j1)*nGauss;
 	for (int j = 0; j < nGauss; j++, ip += nGauss*(nf1-1))
 	  for (int i = 0; i < nGauss; i++, ip++)
 	  {
+	    // Parameter values of current integration point
+	    switch (abs(faceDir)) {
+	    case 1: k2 = i+1; k3 = j+1; k1 = 0; break;
+	    case 2: k1 = i+1; k3 = j+1; k2 = 0; break;
+	    case 3: k1 = i+1; k2 = j+1; k3 = 0; break;
+	    default: k1 = k2 = k3 = 0;
+	    }
+	    if (gpar[0].size() > 1) fe.u = gpar[0](k1,i1-p1+1);
+	    if (gpar[1].size() > 1) fe.v = gpar[1](k2,i2-p2+1);
+	    if (gpar[2].size() > 1) fe.w = gpar[1](k3,i3-p3+1);
+
 	    // Fetch basis function derivatives at current integration point
 	    extractBasis(spline[ip],fe.N,dNdu);
 
@@ -1455,7 +1468,7 @@ bool ASMs3D::integrateEdge (Integrand& integrand, int lEdge,
     else
     {
       int pm1 = svol->order(d) - 1;
-      DoubleIter uit = svol->basis(d).begin() + pm1;
+      RealArray::const_iterator uit = svol->basis(d).begin() + pm1;
       double ucurr, uprev = *(uit++);
       int nCol = svol->numCoefs(d) - pm1;
       gpar[d].resize(nGauss,nCol);
@@ -1484,6 +1497,10 @@ bool ASMs3D::integrateEdge (Integrand& integrand, int lEdge,
   const int p3 = svol->order(2);
 
   FiniteElement fe(p1*p2*p3);
+  fe.u = gpar[0](1,1);
+  fe.v = gpar[1](1,1);
+  fe.w = gpar[2](1,1);
+
   Matrix dNdu, Xnod, Jac;
   Vec4   X;
   Vec3   tang;
@@ -1551,6 +1568,11 @@ bool ASMs3D::integrateEdge (Integrand& integrand, int lEdge,
 	LocalIntegral* elmInt = 0;
 	for (int i = 0; i < nGauss; i++, ip++)
 	{
+	  // Parameter values of current integration point
+	  if (gpar[0].size() > 1) fe.u = gpar[0](i+1,i1-p1+1);
+	  if (gpar[1].size() > 1) fe.v = gpar[1](i+1,i2-p2+1);
+	  if (gpar[2].size() > 1) fe.w = gpar[1](i+1,i3-p3+1);
+
 	  // Fetch basis function derivatives at current integration point
 	  extractBasis(spline[ip],fe.N,dNdu);
 
@@ -1588,7 +1610,7 @@ bool ASMs3D::getGridParameters (RealArray& prm, int dir, int nSegPerSpan) const
     return false;
   }
 
-  DoubleIter uit = svol->basis(dir).begin();
+  RealArray::const_iterator uit = svol->basis(dir).begin();
   double ucurr = 0.0, uprev = *(uit++);
   while (uit != svol->basis(dir).end())
   {
@@ -1627,7 +1649,7 @@ bool ASMs3D::getGrevilleParameters (RealArray& prm, int dir) const
 bool ASMs3D::tesselate (ElementBlock& grid, const int* npe) const
 {
   // Compute parameter values of the nodal points
-  DoubleVec gpar[3];
+  RealArray gpar[3];
   for (int dir = 0; dir < 3; dir++)
     if (!this->getGridParameters(gpar[dir],dir,npe[dir]-1))
       return false;
@@ -1636,7 +1658,7 @@ bool ASMs3D::tesselate (ElementBlock& grid, const int* npe) const
   size_t nx = gpar[0].size();
   size_t ny = gpar[1].size();
   size_t nz = gpar[2].size();
-  DoubleVec XYZ(svol->dimension()*nx*ny*nz);
+  RealArray XYZ(svol->dimension()*nx*ny*nz);
   svol->gridEvaluator(gpar[0],gpar[1],gpar[2],XYZ);
 
   // Establish the block grid coordinates
@@ -1683,17 +1705,52 @@ bool ASMs3D::evalSolution (Matrix& sField, const Vector& locSol,
 			   const int* npe) const
 {
   // Compute parameter values of the result sampling points
-  DoubleVec gpar[3];
+  RealArray gpar[3];
   for (int dir = 0; dir < 3; dir++)
     if (!this->getGridParameters(gpar[dir],dir,npe[dir]-1))
       return false;
 
+  // Evaluate the primary solution at all sampling points
+  return this->evalSolution(sField,locSol,gpar);
+}
+
+
+bool ASMs3D::evalSolution (Matrix& sField, const Vector& locSol,
+			   const RealArray* gpar, bool regular) const
+{
   // Evaluate the basis functions at all points
   std::vector<Go::BasisPts> spline;
+  if (regular)
   {
     PROFILE2("Spline evaluation");
     svol->computeBasisGrid(gpar[0],gpar[1],gpar[2],spline);
   }
+  else if (gpar[0].size() == gpar[1].size() && gpar[0].size() == gpar[2].size())
+  {
+    PROFILE2("Spline evaluation");
+    spline.resize(gpar[0].size());
+    std::vector<Go::BasisPts> tmpSpline(1);
+    for (size_t i = 0; i < spline.size(); i++)
+    {
+      svol->computeBasisGrid(RealArray(1,gpar[0][i]),
+                             RealArray(1,gpar[1][i]),
+                             RealArray(1,gpar[2][i]),
+                             tmpSpline);
+      spline[i] = tmpSpline.front();
+    }
+    // TODO: Request a GoTools method replacing the above:
+    // void SplineVolume::computeBasisGrid(double param_u,
+    //                                     double param_v,
+    //                                     double param_w,
+    //                                     BasisPts& result) const
+    /*
+    spline.resize(gpar[0].size());
+    for (size_t i = 0; i < spline.size(); i++)
+      svol->computeBasis(gpar[0][i],gpar[1][i],gpar[2][i],spline[i]);
+    */
+  }
+  else
+    return false;
 
   const int p1 = svol->order(0);
   const int p2 = svol->order(1);
@@ -1743,9 +1800,11 @@ bool ASMs3D::evalSolution (Matrix& sField, const Integrand& integrand,
   }
   else
   {
+    // Project the secondary solution onto the spline basis
     Go::SplineVolume* v = this->projectSolution(integrand);
     if (v)
     {
+      // Extract control point values from the spline object
       sField.resize(v->dimension(),
 		    v->numCoefs(0)*v->numCoefs(1)*v->numCoefs(2));
       sField.fill(&(*v->coefs_begin()));
@@ -1768,7 +1827,7 @@ Go::GeomObject* ASMs3D::evalSolution (const Integrand& integrand) const
 Go::SplineVolume* ASMs3D::projectSolution (const Integrand& integrand) const
 {
   // Compute parameter values of the result sampling points
-  DoubleVec gpar[3];
+  RealArray gpar[3];
   for (int dir = 0; dir < 3; dir++)
     if (!this->getGrevilleParameters(gpar[dir],dir))
       return 0;
@@ -1785,7 +1844,7 @@ Go::SplineVolume* ASMs3D::projectSolution (const Integrand& integrand) const
   // the result array. Think that is always the case, but beware if trying
   // other projection schemes later.
 
-  DoubleVec weights;
+  RealArray weights;
   if (svol->rational())
     svol->getWeights(weights);
 
@@ -1802,16 +1861,43 @@ Go::SplineVolume* ASMs3D::projectSolution (const Integrand& integrand) const
 
 
 bool ASMs3D::evalSolution (Matrix& sField, const Integrand& integrand,
-			   const RealArray* gpar) const
+			   const RealArray* gpar, bool regular) const
 {
   sField.resize(0,0);
 
   // Evaluate the basis functions and their derivatives at all points
   std::vector<Go::BasisDerivs> spline;
+  if (regular)
   {
     PROFILE2("Spline evaluation");
     svol->computeBasisGrid(gpar[0],gpar[1],gpar[2],spline);
   }
+  else if (gpar[0].size() == gpar[1].size() && gpar[0].size() == gpar[2].size())
+  {
+    PROFILE2("Spline evaluation");
+    spline.resize(gpar[0].size());
+    std::vector<Go::BasisDerivs> tmpSpline(1);
+    for (size_t i = 0; i < spline.size(); i++)
+    {
+      svol->computeBasisGrid(RealArray(1,gpar[0][i]),
+                             RealArray(1,gpar[1][i]),
+                             RealArray(1,gpar[2][i]),
+                             tmpSpline);
+      spline[i] = tmpSpline.front();
+    }
+    // TODO: Request a GoTools method replacing the above:
+    // void SplineVolume::computeBasisGrid(double param_u,
+    //                                     double param_v,
+    //                                     double param_w,
+    //                                     BasisDerivs& result) const
+    /*
+    spline.resize(gpar[0].size());
+    for (size_t i = 0; i < spline.size(); i++)
+      svol->computeBasis(gpar[0][i],gpar[1][i],gpar[2][i],spline[i]);
+    */
+  }
+  else
+    return false;
 
   const int p1 = svol->order(0);
   const int p2 = svol->order(1);
diff --git a/src/ASM/ASMs3D.h b/src/ASM/ASMs3D.h
index cb376da2..ccb3eb74 100644
--- a/src/ASM/ASMs3D.h
+++ b/src/ASM/ASMs3D.h
@@ -1,4 +1,4 @@
-// $Id: ASMs3D.h,v 1.21 2011-01-28 13:34:28 kmo Exp $
+// $Id$
 //==============================================================================
 //!
 //! \file ASMs3D.h
@@ -273,6 +273,20 @@ public:
   virtual bool evalSolution(Matrix& sField, const Vector& locSol,
 			    const int* npe) const;
 
+  //! \brief Evaluates the primary solution field at the given points.
+  //! \param[out] sField Solution field
+  //! \param[in] locSol Solution vector local to current patch
+  //! \param[in] gpar Parameter values of the result sampling points
+  //! \param[in] regular Flag indicating how the sampling points are defined
+  //!
+  //! \details When \a regular is \e true, it is assumed that the parameter
+  //! value array \a gpar forms a regular tensor-product point grid of dimension
+  //! \a gpar[0].size() \a X \a gpar[1].size() \a X \a gpar[2].size().
+  //! Otherwise, we assume that it contains the \a u, \a v and \a w parameters
+  //! directly for each sampling point.
+  virtual bool evalSolution(Matrix& sField, const Vector& locSol,
+			    const RealArray* gpar, bool regular = true) const;
+
   //! \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
@@ -290,17 +304,24 @@ public:
   //! \brief Projects the secondary solution field onto the primary basis.
   virtual Go::GeomObject* evalSolution(const Integrand& integrand) const;
 
-protected:
-
-  // Internal utility methods
-  // ========================
-
   //! \brief Evaluates the secondary solution field at the given points.
   //! \param[out] sField Solution field
   //! \param[in] integrand Object with problem-specific data and methods
   //! \param[in] gpar Parameter values of the result sampling points
+  //! \param[in] regular Flag indicating how the sampling points are defined
+  //!
+  //! \details When \a regular is \e true, it is assumed that the parameter
+  //! value array \a gpar forms a regular tensor-product point grid of dimension
+  //! \a gpar[0].size() \a X \a gpar[1].size() \a X \a gpar[2].size().
+  //! Otherwise, we assume that it contains the \a u, \a v and \a w parameters
+  //! directly for each sampling point.
   bool evalSolution(Matrix& sField, const Integrand& integrand,
-		    const RealArray* gpar) const;
+		    const RealArray* gpar, bool regular = true) const;
+
+protected:
+
+  // Internal utility methods
+  // ========================
 
   //! \brief Calculates parameter values for visualization nodal points.
   //! \param[out] prm Parameter values in given direction for all points
diff --git a/src/ASM/ASMs3DLag.C b/src/ASM/ASMs3DLag.C
index eabb6487..49d92f3c 100644
--- a/src/ASM/ASMs3DLag.C
+++ b/src/ASM/ASMs3DLag.C
@@ -317,7 +317,24 @@ bool ASMs3DLag::integrate (Integrand& integrand, int lIndex,
   const int nely = (ny-1)/(p2-1);
   const int nelz = (nz-1)/(p3-1);
 
+  // Get parametric coordinates of the elements
+  RealArray upar, vpar, wpar;
+  if (t0 == 1)
+    upar.resize(1,faceDir < 0 ? svol->startparam(0) : svol->endparam(0));
+  else if (t0 == 2)
+    vpar.resize(1,faceDir < 0 ? svol->startparam(1) : svol->endparam(1));
+  else if (t0 == 3)
+    wpar.resize(1,faceDir < 0 ? svol->startparam(2) : svol->endparam(2));
+
+  if (upar.empty()) this->getGridParameters(upar,0,1);
+  if (vpar.empty()) this->getGridParameters(vpar,1,1);
+  if (wpar.empty()) this->getGridParameters(wpar,2,1);
+
   FiniteElement fe(p1*p2*p3);
+  fe.u = upar.front();
+  fe.v = vpar.front();
+  fe.w = wpar.front();
+
   Matrix dNdu, Xnod, Jac;
   Vec4   X;
   Vec3   normal;
@@ -359,6 +376,7 @@ bool ASMs3DLag::integrate (Integrand& integrand, int lIndex,
 
 	// --- Integration loop over all Gauss points in each direction --------
 
+	int k1, k2, k3;
 	for (int j = 0; j < nGauss; j++)
 	  for (int i = 0; i < nGauss; i++)
 	  {
@@ -367,6 +385,20 @@ bool ASMs3DLag::integrate (Integrand& integrand, int lIndex,
 	    xi[t1-1] = xg[i];
 	    xi[t2-1] = xg[j];
 
+	    // Parameter values of current integration point
+	    switch (abs(faceDir)) {
+	    case 1: k2 = i; k3 = j; k1 = -1; break;
+	    case 2: k1 = i; k3 = j; k2 = -1; break;
+	    case 3: k1 = i; k2 = j; k3 = -1; break;
+	    default: k1 = k2 = k3 = -1;
+	    }
+	    if (upar.size() > 1)
+	      fe.u = 0.5*(upar[i1]*(1.0-xg[k1]) + upar[i1+1]*(1.0+xg[k1]));
+	    if (vpar.size() > 1)
+	      fe.v = 0.5*(vpar[i2]*(1.0-xg[k2]) + vpar[i2+1]*(1.0+xg[k2]));
+	    if (wpar.size() > 1)
+	      fe.w = 0.5*(wpar[i3]*(1.0-xg[k3]) + wpar[i3+1]*(1.0+xg[k3]));
+
 	    // Compute the basis functions and their derivatives, using
 	    // tensor product of one-dimensional Lagrange polynomials
 	    if (!Lagrange::computeBasis(fe.N,dNdu,p1,xi[0],p2,xi[1],p3,xi[2]))
@@ -549,6 +581,15 @@ bool ASMs3DLag::evalSolution (Matrix& sField, const Vector& locSol,
 }
 
 
+bool ASMs3DLag::evalSolution (Matrix&, const Vector&,
+			      const RealArray*, bool) const
+{
+  std::cerr <<" *** ASMs3DLag::evalSolution(Matrix&,const Vector&,"
+	    <<"const RealArray*,bool): Not implemented."<< std::endl;
+  return false;
+}
+
+
 bool ASMs3DLag::evalSolution (Matrix& sField, const Integrand& integrand,
 			      const int*) const
 {
@@ -608,3 +649,12 @@ bool ASMs3DLag::evalSolution (Matrix& sField, const Integrand& integrand,
 
   return true;
 }
+
+
+bool ASMs3DLag::evalSolution (Matrix&, const Integrand&,
+			      const RealArray*, bool) const
+{
+  std::cerr <<" *** ASMs3DLag::evalSolution(Matrix&,const Integrand&,"
+	    <<"const RealArray*,bool): Not implemented."<< std::endl;
+  return false;
+}
diff --git a/src/ASM/ASMs3DLag.h b/src/ASM/ASMs3DLag.h
index 4f51aa7f..b4c81e4f 100644
--- a/src/ASM/ASMs3DLag.h
+++ b/src/ASM/ASMs3DLag.h
@@ -1,4 +1,4 @@
-// $Id: ASMs3DLag.h,v 1.7 2010-12-29 18:41:38 kmo Exp $
+// $Id$
 //==============================================================================
 //!
 //! \file ASMs3DLag.h
@@ -98,18 +98,30 @@ public:
   //! the Lagrange elements by default.
   //! \param[out] sField Solution field
   //! \param[in] locSol Solution vector in DOF-order
-  //! \param[in] npe Number of visualization nodes over each knot span
   virtual bool evalSolution(Matrix& sField, const Vector& locSol,
-			    const int* npe) const;
+			    const int*) const;
+
+  //! \brief Evaluates the primary solution field at the given points.
+  //! \param[out] sField Solution field
+  //! \param[in] locSol Solution vector local to current patch
+  //! \param[in] gpar Parameter values of the result sampling points
+  virtual bool evalSolution(Matrix& sField, const Vector& locSol,
+			    const RealArray* gpar, bool = true) const;
 
   //! \brief Evaluates the secondary solution field at all visualization points.
   //! \details The number of visualization points is the same as the order of
   //! the Lagrange elements by default.
   //! \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*) const;
+
+  //! \brief Evaluates the secondary solution field at the given points.
+  //! \param[out] sField Solution field
+  //! \param[in] integrand Object with problem-specific data and methods
+  //! \param[in] gpar Parameter values of the result sampling points
+  virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
+			    const RealArray* gpar, bool = true) const;
 
 protected:
 
diff --git a/src/ASM/ASMs3Dmx.C b/src/ASM/ASMs3Dmx.C
index bb60a8e3..9ce85705 100644
--- a/src/ASM/ASMs3Dmx.C
+++ b/src/ASM/ASMs3Dmx.C
@@ -25,9 +25,6 @@
 #include "Vec3Oper.h"
 #include "Vec3.h"
 
-typedef Go::SplineVolume::Dvector DoubleVec;  //!< 1D double array
-typedef DoubleVec::const_iterator DoubleIter; //!< Iterator over DoubleVec
-
 
 ASMs3Dmx::ASMs3Dmx (const char* fileName, bool checkRHS,
 		    char n_f1, unsigned char n_f2)
@@ -312,14 +309,14 @@ bool ASMs3Dmx::getElementCoordinates (Matrix& X, int iel) const
 
   const int n1 = svol->numCoefs(0);
   const int n2 = svol->numCoefs(1);
-  const int ndim = svol->dimension();
-  DoubleIter cit = svol->coefs_begin();
+  const int nd = svol->dimension();
+  RealArray::const_iterator cit = svol->coefs_begin();
   for (size_t n = 0; n < nenod; n++)
   {
     int iI = nodeInd[mnpc[lnod0+n]].I;
     int iJ = nodeInd[mnpc[lnod0+n]].J;
     int iK = nodeInd[mnpc[lnod0+n]].K;
-    int ip = ((iK*n2 + iJ)*n1 + iI)*ndim;
+    int ip = ((iK*n2 + iJ)*n1 + iI)*nd;
     for (size_t i = 0; i < 3; i++)
       X(i+1,n+1) = *(cit+(ip+i));
   }
@@ -342,7 +339,7 @@ Vec3 ASMs3Dmx::getCoord (size_t inod) const
   }
 #endif
 
-  DoubleIter cit;
+  RealArray::const_iterator cit;
   const int I = nodeInd[inod-1].I;
   const int J = nodeInd[inod-1].J;
   const int K = nodeInd[inod-1].K;
@@ -401,7 +398,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
   for (dir = 0; dir < 3; dir++)
   {
     int pm1 = svol->order(dir) - 1;
-    DoubleIter uit = svol->basis(dir).begin() + pm1;
+    RealArray::const_iterator uit = svol->basis(dir).begin() + pm1;
     double ucurr, uprev = *(uit++);
     int nCol = svol->numCoefs(dir) - pm1;
     gpar[dir].resize(nGauss,nCol);
@@ -551,7 +548,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
     else
     {
       int pm1 = svol->order(d) - 1;
-      DoubleIter uit = svol->basis(d).begin() + pm1;
+      RealArray::const_iterator uit = svol->basis(d).begin() + pm1;
       double ucurr, uprev = *(uit++);
       int nCol = svol->numCoefs(d) - pm1;
       gpar[d].resize(nGauss,nCol);
@@ -582,6 +579,9 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
 
   MxFiniteElement fe(basis1->order(0)*basis1->order(1)*basis1->order(2),
 		     basis2->order(0)*basis2->order(1)*basis2->order(2));
+  fe.u = gpar[0](1,1);
+  fe.v = gpar[1](1,1);
+  fe.w = gpar[2](1,1);
 
   Matrix dN1du, dN2du, Xnod, Jac;
   Vec4   X;
@@ -642,10 +642,22 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
 
 	// --- Integration loop over all Gauss points in each direction --------
 
+	int k1, k2, k3;
 	int ip = (j2*nGauss*nf1 + j1)*nGauss;
 	for (int j = 0; j < nGauss; j++, ip += nGauss*(nf1-1))
 	  for (int i = 0; i < nGauss; i++, ip++)
 	  {
+	    // Parameter values of current integration point
+	    switch (abs(faceDir)) {
+	    case 1: k2 = i+1; k3 = j+1; k1 = 0; break;
+	    case 2: k1 = i+1; k3 = j+1; k2 = 0; break;
+	    case 3: k1 = i+1; k2 = j+1; k3 = 0; break;
+	    default: k1 = k2 = k3 = 0;
+	    }
+	    if (gpar[0].size() > 1) fe.u = gpar[0](k1,i1-p1+1);
+	    if (gpar[1].size() > 1) fe.v = gpar[1](k2,i2-p2+1);
+	    if (gpar[2].size() > 1) fe.w = gpar[1](k3,i3-p3+1);
+
 	    // Fetch basis function derivatives at current integration point
 	    extractBasis(spline1[ip],fe.N1,dN1du);
 	    extractBasis(spline2[ip],fe.N2,dN2du);
@@ -686,20 +698,37 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
 
 
 bool ASMs3Dmx::evalSolution (Matrix& sField, const Vector& locSol,
-			     const int* npe) const
+			     const RealArray* gpar, bool regular) const
 {
-  // Compute parameter values of the result sampling points
-  DoubleVec gpar[3];
-  for (int dir = 0; dir < 3; dir++)
-    if (!this->getGridParameters(gpar[dir],dir,npe[dir]-1))
-      return false;
-
   if (!basis1 || !basis2) return false;
 
   // Evaluate the basis functions at all points
   std::vector<Go::BasisPts> spline1, spline2;
-  basis1->computeBasisGrid(gpar[0],gpar[1],gpar[2],spline1);
-  basis2->computeBasisGrid(gpar[0],gpar[1],gpar[2],spline2);
+  if (regular)
+  {
+    basis1->computeBasisGrid(gpar[0],gpar[1],gpar[2],spline1);
+    basis2->computeBasisGrid(gpar[0],gpar[1],gpar[2],spline2);
+  }
+  else if (gpar[0].size() == gpar[1].size() && gpar[0].size() == gpar[2].size())
+  {
+    std::vector<Go::BasisPts> tmpSpline(1);
+    spline1.resize(gpar[0].size());
+    spline2.resize(gpar[0].size());
+    for (size_t i = 0; i < spline1.size(); i++)
+    {
+      basis1->computeBasisGrid(RealArray(1,gpar[0][i]),
+			       RealArray(1,gpar[1][i]),
+			       RealArray(1,gpar[2][i]),
+			       tmpSpline);
+      spline1[i] = tmpSpline.front();
+      basis2->computeBasisGrid(RealArray(1,gpar[0][i]),
+			       RealArray(1,gpar[1][i]),
+			       RealArray(1,gpar[2][i]),
+			       tmpSpline);
+      spline2[i] = tmpSpline.front();
+    }
+    // TODO: Request a GoTools method replacing the above (see ASMs3D)
+  }
 
   const int p1 = basis1->order(0);
   const int p2 = basis1->order(1);
@@ -756,22 +785,39 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const Vector& locSol,
 
 
 bool ASMs3Dmx::evalSolution (Matrix& sField, const Integrand& integrand,
-			     const int* npe) const
+			     const RealArray* gpar, bool regular) const
 {
   sField.resize(0,0);
 
-  // Compute parameter values of the result sampling points
-  DoubleVec gpar[3];
-  for (int dir = 0; dir < 3; dir++)
-    if (!this->getGridParameters(gpar[dir],dir,npe[dir]-1))
-      return false;
-
   if (!basis1 || !basis2) return false;
 
   // Evaluate the basis functions and their derivatives at all points
   std::vector<Go::BasisDerivs> spline1, spline2;
-  basis1->computeBasisGrid(gpar[0],gpar[1],gpar[2],spline1);
-  basis2->computeBasisGrid(gpar[0],gpar[1],gpar[2],spline2);
+  if (regular)
+  {
+    basis1->computeBasisGrid(gpar[0],gpar[1],gpar[2],spline1);
+    basis2->computeBasisGrid(gpar[0],gpar[1],gpar[2],spline2);
+  }
+  else if (gpar[0].size() == gpar[1].size() && gpar[0].size() == gpar[2].size())
+  {
+    std::vector<Go::BasisDerivs> tmpSpline(1);
+    spline1.resize(gpar[0].size());
+    spline2.resize(gpar[0].size());
+    for (size_t i = 0; i < spline1.size(); i++)
+    {
+      basis1->computeBasisGrid(RealArray(1,gpar[0][i]),
+			       RealArray(1,gpar[1][i]),
+			       RealArray(1,gpar[2][i]),
+			       tmpSpline);
+      spline1[i] = tmpSpline.front();
+      basis2->computeBasisGrid(RealArray(1,gpar[0][i]),
+			       RealArray(1,gpar[1][i]),
+			       RealArray(1,gpar[2][i]),
+			       tmpSpline);
+      spline2[i] = tmpSpline.front();
+    }
+    // TODO: Request a GoTools method replacing the above (see ASMs3D)
+  }
 
   const int p1 = basis1->order(0);
   const int p2 = basis1->order(1);
diff --git a/src/ASM/ASMs3Dmx.h b/src/ASM/ASMs3Dmx.h
index 1d454fd0..7e3dd353 100644
--- a/src/ASM/ASMs3Dmx.h
+++ b/src/ASM/ASMs3Dmx.h
@@ -1,4 +1,4 @@
-// $Id: ASMs3Dmx.h,v 1.2 2010-12-30 15:02:02 kmo Exp $
+// $Id$
 //==============================================================================
 //!
 //! \file ASMs3Dmx.h
@@ -95,19 +95,33 @@ public:
   // Post-processing methods
   // =======================
 
-  //! \brief Evaluates the primary solution field at all visualization points.
+  //! \brief Evaluates the primary solution field at the given points.
   //! \param[out] sField Solution field
-  //! \param[in] locSol Solution vector in DOF-order
-  //! \param[in] npe Number of visualization nodes over each knot span
+  //! \param[in] locSol Solution vector local to current patch
+  //! \param[in] gpar Parameter values of the result sampling points
+  //! \param[in] regular Flag indicating how the sampling points are defined
+  //!
+  //! \details When \a regular is \e true, it is assumed that the parameter
+  //! value array \a gpar forms a regular tensor-product point grid of dimension
+  //! \a gpar[0].size() \a X \a gpar[1].size() \a X \a gpar[2].size().
+  //! Otherwise, we assume that it contains the \a u, \a v and \a w parameters
+  //! directly for each sampling point.
   virtual bool evalSolution(Matrix& sField, const Vector& locSol,
-			    const int* npe) const;
+			    const RealArray* gpar, bool regular = true) const;
 
-  //! \brief Evaluates the secondary solution field at all visualization points.
+  //! \brief Evaluates the secondary solution field at the given 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
+  //! \param[in] gpar Parameter values of the result sampling points
+  //! \param[in] regular Flag indicating how the sampling points are defined
+  //!
+  //! \details When \a regular is \e true, it is assumed that the parameter
+  //! value array \a gpar forms a regular tensor-product point grid of dimension
+  //! \a gpar[0].size() \a X \a gpar[1].size() \a X \a gpar[2].size().
+  //! Otherwise, we assume that it contains the \a u, \a v and \a w parameters
+  //! directly for each sampling point.
   virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
-			    const int* npe) const;
+			    const RealArray* gpar, bool regular = true) const;
 
   //! \brief Extracts nodal results for this patch from the global vector.
   //! \param[in] globVec Global solution vector in DOF-order