Representation of scalar field over spline geometry

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@950 e10b68d5-8a6e-419e-a041-bce267b0401d

src/ASM/SplineField.h

@@ -0,0 +1,103 @@
+//==============================================================================
+//!
+//! \file SplineField.h
+//!
+//! \date Mar 28 2011
+//!
+//! \author Runar Holdahl / SINTEF
+//!
+//! \brief Base class for spline based finite element field
+//!
+//==============================================================================
+
+#ifndef _SPLINE_FIELD_H
+#define _SPLINE_FIELD_H
+
+#include "Vec3.h"
+#include "MatVec.h"
+#include "FiniteElement.h"
+
+/* 
+   \brief Base class for spline-based finite element fields.
+
+   \details This class incapsulates the data and methods needed
+   to store and evaluate a spline finite element scalar field.
+   This is an abstract base class and the fields associated with 
+   specific spline objects are implemented as subclasses, for
+   instance 1D, 2D and 3D spline formulations.
+*/
+
+
+class SplineField
+{
+ protected:
+  //! \brief The constructor sets the number of space dimensions and fields.
+  //! \param[in] nsd Number of space dimensions (1, 2 or 3)
+  //! \param[in] name Name of spline field
+ SplineField(unsigned char nsd_, char* name = NULL) 
+   : fieldname(name) {}
+
+ public:
+  //! \brief Empty destructor
+  virtual ~SplineField() {}
+
+  // Returns number of space dimensions
+  unsigned char getNoSpaceDim() const { return nsd; }
+
+  // Returns number of space dimensions
+  int getNoElm() const { return nelm; }
+
+  // Returns number of control points
+  int getNoNodes() const { return nno; }
+
+  // Returns name of spline field
+  const char* getFieldName() const { return fieldname; }
+
+  // Sets the name of the spline field
+  void setFieldName(char* name) { fieldname = name; }
+
+  // Methods to initialize field
+  virtual void fill(Vector& vec) { values = vec; }
+
+  
+  // Methods to compute field values
+  //================================================
+
+  //! \brief Computes the value in a given node/control point
+  //! \param[in] node Node number 
+  virtual double valueNode(int node) const = 0;
+
+  //! \brief Computes the value at a given local coordinate
+  //! \param[in] fe Finite element definition
+  virtual double valueFE(const FiniteElement& fe) const = 0;
+  
+  //! \brief Computed the value at a given global coordinate
+  //! \param[in] x Global/physical coordinate for point
+  virtual double valueCoor(const Vec3 x) const = 0;
+
+  //! \brief Computes the gradient for a given local coordinate
+  //! \param[in] fe Finite element
+  //! \param[out] grad Gradient of solution in a given local coordinate
+  virtual bool gradFE(const FiniteElement& fe, Vector& grad) const = 0;
+
+  //! \brief Computes the gradient for a given global/physical coordinate
+  //! \param[in] x Global coordinate
+  //! \param[out] grad Gradient of solution in a given global coordinate
+  virtual bool gradCoor(const Vec3 x, Vector& grad) const = 0;
+
+ protected:
+  // Dimension of field
+  unsigned char nsd;  //!< Number of space dimensions
+  int nelm;           //!< Number of elements/knot-spans
+  int nno;            //!< Number of nodes/control points
+
+  // Fieldname
+  char* fieldname;    //!< Name of spline element field
+
+  // Field values
+  Vector values;      //!< Field values
+};
+   
+
+
+#endif

src/ASM/SplineField2D.C

