From ac6c1c3a6b864080932aa20b215f5dbf6fea09de Mon Sep 17 00:00:00 2001
From: Knut Morten Okstad <knut.morten.okstad@sintef.no>
Date: Sat, 10 Feb 2018 17:16:22 +0100
Subject: [PATCH] Fixed: Account for that the number of space dimensions in a
 patch can be less than the number of coordinates on file

---
 src/ASM/ASMs1D.C           |  2 +-
 src/ASM/ASMs2D.C           |  2 +-
 src/ASM/ASMs3D.C           |  2 +-
 src/ASM/SplineField2D.C    | 51 ++++++-------------------------
 src/ASM/SplineField2D.h    |  3 ++
 src/ASM/SplineField3D.C    | 53 ++++++++-------------------------
 src/ASM/SplineField3D.h    |  3 ++
 src/ASM/SplineFields2D.C   | 61 +++++++++++---------------------------
 src/ASM/SplineFields2D.h   | 11 ++++---
 src/ASM/SplineFields2Dmx.C | 28 ++++++++---------
 src/ASM/SplineFields3D.C   | 57 ++++++++---------------------------
 src/ASM/SplineFields3D.h   | 11 ++++---
 src/ASM/SplineFields3Dmx.C | 30 +++++++++----------
 13 files changed, 102 insertions(+), 212 deletions(-)

diff --git a/src/ASM/ASMs1D.C b/src/ASM/ASMs1D.C
index 6b1beafe..def15650 100644
--- a/src/ASM/ASMs1D.C
+++ b/src/ASM/ASMs1D.C
@@ -686,7 +686,7 @@ double ASMs1D::getElementEnds (int i, Vec3Vec& XC) const
   XC.reserve(elmCS.empty() ? 2 : 3);
   const double* pt = &XYZ.front();
   for (int j = 0; j < 2; j++, pt += dim)
-    XC.push_back(Vec3(pt,dim));
+    XC.push_back(Vec3(pt,nsd));
 
   // Calculate the characteristic element size
   double h = getElementSize(XC);
diff --git a/src/ASM/ASMs2D.C b/src/ASM/ASMs2D.C
index 6ecd0034..7cc8cfb5 100644
--- a/src/ASM/ASMs2D.C
+++ b/src/ASM/ASMs2D.C
@@ -1456,7 +1456,7 @@ double ASMs2D::getElementCorners (int i1, int i2, Vec3Vec& XC) const
   XC.reserve(4);
   const double* pt = &XYZ.front();
   for (int i = 0; i < 4; i++, pt += dim)
-    XC.push_back(Vec3(pt,dim));
+    XC.push_back(Vec3(pt,nsd));
 
   return getElementSize(XC);
 }
diff --git a/src/ASM/ASMs3D.C b/src/ASM/ASMs3D.C
index 8a27746c..0c8998c3 100644
--- a/src/ASM/ASMs3D.C
+++ b/src/ASM/ASMs3D.C
@@ -1659,7 +1659,7 @@ double ASMs3D::getElementCorners (int i1, int i2, int i3, Vec3Vec& XC) const
   XC.reserve(8);
   const double* pt = &XYZ.front();
   for (int i = 0; i < 8; i++, pt += dim)
-    XC.push_back(Vec3(pt,dim));
+    XC.push_back(Vec3(pt,nsd));
 
   return getElementSize(XC);
 }
diff --git a/src/ASM/SplineField2D.C b/src/ASM/SplineField2D.C
index ba65679d..ab36d08c 100644
--- a/src/ASM/SplineField2D.C
+++ b/src/ASM/SplineField2D.C
@@ -35,6 +35,8 @@ SplineField2D::SplineField2D (const ASMs2D* patch,
   const int p2 = basis->order_v();
   nelm = (n1-p1+1)*(n2-p2+1);
 
+  nsd = patch->getNoSpaceDim();
+
   // Ensure the values array has compatible length, pad with zeros if necessary
   size_t ofs = 0;
   for (char i = 1; i < nbasis; ++i)
@@ -93,7 +95,7 @@ double SplineField2D::valueCoor (const Vec4& x) const
     pt[1] = x[1];
     pt[2] = x[2];
     double clo_u, clo_v, dist;
-  #pragma omp critical
+#pragma omp critical
     surf->closestPoint(pt, clo_u, clo_v, clopt, dist, 1e-10);
     fe.u = clo_u;
     fe.v = clo_v;
@@ -184,9 +186,9 @@ bool SplineField2D::gradFE (const FiniteElement& fe, Vector& grad) const
 		     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);
