diff --git a/Geo/MElement.h b/Geo/MElement.h
index 1b56081c8cc0272dc219b1f8fd6b3285bdf82ba4..2910bd893fc543ab4e410929f8441b400c8e482d 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -223,7 +223,7 @@ class MElement
   virtual const nodalBasis* getFunctionSpace(int order=-1) const { return 0; }
 
   // get the function space for the jacobian of the element
-  virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const { return 0; }
+  virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const { return 0; }
 
   // return parametric coordinates (u,v,w) of a vertex
   virtual void getNode(int num, double &u, double &v, double &w);
diff --git a/Geo/MSubElement.cpp b/Geo/MSubElement.cpp
index 3d8b99cad7058de18a03537d31371b78eec317de..d83eb01ff9fab937f43732908d4fef911338191d 100644
--- a/Geo/MSubElement.cpp
+++ b/Geo/MSubElement.cpp
@@ -12,8 +12,8 @@
 
 MSubTetrahedron::~MSubTetrahedron()
 {
-  if(_owner)
-    delete _orig;
+//   if(_owner)
+//     delete _orig;
 }
 
 const nodalBasis* MSubTetrahedron::getFunctionSpace(int order) const
@@ -28,19 +28,39 @@ const JacobianBasis* MSubTetrahedron::getJacobianFuncSpace(int order) const
   return 0;
 }
 
-void MSubTetrahedron::getShapeFunctions(double u, double v, double w, double s[], int o)
+void MSubTetrahedron::getShapeFunctions(double u, double v, double w, double s[], int order)
 {
-  if(_orig) _orig->getShapeFunctions(u, v, w, s, o);
+  if(_orig) _orig->getShapeFunctions(u, v, w, s, order);
 }
 
-void MSubTetrahedron::getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
+void MSubTetrahedron::getGradShapeFunctions(double u, double v, double w, double s[][3], int order)
 {
-  if(_orig) _orig->getGradShapeFunctions(u, v, w, s, o);
+  if(_orig) _orig->getGradShapeFunctions(u, v, w, s, order);
 }
 
-void MSubTetrahedron::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o)
+void MSubTetrahedron::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order)
 {
-  if(_orig) _orig->getHessShapeFunctions(u, v, w, s, o);
+  if(_orig) _orig->getHessShapeFunctions(u, v, w, s, order);
+}
+
+void MSubTetrahedron::getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order)
+{
+  if(_orig) _orig->getThirdDerivativeShapeFunctions(u, v, w, s, order);
+}
+
+int MSubTetrahedron::getNumShapeFunctions()
+{
+  if(_orig) _orig->getNumShapeFunctions();
+}
+
+int MSubTetrahedron::getNumPrimaryShapeFunctions()
+{
+  if(_orig) _orig->getNumPrimaryShapeFunctions();
+}
+
+MVertex* MSubTetrahedron::getShapeFunctionNode(int i)
+{
+  if(_orig) return _orig->getShapeFunctionNode(i);
 }
 
 bool MSubTetrahedron::isInside(double u, double v, double w)
@@ -110,8 +130,8 @@ void MSubTetrahedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 
 MSubTriangle::~MSubTriangle()
 {
-  if(_owner)
-    delete _orig;
+//   if(_owner)
+//     delete _orig;
 }
 
 const nodalBasis* MSubTriangle::getFunctionSpace(int order) const
@@ -126,19 +146,39 @@ const JacobianBasis* MSubTriangle::getJacobianFuncSpace(int order) const
   return 0;
 }
 
-void MSubTriangle::getShapeFunctions(double u, double v, double w, double s[], int o)
+void MSubTriangle::getShapeFunctions(double u, double v, double w, double s[], int order)
 {
-  if(_orig) _orig->getShapeFunctions(u, v, w, s, o);
+  if(_orig) _orig->getShapeFunctions(u, v, w, s, order);
 }
 
-void MSubTriangle::getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
+void MSubTriangle::getGradShapeFunctions(double u, double v, double w, double s[][3], int order)
 {
-  if(_orig) _orig->getGradShapeFunctions(u, v, w, s, o);
+  if(_orig) _orig->getGradShapeFunctions(u, v, w, s, order);
 }
 