@@ -0,0 +1,308 @@
+//==============================================================================
+//!
+//! \file SplineField2D.C
+//!
+//! \date Mar 28 2011
+//!
+//! \author Runar Holdahl / SINTEF
+//!
+//! \brief Base class for spline based finite element field in 2D
+//!
+//==============================================================================
+
+#include "SplineField2D.h"
+#include <algorithm>
+#include "GoTools/geometry/SplineUtils.h"
+#ifndef __BORLANDC__
+#include "boost/lambda/lambda.hpp"
+#endif
+
+using namespace std;
+#ifndef __BORLANDC__
+using namespace boost::lambda;
+#endif
+
+SplineField2D::SplineField2D(Go::SplineSurface *geometry, char* name)
+  : SplineField(2,name), surf(geometry) 
+{
+  const int n1 = surf->numCoefs_u();
+  const int n2 = surf->numCoefs_v();
+  nno = n1*n2;
+
+  const int p1 = surf->order_u();
+  const int p2 = surf->order_v();
+  nelm = (n1-p1+1)*(n2-p2+1);
+}
+
+
+SplineField2D::~SplineField2D()
+{
+  // Set geometry pointer to NULL, should not be deallocated here
+  surf = NULL;
+}
+
+
+double SplineField2D::valueNode(int node) const
+{
+  // Not implemented yet
+  return 0.0;
+}
+
+
+// double SplineField2D::valueFE(const FiniteElement& fe) const
+// {
+//   Go::Point pt;
+//   surf->pointValue(pt,values,fe.u,fe.v);
+  
+//   return pt[0];
+// }
+
+
+double SplineField2D::valueFE(const FiniteElement& fe) const
+{
+  const int uorder = surf->order_u();
+  const int vorder = surf->order_v();
+  const int unum = surf->numCoefs_u();
+  //const int vnum = surf->numCoefs_v();
+
+  const int dim = surf->dimension();
+  const bool rational = surf->rational();
+  const int kdim = rational ? dim + 1 : dim;
+
+  //const int ncomp = values.size()/(unum*vnum);
+  double val = 0.0;
+  double tempVal;
+ 
+  static Go::ScratchVect<double, 10> Bu(uorder);
+  static Go::ScratchVect<double, 10> Bv(vorder);
+
+  Bu.resize(uorder);
+  Bv.resize(vorder);
+
+  surf->basis_u().computeBasisValues(fe.u,Bu.begin());
+  surf->basis_v().computeBasisValues(fe.v,Bv.begin());
+  
+  const int uleft = surf->basis_u().lastKnotInterval();
+  const int vleft = surf->basis_v().lastKnotInterval();
+
+  // compute the tensor product value
+  const int val_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1));
+
+  register const double* val_ptr = &values[val_start_ix]; 
+
+  if (rational) {
+      double w = 0.0;
+      const int w_start_ix  = (uleft - uorder + 1 + unum * (vleft - vorder + 1)) * kdim + dim;
+      register const double* w_ptr  = &(surf->rcoefs_begin()[w_start_ix]);
+      
+      for (register double* bval_v_ptr = Bv.begin(); bval_v_ptr != Bv.end(); ++bval_v_ptr) {
+        register const double bval_v = *bval_v_ptr;
+	tempVal = 0.0;
+        for (register double* bval_u_ptr = Bu.begin(); bval_u_ptr != Bu.end(); ++bval_u_ptr) {
+          register const double bval_u = *bval_u_ptr;
+	  tempVal += bval_u * (*val_ptr++);
+
+          w += bval_u * (*w_ptr);
+          w_ptr += kdim;
+        }
+
+	val += tempVal * bval_v;
+	val_ptr += unum - uorder;
+
+        w_ptr   += kdim * (unum - uorder);
+      }
+
+      val /= w;
+  }
+  else {
+    for (register double* bval_v_ptr = Bv.begin(); bval_v_ptr != Bv.end(); ++bval_v_ptr) {
+      register const double bval_v = *bval_v_ptr;
+      tempVal = 0.0;
+      for (register double* bval_u_ptr = Bu.begin(); bval_u_ptr != Bu.end(); ++bval_u_ptr) {
+	register const double bval_u = *bval_u_ptr;
+	tempVal += bval_u * (*val_ptr++);
+      }
+
+      val += tempVal * bval_v;
+      val_ptr += unum - uorder;
+    }
+  } 
+  
+  return val;
+}
+
+
+double SplineField2D::valueCoor(const Vec3 x) const
+{
+  // Not implemented yet
+  return 0.0;
+}
+
+
+// bool SplineField2D::gradFE(const FiniteElement& fe, Vector& grad) const
+// {
+//   if (!surf) return false;
+
+//   // Derivatives of solution in parametric space
+//   std::vector<Go::Point> pts(3);
+//   surf->pointValue(pts,values,fe.u,fe.v,1);
+
+//   // Gradient of field wrt parametric coordinates
+//   Vector gradXi(2);
+//   gradXi(1) = pts[1][0];
+//   gradXi(2) = pts[2][0];
+
+//   // Derivatives of coordinate mapping 
+//   std::vector<Go::Point> xpts(3);
+//   surf->point(xpts,fe.u,fe.v,1);
+
+//   // Jacobian matrix
+//   Matrix Jac(2,2);
+//   for (size_t i = 1;i <= 2;i++)
+//     for (size_t n = 0;n < 2;n++)
+//       Jac(n+1,i) = xpts[i][n];
+
+//   // Invert Jacobian matrix
+//   Jac.inverse();
+
+//   // Gradient of solution in given parametric point
+//   // df/dX = Jac^-T * df/dXi = [dX/dXi]^-T * df/dXi
+//   Jac.multiply(gradXi,grad,true,false);
+  
+//   return true;
+// }
+
+
+bool SplineField2D::gradFE(const FiniteElement& fe, Vector& grad) const
+{
+  if (!surf) return false;
+
+  // Gradient of field wrt parametric coordinates
+  Vector gradXi(2);
+
+  const int uorder = surf->order_u();
+  const int vorder = surf->order_v();
+  const int unum = surf->numCoefs_u();
+
+  const int derivs = 1;
+  const int totpts = 3;
+  const bool u_from_right = true;
+  const bool v_from_right = true;
+  const double resolution = 1.0e-12;
+
+  // Dimension
+  const int dim       = surf->dimension();
+  const bool rational = surf->rational();
+  
+  // Take care of the rational case
+  const std::vector<double>::iterator co = rational ? surf->rcoefs_begin() : surf->coefs_begin();
+  int kdim = dim   + (rational ? 1 : 0);
+  int ndim = 1 + (rational ? 1 : 0); 
+  
+  // Make temporary storage for the basis values and a temporary
+  // computation cache.
+  Go::ScratchVect<double, 30> b0(uorder*(derivs+1));
+  Go::ScratchVect<double, 30> b1(vorder*(derivs+1));
+  Go::ScratchVect<double, 30> temp(ndim * totpts);
+  Go::ScratchVect<double, 30> restemp(ndim * totpts);
+  std::fill(restemp.begin(), restemp.end(), 0.0);
+
+  // Compute the basis values and get some data about the spline spaces
+  if (u_from_right) 
+    surf->basis_u().computeBasisValues(fe.u, &b0[0], derivs, resolution);
+  else 
+    surf->basis_u().computeBasisValuesLeft(fe.u, &b0[0], derivs, resolution);  
+  int uleft  = surf->basis_u().lastKnotInterval();
+
+  if (v_from_right) 
+    surf->basis_v().computeBasisValues(fe.v, &b1[0], derivs, resolution);
+  else 
+    surf->basis_v().computeBasisValuesLeft(fe.v, &b1[0], derivs, resolution);  
+  int vleft = surf->basis_v().lastKnotInterval();
+
+  // Compute the tensor product value
+  int coefind = uleft-uorder+1 + unum*(vleft-vorder+1);
+  int derivs_plus1=derivs+1;
+  for (int jj = 0; jj < vorder; ++jj) {
+    int jjd=jj*(derivs_plus1);
+    std::fill(temp.begin(), temp.end(), 0.0);
+    
+    for (int ii = 0; ii < uorder; ++ii) {
+      int iid=ii*(derivs_plus1);
+      const double *val_p = &values[coefind];
+      const double *co_p  = &co[coefind*kdim] + dim;
+      int temp_ind = 0;
+      for (int vder = 0; vder < derivs_plus1; ++vder) 
+	for (int uder = 0; uder < vder+1; ++uder) {
+	  temp[temp_ind] += b0[iid+vder - uder]*(*val_p);
+	  temp_ind += ndim;
+	}
+      
+      for (int dd = dim;dd < kdim;dd++,co_p += kdim) {
+	int temp_ind = 1;
+        for (int vder = 0; vder < derivs_plus1; ++vder) {
+          for (int uder = 0; uder < vder+1; ++uder) {
+            temp[temp_ind] += b0[iid+vder - uder]*(*co_p);
+            temp_ind += ndim;
+          }
+        }
+      }
+
+      coefind += 1;
+    }
+    
+    for (int dd = 0; dd < ndim; ++dd) {
+      int dercount = 0;
+      for (int vder = 0; vder < derivs_plus1; ++vder) {
+        for (int uder = 0; uder < vder + 1; ++uder) {
+          restemp[dercount*ndim + dd] 
+            += temp[dercount*ndim + dd]*b1[uder + jjd];
+          ++dercount;
+        }
+      }
+    }
+    
+    coefind += unum - uorder;
+  }
+
+  // Copy from restemp to result
+  if (rational) {
+    std::vector<double> restemp2(totpts);
+    Go::surface_ratder(&restemp[0],1, derivs, &restemp2[0]);
+    for (int i = 1; i < totpts; ++i) 
+      gradXi(i) = restemp2[i];
+  } else {
+    double* restemp_it=restemp.begin();
+    ++restemp_it;
+    for (int i = 1; i < totpts; ++i) {
+      gradXi(i) = *restemp_it;
+      ++restemp_it;
+    }
+  }  
+
+  // Derivatives of coordinate mapping 
+  std::vector<Go::Point> xpts(3);
+  surf->point(xpts,fe.u,fe.v,1);
+
+  // Jacobian matrix
+  Matrix Jac(2,2);
+  for (size_t i = 1;i <= 2;i++)
+    for (size_t n = 0;n < 2;n++)
+      Jac(n+1,i) = xpts[i][n];
+
+  // Invert Jacobian matrix
+  Jac.inverse();
+
+  // Gradient of solution in given parametric point
+  // df/dX = Jac^-T * df/dXi = [dX/dXi]^-T * df/dXi
+  Jac.multiply(gradXi,grad,true,false);
+  
+  return true;
+}
+
+
+bool SplineField2D::gradCoor(const Vec3 x, Vector& grad) const
+{
+  // Not implemented yet
+  return false;
+}

src/ASM/SplineField2D.h