+  Matrix Xnod(nsd,ip.size()), Jac;
+  for (size_t i = 0; i < ip.size(); i++)
+    Xnod.fillColumn(1+i,&(*surf->coefs_begin())+surf->dimension()*ip[i]);
   utl::Jacobian(Jac,dNdX,Xnod,dNdu);
 
   // Evaluate the gradient of the solution field at the given point
@@ -222,73 +224,38 @@ bool SplineField2D::hessianFE (const FiniteElement& fe, Matrix& H) const
   if (!basis) return false;
   if (!surf)  return false;
 
-  // Order of basis
-  const int uorder = surf->order_u();
-  const int vorder = surf->order_v();
-  const size_t nen = uorder*vorder;
-
-  // Evaluate the basis functions at the given point
-  Go::BasisDerivsSf  spline;
   Go::BasisDerivsSf2 spline2;
   Matrix3D d2Ndu2;
-  Matrix dNdu, dNdX;
   IntVec ip;
   if (surf == basis) {
 #pragma omp critical
     surf->computeBasis(fe.u,fe.v,spline2);
 
-    dNdu.resize(nen,2);
+    const size_t nen = surf->order_u()*surf->order_v();
     d2Ndu2.resize(nen,2,2);
     for (size_t n = 1; n <= nen; n++) {
-      dNdu(n,1) = spline2.basisDerivs_u[n-1];
-      dNdu(n,2) = spline2.basisDerivs_v[n-1];
       d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1];
       d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1];
       d2Ndu2(n,2,2) = spline2.basisDerivs_vv[n-1];
     }
 
     ASMs2D::scatterInd(surf->numCoefs_u(),surf->numCoefs_v(),
-		       uorder,vorder,spline2.left_idx,ip);
+		       surf->order_u(),surf->order_v(),
+		       spline2.left_idx,ip);
   }
   else {
-#pragma omp critical
-    surf->computeBasis(fe.u,fe.v,spline);
-
-    dNdu.resize(nen,2);
-    for (size_t n = 1; n <= nen; n++) {
-      dNdu(n,1) = spline.basisDerivs_u[n-1];
-      dNdu(n,2) = spline.basisDerivs_v[n-1];
-    }
-
-    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
 #pragma omp critical
     basis->computeBasis(fe.u,fe.v,spline2);
 
     const size_t nbf = basis->order_u()*basis->order_v();
-    dNdu.resize(nbf,2);
     d2Ndu2.resize(nbf,2,2);
     for (size_t n = 1; n <= nbf; n++) {
-      dNdu(n,1) = spline2.basisDerivs_u[n-1];
-      dNdu(n,2) = spline2.basisDerivs_v[n-1];
       d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1];
       d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1];
       d2Ndu2(n,2,2) = spline2.basisDerivs_vv[n-1];
     }
 
-    ip.clear();
     ASMs2D::scatterInd(basis->numCoefs_u(),basis->numCoefs_v(),
 		       basis->order_u(),basis->order_v(),
 		       spline2.left_idx,ip);
diff --git a/src/ASM/SplineField2D.h b/src/ASM/SplineField2D.h
index 5f0ce192..4875412d 100644
--- a/src/ASM/SplineField2D.h
+++ b/src/ASM/SplineField2D.h
@@ -77,6 +77,9 @@ public:
 protected:
   const Go::SplineSurface* basis; //!< Spline basis description
   const Go::SplineSurface* surf;  //!< Spline geometry description
+
+private:
+  unsigned char nsd; //!< Number of space dimensions
 };
 
 #endif
diff --git a/src/ASM/SplineField3D.C b/src/ASM/SplineField3D.C
index 3ab7a3d1..a5435d5e 100644
--- a/src/ASM/SplineField3D.C
+++ b/src/ASM/SplineField3D.C
@@ -37,6 +37,8 @@ SplineField3D::SplineField3D (const ASMs3D* patch,
   const int p3 = basis->order(2);
   nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1);
 