-void MSubTriangle::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o)
+void MSubTriangle::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order)
 {
-  if(_orig) _orig->getHessShapeFunctions(u, v, w, s, o);
+  if(_orig) _orig->getHessShapeFunctions(u, v, w, s, order);
+}
+
+void MSubTriangle::getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order)
+{
+  if(_orig) _orig->getThirdDerivativeShapeFunctions(u, v, w, s, order);
+}
+
+int MSubTriangle::getNumShapeFunctions()
+{
+  if(_orig) _orig->getNumShapeFunctions();
+}
+
+int MSubTriangle::getNumPrimaryShapeFunctions()
+{
+  if(_orig) _orig->getNumPrimaryShapeFunctions();
+}
+
+MVertex* MSubTriangle::getShapeFunctionNode(int i)
+{
+  if(_orig) return _orig->getShapeFunctionNode(i);
 }
 
 bool MSubTriangle::isInside(double u, double v, double w)
@@ -205,8 +245,8 @@ void MSubTriangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 
 MSubLine::~MSubLine()
 {
-  if(_owner)
-    delete _orig;
+//   if(_owner)
+//     delete _orig;
 }
 
 const nodalBasis* MSubLine::getFunctionSpace(int order) const
@@ -221,19 +261,39 @@ const JacobianBasis* MSubLine::getJacobianFuncSpace(int order) const
   return 0;
 }
 
-void MSubLine::getShapeFunctions(double u, double v, double w, double s[], int o)
+void MSubLine::getShapeFunctions(double u, double v, double w, double s[], int order)
 {
-  if(_orig) _orig->getShapeFunctions(u, v, w, s, o);
+  if(_orig) _orig->getShapeFunctions(u, v, w, s, order);
 }
 
-void MSubLine::getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
+void MSubLine::getGradShapeFunctions(double u, double v, double w, double s[][3], int order)
 {
-  if(_orig) _orig->getGradShapeFunctions(u, v, w, s, o);
+  if(_orig) _orig->getGradShapeFunctions(u, v, w, s, order);
 }
 
-void MSubLine::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o)
+void MSubLine::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order)
 {
-  if(_orig) _orig->getHessShapeFunctions(u, v, w, s, o);
+  if(_orig) _orig->getHessShapeFunctions(u, v, w, s, order);
+}
+
+void MSubLine::getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order)
+{
+  if(_orig) _orig->getThirdDerivativeShapeFunctions(u, v, w, s, order);
+}
+
+int MSubLine::getNumShapeFunctions()
+{
+  if(_orig) _orig->getNumShapeFunctions();
+}
+
+int MSubLine::getNumPrimaryShapeFunctions()
+{
+  if(_orig) _orig->getNumPrimaryShapeFunctions();
+}
+
+MVertex* MSubLine::getShapeFunctionNode(int i)
+{
+  if(_orig) return _orig->getShapeFunctionNode(i);
 }
 
 bool MSubLine::isInside(double u, double v, double w)
@@ -291,12 +351,13 @@ void MSubLine::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
   *pts = _intpt;
 }
 
+
 // MSubPoint
 
 MSubPoint::~MSubPoint()
 {
-  if(_owner)
-    delete _orig;
+//   if(_owner)
+//     delete _orig;
 }
 
 const nodalBasis* MSubPoint::getFunctionSpace(int order) const
@@ -311,19 +372,39 @@ const JacobianBasis* MSubPoint::getJacobianFuncSpace(int order) const
   return 0;
 }
 
-void MSubPoint::getShapeFunctions(double u, double v, double w, double s[], int o)
+void MSubPoint::getShapeFunctions(double u, double v, double w, double s[], int order)
+{
+  if(_orig) _orig->getShapeFunctions(u, v, w, s, order);
+}
+
+void MSubPoint::getGradShapeFunctions(double u, double v, double w, double s[][3], int order)
+{
+  if(_orig) _orig->getGradShapeFunctions(u, v, w, s, order);
+}
+
+void MSubPoint::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order)
+{
+  if(_orig) _orig->getHessShapeFunctions(u, v, w, s, order);
+}
+
+void MSubPoint::getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order)
+{
+  if(_orig) _orig->getThirdDerivativeShapeFunctions(u, v, w, s, order);
+}
+
+int MSubPoint::getNumShapeFunctions()
 {
-  if(_orig) _orig->getShapeFunctions(u, v, w, s, o);
+  if(_orig) _orig->getNumShapeFunctions();
 }
 