@@ -0,0 +1,76 @@
+//==============================================================================
+//!
+//! \file SplineField2D.h
+//!
+//! \date Mar 28 2011
+//!
+//! \author Runar Holdahl / SINTEF
+//!
+//! \brief Base class for spline based finite element field in 2D
+//!
+//==============================================================================
+
+#ifndef _SPLINE_FIELD_2D_H
+#define _SPLINE_FIELD_2D_H
+
+namespace Go {
+  class SplineSurface;
+  class Point;
+}
+
+#include "SplineField.h"
+#include "GoTools/geometry/SplineSurface.h"
+
+/* 
+   \brief Base class for spline-based finite element fields in 2D.
+
+   \details This class implements the functions required to evaluate
+   a 2D spline field at a given point in parametrical or physical 
+   coordinates.
+*/
+
+
+class SplineField2D : public SplineField
+{
+ public:
+  //! \brief The constructor sets the number of space dimensions and fields.
+  //! \param[in] geometry Spline volume geometry
+  //! \param[in] name Name of spline field
+  SplineField2D(Go::SplineSurface *geometry, char* name = NULL); 
+
+  //! \brief Empty destructor
+  virtual ~SplineField2D();
+
+  // Methods to compute field values
+  //================================================
+
+  //! \brief Computes the value in a given node/control point
+  //! \param[in] node Node number 
+  double valueNode(int 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
+  //! \param[in] x Global/physical coordinate for point
+  double valueCoor(const Vec3 x) const;
+
+  //! \brief Computes the gradient for a given local coordinate
+  //! \param[in] fe Finite element
+  //! \param[out] grad Gradient of solution in a given local coordinate
+  bool gradFE(const FiniteElement& fe, Vector& grad) const;
+
+  //! \brief Computes the gradient for a given global/physical coordinate
+  //! \param[in] x Global coordinate
+  //! \param[out] grad Gradient of solution in a given global coordinate
+  bool gradCoor(const Vec3 x, Vector& grad) const;
+
+ protected:
+  Go::SplineSurface *surf;    //!< Spline surface geometry description
+
+};
+   
+
+
+#endif

src/ASM/SplineField3D.C

@@ -0,0 +1,377 @@
+//==============================================================================
+//!
+//! \file SplineField3D.C
+//!
+//! \date Mar 28 2011
+//!
+//! \author Runar Holdahl / SINTEF
+//!
+//! \brief Base class for spline based finite element field in 3D
+//!
+//==============================================================================
+
+#include "SplineField3D.h"
+#include <algorithm>
+#include "GoTools/geometry/SplineUtils.h"
+#ifndef __BORLANDC__
+#include "boost/lambda/lambda.hpp"
+#endif
+
+using namespace std;
+#ifndef __BORLANDC__
+using namespace boost::lambda;
+#endif
+
+SplineField3D::SplineField3D(Go::SplineVolume *geometry, char* name)
+  : SplineField(3,name), vol(geometry) 
+{
+  const int n1 = vol->numCoefs(0);
+  const int n2 = vol->numCoefs(1);
+  const int n3 = vol->numCoefs(2);
+  nno = n1*n2*n3;
+
+  const int p1 = vol->order(0);
+  const int p2 = vol->order(1);
+  const int p3 = vol->order(2);
+  nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1);
+}
+
+
+SplineField3D::~SplineField3D()
+{
+  // Set geometry pointer to NULL, should not be deallocated here
+  vol = NULL;
+}
+
+
+double SplineField3D::valueNode(int node) const
+{
+  // Not implemented yet
+  return 0.0;
+}
+
+
+// double SplineField3D::valueFE(const FiniteElement& fe) const
+// {
+//   Go::Point pt;
+//   vol->pointValue(pt,values,fe.u,fe.v,fe.w);
+  
+//   return pt[0];
+// }
+
+
+double SplineField3D::valueFE(const FiniteElement& fe) const
+{
+  double val = 0.0;
+
+  const int uorder = vol->order(0);
+  const int vorder = vol->order(1);
+  const int worder = vol->order(2);
+  const int unum = vol->numCoefs(0);
+  const int vnum = vol->numCoefs(1);
+
+  const int dim = vol->dimension();
+  const bool rational = vol->rational();
+  int kdim = rational ? dim + 1 : dim;
+  
+  static Go::ScratchVect<double, 10> Bu(uorder);
+  static Go::ScratchVect<double, 10> Bv(vorder);
+  static Go::ScratchVect<double, 10> Bw(worder);
+  double temp, temp2, tempResult = 0.0;
+  
+  Bu.resize(uorder);
+  Bv.resize(vorder);
+  Bw.resize(worder);
+  
+  // compute tbe basis values and get some data about the spline spaces
+  vol->basis(0).computeBasisValues(fe.u, Bu.begin());
+  vol->basis(1).computeBasisValues(fe.v, Bv.begin());
+  vol->basis(2).computeBasisValues(fe.w, Bw.begin());
+  const int uleft = vol->basis(0).lastKnotInterval();
+  const int vleft = vol->basis(1).lastKnotInterval();
+  const int wleft = vol->basis(2).lastKnotInterval();
+
+  // compute the tensor product value
+  const int val_start_ix = uleft - uorder + 1 + unum * (vleft - vorder + 1 + vnum * (wleft - worder + 1));
+  const double* val_ptr = &values[val_start_ix];
+  
+  double w = 1.0;
+  if (rational) {
+    w = 0.0;
+    const int w_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1 + vnum * (wleft - worder + 1))) * kdim + dim;
+    const double *w_ptr = &(vol->rcoefs_begin()[w_start_ix]);
+    
+    for (double* bval_w_ptr = Bw.begin(); bval_w_ptr != Bw.end(); ++bval_w_ptr) {
+      const double bval_w = *bval_w_ptr;
+      temp = 0.0;
+      for (double* bval_v_ptr = Bv.begin(); bval_v_ptr != Bv.end(); ++bval_v_ptr) {
+	const double bval_v = *bval_v_ptr;
+	temp2 = 0.0;
+	for (double* bval_u_ptr = Bu.begin(); bval_u_ptr != Bu.end(); ++bval_u_ptr) {
+	  const double bval_u = *bval_u_ptr;
+	  temp2 += bval_u * (*val_ptr++);
+
+	  w += bval_u*(*w_ptr);
+	  w_ptr += kdim;
+	}
+
+	temp += temp2 * bval_v;
+
+	val_ptr += unum - uorder;
+	w_ptr   += kdim * (unum - uorder);
+      }
+
+      tempResult += temp * bval_w;
+      w *= bval_w;
+
+      val_ptr += unum * (vnum - vorder);
+      w_ptr   += kdim  * unum * (vnum - vorder);
+    }
+  }
+  else {
+    for (double* bval_w_ptr = Bw.begin(); bval_w_ptr != Bw.end(); ++bval_w_ptr) {
+      const double bval_w = *bval_w_ptr;
+      temp = 0.0;
+      for (double* bval_v_ptr = Bv.begin(); bval_v_ptr != Bv.end(); ++bval_v_ptr) {
+	const double bval_v = *bval_v_ptr;
+	temp2 = 0.0;
+	for (double* bval_u_ptr = Bu.begin(); bval_u_ptr != Bu.end(); ++bval_u_ptr) {
+	  const double bval_u = *bval_u_ptr;
+	  temp2 += bval_u * (*val_ptr++);
+	}
+
+	temp += temp2 * bval_v;
+	val_ptr += unum - uorder;
+      }
+
+      tempResult += temp * bval_w;
+      val_ptr += unum * (vnum - vorder);
+    }
+  }
+  
+  val = tempResult/w;
+  return val;
+}
+
+
+double SplineField3D::valueCoor(const Vec3 x) const
+{
+  // Not implemented yet
+  return 0.0;
+}
+
+
+// bool SplineField3D::gradFE(const FiniteElement& fe, Vector& grad) const
+// {
+//   if (!vol) return false;
+
+//   // Derivatives of solution in parametric space
+//   std::vector<Go::Point> pts;
+//   vol->pointValue(pts,values,fe.u,fe.v,fe.w,1);
+
+//   // Gradient of field wrt parametric coordinates
+//   Vector gradXi(3);
+//   gradXi(1) = pts[1][0];
+//   gradXi(2) = pts[2][0];
+//   gradXi(3) = pts[3][0];
+
+//   // Derivatives of coordinate mapping 
+//   std::vector<Go::Point> xpts(3);
+//   vol->point(xpts,fe.u,fe.v,fe.w,1);
+
+//   // Jacobian matrix
+//   Matrix Jac(3,3);
+//   for (size_t i = 1;i <= 3;i++)
+//     for (size_t n = 0;n < 3;n++)
+//       Jac(n+1,i) = xpts[i][n];
+
+//   // Invert Jacobian matrix
+//   Jac.inverse();
+
+//   // Gradient of solution in given parametric point
+//   // df/dX = Jac^-T * df/dXi = [dX/dXi]^-T * df/dXi
+//   Jac.multiply(gradXi,grad,true,false);
+
+//   return true;
+// }
+
+
+bool SplineField3D::gradFE(const FiniteElement& fe, Vector& grad) const
+{
+  if (!vol) return false;
+
+  // Gradient of field wrt parametric coordinates
+  Vector gradXi(3);
+
+  int derivs = 1;
+  int totpts = (derivs + 1)*(derivs + 2)*(derivs + 3)/6;
+  
+  const int unum = vol->numCoefs(0);
+  const int vnum = vol->numCoefs(1);
+  const int wnum = vol->numCoefs(2);
+
+  const int uorder = vol->order(0);
+  const int vorder = vol->order(1);
+  const int worder = vol->order(2);
+
+  const bool u_from_right = true;
+  const bool v_from_right = true;
+  const bool w_from_right = true;
+
+  const double resolution = 1.0e-12;
+  
+  const int ncomp     = values.size()/(unum*vnum*wnum);
+  const int dim       = vol->dimension();
+  const bool rational = vol->rational();
+  
+  // Take care of the rational case
+  const vector<double>::iterator co = rational ? vol->rcoefs_begin() : vol->coefs_begin();
+  int kdim = dim + (rational ? 1 : 0);
+  int ndim = ncomp + (rational ? 1 : 0);
+  
+  // Make temporary storage for the basis values and a temporary
+  // computation cache.
+  Go::ScratchVect<double, 30> b0(uorder * (derivs+1));
+  Go::ScratchVect<double, 30> b1(vorder * (derivs+1));
+  Go::ScratchVect<double, 30> b2(worder * (derivs+1));
+  Go::ScratchVect<double, 60> temp(ndim * totpts);
+  Go::ScratchVect<double, 60> temp2(ndim * totpts);
+  Go::ScratchVect<double, 60> restemp(ndim * totpts);
+  std::fill(restemp.begin(), restemp.end(), 0.0);
+  
+  // Compute the basis values and get some data about the spline spaces
+  if (u_from_right) 
+    vol->basis(0).computeBasisValues(fe.u, &b0[0], derivs, resolution);
+  else 
+    vol->basis(1).computeBasisValuesLeft(fe.u, &b0[0], derivs, resolution);
+  
+  int uleft = vol->basis(0).lastKnotInterval();
+  if (v_from_right) 
+    vol->basis(1).computeBasisValues(fe.v, &b1[0], derivs, resolution);
+  else 
+    vol->basis(1).computeBasisValuesLeft(fe.v, &b1[0], derivs, resolution);
+  
+  int vleft = vol->basis(1).lastKnotInterval();
+  if (w_from_right) 
+    vol->basis(2).computeBasisValues(fe.w, &b2[0], derivs, resolution);
+  else 
+    vol->basis(2).computeBasisValuesLeft(fe.w, &b2[0], derivs, resolution);
+  
+  int wleft = vol->basis(2).lastKnotInterval();
+  
+  // Compute the tensor product value
+  int coefind = uleft-uorder+1 + unum*(vleft-vorder+1 + vnum*(wleft-worder+1));
+  int derivs_plus1=derivs+1;
+  
+  for (int k = 0; k < worder; ++k) {
+    int kd=k*derivs_plus1;
+    std::fill(temp.begin(), temp.end(), 0.0);
+      
+    for (int j = 0; j < vorder; ++j) {
+      int jd=j*derivs_plus1;
+      std::fill(temp2.begin(), temp2.end(), 0.0);
+      
+      for (int i = 0; i < uorder; ++i) {
+	int id=i*derivs_plus1;
+	const double *val_p = &values[coefind];
+	const double *co_p  = &co[coefind*kdim] + dim;
+	
+	int temp_ind = 0;
+	for (int wder = 0; wder <= derivs; ++wder) {
+	  for (int vder = 0; vder <= wder; ++vder) {
+	    for (int uder = 0; uder <= vder; ++uder) {
+	      temp2[temp_ind]
+		+= b0[id + wder - vder]*(*val_p);
+	      temp_ind+=kdim;
+	    }
+	  }
+	}
+		
+	for (int d = dim;d < kdim;d++, co_p += kdim) {
+	  int temp_ind = 1;
+	  for (int wder = 0; wder <= derivs; ++wder) {
+	    for (int vder = 0; vder <= wder; ++vder) {
+	      for (int uder = 0; uder <= vder; ++uder) {
+		temp2[temp_ind]
+		  += b0[id + wder - vder]*(*co_p);
+		temp_ind+=ndim;
+	      }
+	    }
+	  }
+	}
+	coefind += 1;
+      }
+      
+      for (int d = 0; d < ndim; ++d) {
+	int temp_ind = d;
+	
+	for (int wder = 0; wder <= derivs; ++wder) {
+	  for (int vder = 0; vder <= wder; ++vder) {
+	    for (int uder = 0; uder <= vder; ++uder) {
+	      temp[temp_ind]
+		+= temp2[temp_ind]*b1[jd + vder - uder];
+	      temp_ind += ndim;
+	      }
+	  }
+	}
+      }
+      coefind += unum - uorder;
+    }
+    
+    for (int d = 0; d < ndim; ++d) {
+      int temp_ind = d;
+      
+      for (int wder = 0; wder <= derivs; ++wder) {
+	for (int vder = 0; vder <= wder; ++vder) {
+	  for (int uder = 0; uder <= vder; ++uder) {
+	    restemp[temp_ind]
+	      += temp[temp_ind]*b2[kd + uder];
+	    temp_ind += ndim;
+	  }
+	}
+      }
+    }
+    coefind += unum * (vnum - vorder);
+  }
+  
+  // Copy from restemp to result
+  if (rational) {
+    vector<double> restemp2(totpts);
+    Go::volume_ratder(&restemp[0], dim, derivs, &restemp2[0]);
+    for (int i = 1; i < totpts; ++i) 
+      gradXi(i) = restemp2[i];
+  } 
+  else {
+    double* restemp_it=restemp.begin();
+    for (int i = 1; i < totpts; ++i) {
+      gradXi(i) = *restemp_it;
+      ++restemp_it;
+    }
+  }  
+  
+  // Derivatives of coordinate mapping 
+  std::vector<Go::Point> xpts(3);
+  vol->point(xpts,fe.u,fe.v,fe.w,1);
+
+  // Jacobian matrix
+  Matrix Jac(3,3);
+  for (size_t i = 1;i <= 3;i++)
+    for (size_t n = 0;n < 3;n++)
+      Jac(n+1,i) = xpts[i][n];
+
+  // Invert Jacobian matrix
+  Jac.inverse();
+
+  // Gradient of solution in given parametric point
+  // df/dX = Jac^-T * df/dXi = [dX/dXi]^-T * df/dXi
+  Jac.multiply(gradXi,grad,true,false);
+
+  return true;
+}
+
+
+bool SplineField3D::gradCoor(const Vec3 x, Vector& grad) const
+{
+  // Not implemented yet
+  return false;
+}