+  nsd = patch->getNoSpaceDim();
+
   size_t ofs = 0;
   for (char i = 1; i < nbasis; ++i)
     ofs += patch->getNoNodes(i)*patch->getNoFields(i);
@@ -95,7 +97,7 @@ double SplineField3D::valueCoor (const Vec4& x) const
     pt[1] = x[1];
     pt[2] = x[2];
     double clo_u, clo_v, clo_w, dist;
-  #pragma omp critical
+#pragma omp critical
     vol->closestPoint(pt, clo_u, clo_v, clo_w, clopt, dist, 1e-5);
 
     fe.u = clo_u;
@@ -191,9 +193,9 @@ bool SplineField3D::gradFE (const FiniteElement& fe, Vector& grad) const
 		     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);
+  Matrix Xnod(nsd,ip.size()), Jac;
+  for (size_t i = 0; i < ip.size(); i++)
+    Xnod.fillColumn(1+i,&(*vol->coefs_begin())+vol->dimension()*ip[i]);
   utl::Jacobian(Jac,dNdX,Xnod,dNdu);
 
   // Evaluate the gradient of the solution field at the given point
@@ -230,25 +232,16 @@ bool SplineField3D::hessianFE (const FiniteElement& fe, Matrix& H) const
   if (!basis) return false;
   if (!vol)  return false;
 
-  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;
-
-  // Evaluate the basis functions at the given point
-  Go::BasisDerivs  spline;
   Go::BasisDerivs2 spline2;
   Matrix3D d2Ndu2;
-  Matrix dNdu(nen,3), dNdX;
   IntVec ip;
   if (vol == basis) {
 #pragma omp critical
     vol->computeBasis(fe.u,fe.v,fe.w,spline2);
+
+    const size_t nen = vol->order(0)*vol->order(1)*vol->order(2);
     d2Ndu2.resize(nen,3,3);
     for (size_t n = 1; n <= nen; n++) {
-      dNdu(n,1) = spline2.basisDerivs_u[n-1];
-      dNdu(n,2) = spline2.basisDerivs_v[n-1];
-      dNdu(n,3) = spline2.basisDerivs_w[n-1];
       d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1];
       d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1];
       d2Ndu2(n,1,3) = d2Ndu2(n,3,1) = spline2.basisDerivs_uw[n-1];
@@ -256,40 +249,19 @@ bool SplineField3D::hessianFE (const FiniteElement& fe, Matrix& H) const
       d2Ndu2(n,2,3) = d2Ndu2(n,3,2) = spline2.basisDerivs_vw[n-1];
       d2Ndu2(n,3,3) = spline2.basisDerivs_ww[n-1];
     }
+
     ASMs3D::scatterInd(vol->numCoefs(0),vol->numCoefs(1),vol->numCoefs(2),
-		       uorder,vorder,worder,spline2.left_idx,ip);
+		       vol->order(0),vol->order(1),vol->order(2),
+		       spline2.left_idx,ip);
   }
   else {
-#pragma omp critical
-    vol->computeBasis(fe.u,fe.v,fe.w,spline);
-    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];
-    }
-    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
 #pragma omp critical
     basis->computeBasis(fe.u,fe.v,fe.w,spline2);
 
     const size_t nbf = basis->order(0)*basis->order(1)*basis->order(2);
-    dNdu.resize(nbf,3);
     d2Ndu2.resize(nbf,3,3);
     for (size_t n = 1; n <= nbf; n++) {
-      dNdu(n,1) = spline2.basisDerivs_u[n-1];
-      dNdu(n,2) = spline2.basisDerivs_v[n-1];
-      dNdu(n,3) = spline2.basisDerivs_w[n-1];
       d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1];
       d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1];
       d2Ndu2(n,1,3) = d2Ndu2(n,3,1) = spline2.basisDerivs_uw[n-1];
@@ -298,10 +270,9 @@ bool SplineField3D::hessianFE (const FiniteElement& fe, Matrix& H) const
       d2Ndu2(n,3,3) = spline2.basisDerivs_ww[n-1];
     }
 
-    ip.clear();
     ASMs3D::scatterInd(basis->numCoefs(0),basis->numCoefs(1),basis->numCoefs(2),
 		       basis->order(0),basis->order(1),basis->order(2),
-		       spline.left_idx,ip);
+		       spline2.left_idx,ip);
   }
 
   Vector Vnod;
