Made getElementCoordinates and getNodalCoordinates public. Added function getNoElmDOF. Added some functionality for computation of stabilization parameters

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

src/ASM/ASMbase.h

@@ -114,6 +114,8 @@ public:
   unsigned char getNoParamDim() const { return ndim; }
   //! \brief Returns the number of solution fields.
   virtual unsigned char getNoFields(int b = 0) const { return b > 1 ? 0 : nf; }
+  //! \brief Returns the number of solution fields.
+  virtual size_t getNoElmDOF(int b = 0) const { return b > 1 ? 0 : neldof; }
 
   //! \brief Returns local 1-based index of the node with given global number.
   //! \details If the given node number is not present, 0 is returned.
@@ -136,6 +138,12 @@ public:
   //! \param[out] X nsd\f$\times\f$n-matrix, where \a n is the number of nodes
   //! in the patch
   virtual void getNodalCoordinates(Matrix& X) const = 0;
+  //! \brief Returns a matrix with nodal coordinates for an element.
+  //! \param[in] iel Element index
+  //! \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 = 0;
+
   //! \brief Prints out the nodal coordinates of this patch to the given stream.
   void printNodes(std::ostream& os, const char* heading = 0) const;
 
@@ -440,10 +448,11 @@ public:
 
 protected:
   // Standard finite element data structures
-  unsigned char ndim; //!< Number of parametric dimensions (1, 2 or 3)
+  unsigned char ndim;   //!< Number of parametric dimensions (1, 2 or 3)
-  unsigned char nsd;  //!< Number of space dimensions (ndim <= nsd <= 3)
+  unsigned char nsd;    //!< Number of space dimensions (ndim <= nsd <= 3)
-  unsigned char nf;   //!< Number of primary solution fields (1 or larger)
+  unsigned char nf;     //!< Number of primary solution fields (1 or larger)
-
+  size_t        neldof; //!< Number of degrees of freedom per element
+  
   const IntVec& MLGE; //!< Matrix of Local to Global Element numbers
   const IntVec& MLGN; //!< Matrix of Local to Global Node numbers
   const IntMat& MNPC; //!< Matrix of Nodal Point Correspondance

src/ASM/ASMmxBase.C

@@ -14,8 +14,9 @@
 #include "ASMmxBase.h"
 
 
-bool ASMmxBase::geoUsesBasis1 = false;
+bool ASMmxBase::geoUsesBasis1     = false;
-bool ASMmxBase::useCpminus1   = false;
+bool ASMmxBase::useCpminus1       = false;
+bool ASMmxBase::useLowOrderBasis1 = false;
 
 
 ASMmxBase::ASMmxBase (unsigned char n_f1, unsigned char n_f2)
@@ -25,6 +26,18 @@ ASMmxBase::ASMmxBase (unsigned char n_f1, unsigned char n_f2)
 }
 
 
+size_t ASMmxBase::getNoElmDOF(int basis) const
+{
+  switch(basis) {
+  case 1: return neldof1;
+  case 2: return neldof2;
+  }    
+  
+ 
+ return neldof1 + neldof2;
+}
+
+
 void ASMmxBase::initMx (const std::vector<int>& MLGN, const int* sysMadof)
 {
   MADOF.resize(MLGN.size());

src/ASM/ASMmxBase.h

@@ -49,8 +49,12 @@ protected:
 		     const std::vector<int>& nodes) const;
 
 public:
-  static bool geoUsesBasis1; //!< If \e true, 1st basis represents the geometry
+  static bool geoUsesBasis1;       //!< If \e true, 1st basis represents the geometry
-  static bool useCpminus1;   //!< If \e true, enforce \f$C^{p-1}\f$ continuity
+  static bool useCpminus1;         //!< If \e true, enforce \f$C^{p-1}\f$ continuity
+  static bool useLowOrderBasis1;   //!< If \e true, basis 1 is of lowest order
+
+  //! \brief Returns the number of solution fields.
+  virtual size_t getNoElmDOF(int b = 0) const;
 
 private:
   std::vector<int> MADOF; //!< Matrix of accumulated DOFs for this patch
@@ -61,6 +65,9 @@ protected:
 
   unsigned char nf1; //!< Number of solution fields using first basis
   unsigned char nf2; //!< Number of solution fields using second basis
+
+  size_t neldof1; //!< Number of degrees of freedom per element for basis 1
+  size_t neldof2; //!< Number of degrees of freedom per element for basis 2
 };
 
 #endif

src/ASM/ASMs1D.h