src/ASM/SplineField3D.h

@@ -0,0 +1,77 @@
+//==============================================================================
+//!
+//! \file SplineField3D.h
+//!
+//! \date Mar 28 2011
+//!
+//! \author Runar Holdahl / SINTEF
+//!
+//! \brief Base class for spline based finite element field in 3D
+//!
+//==============================================================================
+
+#ifndef _SPLINE_FIELD_3D_H
+#define _SPLINE_FIELD_3D_H
+
+namespace Go {
+  class SplineVolume;
+  class Point;
+  void volume_ratder(double const eder[],int idim,int ider,double gder[]);
+}
+
+#include "SplineField.h"
+#include "GoTools/trivariate/SplineVolume.h"
+
+/* 
+   \brief Base class for spline-based finite element fields in 3D.
+
+   \details This class implements the functions required to evaluate
+   a 3D spline field at a given point in parametrical or physical 
+   coordinates.
+*/
+
+
+class SplineField3D : public SplineField
+{
+ public:
+  //! \brief The constructor sets the number of space dimensions and fields.
+  //! \param[in] geometry Spline geometry
+  //! \param[in] name Name of spline field
+  SplineField3D(Go::SplineVolume *geometry, char* name = NULL); 
+
+  //! \brief Empty destructor
+  virtual ~SplineField3D();
+
+  // Methods to compute field values
+  //================================================
+
+  //! \brief Computes the value in a given node/control point
+  //! \param[in] node Node number 
+  double valueNode(int 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
+  //! \param[in] x Global/physical coordinate for point
+  double valueCoor(const Vec3 x) const;
+
+  //! \brief Computes the gradient for a given local coordinate
+  //! \param[in] fe Finite element
+  //! \param[out] grad Gradient of solution in a given local coordinate
+  bool gradFE(const FiniteElement& fe, Vector& grad) const;
+
+  //! \brief Computes the gradient for a given global/physical coordinate
+  //! \param[in] x Global coordinate
+  //! \param[out] grad Gradient of solution in a given global coordinate
+  bool gradCoor(const Vec3 x, Vector& grad) const;
+
+ protected:
+  Go::SplineVolume *vol;    //!< Spline volume geometry description
+
+};
+   
+
+
+#endif

src/ASM/SplineFields.h

@@ -0,0 +1,111 @@
+//==============================================================================
+//!
+//! \file SplineFields.h
+//!
+//! \date Mar 28 2011
+//!
+//! \author Runar Holdahl / SINTEF
+//!
+//! \brief Base class for spline based finite element vector field
+//!
+//==============================================================================
+
+#ifndef _SPLINE_FIELDS_H
+#define _SPLINE_FIELDS_H
+
+#include "Vec3.h"
+#include "MatVec.h"
+#include "FiniteElement.h"
+
+/* 
+   \brief Base class for spline-based finite element vector fields.
+
+   \details This class incapsulates the data and methods needed
+   to store and evaluate a spline finite element vector field.
+   This is an abstract base class and the fields associated with 
+   specific spline objects are implemented as subclasses, for
+   instance 1D, 2D and 3D spline formulations.
+*/
+
+
+class SplineFields
+{
+ protected:
+  //! \brief The constructor sets the field name
+  //! \param[in] nsd Number of space dimensions (1, 2 or 3)
+  //! \param[in] name Name of spline field
+ SplineFields(unsigned char nsd_, char* name = NULL) 
+   : nsd(nsd_), fieldname(name) {}
+
+ public:
+  //! \brief Empty destructor
+  virtual ~SplineFields() {}
+
+  // Returns number of space dimensions
+  unsigned char getNoSpaceDim() const { return nsd; }
+
+  // Returns number of fields
+  unsigned char getNoFields() const { return nf; }
+
+  // Returns number of space dimensions
+  int getNoElm() const { return nelm; }
+
+  // Returns number of control points
+  int getNoNodes() const { return nno; }
+
+  // Returns name of spline field
+  const char* getFieldName() const { return fieldname; }
+
+  // Sets the name of the spline field
+  void setFieldName(char* name) { fieldname = name; }
+
+  // Methods to initialize field
+  virtual void fill(Vector& vec) 
+  { values = vec; nf = values.size()/nno; }
+
+  
+  // Methods to compute field values
+  //================================================
+
+  //! \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;
+
+  //! \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
+  virtual bool valueFE(const FiniteElement& fe, Vector& vals) const = 0;
+  
+  //! \brief Computed the value at a given global coordinate
+  //! \param[in] x Global/physical coordinate for point
+  //! \param[in] vals Values in given physical coordinate
+  virtual bool valueCoor(const Vec3 x, Vector& vals) const = 0;
+
+  //! \brief Computes the gradient for a given local coordinate
+  //! \param[in] fe Finite element
+  //! \param[out] grad Gradient of solution in a given local coordinate
+  virtual bool gradFE(const FiniteElement& fe, Matrix& grad) const = 0;
+
+  //! \brief Computes the gradient for a given global/physical coordinate
+  //! \param[in] x Global coordinate
+  //! \param[out] grad Gradient of solution in a given global coordinate
+  virtual bool gradCoor(const Vec3 x, Matrix& grad) const = 0;
+
+ protected:
+  // Dimension of field
+  unsigned char nsd;  //!< Number of space dimensions
+  unsigned char nf;   //!< Number of fields
+  int nelm;           //!< Number of elements/knot-spans
+  int nno;            //!< Number of nodes/control points
+
+  // Fieldname
+  char* fieldname;    //!< Name of spline element field
+
+  // Field values
+  Vector values;      //!< Field values
+};
+   
+
+
+#endif

src/ASM/SplineFields2D.C

@@ -0,0 +1,350 @@
+//==============================================================================
+//!
+//! \file SplineFields2D.C
+//!
+//! \date Mar 28 2011
+//!
+//! \author Runar Holdahl / SINTEF
+//!
+//! \brief Base class for spline based finite element vector field in 2D
+//!
+//==============================================================================
+
+#include "SplineFields2D.h"
+#include <algorithm>
+#include "GoTools/geometry/SplineUtils.h"
+#ifndef __BORLANDC__
+#include "boost/lambda/lambda.hpp"
+#endif
+
+using namespace std;
+#ifndef __BORLANDC__
+using namespace boost::lambda;
+#endif
+
+SplineFields2D::SplineFields2D(Go::SplineSurface *geometry, char* name)
+  : SplineFields(2,name),  surf(geometry) 
+{
+  const int n1 = surf->numCoefs_u();
+  const int n2 = surf->numCoefs_v();
+  nno = n1*n2;
+
+  const int p1 = surf->order_u();
+  const int p2 = surf->order_v();
+  nelm = (n1-p1+1)*(n2-p2+1);
+
+  // Number of fields set in fill
+  nf = 0;
+}
+
+
+SplineFields2D::~SplineFields2D()
+{
+  // Set geometry to NULL; should not be deallocated here
+  surf = NULL;
+} 
+
+
+bool SplineFields2D::valueNode(int node, Vector& vals) const
+{
+  // Not implemented yet
+  return false;
+}
+
+
+// bool SplineFields2D::valueFE(const FiniteElement& fe, Vector& vals) const
+// {
+//   if (!surf) return false;
+
+//   Go::Point pt;
+//   surf->pointValue(pt,values,fe.u,fe.v);
+  
+//   vals.resize(pt.size());
+//   for (int i = 0;i < pt.size();i++)
+//     vals[i] = pt[i];
+
+//   return true;
+// }
+
+
+bool SplineFields2D::valueFE(const FiniteElement& fe, Vector& vals) const
+{
+  if (!surf) return false;
+
+  const int uorder = surf->order_u();
+  const int vorder = surf->order_v();
+  const int unum = surf->numCoefs_u();
+  const int vnum = surf->numCoefs_v();
+
+  const int dim = surf->dimension();
+  const bool rational = surf->rational();
+  int kdim = rational ? dim + 1 : dim;
+
+  const int ncomp = values.size()/(unum*vnum);
+  Go::Point result(ncomp);
+
+  static Go::ScratchVect<double, 10> Bu(uorder);
+  static Go::ScratchVect<double, 10> Bv(vorder);
+  static Go::ScratchVect<double, 4>  tempVal(ncomp);
+  static Go::ScratchVect<double, 4>  tempResult(ncomp);
+
+  Bu.resize(uorder);
+  Bv.resize(vorder);
+  tempVal.resize(ncomp);
+
+  surf->basis_u().computeBasisValues(fe.u,Bu.begin());
+  surf->basis_v().computeBasisValues(fe.v,Bv.begin());
+  
+  const int uleft = surf->basis_u().lastKnotInterval();
+  const int vleft = surf->basis_v().lastKnotInterval();
+
+  //compute the tensor product value
+  const int val_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1)) * ncomp;
+  
+  register double* vtemp;
+  register const double* val_ptr = &values[val_start_ix]; 
+  std::fill(result.begin(), result.end(), double(0));
+
+  if (rational) {
+      double w = 0.0;
+      const int w_start_ix  = (uleft - uorder + 1 + unum * (vleft - vorder + 1)) * kdim + dim;
+      register const double* w_ptr  = &(surf->rcoefs_begin()[w_start_ix]);
+      
+      for (register double* bval_v_ptr = Bv.begin(); bval_v_ptr != Bv.end(); ++bval_v_ptr) {
+        register const double bval_v = *bval_v_ptr;
+	std::fill(tempVal.begin(), tempVal.end(), 0);
+        for (register double* bval_u_ptr = Bu.begin(); bval_u_ptr != Bu.end(); ++bval_u_ptr) {
+          register const double bval_u = *bval_u_ptr;
+          for (vtemp = tempVal.begin(); vtemp != tempVal.end(); ++vtemp) {
+            *vtemp += bval_u * (*val_ptr++);
+          }
+          w += bval_u * (*w_ptr);
+          w_ptr += kdim;
+        }
+        vtemp = tempVal.begin();
+        for (register double* v = result.begin(); v != result.end(); ++v) {
+          *v += (*vtemp++) * bval_v;
+        }
+        val_ptr += ncomp * (unum - uorder);
+        w_ptr   += kdim * (unum - uorder);
+      }
+
+      const double w_inv = double(1) / w;
+#ifdef __BORLANDC__ //C++Builder does not support boost lambda
+      transform(result.begin(), result.end(), result.begin(), ScaleBy(w_inv));
+#else
+      transform(result.begin(), result.end(), result.begin(), _1 * w_inv);
+#endif
+    }
+  else {
+    for (register double* bval_v_ptr = Bv.begin(); bval_v_ptr != Bv.end(); ++bval_v_ptr) {
+      register const double bval_v = *bval_v_ptr;
+      std::fill(tempVal.begin(), tempVal.end(), 0);
+      for (register double* bval_u_ptr = Bu.begin(); bval_u_ptr != Bu.end(); ++bval_u_ptr) {
+	register const double bval_u = *bval_u_ptr;
+	for (vtemp = tempVal.begin(); vtemp != tempVal.end(); ++vtemp) {
+	  *vtemp += bval_u * (*val_ptr++);
+	}
+      }
+      vtemp = tempVal.begin();
+      for (register double* v = result.begin(); v != result.end(); ++v) {
+	*v += (*vtemp++) * bval_v;
+      }
+      val_ptr += ncomp * (unum - uorder);
+    }
+  } 
+
+  vals.resize(result.size());
+  for (int i = 0;i < result.size();i++)
+    vals[i] = result[i];
+
+  return true;
+}
+
+
+
+bool SplineFields2D::valueCoor(const Vec3 x, Vector& vals) const
+{
+  // Not implemented yet
+  return false;
+}
+
+
+// bool SplineFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const
+// {
+//   if (!surf) return false;
+
+//   // Derivatives of solution in parametric space
+//   std::vector<Go::Point> pts(3);
+//   surf->pointValue(pts,values,fe.u,fe.v,1);
+
+//   // Gradient of field wrt parametric coordinates
+//   Matrix gradXi(nf,2);
+//   for (size_t i = 1;i <= 2;i++) 
+//     for (size_t n = 0;n < nf;n++) 
+//       gradXi(n+1,i) = pts[i][n];
+  
+//   // Derivatives of coordinate mapping 
+//   std::vector<Go::Point> xpts(3);
+//   surf->point(xpts,fe.u,fe.v,1);
+
+//   // Jacobian matrix
+//   Matrix Jac(2,2);
+//   for (size_t i = 1;i <= 2;i++)
+//     for (size_t n = 0;n < 2;n++)
+//       Jac(n+1,i) = xpts[i][n];
+  
+//   // Invert Jacobian matrix
+//   Jac.inverse();
+
+//   // Gradient of solution in given parametric point
+//   // df/dX = df/dXi * Jac^-1 = df/dXi * [dX/dXi]^-1
+//   grad.multiply(gradXi,Jac);
+
+//   return true;
+// }
+
+
+bool SplineFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const
+{
+  if (!surf) return false;
+
+  // Gradient of field 
+  Matrix gradXi(nf,2);
+
+  const int uorder = surf->order_u();
+  const int vorder = surf->order_v();
+  const int unum = surf->numCoefs_u();
+  const int vnum = surf->numCoefs_v();
+
+  const int derivs = 1;
+  const bool u_from_right = true;
+  const bool v_from_right = true;
+  const double resolution = 1.0e-12;
+
+  const int totpts = (derivs + 1)*(derivs + 2)/2;
+  const int ncomp = values.size()/(unum*vnum);
+  
+  // Dimension
+  const int dim       = surf->dimension();
+  const bool rational = surf->rational();
+  
+  // Take care of the rational case
+  const std::vector<double>::iterator co = rational ? surf->rcoefs_begin() : surf->coefs_begin();
+  int kdim = dim   + (rational ? 1 : 0);
+  int ndim = ncomp + (rational ? 1 : 0); 
+  
+  // Make temporary storage for the basis values and a temporary
+  // computation cache.
+  Go::ScratchVect<double, 30> b0(uorder*(derivs+1));
+  Go::ScratchVect<double, 30> b1(vorder*(derivs+1));
+  Go::ScratchVect<double, 30> temp(ndim * totpts);
+  Go::ScratchVect<double, 30> restemp(ndim * totpts);
+  std::fill(restemp.begin(), restemp.end(), 0.0);
+
+  // Compute the basis values and get some data about the spline spaces
+  if (u_from_right) {
+    surf->basis_u().computeBasisValues(fe.u, &b0[0], derivs, resolution);
+  } else {
+    surf->basis_u().computeBasisValuesLeft(fe.u, &b0[0], derivs, resolution);
+  }
+  int uleft  = surf->basis_u().lastKnotInterval();
+
+  if (v_from_right) {
+    surf->basis_v().computeBasisValues(fe.v, &b1[0], derivs, resolution);
+  } else {
+    surf->basis_v().computeBasisValuesLeft(fe.v, &b1[0], derivs, resolution);
+  }
+  int vleft = surf->basis_v().lastKnotInterval();
+
+  // Compute the tensor product value
+  int coefind = uleft-uorder+1 + unum*(vleft-vorder+1);
+  int derivs_plus1=derivs+1;
+  for (int jj = 0; jj < vorder; ++jj) {
+    int jjd=jj*(derivs_plus1);
+    std::fill(temp.begin(), temp.end(), 0.0);
+    
+    for (int ii = 0; ii < uorder; ++ii) {
+      int iid=ii*(derivs_plus1);
+      const double *val_p = &values[coefind*ncomp];
+      const double *co_p  = &co[coefind*kdim] + dim;
+      for (int dd = 0; dd < ncomp; ++dd,++val_p) {
+        int temp_ind=dd;
+        for (int vder = 0; vder < derivs_plus1; ++vder) {
+          for (int uder = 0; uder < vder+1; ++uder) {
+            temp[temp_ind] += b0[iid+vder - uder]*(*val_p);
+            temp_ind += ndim;
+          }
+        }
+      }
+      for (int dd = dim;dd < kdim;dd++,co_p += kdim) {
+        int temp_ind = ncomp;
+        for (int vder = 0; vder < derivs_plus1; ++vder) {
+          for (int uder = 0; uder < vder+1; ++uder) {
+            temp[temp_ind] += b0[iid+vder - uder]*(*co_p);
+            temp_ind += ndim;
+          }
+        }
+      }
+
+      coefind += 1;
+    }
+    
+    for (int dd = 0; dd < ndim; ++dd) {
+      int dercount = 0;
+      for (int vder = 0; vder < derivs_plus1; ++vder) {
+        for (int uder = 0; uder < vder + 1; ++uder) {
+          restemp[dercount*ndim + dd] 
+            += temp[dercount*ndim + dd]*b1[uder + jjd];
+          ++dercount;
+        }
+      }
+    }
+    
+    coefind += unum - uorder;
+  }
+
+  // Copy from restemp to result
+  if (rational) {
+    std::vector<double> restemp2(totpts*ncomp);
+    Go::surface_ratder(&restemp[0], ncomp, derivs, &restemp2[0]);
+    for (int i = 1; i < totpts; ++i) {
+      for (int dd = 0; dd < ncomp; ++dd) 
+	gradXi(dd+1,i) = restemp2[i*ncomp + dd];
+    }
+  } else {
+    double* restemp_it=restemp.begin() + ncomp;
+    for (int i = 1; i < totpts; ++i) {
+      for (int dd = 0; dd < ncomp; ++dd) {
+	gradXi(dd+1,i) = *restemp_it;
+        ++restemp_it;
+      }
+    }
+  }  
+
+  // Derivatives of coordinate mapping 
+  std::vector<Go::Point> xpts(3);
+  surf->point(xpts,fe.u,fe.v,1);
+
+  // Jacobian matrix
+  Matrix Jac(2,2);
+  for (size_t i = 1;i <= 2;i++)
+    for (size_t n = 0;n < 2;n++)
+      Jac(n+1,i) = xpts[i][n];
+  
+  // Invert Jacobian matrix
+  Jac.inverse();
+
+  // Gradient of solution in given parametric point
+  // df/dX = df/dXi * Jac^-1 = df/dXi * [dX/dXi]^-1
+  grad.multiply(gradXi,Jac);
+
+  return true;
+}
+
+
+bool SplineFields2D::gradCoor(const Vec3 x, Matrix& grad) const
+{
+  // Not implemented yet
+  return false;
+}