diff --git a/src/ASM/SplineField3D.h b/src/ASM/SplineField3D.h
index 963ce497..844f9ec1 100644
--- a/src/ASM/SplineField3D.h
+++ b/src/ASM/SplineField3D.h
@@ -77,6 +77,9 @@ public:
 protected:
   const Go::SplineVolume* basis; //!< Spline basis description
   const Go::SplineVolume* vol;   //!< Spline geometry description
+
+private:
+  unsigned char nsd; //!< Number of space dimensions
 };
 
 #endif
diff --git a/src/ASM/SplineFields2D.C b/src/ASM/SplineFields2D.C
index 1f8d40f0..919357e9 100644
--- a/src/ASM/SplineFields2D.C
+++ b/src/ASM/SplineFields2D.C
@@ -11,6 +11,8 @@
 //!
 //==============================================================================
 
+#include "GoTools/geometry/SplineSurface.h"
+
 #include "SplineFields2D.h"
 #include "ASMs2D.h"
 #include "FiniteElement.h"
@@ -18,8 +20,6 @@
 #include "Utilities.h"
 #include "Vec3.h"
 
-#include "GoTools/geometry/SplineSurface.h"
-
 
 SplineFields2D::SplineFields2D (const ASMs2D* patch,
                                 const RealArray& v, char nbasis,
@@ -34,6 +34,8 @@ SplineFields2D::SplineFields2D (const ASMs2D* patch,
   const int p2 = basis->order_v();
   nelm = (n1-p1+1)*(n2-p2+1);
 
+  nsd = patch->getNoSpaceDim();
+
   size_t ofs = 0;
   for (char i = 1; i < nbasis; ++i)
     ofs += patch->getNoNodes(i)*patch->getNoFields(i);
@@ -61,8 +63,9 @@ SplineFields2D::SplineFields2D (const Go::SplineSurface* srf,
                                 const RealArray& v, int cmp, const char* name)
   : Fields(name), basis(srf), surf(srf)
 {
-  values = v;
+  nsd = 2; // That is, this constructor can not be used for shells (nsd=3)
   nf = cmp;
+  values = v;
 }
 
 
@@ -113,7 +116,7 @@ bool SplineFields2D::valueCoor (const Vec4& x, Vector& vals) const
     pt[1] = x[1];
     pt[2] = x[2];
     double clo_u, clo_v, dist;
-  #pragma omp critical
+#pragma omp critical
     surf->closestPoint(pt, clo_u, clo_v, clopt, dist, 1e-5);
 
     fe.u = clo_u;
@@ -150,9 +153,9 @@ bool SplineFields2D::gradFE (const FiniteElement& fe, Matrix& grad) const
 		     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);
+  Matrix Xnod(nsd,ip.size()), Jac;
+  for (size_t i = 0; i < ip.size(); i++)
+    Xnod.fillColumn(1+i,&(*surf->coefs_begin())+surf->dimension()*ip[i]);
   utl::Jacobian(Jac,dNdX,Xnod,dNdu);
 
   // Evaluate the gradient of the solution field at the given point
@@ -189,72 +192,42 @@ bool SplineFields2D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
   if (!basis) return false;
   if (!surf)  return false;
 
-  // Order of basis
-  const int uorder = surf->order_u();
-  const int vorder = surf->order_v();
-  const size_t nen = uorder*vorder;
-
-  // Evaluate the basis functions at the given point
-  Go::BasisDerivsSf  spline;
   Go::BasisDerivsSf2 spline2;
   Matrix3D d2Ndu2;
-  Matrix dNdu(nen,2), dNdX;
   IntVec ip;
+
   if (surf == basis) {
 #pragma omp critical
     surf->computeBasis(fe.u,fe.v,spline2);
+
+    const size_t nen = surf->order_u()*surf->order_v();
     d2Ndu2.resize(nen,2,2);
     for (size_t n = 1; n <= nen; n++) {
-      dNdu(n,1) = spline2.basisDerivs_u[n-1];
-      dNdu(n,2) = spline2.basisDerivs_v[n-1];
       d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1];
       d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1];
       d2Ndu2(n,2,2) = spline2.basisDerivs_vv[n-1];
     }
 
     ASMs2D::scatterInd(surf->numCoefs_u(),surf->numCoefs_v(),
-		       uorder,vorder,spline2.left_idx,ip);
+                       surf->order_u(),surf->order_v(),
+                       spline2.left_idx,ip);
   }
   else {
-#pragma omp critical
-    surf->computeBasis(fe.u,fe.v,spline);
-    for (size_t n = 1; n <= nen; n++) {
-      dNdu(n,1) = spline.basisDerivs_u[n-1];
-      dNdu(n,2) = spline.basisDerivs_v[n-1];
-    }
-
-    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
 #pragma omp critical
     basis->computeBasis(fe.u,fe.v,spline2);
 
     const size_t nbf = basis->order_u()*basis->order_v();
-    dNdu.resize(nbf,2);
     d2Ndu2.resize(nbf,2,2);
     for (size_t n = 1; n <= nbf; n++) {
-      dNdu(n,1) = spline2.basisDerivs_u[n-1];
-      dNdu(n,2) = spline2.basisDerivs_v[n-1];
       d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1];
       d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1];
       d2Ndu2(n,2,2) = spline2.basisDerivs_vv[n-1];
     }
 
