diff --git a/src/ASM/ASMs2D.h b/src/ASM/ASMs2D.h
index b9b9d9a6..6dcc5d1d 100644
--- a/src/ASM/ASMs2D.h
+++ b/src/ASM/ASMs2D.h
@@ -76,8 +76,10 @@ public:
   //! \brief Empty destructor.
   virtual ~ASMs2D() {}
 
-  //! \brief Returns the spline surface representing this patch.
+  //! \brief Returns the spline surface representing the geometry of this patch.
   Go::SplineSurface* getSurface() const { return surf; }
+  //! \brief Returns the spline surface representing the basis of this patch.
+  virtual Go::SplineSurface* getBasis() const { return surf; }
 
 
   // Methods for model generation
diff --git a/src/ASM/ASMs2DLag.h b/src/ASM/ASMs2DLag.h
index 35f23fcc..96dca6f4 100644
--- a/src/ASM/ASMs2DLag.h
+++ b/src/ASM/ASMs2DLag.h
@@ -138,6 +138,8 @@ protected:
   //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
   //! in one element
   virtual bool getElementCoordinates(Matrix& X, int iel) const;
+
+public:
   //! \brief Returns a matrix with all nodal coordinates within the patch.
   //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
   //! in the patch
diff --git a/src/ASM/ASMs2Dmx.h b/src/ASM/ASMs2Dmx.h
index 8836bf48..143c4789 100644
--- a/src/ASM/ASMs2Dmx.h
+++ b/src/ASM/ASMs2Dmx.h
@@ -39,6 +39,9 @@ public:
   //! \brief Empty destructor.
   virtual ~ASMs2Dmx() {}
 
+  //! \brief Returns the spline surface representing the basis of this patch.
+  virtual Go::SplineSurface* getBasis() const { return basis1; }
+
 
   //! \brief Returns the spline surface representing this patch.
   Go::SplineSurface* getSurface() const { return surf; }
diff --git a/src/ASM/ASMs3D.h b/src/ASM/ASMs3D.h
index 55643bd4..60fb4765 100644
--- a/src/ASM/ASMs3D.h
+++ b/src/ASM/ASMs3D.h
@@ -98,6 +98,8 @@ public:
 
   //! \brief Returns the spline volume representing this patch.
   Go::SplineVolume* getVolume() const { return svol; }
+  //! \brief Returns the spline volume representing the basis of this patch.
+  virtual Go::SplineVolume* getBasis() const { return svol; }
 
 
   // Methods for model generation
diff --git a/src/ASM/ASMs3DLag.h b/src/ASM/ASMs3DLag.h
index d8a5c651..d78d2168 100644
--- a/src/ASM/ASMs3DLag.h
+++ b/src/ASM/ASMs3DLag.h
@@ -148,6 +148,8 @@ protected:
   //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
   //! in one element
   virtual bool getElementCoordinates(Matrix& X, int iel) const;
+
+public:
   //! \brief Returns a matrix with all nodal coordinates within the patch.
   //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
   //! in the patch
diff --git a/src/ASM/ASMs3Dmx.h b/src/ASM/ASMs3Dmx.h
index 9280ebb1..67736b69 100644
--- a/src/ASM/ASMs3Dmx.h
+++ b/src/ASM/ASMs3Dmx.h
@@ -42,6 +42,9 @@ public:
   //! \brief Empty destructor.
   virtual ~ASMs3Dmx() {}
 
+  //! \brief Returns the spline volume representing the basis of this patch.
+  virtual Go::SplineVolume* getBasis() const { return basis1; }
+
 
   // Methods for model generation
   // ============================
diff --git a/src/ASM/Field.C b/src/ASM/Field.C
new file mode 100644
index 00000000..c7d045b9
--- /dev/null
+++ b/src/ASM/Field.C
@@ -0,0 +1,37 @@
+// $Id$
+//==============================================================================
+//!
+//! \file Field.C
+//!
+//! \date Mar 28 2011
+//!
+//! \author Runar Holdahl / SINTEF
+//!
+//! \brief Base class for scalar fields.
+//!
+//==============================================================================
+
+#include "SplineField2D.h"
+#include "SplineField3D.h"
+#include "LagrangeField2D.h"
+#include "LagrangeField3D.h"
+#include "ASMs2DLag.h"
+#include "ASMs3DLag.h"
+
+
+Field* Field::create (const ASMbase* pch, const RealArray& v, const char* name)
+{
+  const ASMs2DLag* pl2 = dynamic_cast<const ASMs2DLag*>(pch);
+  if (pl2) return new LagrangeField2D(pl2,v,name);
+
+  const ASMs2D* ps2 = dynamic_cast<const ASMs2D*>(pch);
+  if (ps2) return new SplineField2D(ps2,v,name);
+
+  const ASMs3DLag* pl3 = dynamic_cast<const ASMs3DLag*>(pch);
+  if (pl3) return new LagrangeField3D(pl3,v,name);
+
+  const ASMs3D* ps3 = dynamic_cast<const ASMs3D*>(pch);
+  if (ps3) return new SplineField3D(ps3,v,name);
+
+  return NULL;
+}
diff --git a/src/ASM/Field.h b/src/ASM/Field.h
index edb86d27..79073d89 100644
--- a/src/ASM/Field.h
+++ b/src/ASM/Field.h
@@ -17,6 +17,7 @@
 #include "MatVec.h"
 #include <string>
 
+class ASMbase;
 class FiniteElement;
 class Vec3;
 
@@ -36,7 +37,8 @@ protected:
   //! \brief The constructor sets the number of space dimensions and field name.
   //! \param[in] n Number of space dimensions (1, 2 or 3)
   //! \param[in] name Name of field
-  Field(unsigned char n, char* name = 0) : nsd(n) { if (name) fname = name; }
+  Field(unsigned char n, const char* name = 0) : nsd(n), nelm(0), nno(0)
+  { if (name) fname = name; }
 
 public:
   //! \brief Empty destructor.
@@ -46,27 +48,27 @@ public:
   unsigned char getNoSpaceDim() const { return nsd; }
 
   //! \brief Returns number of elements.
-  int getNoElm() const { return nelm; }
+  size_t getNoElm() const { return nelm; }
 
-  //! \brief Returns number of control points.
-  int getNoNodes() const { return nno; }
+  //! \brief Returns number of nodal/control points.
+  size_t getNoNodes() const { return nno; }
 
-  //! \brief Returns name of field.
+  //! \brief Returns the name of field.
   const char* getFieldName() const { return fname.c_str(); }
 
-  //! \brief Sets the name of the field.
-  void setFieldName(const char* name) { fname = name; }
+  //! \brief Creates a dynamically allocated field object.
+  //! \param[in] pch The spline patch on which the field is to be defined on
+  //! \param[in] v Array of nodal/control point field values
+  //! \param[in] name Name of field
+  static Field* create(const ASMbase* pch, const RealArray& v,
+		       const char* name = NULL);
 
-  //! \brief Initializes the field values.
-  void fill(const Vector& vec) { values = vec; }
-
-
-  // Methods to compute field values
-  //================================
+  // Methods to evaluate the field
+  //==============================
 
   //! \brief Computes the value in a given node/control point.
   //! \param[in] node Node number
-  virtual double valueNode(int node) const = 0;
+  virtual double valueNode(size_t node) const = 0;
 
   //! \brief Computes the value at a given local coordinate.
   //! \param[in] fe Finite element definition
@@ -88,8 +90,8 @@ public:
 
 protected:
   unsigned char nsd; //!< Number of space dimensions
-  int nelm;          //!< Number of elements/knot-spans
-  int nno;           //!< Number of nodes/control points
+  size_t nelm;       //!< Number of elements/knot-spans
+  size_t nno;        //!< Number of nodes/control points
   std::string fname; //!< Name of the field
   Vector values;     //!< Field values
 };