@@ -59,6 +59,17 @@ public:
   //! \param[in] inod 1-based node index local to current patch
   virtual Vec3 getCoord(size_t inod) const;
 
+  //! \brief Returns a matrix with nodal coordinates for an element.
+  //! \param[in] iel Element index
+  //! \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;
+
+  //! \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
+  virtual void getNodalCoordinates(Matrix& X) const;
+
   //! \brief Updates the nodal coordinates for this patch.
   //! \param[in] displ Incremental displacements to update the coordinates with
   virtual bool updateCoords(const Vector& displ);
@@ -223,16 +234,6 @@ protected:
   //! \param[in] iel Element index
   double getParametricLength(int iel) const;
 
-  //! \brief Returns a matrix with nodal coordinates for an element.
-  //! \param[in] iel Element index
-  //! \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;
-  //! \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
-  virtual void getNodalCoordinates(Matrix& X) const;
-
   //! \brief Returns the number of nodal points in the patch.
   //! \param[in] basis Which basis to return size parameters for (mixed methods)
   virtual int getSize(int basis = 0) const;

src/ASM/ASMs1DLag.h

@@ -47,6 +47,17 @@ public:
   //! This is used to reinitialize the patch after it has been refined.
   virtual void clear(bool retainGeometry = false);
 
+  //! \brief Returns a matrix with nodal coordinates for an element.
+  //! \param[in] iel Element index
+  //! \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;
+
+  //! \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
+  virtual void getNodalCoordinates(Matrix& X) const;
+
   //! \brief Returns the global coordinates for the given node.
   //! \param[in] inod 1-based node index local to current patch
   virtual Vec3 getCoord(size_t inod) const;
@@ -131,16 +142,6 @@ protected:
   // Internal utility methods
   // ========================
 
-  //! \brief Returns a matrix with nodal coordinates for an element.
-  //! \param[in] iel Element index
-  //! \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;
-  //! \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
-  virtual void getNodalCoordinates(Matrix& X) const;
-
   //! \brief Returns the number of nodal points in the patch.
   virtual int getSize(int = 0) const { return nx; }
 

src/ASM/ASMs2D.C

@@ -86,6 +86,9 @@ bool ASMs2D::read (std::istream& is)
     nsd = surf->dimension();
   }
 
+  // Define number of DOFs pr element
+  neldof = surf->order_u()*surf->order_v()*nf;
+
   geo = surf;
   return true;
 }
@@ -668,6 +671,7 @@ double ASMs2D::getParametricArea (int iel) const
 
   double du = surf->knotSpan(0,nodeInd[inod1].I);
   double dv = surf->knotSpan(1,nodeInd[inod1].J);
+
   return du*dv;
 }
 
@@ -889,6 +893,25 @@ static double getElmSize (int p1, int p2, const Matrix& X)
 }
 
 