-    ip.clear();
     ASMs2D::scatterInd(basis->numCoefs_u(),basis->numCoefs_v(),
-		       basis->order_u(),basis->order_v(),
-		       spline2.left_idx,ip);
+                       basis->order_u(),basis->order_v(),
+                       spline2.left_idx,ip);
   }
 
   Matrix Vnod;
diff --git a/src/ASM/SplineFields2D.h b/src/ASM/SplineFields2D.h
index 2e592e7c..6f3ddc3e 100644
--- a/src/ASM/SplineFields2D.h
+++ b/src/ASM/SplineFields2D.h
@@ -59,22 +59,22 @@ public:
   //! \brief Computes the value in a given node/control point.
   //! \param[in] node Node number
   //! \param[out] vals Node values
-  bool valueNode(size_t node, Vector& vals) const;
+  virtual 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;
+  virtual bool valueFE(const FiniteElement& fe, Vector& vals) const;
 
   //! \brief Computes the value at a given global coordinate.
   //! \param[in] x Global/physical coordinate for point
   //! \param[out] vals Values in given physical coordinate
-  bool valueCoor(const Vec4& x, Vector& vals) const;
+  virtual bool valueCoor(const Vec4& 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;
+  virtual bool gradFE(const FiniteElement& fe, Matrix& grad) const;
 
   //! \brief Computes the hessian for a given local coordinate.
   //! \param[in] fe Finite element quantities
@@ -84,6 +84,9 @@ public:
 protected:
   const Go::SplineSurface* basis; //!< Spline basis description
   const Go::SplineSurface* surf;  //!< Spline geometry description
+
+private:
+  unsigned char nsd; //!< Number of space dimensions
 };
 
 #endif