diff --git a/src/ASM/Fields.C b/src/ASM/Fields.C
new file mode 100644
index 00000000..461fe82b
--- /dev/null
+++ b/src/ASM/Fields.C
@@ -0,0 +1,38 @@
+// $Id$
+//==============================================================================
+//!
+//! \file Fields.C
+//!
+//! \date Mar 28 2011
+//!
+//! \author Runar Holdahl / SINTEF
+//!
+//! \brief Base class for vector fields.
+//!
+//==============================================================================
+
+#include "SplineFields2D.h"
+#include "SplineFields3D.h"
+#include "LagrangeFields2D.h"
+#include "LagrangeFields3D.h"
+#include "ASMs2DLag.h"
+#include "ASMs3DLag.h"
+
+
+Fields* Fields::create (const ASMbase* pch, const RealArray& v,
+			const char* name)
+{
+  const ASMs2DLag* pl2 = dynamic_cast<const ASMs2DLag*>(pch);
+  if (pl2) return new LagrangeFields2D(pl2,v,name);
+
+  const ASMs2D* ps2 = dynamic_cast<const ASMs2D*>(pch);
+  if (ps2) return new SplineFields2D(ps2,v,name);
+
+  const ASMs3DLag* pl3 = dynamic_cast<const ASMs3DLag*>(pch);
+  if (pl3) return new LagrangeFields3D(pl3,v,name);
+
+  const ASMs3D* ps3 = dynamic_cast<const ASMs3D*>(pch);
+  if (ps3) return new SplineFields3D(ps3,v,name);
+
+  return NULL;
+}
diff --git a/src/ASM/Fields.h b/src/ASM/Fields.h
index 8a9189c3..d10ea140 100644
--- a/src/ASM/Fields.h
+++ b/src/ASM/Fields.h
@@ -17,6 +17,7 @@
 #include "MatVec.h"
 #include <string>
 
+class ASMbase;
 class FiniteElement;
 class Vec3;
 
@@ -58,20 +59,20 @@ public:
   //! \brief Returns name of field.
   const char* getFieldName() const { return fname.c_str(); }
 
-  //! \brief Sets the name of the field.
-  void setFieldName(const char* name) { fname = name; }
+  //! \brief Creates a dynamically allocated field object.
+  //! \param[in] pch The spline patch on which the field is to be defined on
+  //! \param[in] v Array of nodal/control point field values
+  //! \param[in] name Name of field
+  static Fields* create(const ASMbase* pch, const RealArray& v,
+                        const char* name = NULL);
 
-  //! \brief Initializes the field values.
-  void fill(const Vector& vec) { values = vec; nf = values.size()/nno; }
-
-
-  // Methods to compute field values
-  //================================
+  // Methods to evaluate the field
+  //==============================
 
   //! \brief Computes the value in a given node/control point.
   //! \param[in] node Node number
   //! \param[out] vals Node values
-  virtual bool valueNode(int node, Vector& vals) const = 0;
+  virtual bool valueNode(size_t node, Vector& vals) const = 0;
 
   //! \brief Computes the value at a given local coordinate.
   //! \param[in] fe Finite element definition
diff --git a/src/ASM/LagrangeField2D.C b/src/ASM/LagrangeField2D.C
index 5b35dae9..58e05c79 100644
--- a/src/ASM/LagrangeField2D.C
+++ b/src/ASM/LagrangeField2D.C
@@ -12,24 +12,32 @@
 //==============================================================================
 
 #include "LagrangeField2D.h"
+#include "ASMs2DLag.h"
 #include "FiniteElement.h"
 #include "Lagrange.h"
 #include "CoordinateMapping.h"
 #include "Vec3.h"
 
 
-LagrangeField2D::LagrangeField2D (const Matrix& X, int nx, int ny,
-				  int px, int py, char* name)
-  : Field(2,name), coord(X), n1(nx), n2(ny), p1(px), p2(py)
+LagrangeField2D::LagrangeField2D (const ASMs2DLag* patch, const RealArray& v,
+				  const char* name) : Field(2,name)
 {
+  patch->getNodalCoordinates(coord);
+  patch->getSize(n1,n2);
+  patch->getOrder(p1,p2);
   nno = n1*n2;
   nelm = (n1-1)*(n2-1)/(p1*p2);
+
+  // Ensure the values array has compatible length, pad with zeros if necessary
+  values.resize(nno);
+  RealArray::const_iterator end = v.size() > nno ? v.begin()+nno : v.end();
+  std::copy(v.begin(),end,values.begin());
 }
 
 
-double LagrangeField2D::valueNode (int node) const
+double LagrangeField2D::valueNode (size_t node) const
 {
-  return values(node);
+  return node > 0 && node <= nno ? values(node) : 0.0;
 }
 
 
diff --git a/src/ASM/LagrangeField2D.h b/src/ASM/LagrangeField2D.h
index 56a40f70..3ca2d844 100644
--- a/src/ASM/LagrangeField2D.h
+++ b/src/ASM/LagrangeField2D.h
@@ -16,6 +16,8 @@
 
 #include "Field.h"
 