+/*!
+  \brief Computes the characteristic element length from inverse Jacobian.
+*/
+
+static bool getG (const Matrix& Jinv, const Vector& du, Matrix& G)
+{
+  const size_t nsd = Jinv.cols();
+
+  G.resize(nsd,nsd,true);
+
+  for (size_t k = 1;k <= nsd;k++) 
+    for (size_t l = 1;l <= nsd;l++)
+      for (size_t m = 1;m <= nsd;m++)
+	G(k,l) += Jinv(m,k)*Jinv(m,l)/(du(k)*du(l));
+
+  return true;
+}
+
+
 void ASMs2D::extractBasis (const Go::BasisDerivsSf& spline,
 			   Vector& N, Matrix& dNdu)
 {
@@ -994,7 +1017,6 @@ bool ASMs2D::integrate (Integrand& integrand,
   Matrix3D d2Ndu2, Hess;
   Vec4     X;
 
-
   // === Assembly loop over all elements in the patch ==========================
 
   int iel = 1;
@@ -1012,9 +1034,8 @@ bool ASMs2D::integrate (Integrand& integrand,
       if (!this->getElementCoordinates(Xnod,iel)) return false;
 
       // Compute characteristic element length, if needed
-      if (integrand.getIntegrandType() == 2)
+      if (integrand.getIntegrandType() == 2) 
-	fe.h = getElmSize(p1,p2,Xnod);
+	fe.h   = getElmSize(p1,p2,Xnod);
-
       else if (integrand.getIntegrandType() == 3)
       {
 	// --- Compute average value of basis functions over the element -------
@@ -1062,7 +1083,6 @@ bool ASMs2D::integrate (Integrand& integrand,
       // If the array is shorter than this, expect a segmentation fault.
       LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
 
-
       if (integrand.getIntegrandType() > 10)
       {
 	// --- Selective reduced integration loop ------------------------------
@@ -1121,10 +1141,18 @@ bool ASMs2D::integrate (Integrand& integrand,
 	  if (fe.detJxW == 0.0) continue; // skip singular points
 
 	  // Compute Hessian of coordinate mapping and 2nd order derivatives
-	  if (integrand.getIntegrandType() == 2)
+	  if (integrand.getIntegrandType() == 2) {
 	    if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu))
 	      return false;
 
+	    // RUNAR
+	    Vector dXidu(2);
+	    dXidu(1) = this->getParametricLength(iel,1);
+	    dXidu(2) = this->getParametricLength(iel,2);
+	    if (!getG(Jac,dXidu,fe.G))
+	      return false;
+	  }
+
 #if SP_DEBUG > 4
 	  std::cout <<"\niel, ip = "<< iel <<" "<< ip
 		    <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX << std::endl;
@@ -1146,7 +1174,7 @@ bool ASMs2D::integrate (Integrand& integrand,
 
       // Assembly of global system integral
       if (!glInt.assemble(elmInt,fe.iel))
-	return false;
+      	return false;
     }
 
   return true;

src/ASM/ASMs2D.h

@@ -100,6 +100,17 @@ public:
   //! This is used to reinitialize the patch after it has been refined.
   virtual void clear(bool retainGeometry = false);
 
+    //! \brief Returns a matrix with nodal coordinates for an element.
+  //! \param[in] iel Element index
+  //! \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;
+
+  //! \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
+  virtual void getNodalCoordinates(Matrix& X) const;
+
   //! \brief Returns the global coordinates for the given node.
   //! \param[in] inod 1-based node index local to current patch
   virtual Vec3 getCoord(size_t inod) const;
@@ -331,15 +342,6 @@ protected:
   //! \param[in] iel Element index
   //! \param[in] dir Local index of the boundary edge
   double getParametricLength(int iel, int dir) const;
-  //! \brief Returns a matrix with nodal coordinates for an element.
-  //! \param[in] iel Element index
-  //! \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;
-  //! \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
-  virtual void getNodalCoordinates(Matrix& X) const;
 
   //! \brief Returns the number of nodal points in each parameter direction.
   //! \param[out] n1 Number of nodes in first (u) direction

src/ASM/ASMs2DLag.h

@@ -47,6 +47,17 @@ public:
   //! This is used to reinitialize the patch after it has been refined.
   virtual void clear(bool retainGeometry = false);
 
+  //! \brief Returns a matrix with nodal coordinates for an element.
+  //! \param[in] iel Element index
+  //! \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;
+
+  //! \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
+  virtual void getNodalCoordinates(Matrix& X) const;
+
   //! \brief Returns the global coordinates for the given node.
   //! \param[in] inod 1-based node index local to current patch
   virtual Vec3 getCoord(size_t inod) const;
@@ -130,23 +141,6 @@ public:
   virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
 			    const RealArray* gpar, bool regular = true) const;
 
-protected:
-
-  // Internal utility methods
-  // ========================
-
-  //! \brief Returns a matrix with nodal coordinates for an element.
-  //! \param[in] iel Element index
-  //! \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
-  virtual void getNodalCoordinates(Matrix& X) const;
-
   //! \brief Returns the number of nodal points in each parameter direction.
   //! \param[out] n1 Number of nodes in first (u) direction
   //! \param[out] n2 Number of nodes in second (v) direction

src/ASM/ASMs2Dmx.C

@@ -189,7 +189,7 @@ bool ASMs2Dmx::generateFEMTopology ()
 							     ug,vg,XYZ,ndim,
 							     false,XYZ);
     }
-    else
+    else 
     {
       // Order-elevate basis1 such that it is of one degree higher than basis2
       // but only C^p-2 continuous
@@ -197,6 +197,12 @@ bool ASMs2Dmx::generateFEMTopology ()
       basis1->raiseOrder(1,1);
     }
     basis2 = surf;
+    
+    if (useLowOrderBasis1) {
+      basis2 = basis1;
+      basis1 = surf;
+      surf = basis2;
+    }
 
     // Define which basis that should be used to represent the geometry
     if (geoUsesBasis1) surf = basis1;
@@ -249,7 +255,7 @@ bool ASMs2Dmx::generateFEMTopology ()
   if (p1 > n1 || p2 > n2) return false;
   if (q1 > m1 || q2 > m2) return false;
 
-  myMLGE.resize((n1-p1+1)*(n2-p2+1),0);
+  myMLGE.resize(geoUsesBasis1 ? (n1-p1+1)*(n2-p2+1) : (m1-q1+1)*(m2-q2+1),0);
   myMLGN.resize(nb1 + nb2);
   myMNPC.resize(myMLGE.size());
   myNodeInd.resize(myMLGN.size());
@@ -340,7 +346,7 @@ bool ASMs2Dmx::generateFEMTopology ()
     iel = inod = 0;
     for (i2 = 1; i2 <= n2; i2++)
       for (i1 = 1; i1 <= n1; i1++, inod++)
-	if (i1 >= p1 && i2 >= p2)
+	if (i1 >= p1 && i2 >= p2) {
 	  if (basis1->knotSpan(0,i1-1) > 0.0)
 	    if (basis1->knotSpan(1,i2-1) > 0.0)
 	    {
@@ -352,9 +358,15 @@ bool ASMs2Dmx::generateFEMTopology ()
 		  myMNPC[iel][lnod++] = inod - n1*j2 - j1;
 
 	      iel++;
-	    }
+	    }	
+	}
   }
 