-void MSubPoint::getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
+int MSubPoint::getNumPrimaryShapeFunctions()
 {
-  if(_orig) _orig->getGradShapeFunctions(u, v, w, s, o);
+  if(_orig) _orig->getNumPrimaryShapeFunctions();
 }
 
-void MSubPoint::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o)
+MVertex* MSubPoint::getShapeFunctionNode(int i)
 {
-  if(_orig) _orig->getHessShapeFunctions(u, v, w, s, o);
+  if(_orig) return _orig->getShapeFunctionNode(i);
 }
 
 bool MSubPoint::isInside(double u, double v, double w)
diff --git a/Geo/MSubElement.h b/Geo/MSubElement.h
index c2cd0ec5a1f6503925b924dc173cc7314b21b959..33ebf67efc401fefeeaa455e83e6fffe26910f3e 100644
--- a/Geo/MSubElement.h
+++ b/Geo/MSubElement.h
@@ -38,9 +38,13 @@ class MSubTetrahedron : public MTetrahedron
   virtual int getTypeForMSH() const { return MSH_TET_SUB; }
   virtual const nodalBasis* getFunctionSpace(int order=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
-  virtual void getShapeFunctions(double u, double v, double w, double s[], int o);
-  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o);
-  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o);
+  virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1);
+  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int order=-1);
+  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order=-1);
+  virtual void getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order=-1);
+  virtual int getNumShapeFunctions();
+  virtual int getNumPrimaryShapeFunctions();
+  virtual MVertex* getShapeFunctionNode(int i);
   // the parametric coordinates are the coordinates in the local parent element
   virtual bool isInside(double u, double v, double w);
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
@@ -75,9 +79,13 @@ class MSubTriangle : public MTriangle
   virtual int getTypeForMSH() const { return MSH_TRI_SUB; }
   virtual const nodalBasis* getFunctionSpace(int order=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
-  virtual void getShapeFunctions(double u, double v, double w, double s[], int o);
-  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o);
-  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o);
+  virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1);
+  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int order=-1);
+  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order=-1);
+  virtual void getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order=-1);
+  virtual int getNumShapeFunctions();
+  virtual int getNumPrimaryShapeFunctions();
+  virtual MVertex* getShapeFunctionNode(int i);
   // the parametric coordinates are the coordinates in the local parent element
   virtual bool isInside(double u, double v, double w);
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
@@ -113,9 +121,13 @@ class MSubLine : public MLine
   virtual int getTypeForMSH() const { return MSH_LIN_SUB; }
   virtual const nodalBasis* getFunctionSpace(int order=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
-  virtual void getShapeFunctions(double u, double v, double w, double s[], int o);
-  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o);
-  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o);
+  virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1);
+  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int order=-1);
+  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order=-1);
+  virtual void getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order=-1);
+  virtual int getNumShapeFunctions();
+  virtual int getNumPrimaryShapeFunctions();
+  virtual MVertex* getShapeFunctionNode(int i);
   // the parametric coordinates are the coordinates in the local parent element
   virtual bool isInside(double u, double v, double w);
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
@@ -150,9 +162,13 @@ class MSubPoint : public MPoint
   virtual int getTypeForMSH() const { return MSH_PNT_SUB; }
   virtual const nodalBasis* getFunctionSpace(int order=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
-  virtual void getShapeFunctions(double u, double v, double w, double s[], int o);
-  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o);
-  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o);
+  virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1);
+  virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int order=-1);
+  virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order=-1);
+  virtual void getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order=-1);
+  virtual int getNumShapeFunctions();
+  virtual int getNumPrimaryShapeFunctions();
+  virtual MVertex* getShapeFunctionNode(int i);
   // the parametric coordinates are the coordinates in the local parent element
   virtual bool isInside(double u, double v, double w);
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);