+class ASMs2DLag;
+
 
 /*!
   \brief Class for Lagrange-based finite element scalar fields in 2D.
@@ -28,23 +30,20 @@ class LagrangeField2D : public Field
 {
 public:
   //! \brief The constructor sets the number of space dimensions and fields.
-  //! \param[in] X  Matrix of nodal coordinates
-  //! \param[in] nx Number of nodes in first parameter direction
-  //! \param[in] ny Number of nodes in second parameter direction
-  //! \param[in] px Element order in first parameter direction
-  //! \param[in] py Element order in second parameter direction
+  //! \param[in] patch The spline patch on which the field is to be defined
+  //! \param[in] v Array of control point field values
   //! \param[in] name Name of field
-  LagrangeField2D(const Matrix& X,
-		  int nx, int ny, int px, int py, char* name = NULL);
+  LagrangeField2D(const ASMs2DLag* patch, const RealArray& v,
+		  const char* name = NULL);
   //! \brief Empty destructor.
   virtual ~LagrangeField2D() {}
 
-  // Methods to compute field values
-  //================================
+  // Methods to evaluate the field
+  //==============================
 
   //! \brief Computes the value in a given node/control point.
   //! \param[in] node Node number
-  double valueNode(int node) const;
+  double valueNode(size_t node) const;
 
   //! \brief Computes the value at a given local coordinate.
   //! \param[in] fe Finite element definition
diff --git a/src/ASM/LagrangeField3D.C b/src/ASM/LagrangeField3D.C
index 43b8e468..524c4581 100644
--- a/src/ASM/LagrangeField3D.C
+++ b/src/ASM/LagrangeField3D.C
@@ -12,24 +12,32 @@
 //==============================================================================
 
 #include "LagrangeField3D.h"
+#include "ASMs3DLag.h"
 #include "FiniteElement.h"
 #include "Lagrange.h"
 #include "CoordinateMapping.h"
 #include "Vec3.h"
 
 
-LagrangeField3D::LagrangeField3D (const Matrix& X, int nx, int ny, int nz,
-				  int px, int py, int pz, char* name)
-  : Field(3,name), coord(X), n1(nx), n2(ny), n3(nz), p1(px), p2(py), p3(pz)
+LagrangeField3D::LagrangeField3D (const ASMs3DLag* patch, const RealArray& v,
+				  const char* name) : Field(3,name)
 {
+  patch->getNodalCoordinates(coord);
+  patch->getSize(n1,n2,n3);
+  patch->getOrder(p1,p2,p3);
   nno = n1*n2*n3;
   nelm = (n1-1)*(n2-1)*(n3-1)/(p1*p2*p3);
+
+  // Ensure the values array has compatible length, pad with zeros if necessary
+  values.resize(nno);
+  RealArray::const_iterator end = v.size() > nno ? v.begin()+nno : v.end();
+  std::copy(v.begin(),end,values.begin());
 }
 
 
-double LagrangeField3D::valueNode (int node) const
+double LagrangeField3D::valueNode (size_t node) const
 {
-  return values(node);
+  return node > 0 && node <= nno ? values(node) : 0.0;
 }
 
 
diff --git a/src/ASM/LagrangeField3D.h b/src/ASM/LagrangeField3D.h
index a8f276a8..130bc208 100644
--- a/src/ASM/LagrangeField3D.h
+++ b/src/ASM/LagrangeField3D.h
@@ -16,6 +16,8 @@
 
 #include "Field.h"
 
+class ASMs3DLag;
+
 
 /*!
   \brief Class for Lagrange-based finite element scalar fields in 3D.
@@ -28,31 +30,26 @@ class LagrangeField3D : public Field
 {
 public:
   //! \brief The constructor sets the number of space dimensions and fields.
-  //! \param[in] X  Matrix of nodal coordinates
-  //! \param[in] nx Number of nodes in first parameter direction
-  //! \param[in] ny Number of nodes in second parameter direction
-  //! \param[in] nz Number of nodes in third parameter direction
-  //! \param[in] px Element order in first parameter direction
-  //! \param[in] py Element order in second parameter direction
-  //! \param[in] pz Element order in third parameter direction
+  //! \param[in] patch The spline patch on which the field is to be defined
+  //! \param[in] v Array of control point field values
   //! \param[in] name Name of field
-  LagrangeField3D(const Matrix& X, int nx, int ny, int nz,
-		  int px, int py, int pz, char* name = NULL);
+  LagrangeField3D(const ASMs3DLag* patch, const RealArray& v,
+		  const char* name = NULL);
   //! \brief Empty destructor.
   virtual ~LagrangeField3D() {}
 
-  // Methods to compute field values
-  //================================
+  // Methods to evaluate the field
+  //==============================
 
   //! \brief Computes the value in a given node/control point.
   //! \param[in] node Node number
-  double valueNode(int node) const;
+  double valueNode(size_t node) const;
 
   //! \brief Computes the value at a given local coordinate.
   //! \param[in] fe Finite element definition
   double valueFE(const FiniteElement& fe) const;
 
-  //! \brief Computed the value at a given global coordinate.
+  //! \brief Computes the value at a given global coordinate.
   //! \param[in] x Global/physical coordinate for point
   double valueCoor(const Vec3& x) const;
 
diff --git a/src/ASM/LagrangeFields2D.C b/src/ASM/LagrangeFields2D.C
index 0fcb7293..89f829dd 100644
--- a/src/ASM/LagrangeFields2D.C
+++ b/src/ASM/LagrangeFields2D.C
@@ -12,24 +12,33 @@
 //==============================================================================
 
 #include "LagrangeFields2D.h"
+#include "ASMs2DLag.h"
 #include "FiniteElement.h"
 #include "Lagrange.h"
 #include "CoordinateMapping.h"
 #include "Vec3.h"
 
 
-LagrangeFields2D::LagrangeFields2D (const Matrix& X, int nx, int ny,
-				    int px, int py, char* name)
-  : Fields(2,name), coord(X), n1(nx), n2(ny), p1(px), p2(py)
+LagrangeFields2D::LagrangeFields2D (const ASMs2DLag* patch, const RealArray& v,
+				    const char* name) : Fields(2,name)
 {
+  patch->getNodalCoordinates(coord);
+  patch->getSize(n1,n2);
+  patch->getOrder(p1,p2);
   nno = n1*n2;
   nelm = (n1-1)*(n2-1)/(p1*p2);
+  nf = v.size()/nno;
+
+  // Ensure the values array has compatible length, pad with zeros if necessary
+  values.resize(nf*nno);
+  RealArray::const_iterator end = v.size() > nf*nno ? v.begin()+nf*nno:v.end();
+  std::copy(v.begin(),end,values.begin());
 }
 
 
-bool LagrangeFields2D::valueNode (int node, Vector& vals) const
+bool LagrangeFields2D::valueNode (size_t node, Vector& vals) const
 {
-  if (node < 1 || (size_t)node > nno) return false;
+  if (node < 1 || node > nno) return false;
 
   vals.resize(nf);
   vals.fill(values.ptr()+nf*(node-1));
diff --git a/src/ASM/LagrangeFields2D.h b/src/ASM/LagrangeFields2D.h
index c85a916f..78218d11 100644
--- a/src/ASM/LagrangeFields2D.h
+++ b/src/ASM/LagrangeFields2D.h
@@ -16,6 +16,8 @@
 
 #include "Fields.h"
 
+class ASMs2DLag;
+
 
 /*!
   \brief Class for Lagrange-based finite element vector fields in 2D.
@@ -27,25 +29,22 @@
 class LagrangeFields2D : public Fields
 {
 public:
-  //! \brief The constructor sets the field name.
-  //! \param[in] X  Matrix of nodal coordinates
-  //! \param[in] nx Number of nodes in first parameter direction
-  //! \param[in] ny Number of nodes in second parameter direction
-  //! \param[in] px Element order in first parameter direction
-  //! \param[in] py Element order in second parameter direction
+  //! \brief The constructor sets the number of space dimensions and fields.
+  //! \param[in] patch The spline patch on which the field is to be defined
+  //! \param[in] v Array of control point field values
   //! \param[in] name Name of field
-  LagrangeFields2D(const Matrix& X,
-		   int nx, int ny, int px, int py, char* name = NULL);
+  LagrangeFields2D(const ASMs2DLag* patch, const RealArray& v,
+                   const char* name = NULL);
   //! \brief Empty destructor.
   virtual ~LagrangeFields2D() {}
 
-  // Methods to compute field values
-  //================================
+  // Methods to evaluate the field
+  //==============================
 
   //! \brief Computes the value in a given node/control point.
   //! \param[in] node Node number
   //! \param[out] vals Node values
-  bool valueNode(int node, Vector& vals) const;
+  bool valueNode(size_t node, Vector& vals) const;
 
   //! \brief Computes the value at a given local coordinate.
   //! \param[in] fe Finite element definition
diff --git a/src/ASM/LagrangeFields3D.C b/src/ASM/LagrangeFields3D.C
index 52533d21..6ad04cb3 100644
--- a/src/ASM/LagrangeFields3D.C
+++ b/src/ASM/LagrangeFields3D.C
@@ -12,24 +12,33 @@
 //==============================================================================
 
 #include "LagrangeFields3D.h"
+#include "ASMs3DLag.h"
 #include "FiniteElement.h"
 #include "Lagrange.h"
 #include "CoordinateMapping.h"
 #include "Vec3.h"
 
 
-LagrangeFields3D::LagrangeFields3D (const Matrix& X, int nx, int ny, int nz,
-				    int px, int py, int pz, char* name)
-  : Fields(3,name), coord(X), n1(nx), n2(ny), n3(nz), p1(px), p2(py), p3(pz)
+LagrangeFields3D::LagrangeFields3D (const ASMs3DLag* patch, const RealArray& v,
+				    const char* name) : Fields(3,name)
 {
+  patch->getNodalCoordinates(coord);
+  patch->getSize(n1,n2,n3);
+  patch->getOrder(p1,p2,p3);
   nno = n1*n2*n3;
   nelm = (n1-1)*(n2-1)*(n3-1)/(p1*p2*p3);
+  nf = v.size()/nno;
+
+  // Ensure the values array has compatible length, pad with zeros if necessary
+  values.resize(nf*nno);
+  RealArray::const_iterator end = v.size() > nf*nno ? v.begin()+nf*nno:v.end();
+  std::copy(v.begin(),end,values.begin());
 }
 
 
-bool LagrangeFields3D::valueNode (int node, Vector& vals) const
+bool LagrangeFields3D::valueNode (size_t node, Vector& vals) const
 {
-  if (node < 1 || (size_t)node > nno) return false;
+  if (node < 1 || node > nno) return false;
 
   vals.resize(nf);
   vals.fill(values.ptr()+nf*(node-1));
diff --git a/src/ASM/LagrangeFields3D.h b/src/ASM/LagrangeFields3D.h
index 94606430..a30dd70f 100644
--- a/src/ASM/LagrangeFields3D.h
+++ b/src/ASM/LagrangeFields3D.h
@@ -16,6 +16,8 @@
 
 #include "Fields.h"
 
+class ASMs3DLag;
+
 
 /*!
   \brief Class for Lagrange-based finite element vector fields in 3D.
@@ -27,17 +29,12 @@
 class LagrangeFields3D : public Fields
 {
 public:
-  //! \brief The constructor sets the field name.
-  //! \param[in] X  Matrix of nodal coordinates
-  //! \param[in] nx Number of nodes in first parameter direction
-  //! \param[in] ny Number of nodes in second parameter direction
-  //! \param[in] nz Number of nodes in third parameter direction
-  //! \param[in] px Element order in first parameter direction
-  //! \param[in] py Element order in second parameter direction
-  //! \param[in] pz Element order in third parameter direction
+  //! \brief The constructor sets the number of space dimensions and fields.
+  //! \param[in] patch The spline patch on which the field is to be defined
+  //! \param[in] v Array of control point field values
   //! \param[in] name Name of field
-  LagrangeFields3D(const Matrix& X, int nx, int ny, int nz,
-		   int px, int py, int pz, char* name = NULL);
+  LagrangeFields3D(const ASMs3DLag* patch, const RealArray& v,
+                   const char* name = NULL);
   //! \brief Empty destructor.
   virtual ~LagrangeFields3D() {}
 
@@ -47,7 +44,7 @@ public:
   //! \brief Computes the value in a given node/control point.
   //! \param[in] node Node number
   //! \param[out] vals Node values
-  bool valueNode(int node, Vector& vals) const;
+  bool valueNode(size_t node, Vector& vals) const;
 
   //! \brief Computes the value at a given local coordinate.
   //! \param[in] fe Finite element definition
diff --git a/src/ASM/SplineField2D.C b/src/ASM/SplineField2D.C
index 4824e021..7fee8c94 100644
--- a/src/ASM/SplineField2D.C
+++ b/src/ASM/SplineField2D.C
@@ -12,15 +12,18 @@
 //==============================================================================
 
 #include "SplineField2D.h"
+#include "ASMs2D.h"
 #include "FiniteElement.h"
+#include "CoordinateMapping.h"
+#include "Utilities.h"
 #include "Vec3.h"
 
 #include "GoTools/geometry/SplineSurface.h"
-#include "GoTools/geometry/SplineUtils.h"
 
 
-SplineField2D::SplineField2D(Go::SplineSurface *bf, char* name)
-  : Field(2,name), basis(bf), surf(bf) 
+SplineField2D::SplineField2D (const ASMs2D* patch,
+                              const RealArray& v, const char* name)
+  : Field(2,name), basis(patch->getBasis()), surf(patch->getSurface())
 {
   const int n1 = basis->numCoefs_u();
   const int n2 = basis->numCoefs_v();
@@ -29,44 +32,41 @@ SplineField2D::SplineField2D(Go::SplineSurface *bf, char* name)
   const int p1 = basis->order_u();
   const int p2 = basis->order_v();
   nelm = (n1-p1+1)*(n2-p2+1);
+
+  // Ensure the values array has compatible length, pad with zeros if necessary
+  values.resize(nno);
+  RealArray::const_iterator end = v.size() > nno ? v.begin()+nno : v.end();
+  std::copy(v.begin(),end,values.begin());
 }
 
 
-SplineField2D::SplineField2D(Go::SplineSurface *bf, 
-			     Go::SplineSurface *geometry, 
-			     char* name)
-  : Field(2,name), basis(bf), surf(geometry) 
+double SplineField2D::valueNode (size_t node) const
 {
-  const int n1 = basis->numCoefs_u();
-  const int n2 = basis->numCoefs_v();
-  nno = n1*n2;
-
-  const int p1 = basis->order_u();
-  const int p2 = basis->order_v();
-  nelm = (n1-p1+1)*(n2-p2+1);
+  return node > 0 && node <= nno ? values(node) : 0.0;
 }
 
 
-SplineField2D::~SplineField2D()
+double SplineField2D::valueFE (const FiniteElement& fe) const
 {
-  // Set basis and geometry pointers to NULL, should not be deallocated here
-  basis = surf = NULL;
+  if (!basis) return false;
+
+  // Evaluate the basis functions at the given point
+  Go::BasisPtsSf spline;
+  basis->computeBasis(fe.u,fe.v,spline);
+
+  // Evaluate the solution field at the given point
+  std::vector<int> ip;
+  ASMs2D::scatterInd(basis->numCoefs_u(),basis->numCoefs_v(),
+		     basis->order_u(),basis->order_v(),
+		     spline.left_idx,ip);
+
+  Vector Vnod;
+  utl::gather(ip,1,values,Vnod);
+  return Vnod.dot(spline.basisValues);
 }
 
-
-double SplineField2D::valueNode(int node) const
-{
-  // Not implemented yet
-  return 0.0;
-}
-
-
-double SplineField2D::valueFE(const FiniteElement& fe) const
-{
-//   Go::Point pt;
-//   basis->pointValue(pt,values,fe.u,fe.v);
-//   return pt[0];
-
+/* Old procedure, way too complex...
+   The above code is more in line with how we do it in ASMs2D.
   const int uorder = basis->order_u();
   const int vorder = basis->order_v();
   const int unum = basis->numCoefs_u();
@@ -140,19 +140,82 @@ double SplineField2D::valueFE(const FiniteElement& fe) const
   
   return val;
 }
+*/
 