+  // Number of DOFs pr element for each basis
+  neldof1 = basis1->order_u()*basis1->order_v()*nf1;
+  neldof2 = basis2->order_u()*basis2->order_v()*nf2;
+  neldof = neldof1 + neldof2;
+
 #ifdef SP_DEBUG
   std::cout <<"NEL = "<< MLGE.size() <<" NNOD = "<< MLGN.size() << std::endl;
 #endif
@@ -499,7 +511,6 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
   Matrix dN1du, dN2du, Xnod, Jac;
   Vec4   X;
 
-
   // === Assembly loop over all elements in the patch ==========================
 
   int iel = 1;

src/ASM/ASMs2Dmx.h

@@ -60,6 +60,12 @@ public:
   //! This is used to reinitialize the patch after it has been refined.
   virtual void clear(bool retainGeometry = false);
 
+  //! \brief Returns a matrix with nodal coordinates for an element.
+  //! \param[in] iel Element index
+  //! \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;
+
   //! \brief Returns the global coordinates for the given node.
   //! \param[in] inod 1-based node index local to current patch
   virtual Vec3 getCoord(size_t inod) const;
@@ -175,12 +181,6 @@ protected:
   // Internal utility methods
   // ========================
 
-  //! \brief Returns a matrix with nodal coordinates for an element.
-  //! \param[in] iel Element index
-  //! \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;
-
   //! \brief Returns the number of nodal points in each parameter direction.
   //! \param[out] n1 Number of nodes in first (u) direction
   //! \param[out] n2 Number of nodes in second (v) direction

src/ASM/ASMs3D.C

@@ -75,6 +75,9 @@ bool ASMs3D::read (std::istream& is)
     return false;
   }
 
+  // Number of DOFs pr element
+  neldof = svol->order(0)*svol->order(1)*svol->order(2)*nf;
+
   geo = svol;
   return true;
 }
@@ -1097,6 +1100,25 @@ static double getElmSize (int p1, int p2, int p3, const Matrix& X)
 }
 
 
+/*!
+  \brief Computes the characteristic element length from inverse Jacobian.
+*/
+
+static bool getG (const Matrix& Jinv, const Vector& du, Matrix& G)
+{
+  const size_t nsd = Jinv.cols();
+
+  G.resize(nsd,nsd,true);
+
+  for (size_t k = 1;k <= nsd;k++) 
+    for (size_t l = 1;l <= nsd;l++)
+      for (size_t m = 1;m <= nsd;m++)
+	G(k,l) += Jinv(m,k)*Jinv(m,l)/(du(k)*du(l));
+
+  return true;
+}
+
+
 void ASMs3D::extractBasis (const Go::BasisDerivs& spline,
 			   Vector& N, Matrix& dNdu)
 {
@@ -1335,10 +1357,19 @@ bool ASMs3D::integrate (Integrand& integrand,
 	      if (fe.detJxW == 0.0) continue; // skip singular points
 
  	      // Compute Hessian of coordinate mapping and 2nd order derivatives
-	      if (integrand.getIntegrandType() == 2)
+	      if (integrand.getIntegrandType() == 2) {
 		if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu))
 		  return false;
 
+		// RUNAR
+		// Vector dXidu(3);
+// 		dXidu(1) = this->getParametricLength(iel,1);
+// 		dXidu(2) = this->getParametricLength(iel,2);
+// 		dXidu(3) = this->getParametricLength(iel,3);
+// 		if (!getG(Jac,dXidu,fe.G))
+// 		  return false;
+	      }
+
 #if SP_DEBUG > 4
 	      std::cout <<"\niel, ip = "<< iel <<" "<< ip
 			<<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX << std::endl;