diff --git a/src/ASM/SplineFields2Dmx.C b/src/ASM/SplineFields2Dmx.C
index 20c0e145..02162b82 100644
--- a/src/ASM/SplineFields2Dmx.C
+++ b/src/ASM/SplineFields2Dmx.C
@@ -32,8 +32,8 @@ SplineFields2Dmx::SplineFields2Dmx (const ASMs2Dmx* patch,
   for (int i = 1; i < *bases.begin(); ++i)
     ofs += patch->getNoNodes(i)*patch->getNoFields(i);
   vit += ofs;
-  for (const auto& it : bases) {
-    size_t nno = patch->getNoNodes(it)*patch->getNoFields(it);
+  for (int b : bases) {
+    size_t nno = patch->getNoNodes(b)*patch->getNoFields(b);
     RealArray::const_iterator end = v.size() > nno+ofs ? vit+nno : v.end();
     std::copy(vit,end,std::back_inserter(values));
     vit += nno;
@@ -58,8 +58,8 @@ bool SplineFields2Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
   // Evaluate the basis functions at the given point
   auto vit = values.begin();
   auto rit = vals.begin();
-  for (const auto& it : bases) {
-    Go::SplineSurface* basis = surf->getBasis(it);
+  for (int b : bases) {
+    Go::SplineSurface* basis = surf->getBasis(b);
     Go::BasisPtsSf spline;
 #pragma omp critical
     basis->computeBasis(fe.u,fe.v,spline);
@@ -71,11 +71,11 @@ bool SplineFields2Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
                        spline.left_idx,ip);
 
     Matrix Vnod;
-    utl::gather(ip,1,Vector(&*vit,surf->getNoNodes(it)),Vnod);
+    utl::gather(ip,1,Vector(&*vit,surf->getNoNodes(b)),Vnod);
     Vector val2;
     Vnod.multiply(spline.basisValues,val2); // vals = Vnod * basisValues
     *rit++ = val2.front();
-    vit += surf->getNoNodes(it);
+    vit += surf->getNoNodes(b);
   }
 
   return true;
@@ -108,16 +108,16 @@ bool SplineFields2Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
 		     uorder,vorder,spline.left_idx,ip);
 
   // Evaluate the Jacobian inverse
-  Matrix Xnod, Jac;
-  Vector Xctrl(&(*gsurf->coefs_begin()),gsurf->coefs_end()-gsurf->coefs_begin());
-  utl::gather(ip,gsurf->dimension(),Xctrl,Xnod);
+  Matrix Xnod(surf->getNoSpaceDim(),ip.size()), Jac;
+  for (size_t i = 0; i < ip.size(); i++)
+    Xnod.fillColumn(1+i,&(*gsurf->coefs_begin())+gsurf->dimension()*ip[i]);
   utl::Jacobian(Jac,dNdX,Xnod,dNdu);
 
   // Evaluate the gradient of the solution field at the given point
   auto vit = values.begin();
   size_t row = 1;
-  for (const auto& it : bases) {
-    const Go::SplineSurface* basis = surf->getBasis(it);
+  for (int b : bases) {
+    const Go::SplineSurface* basis = surf->getBasis(b);
 #pragma omp critical
     basis->computeBasis(fe.u,fe.v,spline);
 
@@ -135,13 +135,13 @@ bool SplineFields2Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
 		       basis->order_u(),basis->order_v(),
 		       spline.left_idx,ip);
 
-    utl::gather(ip,1,Vector(&*vit, surf->getNoNodes(it)*
-                                   surf->getNoFields(it)),Xnod);
+    const size_t nval = surf->getNoNodes(b)*surf->getNoFields(b);
+    utl::gather(ip,1,Vector(&*vit,nval),Xnod);
     Matrix grad2;
     grad2.multiply(Xnod,dNdX); // grad = Xnod * dNdX
     grad(row,1) = grad2(1,1);
     grad(row++,2) = grad2(1,2);
-    vit += surf->getNoNodes(it)*surf->getNoFields(it);
+    vit += nval;
   }
 
   return true;
diff --git a/src/ASM/SplineFields3D.C b/src/ASM/SplineFields3D.C
index 63c3389f..db35af32 100644
--- a/src/ASM/SplineFields3D.C
+++ b/src/ASM/SplineFields3D.C
@@ -36,6 +36,8 @@ SplineFields3D::SplineFields3D (const ASMs3D* patch,
   const int p3 = basis->order(2);
   nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1);
 
+  nsd = patch->getNoSpaceDim();
+
   size_t ofs = 0;
   for (char i = 1; i < nbasis; ++i)
     ofs += patch->getNoNodes(i)*patch->getNoFields(i);
@@ -62,8 +64,9 @@ SplineFields3D::SplineFields3D (const Go::SplineVolume* svol,
                                 const RealArray& v, int cmp, const char* name)
   : Fields(name), basis(svol), vol(svol)
 {
-  values = v;
+  nsd = 3;
   nf = cmp;
+  values = v;
 }
 
 
@@ -115,7 +118,7 @@ bool SplineFields3D::valueCoor (const Vec4& x, Vector& vals) const
     pt[1] = x[1];
     pt[2] = x[2];
     double clo_u, clo_v, clo_w, dist;
-  #pragma omp critical
+#pragma omp critical
     vol->closestPoint(pt, clo_u, clo_v, clo_w, clopt, dist, 1e-5);
 
     fe.u = clo_u;
@@ -155,9 +158,9 @@ bool SplineFields3D::gradFE (const FiniteElement& fe, Matrix& grad) const
 		     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);
+  Matrix Xnod(nsd,ip.size()), Jac;
+  for (size_t i = 0; i < ip.size(); i++)
+    Xnod.fillColumn(1+i,&(*vol->coefs_begin())+vol->dimension()*ip[i]);
   utl::Jacobian(Jac,dNdX,Xnod,dNdu);
 
   // Evaluate the gradient of the solution field at the given point
@@ -195,27 +198,17 @@ bool SplineFields3D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
   if (!basis) return false;
   if (!vol)  return false;
 
-  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;
-
   // Evaluate the basis functions at the given point
-  Go::BasisDerivs  spline;
   Go::BasisDerivs2 spline2;
   Matrix3D d2Ndu2;
-  Matrix dNdu, dNdX;
   IntVec ip;
   if (vol == basis) {
 #pragma omp critical
     vol->computeBasis(fe.u,fe.v,fe.w,spline2);
 
-    dNdu.resize(nen,3);
+    const size_t nen = vol->order(0)*vol->order(1)*vol->order(2);
     d2Ndu2.resize(nen,3,3);
     for (size_t n = 1; n <= nen; n++) {
-      dNdu(n,1) = spline2.basisDerivs_u[n-1];
-      dNdu(n,2) = spline2.basisDerivs_v[n-1];
-      dNdu(n,3) = spline2.basisDerivs_w[n-1];
       d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1];
       d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1];
       d2Ndu2(n,1,3) = d2Ndu2(n,3,1) = spline2.basisDerivs_uw[n-1];
@@ -225,42 +218,17 @@ bool SplineFields3D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
     }
 
     ASMs3D::scatterInd(vol->numCoefs(0),vol->numCoefs(1),vol->numCoefs(2),
-		       uorder,vorder,worder,spline2.left_idx,ip);
+		       vol->order(0),vol->order(1),vol->order(2),
+		       spline2.left_idx,ip);
   }
   else {
-#pragma omp critical
-    vol->computeBasis(fe.u,fe.v,fe.w,spline);
-
-    dNdu.resize(nen,3);
-    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];
-    }
-
-    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
 #pragma omp critical
     basis->computeBasis(fe.u,fe.v,fe.w,spline2);
 
     const size_t nbf = basis->order(0)*basis->order(1)*basis->order(2);