-
-double SplineField2D::valueCoor(const Vec3& x) const
+double SplineField2D::valueCoor (const Vec3& x) const
 {
   // Not implemented yet
   return 0.0;
 }
 
 
-bool SplineField2D::gradFE(const FiniteElement& fe, Vector& grad) const
+bool SplineField2D::gradFE (const FiniteElement& fe, Vector& grad) const
 {
   if (!basis) return false;
+  if (!surf)  return false;
 
+  // Evaluate the basis functions at the given point
+/*
+  Go::BasisDerivsSf spline;
+  surf->computeBasis(fe.u,fe.v,spline);
+  TODO: The above is not available yet, the below is temporary workaround.
+*/
+  std::vector<Go::BasisDerivsSf> tmpSpline(1);
+  surf->computeBasisGrid(RealArray(1,fe.u),RealArray(1,fe.v),tmpSpline);
+  Go::BasisDerivsSf& spline = tmpSpline.front();
+
+  const int uorder = surf->order_u();
+  const int vorder = surf->order_v();
+  const size_t nen = uorder*vorder;
+
+  Matrix dNdu(nen,2), dNdX;
+  for (size_t n = 1; n <= nen; n++)
+  {
+    dNdu(n,1) = spline.basisDerivs_u[n-1];
+    dNdu(n,2) = spline.basisDerivs_v[n-1];
+  }
+
+  std::vector<int> ip;
+  ASMs2D::scatterInd(surf->numCoefs_u(),surf->numCoefs_v(),
+		     uorder,vorder,spline.left_idx,ip);
+
+  // Evaluate the Jacobian inverse
+  Matrix Xnod, Jac;
+  Vector Xctrl(&(*surf->coefs_begin()),surf->coefs_end()-surf->coefs_begin());
+  utl::gather(ip,surf->dimension(),Xctrl,Xnod);
+  utl::Jacobian(Jac,dNdX,Xnod,dNdu);
+
+  // Evaluate the gradient of the solution field at the given point
+  if (basis != surf)
+  {
+    // Mixed formulation, the solution uses a different basis than the geometry
+    /*
+    basis->computeBasis(fe.u,fe.v,spline);
+    TODO: The above is not available yet, the below is temporary workaround.
+    */
+    basis->computeBasisGrid(RealArray(1,fe.u),RealArray(1,fe.v),tmpSpline);
+    Go::BasisDerivsSf& spline = tmpSpline.front();
+
+    const size_t nbf = basis->order_u()*basis->order_v();
+    dNdu.resize(nbf,2);
+    for (size_t n = 1; n <= nbf; n++)
+    {
+      dNdu(n,1) = spline.basisDerivs_u[n-1];
+      dNdu(n,2) = spline.basisDerivs_v[n-1];
+    }
+    dNdX.multiply(dNdu,Jac); // dNdX = dNdu * Jac
+
+    ASMs2D::scatterInd(basis->numCoefs_u(),basis->numCoefs_v(),
+		       basis->order_u(),basis->order_v(),
+		       spline.left_idx,ip);
+  }
+
+  Vector Vnod;
+  utl::gather(ip,1,values,Vnod);
+  return dNdX.multiply(Vnod,grad,true);
+}
+
+/*
 //   // Derivatives of solution in parametric space
 //   std::vector<Go::Point> pts(3);
 //   basis->pointValue(pts,values,fe.u,fe.v,1);
@@ -303,9 +366,10 @@ bool SplineField2D::gradFE(const FiniteElement& fe, Vector& grad) const
   
   return true;
 }
+*/
 
 
-bool SplineField2D::gradCoor(const Vec3& x, Vector& grad) const
+bool SplineField2D::gradCoor (const Vec3& x, Vector& grad) const
 {
   // Not implemented yet
   return false;
diff --git a/src/ASM/SplineField2D.h b/src/ASM/SplineField2D.h
index 03ebe0f1..f38bf48f 100644
--- a/src/ASM/SplineField2D.h
+++ b/src/ASM/SplineField2D.h
@@ -16,6 +16,8 @@
 
 #include "Field.h"
 
+class ASMs2D;
+
 namespace Go {
   class SplineSurface;
 }
@@ -24,40 +26,34 @@ namespace Go {
 /*!
   \brief Class for spline-based finite element scalar fields in 2D.
 
-  \details This class implements the functions required to evaluate a 2D
+  \details This class implements the methods required to evaluate a 2D
   spline scalar field at a given point in parametrical or physical coordinates.
 */
 
-
 class SplineField2D : public Field
 {
 public:
   //! \brief The constructor sets the number of space dimensions and fields.
-  //! \param[in] bf Spline basis description
+  //! \param[in] patch The spline patch on which the field is to be defined
+  //! \param[in] v Array of control point field values
   //! \param[in] name Name of spline field
-  SplineField2D(Go::SplineSurface *bf, char* name = NULL);
-  //! \brief The constructor sets the number of space dimensions and fields.
-  //! \param[in] bf Spline basis description
-  //! \param[in] geometry Spline geometry description
-  //! \param[in] name Name of spline field
-  SplineField2D(Go::SplineSurface *bf, 
-		Go::SplineSurface *geometry, 
-		char* name = NULL);
+  SplineField2D(const ASMs2D* patch, const RealArray& v,
+                const char* name = NULL);
   //! \brief Empty destructor.
-  virtual ~SplineField2D();
+  virtual ~SplineField2D() {}
 
-  // Methods to compute field values
-  //================================
+  // Methods to evaluate the field
+  //==============================
 
   //! \brief Computes the value in a given node/control point.
   //! \param[in] node Node number
-  double valueNode(int node) const;
+  double valueNode(size_t node) const;
 
   //! \brief Computes the value at a given local coordinate.
   //! \param[in] fe Finite element definition
   double valueFE(const FiniteElement& fe) const;
 
-  //! \brief Computed the value at a given global coordinate.
+  //! \brief Computes the value at a given global coordinate.
   //! \param[in] x Global/physical coordinate for point
   double valueCoor(const Vec3& x) const;
 
@@ -72,8 +68,8 @@ public:
   bool gradCoor(const Vec3& x, Vector& grad) const;
 
 protected:
-  Go::SplineSurface* basis; //!< Spline basis description
-  Go::SplineSurface* surf;  //!< Spline geometry description
+  const Go::SplineSurface* basis; //!< Spline basis description
+  const Go::SplineSurface* surf;  //!< Spline geometry description
 };
 
 #endif
diff --git a/src/ASM/SplineField3D.C b/src/ASM/SplineField3D.C
index 8ea9d665..d977be58 100644
--- a/src/ASM/SplineField3D.C
+++ b/src/ASM/SplineField3D.C
@@ -12,19 +12,18 @@
 //==============================================================================
 
 #include "SplineField3D.h"
+#include "ASMs3D.h"
 #include "FiniteElement.h"
+#include "CoordinateMapping.h"
+#include "Utilities.h"
 #include "Vec3.h"
 
 #include "GoTools/trivariate/SplineVolume.h"
-#include "GoTools/geometry/SplineUtils.h"
-
-namespace Go {
-  void volume_ratder(double const eder[],int idim,int ider,double gder[]);
-}
 
 
-SplineField3D::SplineField3D(Go::SplineVolume *bf, char* name)
-  : Field(3,name), basis(bf), vol(bf) 
+SplineField3D::SplineField3D (const ASMs3D* patch,
+                              const RealArray& v, const char* name)
+  : Field(3,name), basis(patch->getBasis()), vol(patch->getVolume())
 {
   const int n1 = basis->numCoefs(0);
   const int n2 = basis->numCoefs(1);
@@ -35,46 +34,41 @@ SplineField3D::SplineField3D(Go::SplineVolume *bf, char* name)
   const int p2 = basis->order(1);
   const int p3 = basis->order(2);
   nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1);
+
+  // Ensure the values array has compatible length, pad with zeros if necessary
+  values.resize(nno);
+  RealArray::const_iterator end = v.size() > nno ? v.begin()+nno : v.end();
+  std::copy(v.begin(),end,values.begin());
 }
 
 
-SplineField3D::SplineField3D(Go::SplineVolume *bf, 
-			     Go::SplineVolume *geometry, 
-			     char* name)
-  : Field(3,name), basis(bf), vol(geometry) 
+double SplineField3D::valueNode (size_t node) const
 {
-  const int n1 = basis->numCoefs(0);
-  const int n2 = basis->numCoefs(1);
-  const int n3 = basis->numCoefs(2);
-  nno = n1*n2*n3;
-
-  const int p1 = basis->order(0);
-  const int p2 = basis->order(1);
-  const int p3 = basis->order(2);
-  nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1);
+  return node > 0 && node <= nno ? values(node) : 0.0;
 }
 
 