src/ASM/SplineFields2D.h

@@ -0,0 +1,78 @@
+//==============================================================================
+//!
+//! \file SplineFields2D.h
+//!
+//! \date Mar 28 2011
+//!
+//! \author Runar Holdahl / SINTEF
+//!
+//! \brief Base class for spline based finite element vector field in 2D
+//!
+//==============================================================================
+
+#ifndef _SPLINE_FIELDS_2D_H
+#define _SPLINE_FIELDS_2D_H
+
+namespace Go {
+  class SplineSurface;
+}
+
+#include "SplineFields.h"
+#include "GoTools/geometry/SplineSurface.h"
+
+/* 
+   \brief Class for spline-based finite element vector fields in 2D.
+
+   \details This class implements the functions required to evaluate
+   a 2D spline vector field at a given point in parametrical or physical 
+   coordinates.
+
+*/
+
+
+class SplineFields2D : public SplineFields
+{
+ public:
+  //! \brief The constructor sets the field name
+  //! \param[in] geometry Spline volume geometry
+  //! \param[in] name Name of spline field
+  SplineFields2D(Go::SplineSurface *geometry, char* name = NULL); 
+ 
+//! \brief Empty destructor
+  virtual ~SplineFields2D();
+
+  // Methods to compute field values
+  //================================================
+
+  //! \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;
+
+  //! \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
+  //! \param[in] x Global/physical coordinate for point
+  //! \param[in] vals Values in given physical coordinate
+  bool valueCoor(const Vec3 x, Vector& vals) const;
+
+  //! \brief Computes the gradient for a given local coordinate
+  //! \param[in] fe Finite element
+  //! \param[out] grad Gradient of solution in a given local coordinate
+  bool gradFE(const FiniteElement& fe, Matrix& grad) const;
+
+  //! \brief Computes the gradient for a given global/physical coordinate
+  //! \param[in] x Global coordinate
+  //! \param[out] grad Gradient of solution in a given global coordinate
+  bool gradCoor(const Vec3 x, Matrix& grad) const;
+
+ protected:
+  Go::SplineSurface *surf;    //!< Spline surface geometry description
+};
+   
+
+
+#endif