-    dNdu.resize(nbf,3);
     d2Ndu2.resize(nbf,3,3);
     for (size_t n = 1; n <= nbf; n++) {
-      dNdu(n,1) = spline2.basisDerivs_u[n-1];
-      dNdu(n,2) = spline2.basisDerivs_v[n-1];
-      dNdu(n,3) = spline2.basisDerivs_w[n-1];
       d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1];
       d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1];
       d2Ndu2(n,1,3) = d2Ndu2(n,3,1) = spline2.basisDerivs_uw[n-1];
@@ -269,10 +237,9 @@ bool SplineFields3D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
       d2Ndu2(n,3,3) = spline2.basisDerivs_ww[n-1];
     }
 
-    ip.clear();
     ASMs3D::scatterInd(basis->numCoefs(0),basis->numCoefs(1),basis->numCoefs(2),
 		       basis->order(0),basis->order(1),basis->order(2),
-		       spline.left_idx,ip);
+		       spline2.left_idx,ip);
   }
 
   Matrix Vnod;
diff --git a/src/ASM/SplineFields3D.h b/src/ASM/SplineFields3D.h
index f7a4890b..9f5b82fb 100644
--- a/src/ASM/SplineFields3D.h
+++ b/src/ASM/SplineFields3D.h
@@ -59,22 +59,22 @@ public:
   //! \brief Computes the value in a given node/control point.
   //! \param[in] node Node number
   //! \param[out] vals Node values
-  bool valueNode(size_t node, Vector& vals) const;
+  virtual 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;
+  virtual bool valueFE(const FiniteElement& fe, Vector& vals) const;
 
   //! \brief Computes the value at a given global coordinate.
   //! \param[in] x Global/physical coordinate for point
   //! \param[out] vals Values in given physical coordinate
-  bool valueCoor(const Vec4& x, Vector& vals) const;
+  virtual bool valueCoor(const Vec4& 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;
+  virtual bool gradFE(const FiniteElement& fe, Matrix& grad) const;
 
   //! \brief Computes the hessian for a given local coordinate.
   //! \param[in] fe Finite element quantities
@@ -84,6 +84,9 @@ public:
 protected:
   const Go::SplineVolume* basis; //!< Spline basis description
   const Go::SplineVolume* vol;   //!< Spline geometry description
+
+private:
+  unsigned char nsd; //!< Number of space dimensions
 };
 
 #endif