-SplineField3D::~SplineField3D()
+double SplineField3D::valueFE (const FiniteElement& fe) const
 {
-  // Set basis and geometry pointers to NULL, should not be deallocated here
-  basis = vol = NULL;
+  if (!basis) return false;
+
+  // Evaluate the basis functions at the given point
+  Go::BasisPts spline;
+  basis->computeBasis(fe.u,fe.v,fe.w,spline);
+
+  // Evaluate the solution field at the given point
+  std::vector<int> ip;
+  ASMs3D::scatterInd(basis->numCoefs(0),basis->numCoefs(1),basis->numCoefs(2),
+		     basis->order(0),basis->order(1),basis->order(2),
+		     spline.left_idx,ip);
+
+  Vector Vnod;
+  utl::gather(ip,1,values,Vnod);
+  return Vnod.dot(spline.basisValues);
 }
 
-
-double SplineField3D::valueNode(int node) const
-{
-  // Not implemented yet
-  return 0.0;
-}
-
-
-double SplineField3D::valueFE(const FiniteElement& fe) const
-{
-//   Go::Point pt;
-//   basis->pointValue(pt,values,fe.u,fe.v,fe.w);
-//   return pt[0];
-
+/* Old procedure, way too complex...
+   The above code is more in line with how we do it in ASMs3D.
   double val = 0.0;
 
   const int uorder = basis->order(0);
@@ -168,19 +162,85 @@ double SplineField3D::valueFE(const FiniteElement& fe) const
   val = tempResult/w;
   return val;
 }
+*/
 
-
-double SplineField3D::valueCoor(const Vec3& x) const
+double SplineField3D::valueCoor (const Vec3& x) const
 {
   // Not implemented yet
   return 0.0;
 }
 
 
-bool SplineField3D::gradFE(const FiniteElement& fe, Vector& grad) const
+bool SplineField3D::gradFE (const FiniteElement& fe, Vector& grad) const
 {
   if (!basis) return false;
+  if (!vol)   return false;
 
+  // Evaluate the basis functions at the given point
+/*
+  Go::BasisDerivs spline;
+  vol->computeBasis(fe.u,fe.v,fe.w,spline);
+  TODO: The above is not available yet, the below is temporary workaround.
+*/
+  std::vector<Go::BasisDerivs> tmpSpline(1);
+  vol->computeBasisGrid(RealArray(1,fe.u),RealArray(1,fe.v),RealArray(1,fe.w),tmpSpline);
+  Go::BasisDerivs& spline = tmpSpline.front();
+
+  const int uorder = vol->order(0);
+  const int vorder = vol->order(1);
+  const int worder = vol->order(2);
+  const size_t nen = uorder*vorder*worder;
+
+  Matrix dNdu(nen,3), dNdX;
+  for (size_t n = 1; n <= nen; n++)
+  {
+    dNdu(n,1) = spline.basisDerivs_u[n-1];
+    dNdu(n,2) = spline.basisDerivs_v[n-1];
+    dNdu(n,3) = spline.basisDerivs_w[n-1];
+  }
+
+  std::vector<int> ip;
+  ASMs3D::scatterInd(vol->numCoefs(0),vol->numCoefs(1),vol->numCoefs(2),
+		     uorder,vorder,worder,spline.left_idx,ip);
+
+  // Evaluate the Jacobian inverse
+  Matrix Xnod, Jac;
+  Vector Xctrl(&(*vol->coefs_begin()),vol->coefs_end()-vol->coefs_begin());
+  utl::gather(ip,vol->dimension(),Xctrl,Xnod);
+  utl::Jacobian(Jac,dNdX,Xnod,dNdu);
+
+  // Evaluate the gradient of the solution field at the given point
+  if (basis != vol)
+  {
+    // Mixed formulation, the solution uses a different basis than the geometry
+    /*
+    basis->computeBasis(fe.u,fe.v,fe.w,spline);
+    TODO: The above is not available yet, the below is temporary workaround.
+    */
+    basis->computeBasisGrid(RealArray(1,fe.u),RealArray(1,fe.v),RealArray(1,fe.w),tmpSpline);
+    Go::BasisDerivs& spline = tmpSpline.front();
+
+    const size_t nbf = basis->order(0)*basis->order(1)*basis->order(2);
+    dNdu.resize(nbf,3);
+    for (size_t n = 1; n <= nbf; n++)
+    {
+      dNdu(n,1) = spline.basisDerivs_u[n-1];
+      dNdu(n,2) = spline.basisDerivs_v[n-1];
+      dNdu(n,3) = spline.basisDerivs_w[n-1];
+    }
+    dNdX.multiply(dNdu,Jac);
+
+    ASMs3D::scatterInd(basis->numCoefs(0),basis->numCoefs(1),basis->numCoefs(2),
+		       basis->order(0),basis->order(1),basis->order(2),
+		       spline.left_idx,ip);
+  }
+
+  Vector Vnod;
+  utl::gather(ip,1,values,Vnod);
+  return dNdX.multiply(Vnod,grad,true);
+}
+
+/*
 //   // Derivatives of solution in parametric space
 //   std::vector<Go::Point> pts;
 //   basis->pointValue(pts,values,fe.u,fe.v,fe.w,1);
@@ -376,9 +436,9 @@ bool SplineField3D::gradFE(const FiniteElement& fe, Vector& grad) const
 
   return true;
 }
+*/
 
-
-bool SplineField3D::gradCoor(const Vec3& x, Vector& grad) const
+bool SplineField3D::gradCoor (const Vec3& x, Vector& grad) const
 {
   // Not implemented yet
   return false;
diff --git a/src/ASM/SplineField3D.h b/src/ASM/SplineField3D.h
index c4bd570a..adc02aa4 100644
--- a/src/ASM/SplineField3D.h
+++ b/src/ASM/SplineField3D.h
@@ -16,6 +16,8 @@
 
 #include "Field.h"
 
+class ASMs3D;
+
 namespace Go {
   class SplineVolume;
 }
@@ -28,35 +30,30 @@ namespace Go {
   spline scalar field at a given point in parametrical or physical coordinates.
 */
 
-
 class SplineField3D : public Field
 {
 public:
   //! \brief The constructor sets the number of space dimensions and fields.
-  //! \param[in] bf Spline basis description
+  //! \param[in] patch The spline patch on which the field is to be defined
+  //! \param[in] v Array of control point field values
   //! \param[in] name Name of spline field
-  SplineField3D(Go::SplineVolume *bf, char* name = NULL);
-  //! \param[in] bf Spline basis description
-  //! \param[in] geometry Spline geometry description
-  //! \param[in] name Name of spline field
-  SplineField3D(Go::SplineVolume *bf, 
-		Go::SplineVolume *geometry, 
-		char* name = NULL);
+  SplineField3D(const ASMs3D* patch, const RealArray& v,
+                const char* name = NULL);
   //! \brief Empty destructor.
-  virtual ~SplineField3D();
+  virtual ~SplineField3D() {}
 
-  // Methods to compute field values
-  //================================
+  // Methods to evaluate the field
+  //==============================
 
   //! \brief Computes the value in a given node/control point.
   //! \param[in] node Node number
-  double valueNode(int node) const;
+  double valueNode(size_t node) const;
 
   //! \brief Computes the value at a given local coordinate.
   //! \param[in] fe Finite element definition
   double valueFE(const FiniteElement& fe) const;
 
-  //! \brief Computed the value at a given global coordinate.
+  //! \brief Computes the value at a given global coordinate.
   //! \param[in] x Global/physical coordinate for point
   double valueCoor(const Vec3& x) const;
 
@@ -71,8 +68,8 @@ public:
   bool gradCoor(const Vec3& x, Vector& grad) const;
 
 protected:
-  Go::SplineVolume* basis; //!< Spline basis description
-  Go::SplineVolume* vol;   //!< Spline volume description
+  const Go::SplineVolume* basis; //!< Spline basis description
+  const Go::SplineVolume* vol;   //!< Spline geometry description
 };
 
 #endif
diff --git a/src/ASM/SplineFields2D.C b/src/ASM/SplineFields2D.C
index 23d4e8c0..f0dc535c 100644
--- a/src/ASM/SplineFields2D.C
+++ b/src/ASM/SplineFields2D.C
@@ -12,18 +12,18 @@
 //==============================================================================
 
 #include "SplineFields2D.h"
+#include "ASMs2D.h"
 #include "FiniteElement.h"
+#include "CoordinateMapping.h"
+#include "Utilities.h"
 #include "Vec3.h"
 
 #include "GoTools/geometry/SplineSurface.h"
-#include "GoTools/geometry/SplineUtils.h"
-#ifndef __BORLANDC__
-#include "boost/lambda/lambda.hpp"
-#endif
 
 
-SplineFields2D::SplineFields2D(Go::SplineSurface *bf, char* name)
-  : Fields(2,name),  basis(bf), surf(bf) 
+SplineFields2D::SplineFields2D (const ASMs2D* patch,
+                                const RealArray& v, const char* name)
+  : Fields(2,name), basis(patch->getBasis()), surf(patch->getSurface())
 {
   const int n1 = basis->numCoefs_u();
   const int n2 = basis->numCoefs_v();
@@ -33,53 +33,47 @@ SplineFields2D::SplineFields2D(Go::SplineSurface *bf, char* name)
   const int p2 = basis->order_v();
   nelm = (n1-p1+1)*(n2-p2+1);
 
-  // Number of fields set in fill
-  nf = 0;
+  // Ensure the values array has compatible length, pad with zeros if necessary
+  nf = v.size()/nno;
+  values.resize(nf*nno);
+  RealArray::const_iterator end = v.size() > nf*nno ? v.begin()+nf*nno:v.end();
+  std::copy(v.begin(),end,values.begin());
 }
 
 
-SplineFields2D::SplineFields2D(Go::SplineSurface *bf, 
-			       Go::SplineSurface *geometry, 
-			       char* name)
-  : Fields(2,name),  basis(bf), surf(geometry) 
+bool SplineFields2D::valueNode (size_t node, Vector& vals) const
 {
-  const int n1 = basis->numCoefs_u();
-  const int n2 = basis->numCoefs_v();
-  nno = n1*n2;
+  if (node < 1 || node > nno) return false;
 
-  const int p1 = basis->order_u();
-  const int p2 = basis->order_v();
-  nelm = (n1-p1+1)*(n2-p2+1);
-
-  // Number of fields set in fill
-  nf = 0;
+  vals.resize(nf);
+  vals.fill(values.ptr()+(node-1)*nf);
+  return true;
 }
 
 
-SplineFields2D::~SplineFields2D()
-{
-  // Set geometry to NULL; should not be deallocated here
-  basis = surf = NULL;
-} 
-
-
-bool SplineFields2D::valueNode(int node, Vector& vals) const
-{
-  // Not implemented yet
-  return false;
-}
-
-
-bool SplineFields2D::valueFE(const FiniteElement& fe, Vector& vals) const
+bool SplineFields2D::valueFE (const FiniteElement& fe, Vector& vals) const
 {
   if (!basis) return false;
 
-//   Go::Point pt;
-//   basis->pointValue(pt,values,fe.u,fe.v);
-//   vals.resize(pt.size());
-//   for (int i = 0;i < pt.size();i++)
-//     vals[i] = pt[i];
+  // Evaluate the basis functions at the given point
+  Go::BasisPtsSf spline;
+  basis->computeBasis(fe.u,fe.v,spline);
 
+  // Evaluate the solution field at the given point
+  std::vector<int> ip;
+  ASMs2D::scatterInd(basis->numCoefs_u(),basis->numCoefs_v(),
+		     basis->order_u(),basis->order_v(),
+		     spline.left_idx,ip);
+
+  Matrix Vnod;
+  utl::gather(ip,nf,values,Vnod);
+  Vnod.multiply(spline.basisValues,vals); // vals = Vnod * basisValues
+
+  return true;
+}
+
+/* Old procedure, way too complex...
+   The above code is more in line with how we do it in ASMs2D.
   const int uorder = basis->order_u();
   const int vorder = basis->order_v();
   const int unum = basis->numCoefs_u();
@@ -175,19 +169,82 @@ bool SplineFields2D::valueFE(const FiniteElement& fe, Vector& vals) const
 
   return true;
 }
+*/
 
 
-
-bool SplineFields2D::valueCoor(const Vec3& x, Vector& vals) const
+bool SplineFields2D::valueCoor (const Vec3& x, Vector& vals) const
 {
   // Not implemented yet
   return false;
 }
 
 
-bool SplineFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const
+bool SplineFields2D::gradFE (const FiniteElement& fe, Matrix& grad) const
 {
   if (!basis) return false;
+  if (!surf)  return false;
+
+  // Evaluate the basis functions at the given point
+/*
+  Go::BasisDerivsSf spline;
+  surf->computeBasis(fe.u,fe.v,spline);
+  TODO: The above is not available yet, the below is temporary workaround.
+*/
+  std::vector<Go::BasisDerivsSf> tmpSpline(1);
+  surf->computeBasisGrid(RealArray(1,fe.u),RealArray(1,fe.v),tmpSpline);
+  Go::BasisDerivsSf& spline = tmpSpline.front();
+
+  const int uorder = surf->order_u();
+  const int vorder = surf->order_v();
+  const size_t nen = uorder*vorder;
+
+  Matrix dNdu(nen,2), dNdX;
+  for (size_t n = 1; n <= nen; n++)
+  {
+    dNdu(n,1) = spline.basisDerivs_u[n-1];
+    dNdu(n,2) = spline.basisDerivs_v[n-1];
+  }
+
+  std::vector<int> ip;
+  ASMs2D::scatterInd(surf->numCoefs_u(),surf->numCoefs_v(),
+		     uorder,vorder,spline.left_idx,ip);
+
+  // Evaluate the Jacobian inverse
+  Matrix Xnod, Jac;
+  Vector Xctrl(&(*surf->coefs_begin()),surf->coefs_end()-surf->coefs_begin());
+  utl::gather(ip,surf->dimension(),Xctrl,Xnod);
+  utl::Jacobian(Jac,dNdX,Xnod,dNdu);
+
+  // Evaluate the gradient of the solution field at the given point
+  if (basis != surf)
+  {
+    // Mixed formulation, the solution uses a different basis than the geometry
+    /*
+    basis->computeBasis(fe.u,fe.v,spline);
+    TODO: The above is not available yet, the below is temporary workaround.
+    */
+    basis->computeBasisGrid(RealArray(1,fe.u),RealArray(1,fe.v),tmpSpline);
+    Go::BasisDerivsSf& spline = tmpSpline.front();
+
+    const size_t nbf = basis->order_u()*basis->order_v();
+    dNdu.resize(nbf,2);
+    for (size_t n = 1; n <= nbf; n++)
+    {
+      dNdu(n,1) = spline.basisDerivs_u[n-1];
+      dNdu(n,2) = spline.basisDerivs_v[n-1];
+    }
+    dNdX.multiply(dNdu,Jac); // dNdX = dNdu * Jac
+
+    ASMs2D::scatterInd(basis->numCoefs_u(),basis->numCoefs_v(),
+		       basis->order_u(),basis->order_v(),
+		       spline.left_idx,ip);
+  }
+
+  utl::gather(ip,nf,values,Xnod);
+  grad.multiply(Xnod,dNdX); // grad = Xnod * dNdX
+
+  return true;
+}
 
 //   // Derivatives of solution in parametric space
 //   std::vector<Go::Point> pts(3);
@@ -215,7 +272,7 @@ bool SplineFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const
 //   // Gradient of solution in given parametric point
 //   // df/dX = df/dXi * Jac^-1 = df/dXi * [dX/dXi]^-1
 //   grad.multiply(gradXi,Jac);
-
+/*
   // Gradient of field 
   Matrix gradXi(nf,2);
 
@@ -351,9 +408,10 @@ bool SplineFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const
 
   return true;
 }
+*/
 
 
-bool SplineFields2D::gradCoor(const Vec3& x, Matrix& grad) const
+bool SplineFields2D::gradCoor (const Vec3& x, Matrix& grad) const
 {
   // Not implemented yet
   return false;
diff --git a/src/ASM/SplineFields2D.h b/src/ASM/SplineFields2D.h
index 83b1d4cf..76c4fb54 100644
--- a/src/ASM/SplineFields2D.h
+++ b/src/ASM/SplineFields2D.h
@@ -16,6 +16,8 @@
 
 #include "Fields.h"
 
+class ASMs2D;
+
 namespace Go {
   class SplineSurface;
 }
@@ -24,44 +26,38 @@ namespace Go {
 /*!
   \brief Class for spline-based finite element vector fields in 2D.
 
-  \details This class implements the functions required to evaluate a 2D
+  \details This class implements the methods required to evaluate a 2D
   spline vector field at a given point in parametrical or physical coordinates.
 */
 
-
 class SplineFields2D : public Fields
 {
 public:
-  //! \brief The constructor sets the field name.
-  //! \param[in] bf Spline basis description
+  //! \brief The constructor sets the number of space dimensions and fields.
+  //! \param[in] patch The spline patch on which the field is to be defined
+  //! \param[in] v Array of control point field values
   //! \param[in] name Name of spline field
-  SplineFields2D(Go::SplineSurface *bf, char* name = NULL);
-  //! \brief The constructor sets the field name.
-  //! \param[in] bf Spline basis description
-  //! \param[in] geometry Spline geometry description
-  //! \param[in] name Name of spline field
-  SplineFields2D(Go::SplineSurface *bf, 
-		 Go::SplineSurface *geoemtry, 
-		 char* name = NULL);
+  SplineFields2D(const ASMs2D* patch, const RealArray& v,
+                 const char* name = NULL);
   //! \brief Empty destructor.
-  virtual ~SplineFields2D();
+  virtual ~SplineFields2D() {}
 
-  // Methods to compute field values
-  //================================
+  // Methods to evaluate the field
+  //==============================
 
   //! \brief Computes the value in a given node/control point.
   //! \param[in] node Node number
   //! \param[out] vals Node values
-  bool valueNode(int node, Vector& vals) const;
+  bool valueNode(size_t node, Vector& vals) const;
 
   //! \brief Computes the value at a given local coordinate.
   //! \param[in] fe Finite element definition
   //! \param[out] vals Values in local point in given element
   bool valueFE(const FiniteElement& fe, Vector& vals) const;
 
-  //! \brief Computed the value at a given global coordinate.
+  //! \brief Computes the value at a given global coordinate.
   //! \param[in] x Global/physical coordinate for point
-  //! \param[in] vals Values in given physical coordinate
+  //! \param[out] vals Values in given physical coordinate
   bool valueCoor(const Vec3& x, Vector& vals) const;
 
   //! \brief Computes the gradient for a given local coordinate.
@@ -75,8 +71,8 @@ public:
   bool gradCoor(const Vec3& x, Matrix& grad) const;
 
 protected:
-  Go::SplineSurface* basis; //!< Spline basis description
-  Go::SplineSurface* surf;  //!< Spline geoemtry description
+  const Go::SplineSurface* basis; //!< Spline basis description
+  const Go::SplineSurface* surf;  //!< Spline geometry description
 };
 
 #endif
diff --git a/src/ASM/SplineFields3D.C b/src/ASM/SplineFields3D.C
index dc039edb..c5a6e1f9 100644
--- a/src/ASM/SplineFields3D.C
+++ b/src/ASM/SplineFields3D.C
@@ -12,20 +12,18 @@
 //==============================================================================
 
 #include "SplineFields3D.h"
+#include "ASMs3D.h"
 #include "FiniteElement.h"
+#include "CoordinateMapping.h"
+#include "Utilities.h"
 #include "Vec3.h"
 
 #include "GoTools/trivariate/SplineVolume.h"
-#include "GoTools/geometry/SplineUtils.h"
-#include "boost/lambda/lambda.hpp"
-
-namespace Go {
-  void volume_ratder(double const eder[],int idim,int ider,double gder[]);
-}
 
 
-SplineFields3D::SplineFields3D(Go::SplineVolume *bf, char* name)
-  : Fields(3,name), basis(bf), vol(bf) 
+SplineFields3D::SplineFields3D (const ASMs3D* patch,
+                                const RealArray& v, const char* name)
+  : Fields(3,name), basis(patch->getBasis()), vol(patch->getVolume())
 {
   const int n1 = basis->numCoefs(0);
   const int n2 = basis->numCoefs(1);
@@ -37,56 +35,48 @@ SplineFields3D::SplineFields3D(Go::SplineVolume *bf, char* name)
   const int p3 = basis->order(2);
   nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1);
 
-  // Number of fields set in fill
-  nf = 0;
+  // Ensure the values array has compatible length, pad with zeros if necessary
+  nf = v.size()/nno;
+  values.resize(nf*nno);
+  RealArray::const_iterator end = v.size() > nf*nno ? v.begin()+nf*nno:v.end();
+  std::copy(v.begin(),end,values.begin());
 }
 
 
-SplineFields3D::SplineFields3D(Go::SplineVolume *bf, 
-			       Go::SplineVolume *geometry, 
-			       char* name)
-  : Fields(3,name), basis(bf), vol(geometry) 
+bool SplineFields3D::valueNode (size_t node, Vector& vals) const
 {
-  const int n1 = basis->numCoefs(0);
-  const int n2 = basis->numCoefs(1);
-  const int n3 = basis->numCoefs(2);
-  nno = n1*n2*n3;
+  if (node < 1 || node > nno) return false;
 
-  const int p1 = basis->order(0);
-  const int p2 = basis->order(1);
-  const int p3 = basis->order(2);
-  nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1);
-
-  // Number of fields set in fill
-  nf = 0;
+  vals.resize(nf);
+  vals.fill(values.ptr()+(node-1)*nf);
+  return true;
 }
 
 
-
-SplineFields3D::~SplineFields3D()
-{
-  // Set geometry to NULL; should not be deallocated here
-  basis = NULL;
-} 
-
-
-bool SplineFields3D::valueNode(int node, Vector& vals) const
-{
-  // Not implemented yet
-  return false;
-}
-
-
-bool SplineFields3D::valueFE(const FiniteElement& fe, Vector& vals) const
+bool SplineFields3D::valueFE (const FiniteElement& fe, Vector& vals) const
 {
   if (!basis) return false;
 
-  Go::Point pt;
+  // Evaluate the basis functions at the given point
+  Go::BasisPts spline;
+  basis->computeBasis(fe.u,fe.v,fe.w,spline);
 
-//   basis->pointValue(pt,values,fe.u,fe.v,fe.w);
-//   vals.resize(pt.size());
-//   for (int i = 0;i < pt.size();i++)
-//     vals[i] = pt[i];
+  // Evaluate the solution field at the given point
+  std::vector<int> ip;
+  ASMs3D::scatterInd(basis->numCoefs(0),basis->numCoefs(1),basis->numCoefs(2),
+		     basis->order(0),basis->order(1),basis->order(2),
+		     spline.left_idx,ip);
+
+  Matrix Vnod;
+  utl::gather(ip,nf,values,Vnod);
+  Vnod.multiply(spline.basisValues,vals); // vals = Vnod * basisValues
+
+  return true;
+}
+
+/* Old procedure, way too complex...
+   The above code is more in line with how we do it in ASMs3D.
+  Go::Point pt;
 
   const int uorder = basis->order(0);
   const int vorder = basis->order(1);
@@ -213,10 +203,10 @@ bool SplineFields3D::valueFE(const FiniteElement& fe, Vector& vals) const
 
   return true;
 }
+*/
 
 
-
-bool SplineFields3D::valueCoor(const Vec3& x, Vector& vals) const
+bool SplineFields3D::valueCoor (const Vec3& x, Vector& vals) const
 {
   // Not implemented yet
   return false;
@@ -226,7 +216,73 @@ bool SplineFields3D::valueCoor(const Vec3& x, Vector& vals) const
 bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const
 {
   if (!basis) return false;
-  
+  if (!vol)   return false;
+
+  // Evaluate the basis functions at the given point
+/*
+  Go::BasisDerivs spline;
+  vol->computeBasis(fe.u,fe.v,fe.w,spline);
+  TODO: The above is not available yet, the below is temporary workaround.
+*/
+  std::vector<Go::BasisDerivs> tmpSpline(1);
+  vol->computeBasisGrid(RealArray(1,fe.u),RealArray(1,fe.v),RealArray(1,fe.w),tmpSpline);
+  Go::BasisDerivs& spline = tmpSpline.front();
+
+  const int uorder = vol->order(0);
+  const int vorder = vol->order(1);
+  const int worder = vol->order(2);
+  const size_t nen = uorder*vorder*worder;
+
+  Matrix dNdu(nen,3), dNdX;
+  for (size_t n = 1; n <= nen; n++)
+  {
+    dNdu(n,1) = spline.basisDerivs_u[n-1];
+    dNdu(n,2) = spline.basisDerivs_v[n-1];
+    dNdu(n,3) = spline.basisDerivs_w[n-1];
+  }
+
+  std::vector<int> ip;
+  ASMs3D::scatterInd(vol->numCoefs(0),vol->numCoefs(1),vol->numCoefs(2),
+		     uorder,vorder,worder,spline.left_idx,ip);
+
+  // Evaluate the Jacobian inverse
+  Matrix Xnod, Jac;
+  Vector Xctrl(&(*vol->coefs_begin()),vol->coefs_end()-vol->coefs_begin());
+  utl::gather(ip,vol->dimension(),Xctrl,Xnod);
+  utl::Jacobian(Jac,dNdX,Xnod,dNdu);
+
+  // Evaluate the gradient of the solution field at the given point
+  if (basis != vol)
+  {
+    // Mixed formulation, the solution uses a different basis than the geometry
+    /*
+    basis->computeBasis(fe.u,fe.v,fe.w,spline);
+    TODO: The above is not available yet, the below is temporary workaround.
+    */
+    basis->computeBasisGrid(RealArray(1,fe.u),RealArray(1,fe.v),RealArray(1,fe.w),tmpSpline);
+    Go::BasisDerivs& spline = tmpSpline.front();
+
+    const size_t nbf = basis->order(0)*basis->order(1)*basis->order(2);
+    dNdu.resize(nbf,3);
+    for (size_t n = 1; n <= nbf; n++)
+    {
+      dNdu(n,1) = spline.basisDerivs_u[n-1];
+      dNdu(n,2) = spline.basisDerivs_v[n-1];
+      dNdu(n,3) = spline.basisDerivs_w[n-1];
+    }
+    dNdX.multiply(dNdu,Jac); // dNdX = dNdu * Jac
+
+    ASMs3D::scatterInd(basis->numCoefs(0),basis->numCoefs(1),basis->numCoefs(2),
+		       basis->order(0),basis->order(1),basis->order(2),
+		       spline.left_idx,ip);
+  }
+
+  utl::gather(ip,nf,values,Xnod);
+  grad.multiply(Xnod,dNdX); // grad = Xnod * dNdX
+
+  return true;
+}
+
 //   // Derivatives of solution in parametric space
 //   std::vector<Go::Point> pts;
 //   basis->pointValue(pts,values,fe.u,fe.v,fe.w,1);
@@ -252,6 +308,7 @@ bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const
 //   // df/dX = df/dXi * Jac^-1 = df/dXi * [dX/dXi]^-1
 //   grad.multiply(gradXi,Jac);
 
+/*
   // Derivatives of solution in parametric space
   std::vector<Go::Point> pts(4);
 
@@ -438,9 +495,9 @@ bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const
 
   return true;
 }
+*/
 
-
-bool SplineFields3D::gradCoor(const Vec3& x, Matrix& grad) const
+bool SplineFields3D::gradCoor (const Vec3& x, Matrix& grad) const
 {
   // Not implemented yet
   return false;
diff --git a/src/ASM/SplineFields3D.h b/src/ASM/SplineFields3D.h
index e8466659..b0a9bc66 100644
--- a/src/ASM/SplineFields3D.h
+++ b/src/ASM/SplineFields3D.h
@@ -16,6 +16,8 @@
 
 #include "Fields.h"
 
+class ASMs3D;
+
 namespace Go {
   class SplineVolume;
 }
@@ -24,44 +26,38 @@ namespace Go {
 /*!
   \brief Class for spline-based finite element vector fields in 3D.
 
-  \details This class implements the functions required to evaluate a 3D
+  \details This class implements the methods required to evaluate a 3D
   spline vector field at a given point in parametrical or physical coordinates.
 */
 
-
 class SplineFields3D : public Fields
 {
 public:
-  //! \brief The constructor sets the field name.
-  //! \param[in] bf Spline basis description
+  //! \brief The constructor sets the number of space dimensions and fields.
+  //! \param[in] patch The spline patch on which the field is to be defined
+  //! \param[in] v Array of control point field values
   //! \param[in] name Name of spline field
-  SplineFields3D(Go::SplineVolume *bf, char* name = NULL);
-  //! \brief The constructor sets the field name.
-  //! \param[in] bf Spline basis description
-  //! \param[in] geometry Spline geometry description
-  //! \param[in] name Name of spline field
-  SplineFields3D(Go::SplineVolume *bf, 
-		 Go::SplineVolume *geometry, 
-		 char* name = NULL);
+  SplineFields3D(const ASMs3D* patch, const RealArray& v,
+                 const char* name = NULL);
   //! \brief Empty destructor.
-  virtual ~SplineFields3D();
+  virtual ~SplineFields3D() {}
 
-  // Methods to compute field values
-  //================================
+  // Methods to evaluate the field
+  //==============================
 
   //! \brief Computes the value in a given node/control point.
   //! \param[in] node Node number
   //! \param[out] vals Node values
-  bool valueNode(int node, Vector& vals) const;
+  bool valueNode(size_t node, Vector& vals) const;
 
   //! \brief Computes the value at a given local coordinate.
   //! \param[in] fe Finite element definition
   //! \param[out] vals Values in local point in given element
   bool valueFE(const FiniteElement& fe, Vector& vals) const;
 
-  //! \brief Computed the value at a given global coordinate.
+  //! \brief Computes the value at a given global coordinate.
   //! \param[in] x Global/physical coordinate for point
-  //! \param[in] vals Values in given physical coordinate
+  //! \param[out] vals Values in given physical coordinate
   bool valueCoor(const Vec3& x, Vector& vals) const;
 
   //! \brief Computes the gradient for a given local coordinate.
@@ -75,8 +71,8 @@ public:
   bool gradCoor(const Vec3& x, Matrix& grad) const;
 
 protected:
-  Go::SplineVolume* basis;    //!< Spline basis description
-  Go::SplineVolume* vol;      //!< Spline geometry description
+  const Go::SplineVolume* basis; //!< Spline basis description
+  const Go::SplineVolume* vol;   //!< Spline geometry description
 };
 
 #endif