src/ASM/SplineFields3D.C

@@ -0,0 +1,435 @@
+//==============================================================================
+//!
+//! \file SplineFields3D.C
+//!
+//! \date Mar 28 2011
+//!
+//! \author Runar Holdahl / SINTEF
+//!
+//! \brief Base class for spline based finite element vector field in 3D
+//!
+//==============================================================================
+
+
+#include "SplineFields3D.h"
+#include <algorithm>
+#include "GoTools/geometry/SplineUtils.h"
+#ifndef __BORLANDC__
+#include "boost/lambda/lambda.hpp"
+#endif
+
+using namespace std;
+#ifndef __BORLANDC__
+using namespace boost::lambda;
+#endif
+
+SplineFields3D::SplineFields3D(Go::SplineVolume *geometry, char* name)
+  : SplineFields(3,name), vol(geometry) 
+{
+  const int n1 = vol->numCoefs(0);
+  const int n2 = vol->numCoefs(1);
+  const int n3 = vol->numCoefs(2);
+  nno = n1*n2*n3;
+
+  const int p1 = vol->order(0);
+  const int p2 = vol->order(1);
+  const int p3 = vol->order(2);
+  nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1);
+
+  // Number of fields set in fill
+  nf = 0;
+}
+
+
+SplineFields3D::~SplineFields3D()
+{
+  // Set geometry to NULL; should not be deallocated here
+  vol = NULL;
+} 
+
+
+bool SplineFields3D::valueNode(int node, Vector& vals) const
+{
+  // Not implemented yet
+  return false;
+}
+
+
+// bool SplineFields3D::valueFE(const FiniteElement& fe, Vector& vals) const
+// {
+//   if (!vol) return false;
+
+//   Go::Point pt;
+//   vol->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];
+
+//   return true;
+// }
+
+
+bool SplineFields3D::valueFE(const FiniteElement& fe, Vector& vals) const
+{
+  if (!vol) return false;
+
+  Go::Point pt;
+
+  const int uorder = vol->order(0);
+  const int vorder = vol->order(1);
+  const int worder = vol->order(2);
+  const int unum = vol->numCoefs(0);
+  const int vnum = vol->numCoefs(1);
+  const int wnum = vol->numCoefs(2);
+
+  const int dim = vol->dimension();
+  const bool rational = vol->rational();
+  int kdim = rational ? dim + 1 : dim;
+  
+  const int ncomp = values.size()/(unum*vnum*wnum);
+  
+  static Go::ScratchVect<double, 10> Bu(uorder);
+  static Go::ScratchVect<double, 10> Bv(vorder);
+  static Go::ScratchVect<double, 10> Bw(worder);
+  static Go::ScratchVect<double, 4> tempPt(ncomp);
+  static Go::ScratchVect<double, 4> tempPt2(ncomp);
+  static Go::ScratchVect<double, 4> tempResult(ncomp);
+  
+  Bu.resize(uorder);
+  Bv.resize(vorder);
+  Bw.resize(worder);
+  tempPt.resize(ncomp);
+  tempPt2.resize(ncomp);
+  tempResult.resize(ncomp);
+  
+  // compute tbe basis values and get some data about the spline spaces
+  vol->basis(0).computeBasisValues(fe.u, Bu.begin());
+  vol->basis(1).computeBasisValues(fe.v, Bv.begin());
+  vol->basis(2).computeBasisValues(fe.w, Bw.begin());
+  const int uleft = vol->basis(0).lastKnotInterval();
+  const int vleft = vol->basis(1).lastKnotInterval();
+  const int wleft = vol->basis(2).lastKnotInterval();
+  
+  // compute the tensor product value
+  const int val_start_ix =  (uleft - uorder + 1 + unum * (vleft - vorder + 1 + vnum * (wleft - worder + 1))) * ncomp;
+  
+  double* vtemp;
+  const double* val_ptr = &values[val_start_ix];
+  std::fill(pt.begin(), pt.end(), double(0));
+  
+  double w = 1.0;
+  if (rational) {
+    w = 0.0;
+    const int w_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1 + vnum * (wleft - worder + 1))) * kdim + dim;
+    const double *w_ptr = &(vol->rcoefs_begin()[w_start_ix]);
+    
+    for (double* bval_w_ptr = Bw.begin(); bval_w_ptr != Bw.end(); ++bval_w_ptr) {
+      const double bval_w = *bval_w_ptr;
+      std::fill(tempPt.begin(), tempPt.end(), 0);
+      for (double* bval_v_ptr = Bv.begin(); bval_v_ptr != Bv.end(); ++bval_v_ptr) {
+	const double bval_v = *bval_v_ptr;
+	std::fill(tempPt2.begin(), tempPt2.end(), 0);
+	for (double* bval_u_ptr = Bu.begin(); bval_u_ptr != Bu.end(); ++bval_u_ptr) {
+	  const double bval_u = *bval_u_ptr;
+	  for (vtemp = tempPt2.begin(); vtemp != tempPt2.end(); ++vtemp) {
+	    *vtemp += bval_u * (*val_ptr++);
+	  }
+	  w += bval_u*(*w_ptr);
+	  w_ptr += kdim;
+	}
+	vtemp = tempPt2.begin();
+	for (double* v = tempPt.begin(); v != tempPt.end(); ++v) {
+	  *v += (*vtemp++) * bval_v;
+	}
+	val_ptr += ncomp * (unum - uorder);
+	w_ptr   += kdim * (unum - uorder);
+      }
+      vtemp = tempPt.begin();
+      for (double* v = tempResult.begin(); v != tempResult.end(); ++v) {
+	*v += (*vtemp++) * bval_w;
+      }
+      w *= bval_w;
+      val_ptr += ncomp * unum * (vnum - vorder);
+      w_ptr   += kdim  * unum * (vnum - vorder);
+    }
+  }
+  else {
+    for (double* bval_w_ptr = Bw.begin(); bval_w_ptr != Bw.end(); ++bval_w_ptr) {
+      const double bval_w = *bval_w_ptr;
+      std::fill(tempPt.begin(), tempPt.end(), 0);
+      for (double* bval_v_ptr = Bv.begin(); bval_v_ptr != Bv.end(); ++bval_v_ptr) {
+	const double bval_v = *bval_v_ptr;
+	std::fill(tempPt2.begin(), tempPt2.end(), 0);
+	for (double* bval_u_ptr = Bu.begin(); bval_u_ptr != Bu.end(); ++bval_u_ptr) {
+	  const double bval_u = *bval_u_ptr;
+	  for (vtemp = tempPt2.begin(); vtemp != tempPt2.end(); ++vtemp) {
+	    *vtemp += bval_u * (*val_ptr++);
+	  }
+	}
+	vtemp = tempPt2.begin();
+	for (double* v = tempPt.begin(); v != tempPt.end(); ++v) {
+	  *v += (*vtemp++) * bval_v;
+	}
+	val_ptr += ncomp * (unum - uorder);
+      }
+      vtemp = tempPt.begin();
+      for (double* v = tempResult.begin(); v != tempResult.end(); ++v) {
+	*v += (*vtemp++) * bval_w;
+      }
+      val_ptr += ncomp * unum * (vnum - vorder);
+    }
+  }
+  
+  copy(tempResult.begin(), tempResult.begin() + dim, pt.begin());
+  if (rational) {
+    const double w_inv = double(1) / w;
+    transform(pt.begin(), pt.end(), pt.begin(), _1 * w_inv);
+  }    
+  
+  vals.resize(pt.size());
+  for (int i = 0;i < pt.size();i++)
+    vals[i] = pt[i];
+
+  return true;
+}
+
+
+
+bool SplineFields3D::valueCoor(const Vec3 x, Vector& vals) const
+{
+  // Not implemented yet
+  return false;
+}
+
+
+// bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const
+// {
+//   if (!vol) return false;
+  
+//   // Derivatives of solution in parametric space
+//   std::vector<Go::Point> pts;
+//   vol->pointValue(pts,values,fe.u,fe.v,fe.w,1);
+
+//  // Gradient of field wrt parametric coordinates 
+//   Matrix gradXi(nf,3);
+//   for (size_t i = 1;i <= 3;i++) 
+//     for (size_t n = 0;n < nf;n++) 
+//       gradXi(n+1,i) = pts[i][n];
+
+//   // Derivatives of coordinate mapping 
+//   std::vector<Go::Point> xpts(3);
+//   vol->point(xpts,fe.u,fe.v,fe.w,1);
+
+//   // Jacobian matrix
+//   Matrix Jac(3,3);
+//   for (size_t i = 1;i <= 3;i++)
+//     for (size_t n = 0;n < 3;n++)
+//       Jac(n+1,i) = xpts[i][n];
+//   Jac.inverse();
+  
+//   // Gradient of solution in given parametric point
+//   // df/dX = df/dXi * Jac^-1 = df/dXi * [dX/dXi]^-1
+//   grad.multiply(gradXi,Jac);
+
+//   return true;
+// }
+
+
+bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const
+{
+  if (!vol) return false;
+
+  // Derivatives of solution in parametric space
+  std::vector<Go::Point> pts(4);
+
+  int derivs = 1;
+  int totpts = (derivs + 1)*(derivs + 2)*(derivs + 3)/6;
+  
+  const int unum = vol->numCoefs(0);
+  const int vnum = vol->numCoefs(1);
+  const int wnum = vol->numCoefs(2);
+
+  const int uorder = vol->order(0);
+  const int vorder = vol->order(1);
+  const int worder = vol->order(2);
+
+  const bool u_from_right = true;
+  const bool v_from_right = true;
+  const bool w_from_right = true;
+
+  const double resolution = 1.0e-12;
+  
+  const int ncomp     = values.size()/(unum*vnum*wnum);
+  const int dim       = vol->dimension();
+  const bool rational = vol->rational();
+  
+  for (int i = 0; i < totpts; ++i) {
+    if (pts[i].dimension() != ncomp) {
+      pts[i].resize(ncomp);
+    }
+  }
+  
+  // Take care of the rational case
+  const vector<double>::iterator co = rational ? vol->rcoefs_begin() : vol->coefs_begin();
+  int kdim = dim + (rational ? 1 : 0);
+  int ndim = ncomp + (rational ? 1 : 0);
+  
+  // Make temporary storage for the basis values and a temporary
+  // computation cache.
+  Go::ScratchVect<double, 30> b0(uorder * (derivs+1));
+  Go::ScratchVect<double, 30> b1(vorder * (derivs+1));
+  Go::ScratchVect<double, 30> b2(worder * (derivs+1));
+  Go::ScratchVect<double, 60> temp(ndim * totpts);
+  Go::ScratchVect<double, 60> temp2(ndim * totpts);
+  Go::ScratchVect<double, 60> restemp(ndim * totpts);
+  std::fill(restemp.begin(), restemp.end(), 0.0);
+  
+  // Compute the basis values and get some data about the spline spaces
+  if (u_from_right) 
+    vol->basis(0).computeBasisValues(fe.u, &b0[0], derivs, resolution);
+  else 
+    vol->basis(1).computeBasisValuesLeft(fe.u, &b0[0], derivs, resolution);
+  
+  int uleft = vol->basis(0).lastKnotInterval();
+  if (v_from_right) 
+    vol->basis(1).computeBasisValues(fe.v, &b1[0], derivs, resolution);
+  else 
+    vol->basis(1).computeBasisValuesLeft(fe.v, &b1[0], derivs, resolution);
+  
+  int vleft = vol->basis(1).lastKnotInterval();
+  if (w_from_right) 
+    vol->basis(2).computeBasisValues(fe.w, &b2[0], derivs, resolution);
+  else 
+    vol->basis(2).computeBasisValuesLeft(fe.w, &b2[0], derivs, resolution);
+  
+  int wleft = vol->basis(2).lastKnotInterval();
+  
+  // Compute the tensor product value
+  int coefind = uleft-uorder+1 + unum*(vleft-vorder+1 + vnum*(wleft-worder+1));
+  int derivs_plus1=derivs+1;
+  
+  for (int k = 0; k < worder; ++k) {
+    int kd=k*derivs_plus1;
+    std::fill(temp.begin(), temp.end(), 0.0);
+      
+    for (int j = 0; j < vorder; ++j) {
+      int jd=j*derivs_plus1;
+      std::fill(temp2.begin(), temp2.end(), 0.0);
+      
+      for (int i = 0; i < uorder; ++i) {
+	int id=i*derivs_plus1;
+	const double *val_p = &values[coefind*ncomp]; 
+	const double *co_p  = &co[coefind*kdim] + dim;
+	
+	for (int d = 0; d < ncomp; ++d,++val_p) {
+	  int temp_ind = d;
+	  
+	  for (int wder = 0; wder <= derivs; ++wder) {
+	    for (int vder = 0; vder <= wder; ++vder) {
+	      for (int uder = 0; uder <= vder; ++uder) {
+		temp2[temp_ind]
+		  += b0[id + wder - vder]*(*val_p);
+		temp_ind+=kdim;
+	      }
+	    }
+	  }
+	}
+	
+	for (int d = dim;d < kdim;d++, co_p += kdim) {
+	  int temp_ind = ncomp;
+	  
+	  for (int wder = 0; wder <= derivs; ++wder) {
+	    for (int vder = 0; vder <= wder; ++vder) {
+	      for (int uder = 0; uder <= vder; ++uder) {
+		temp2[temp_ind]
+		  += b0[id + wder - vder]*(*co_p);
+		temp_ind+=ndim;
+	      }
+	    }
+	  }
+	}
+	coefind += 1;
+      }
+      
+      for (int d = 0; d < ndim; ++d) {
+	int temp_ind = d;
+	
+	for (int wder = 0; wder <= derivs; ++wder) {
+	  for (int vder = 0; vder <= wder; ++vder) {
+	    for (int uder = 0; uder <= vder; ++uder) {
+	      temp[temp_ind]
+		+= temp2[temp_ind]*b1[jd + vder - uder];
+	      temp_ind += ndim;
+	      }
+	  }
+	}
+      }
+      coefind += unum - uorder;
+    }
+    
+    for (int d = 0; d < ndim; ++d) {
+      int temp_ind = d;
+      
+      for (int wder = 0; wder <= derivs; ++wder) {
+	for (int vder = 0; vder <= wder; ++vder) {
+	  for (int uder = 0; uder <= vder; ++uder) {
+	    restemp[temp_ind]
+	      += temp[temp_ind]*b2[kd + uder];
+	    temp_ind += ndim;
+	  }
+	}
+      }
+    }
+    coefind += unum * (vnum - vorder);
+  }
+  
+  // Copy from restemp to result
+  if (rational) {
+    vector<double> restemp2(totpts*ncomp);
+    Go::volume_ratder(&restemp[0], dim, derivs, &restemp2[0]);
+    for (int i = 0; i < totpts; ++i) {
+      for (int d = 0; d < ncomp; ++d) {
+	pts[i][d] = restemp2[i*ncomp + d];
+      }
+    }
+  } else {
+    double* restemp_it=restemp.begin();
+    for (int i = 0; i < totpts; ++i) {
+      for (int d = 0; d < ncomp; ++d) {
+	pts[i][d] = *restemp_it;
+	++restemp_it;
+      }
+    }
+  }  
+
+ // Gradient of field wrt parametric coordinates 
+  Matrix gradXi(nf,3);
+  for (size_t i = 1;i <= 3;i++) 
+    for (size_t n = 0;n < nf;n++) 
+      gradXi(n+1,i) = pts[i][n];
+
+  // Derivatives of coordinate mapping 
+  std::vector<Go::Point> xpts(3);
+  vol->point(xpts,fe.u,fe.v,fe.w,1);
+
+  // Jacobian matrix
+  Matrix Jac(3,3);
+  for (size_t i = 1;i <= 3;i++)
+    for (size_t n = 0;n < 3;n++)
+      Jac(n+1,i) = xpts[i][n];
+  Jac.inverse();
+  
+  // Gradient of solution in given parametric point
+  // df/dX = df/dXi * Jac^-1 = df/dXi * [dX/dXi]^-1
+  grad.multiply(gradXi,Jac);
+
+  return true;
+}
+
+
+bool SplineFields3D::gradCoor(const Vec3 x, Matrix& grad) const
+{
+  // Not implemented yet
+  return false;
+}