src/ASM/ASMs3D.h

@@ -118,6 +118,17 @@ public:
   //! This is used to reinitialize the patch after it has been refined.
   virtual void clear(bool retainGeometry = false);
 
+  //! \brief Returns a matrix with nodal coordinates for an element.
+  //! \param[in] iel Element index
+  //! \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;
+
+  //! \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
+  virtual void getNodalCoordinates(Matrix& X) const;
+
   //! \brief Returns the global coordinates for the given node.
   //! \param[in] inod 1-based node index local to current patch
   virtual Vec3 getCoord(size_t inod) const;
@@ -390,16 +401,6 @@ protected:
   //! \param[in] dir Local face index of the boundary face
   double getParametricArea(int iel, int dir) const;
 
-  //! \brief Returns a matrix with nodal coordinates for an element.
-  //! \param[in] iel Element index
-  //! \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;
-  //! \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
-  virtual void getNodalCoordinates(Matrix& X) const;
-
   //! \brief Returns the number of nodal points in each parameter direction.
   //! \param[out] n1 Number of nodes in first (u) direction
   //! \param[out] n2 Number of nodes in second (v) direction

src/ASM/ASMs3DLag.h

@@ -138,18 +138,12 @@ public:
   virtual bool evalSolution(Matrix& sField, const Integrand& integrand,
 			    const RealArray* gpar, bool regular = true) const;
 
-protected:
-
-  // Internal utility methods
-  // ========================
-
   //! \brief Returns a matrix with nodal coordinates for an element.
   //! \param[in] iel Element index
   //! \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

src/ASM/ASMs3Dmx.C

@@ -153,6 +153,17 @@ bool ASMs3Dmx::generateFEMTopology ()
       Go::BsplineBasis b1 = svol->basis(0).extendedBasis(svol->order(0)+1);
       Go::BsplineBasis b2 = svol->basis(1).extendedBasis(svol->order(1)+1);
       Go::BsplineBasis b3 = svol->basis(2).extendedBasis(svol->order(2)+1);
+      /*  To lower order and regularity this can be used instead     
+      std::vector<double>::const_iterator first =  ++svol->basis(0).begin();
+      std::vector<double>::const_iterator last  =  --svol->basis(0).end();
+      Go::BsplineBasis b1 = Go::BsplineBasis(svol->order(0)-1,first,last);
+      first =  ++svol->basis(1).begin();
+      last  =  --svol->basis(1).end();
+      Go::BsplineBasis b2 = Go::BsplineBasis(svol->order(1)-1,first,last);
+      first =  ++svol->basis(2).begin();
+      last  =  --svol->basis(2).end();
+      Go::BsplineBasis b3 = Go::BsplineBasis(svol->order(2)-1,first,last);
+      */
 
       // Note: Currently this is implemented for non-rational splines only.
       // TODO: Ask the splines people how to fix this properly, that is, how
@@ -367,6 +378,11 @@ bool ASMs3Dmx::generateFEMTopology ()
 		}
   }
 
+  // Number of DOFs pr element for each basis
+  neldof1 = basis1->order(0)*basis1->order(1)*basis1->order(2)*nf1;
+  neldof2 = basis2->order(0)*basis2->order(1)*basis2->order(2)*nf2;
+  neldof = neldof1+neldof2;
+
 #ifdef SP_DEBUG
   std::cout <<"NEL = "<< MLGE.size() <<" NNOD = "<< MLGN.size() << std::endl;
 #endif

src/ASM/ASMs3Dmx.h

@@ -55,6 +55,12 @@ public:
   //! \param[in] retainGeometry If \e true, the spline geometry is not cleared.
   //! This is used to reinitialize the patch after it has been refined.
   virtual void clear(bool retainGeometry = false);
+  
+  //! \brief Returns a matrix with nodal coordinates for an element.
+  //! \param[in] iel Element index
+  //! \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;
 
   //! \brief Returns the global coordinates for the given node.
   //! \param[in] inod 1-based node index local to current patch
@@ -174,12 +180,6 @@ protected:
   // Internal utility methods
   // ========================
 
-  //! \brief Returns a matrix with nodal coordinates for an element.
-  //! \param[in] iel Element index
-  //! \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;
-
   //! \brief Returns the number of nodal points in each parameter direction.
   //! \param[out] n1 Number of nodes in first (u) direction
   //! \param[out] n2 Number of nodes in second (v) direction