diff --git a/src/ASM/SplineFields3Dmx.C b/src/ASM/SplineFields3Dmx.C
index a611de54..44d2f3c4 100644
--- a/src/ASM/SplineFields3Dmx.C
+++ b/src/ASM/SplineFields3Dmx.C
@@ -32,8 +32,8 @@ SplineFields3Dmx::SplineFields3Dmx (const ASMs3Dmx* patch,
   for (int i = 1; i < *bases.begin(); ++i)
     ofs += patch->getNoNodes(i)*patch->getNoFields(i);
   vit += ofs;
-  for (const auto& it : bases) {
-    size_t nno = patch->getNoFields(it)*patch->getNoNodes(it);
+  for (int b : bases) {
+    size_t nno = patch->getNoFields(b)*patch->getNoNodes(b);
     RealArray::const_iterator end = v.size() > nno+ofs ? vit+nno : v.end();
     std::copy(vit,end,std::back_inserter(values));
     vit += nno;
@@ -58,8 +58,8 @@ bool SplineFields3Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
   // Evaluate the basis functions at the given point
   auto vit = values.begin();
   auto rit = vals.begin();
-  for (const auto& it : bases) {
-    Go::SplineVolume* basis = svol->getBasis(it);
+  for (int b : bases) {
+    Go::SplineVolume* basis = svol->getBasis(b);
     Go::BasisPts spline;
 #pragma omp critical
     basis->computeBasis(fe.u,fe.v,fe.w,spline);
@@ -72,11 +72,11 @@ bool SplineFields3Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
                        spline.left_idx,ip);
 
     Matrix Vnod;
-    utl::gather(ip,1,Vector(&*vit,svol->getNoNodes(it)),Vnod);
+    utl::gather(ip,1,Vector(&*vit,svol->getNoNodes(b)),Vnod);
     Vector val2;
     Vnod.multiply(spline.basisValues,val2); // vals = Vnod * basisValues
     *rit++ = val2.front();
-    vit += svol->getNoNodes(it);
+    vit += svol->getNoNodes(b);
   }
 
   return true;
@@ -111,16 +111,16 @@ bool SplineFields3Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
 		     uorder,vorder,worder,spline.left_idx,ip);
 
   // Evaluate the Jacobian inverse
-  Matrix Xnod, Jac;
-  Vector Xctrl(&(*gvol->coefs_begin()),gvol->coefs_end()-gvol->coefs_begin());
-  utl::gather(ip,gvol->dimension(),Xctrl,Xnod);
-  utl::Jacobian(Jac,dNdX,Xnod,dNdu);
+  Matrix Xnod(svol->getNoSpaceDim(),ip.size()), Jac;
+  for (size_t i = 0; i < ip.size(); i++)
+    Xnod.fillColumn(1+i,&(*gvol->coefs_begin())+gvol->dimension()*ip[i]);
+   utl::Jacobian(Jac,dNdX,Xnod,dNdu);
 
   // Evaluate the gradient of the solution field at the given point
   auto vit = values.begin();
   size_t row = 1;
-  for (const auto& it : bases) {
-    const Go::SplineVolume* basis = svol->getBasis(it);
+  for (int b : bases) {
+    const Go::SplineVolume* basis = svol->getBasis(b);
 #pragma omp critical
     basis->computeBasis(fe.u,fe.v,fe.w,spline);
 
@@ -139,14 +139,14 @@ bool SplineFields3Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
 		       basis->order(0),basis->order(1),basis->order(2),
 		       spline.left_idx,ip);
 
-    utl::gather(ip,1,Vector(&*vit, svol->getNoNodes(it)*
-                                   svol->getNoFields(it)),Xnod);
+    const size_t nval = svol->getNoNodes(b)*svol->getNoFields(b);
+    utl::gather(ip,1,Vector(&*vit,nval),Xnod);
     Matrix grad2;
     grad2.multiply(Xnod,dNdX); // grad = Xnod * dNdX
     grad(row,1) = grad2(1,1);
     grad(row,2) = grad(1,2);
     grad(row++,3) = grad(1,3);
-    vit += svol->getNoNodes(it)*svol->getNoFields(it);
+    vit += nval;
   }
 
   return true;