src/ASM/SplineFields3D.h

@@ -0,0 +1,80 @@
+//==============================================================================
+//!
+//! \file SplineFields3D.h
+//!
+//! \date Mar 28 2011
+//!
+//! \author Runar Holdahl / SINTEF
+//!
+//! \brief Base class for spline based finite element vector field in 3D
+//!
+//==============================================================================
+
+#ifndef _SPLINE_FIELDS_3D_H
+#define _SPLINE_FIELDS_3D_H
+
+namespace Go {
+  class SplineVolume;
+  class Point;
+  void volume_ratder(double const eder[],int idim,int ider,double gder[]);
+}
+
+#include "SplineFields.h"
+#include "GoTools/trivariate/SplineVolume.h"
+
+/* 
+   \brief Class for spline-based finite element vector fields in 3D.
+
+   \details This class implements the functions required to evaluate
+   a 3D spline vector field at a given point in parametrical or physical 
+   coordinates.
+
+*/
+
+
+class SplineFields3D : public SplineFields
+{
+ public:
+  //! \brief The constructor sets the field name
+  //! \param[in] geometry Spline volume geometry
+  //! \param[in] name Name of spline field
+  SplineFields3D(Go::SplineVolume *geometry, char* name = NULL); 
+ 
+//! \brief Empty destructor
+  virtual ~SplineFields3D();
+
+  // Methods to compute field values
+  //================================================
+
+  //! \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;
+
+  //! \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
+  //! \param[in] x Global/physical coordinate for point
+  //! \param[in] vals Values in given physical coordinate
+  bool valueCoor(const Vec3 x, Vector& vals) const;
+
+  //! \brief Computes the gradient for a given local coordinate
+  //! \param[in] fe Finite element
+  //! \param[out] grad Gradient of solution in a given local coordinate
+  bool gradFE(const FiniteElement& fe, Matrix& grad) const;
+
+  //! \brief Computes the gradient for a given global/physical coordinate
+  //! \param[in] x Global coordinate
+  //! \param[out] grad Gradient of solution in a given global coordinate
+  bool gradCoor(const Vec3 x, Matrix& grad) const;
+
+ protected:
+  Go::SplineVolume *vol;    //!< Spline volume geometry description
+};
+   
+
+
+#endif