diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index adcd2ec447bb3843c20f3595e152292768ca6cb6..f277c07b45f2597e822e0d7c41ec874d76693d41 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -309,6 +309,13 @@ std::string MElement::getInfoString()
   return std::string(tmp);
 }
 
+const nodalBasis* MElement::getFunctionSpace(int order, bool serendip) const
+{
+  if (order == -1) return BasisFactory::getNodalBasis(getTypeForMSH());
+  int tag = ElementType::getTag(getType(), order, serendip);
+  return tag ? BasisFactory::getNodalBasis(tag) : NULL;
+}
+
 static double _computeDeterminantAndRegularize(const MElement *ele, double jac[3][3])
 {
   double dJ = 0;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 99637b1b4895611780962bb3ae1ca3f47795dacb..7df20c6ae0d1fa867bb22ed7746cd698ddcf2db6 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -225,7 +225,7 @@ class MElement
   virtual std::string getInfoString();
 
   // get the function space for the element
-  virtual const nodalBasis* getFunctionSpace(int order=-1) const { return 0; }
+  virtual const nodalBasis* getFunctionSpace(int order=-1, bool serendip=false) const;
 
   // get the function space for the jacobian of the element
   virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const { return 0; }
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index f1a011cb7456378bcbbb927dd5037097b78d56ee..a8229ea0219f4d4c73e2d4110818457d61836a1c 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -121,9 +121,9 @@ class MPolyhedron : public MElement {
       vol += _parts[i]->getVolume();
     return vol;
   }
-  virtual const nodalBasis* getFunctionSpace(int order=-1) const
+  virtual const nodalBasis* getFunctionSpace(int order=-1, bool serendip=false) const
   {
-    return (_orig ? _orig->getFunctionSpace(order) : 0);
+    return (_orig ? _orig->getFunctionSpace(order, serendip) : 0);
   }
   virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const
   {
@@ -271,9 +271,9 @@ class MPolygon : public MElement {
   virtual int getNumChildren() const { return _parts.size(); }
   virtual MElement *getChild(int i) const { return _parts[i]; }
   virtual bool ownsParent() const { return _owner; }
-  virtual const nodalBasis* getFunctionSpace(int order=-1) const
+  virtual const nodalBasis* getFunctionSpace(int order=-1, bool serendip=false) const
   {
-    return (_orig ? _orig->getFunctionSpace(order) : 0);
+    return (_orig ? _orig->getFunctionSpace(order, serendip) : 0);
   }
   virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const
   {
@@ -341,9 +341,9 @@ class MLineChild : public MLine {
       delete _orig;
   }
   virtual int getTypeForMSH() const { return MSH_LIN_C; }
-  virtual const nodalBasis* getFunctionSpace(int order=-1) const
+  virtual const nodalBasis* getFunctionSpace(int order=-1, bool serendip=false) const
   {
-    if(_orig) return _orig->getFunctionSpace(order);
+    if(_orig) return _orig->getFunctionSpace(order, serendip);
     return 0;
   }
   virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const
diff --git a/Geo/MHexahedron.cpp b/Geo/MHexahedron.cpp
index f862e644f496da580bcfc93f25bf20dff7747d60..31f5e74d7184e8778d4aa5280d7c06a7a0f97662 100644
--- a/Geo/MHexahedron.cpp
+++ b/Geo/MHexahedron.cpp
@@ -193,14 +193,6 @@ int MHexahedronN::getNumEdgesRep(bool curved)
   return curved ? 12 * CTX::instance()->mesh.numSubEdges : 12;
 }
 
-const nodalBasis* MHexahedron::getFunctionSpace(int order) const
-{
-  if (order == -1) return BasisFactory::getNodalBasis(getTypeForMSH());
-
-  int tag = ElementType::getTag(TYPE_HEX, order);
-  return tag ? BasisFactory::getNodalBasis(tag) : NULL;
-}
-
 const JacobianBasis* MHexahedron::getJacobianFuncSpace(int order) const
 {
   if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
diff --git a/Geo/MHexahedron.h b/Geo/MHexahedron.h
index 22ddc41034d1922e978e53faa580874fc4bae585..2b29dea0a8bd330f9c437a434e5e2eba784b1e42 100644
--- a/Geo/MHexahedron.h
+++ b/Geo/MHexahedron.h
@@ -59,7 +59,6 @@ class MHexahedron : public MElement {
   virtual MVertex *getVertex(int num){ return _v[num]; }
   virtual const MVertex *getVertex(int num) const { return _v[num]; }
   virtual void setVertex(int num,  MVertex *v){ _v[num] = v; }
-  virtual const nodalBasis* getFunctionSpace(int o=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
   virtual MVertex *getVertexDIFF(int num)
   {
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index 4dfc65ff85bd7295c09b72c207a2c32f64166390..c0859e53588ff2467dca0f8b6c3902f83b905b9e 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -11,14 +11,6 @@
 #include "Context.h"
 #include "qualityMeasures.h"
 
-const nodalBasis* MLine::getFunctionSpace(int order) const
-{
-  if (order == -1) return BasisFactory::getNodalBasis(getTypeForMSH());
-
-  int tag = ElementType::getTag(TYPE_LIN, order);
-  return tag ? BasisFactory::getNodalBasis(tag) : NULL;
-}
-
 const JacobianBasis* MLine::getJacobianFuncSpace(int order) const
 {
   if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
diff --git a/Geo/MLine.h b/Geo/MLine.h
index d6e47489c401f9cc3d1d2677000b304c2f5b3dbb..98d6d240a7ac0dd22df33a965ff39a4a48f91b14 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -74,7 +74,6 @@ class MLine : public MElement {
   {
     MVertex *tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
   }
-  virtual const nodalBasis* getFunctionSpace(int o=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
   virtual bool isInside(double u, double v, double w) const
   {
diff --git a/Geo/MPoint.h b/Geo/MPoint.h
index 0ce0bee58e6f32dfe7ee7c20aa800d104eebea86..df733f19100f652fee8aea3e571708db6cabb5cf 100644
--- a/Geo/MPoint.h
+++ b/Geo/MPoint.h
@@ -61,10 +61,6 @@ class MPoint : public MElement {
   {
     s[0][0] = s[0][1] = s[0][2] = 0.;
   }
-  virtual const nodalBasis* getFunctionSpace(int o) const
-  {
-    return BasisFactory::getNodalBasis(MSH_PNT);
-  }
   virtual const JacobianBasis* getJacobianFuncSpace(int o) const
   {
     return BasisFactory::getJacobianBasis(MSH_PNT);
diff --git a/Geo/MPrism.cpp b/Geo/MPrism.cpp
index b83af8b1a6291a2007ffa0658e0cd4f9484e0238..36dcfc1d100a6f8101aaa068cb8a3fbb4412db98 100644
--- a/Geo/MPrism.cpp
+++ b/Geo/MPrism.cpp
@@ -32,14 +32,6 @@ void MPrism::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
   *pts = getGQPriPts(pOrder);
 }
 
-const nodalBasis* MPrism::getFunctionSpace(int order) const
-{
-  if (order == -1) return BasisFactory::getNodalBasis(getTypeForMSH());
-
-  int tag = ElementType::getTag(TYPE_PRI, order);
-  return tag ? BasisFactory::getNodalBasis(tag) : NULL;
-}
-
 const JacobianBasis* MPrism::getJacobianFuncSpace(int order) const
 {
   if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
diff --git a/Geo/MPrism.h b/Geo/MPrism.h
index 2a2075f1743d54b16ae15da317e8411889d42098..1ac48c5d0d4656076acf72570321cd46267059dd 100644
--- a/Geo/MPrism.h
+++ b/Geo/MPrism.h
@@ -128,7 +128,6 @@ class MPrism : public MElement {
     tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
     tmp = _v[3]; _v[3] = _v[4]; _v[4] = tmp;
   }
-  virtual const nodalBasis* getFunctionSpace(int o=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
   virtual int getVolumeSign();
   virtual void getNode(int num, double &u, double &v, double &w) const
diff --git a/Geo/MPyramid.cpp b/Geo/MPyramid.cpp
index 4bf91037f721b966fec93a87a21850b70ba643fc..1d2af2da9164281b66d813507722023653676f36 100644
--- a/Geo/MPyramid.cpp
+++ b/Geo/MPyramid.cpp
@@ -33,14 +33,6 @@ void MPyramid::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
   *pts = getGQPyrPts(pOrder);
 }
 
-const nodalBasis* MPyramid::getFunctionSpace(int order) const
-{
-  if (order == -1) return BasisFactory::getNodalBasis(getTypeForMSH());
-
-  int tag = ElementType::getTag(TYPE_PYR, order);
-  return tag ? BasisFactory::getNodalBasis(tag) : NULL;
-}
-
 const JacobianBasis* MPyramid::getJacobianFuncSpace(int o) const
 {
   // FIXME add other order and see MPyramid::getFunctionSpace for 'design'
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index 423a7d0ae1e86331b01066b9bad4443ffdfcc3b1..7a3623ffabb9f4da3a97a31af631d30dd08f43a9 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -69,7 +69,6 @@ class MPyramid : public MElement {
   virtual MVertex *getVertex(int num){ return _v[num]; }
   virtual const MVertex *getVertex(int num) const{ return _v[num]; }
   virtual void setVertex(int num,  MVertex *v){ _v[num] = v; }
-  virtual const nodalBasis* getFunctionSpace(int o=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
   virtual int getNumEdges(){ return 8; }
   virtual MEdge getEdge(int num) const
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index 271c6d70082be4988bb5eb126dcc0880b44efdde..95c6313636a40bc0815a3fee083344b9722d4ee8 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -17,14 +17,6 @@
 
 #define SQU(a)      ((a)*(a))
 
-const nodalBasis* MQuadrangle::getFunctionSpace(int order) const
-{
-  if (order == -1) return BasisFactory::getNodalBasis(getTypeForMSH());
-
-  int tag = ElementType::getTag(TYPE_QUA, order);
-  return tag ? BasisFactory::getNodalBasis(tag) : NULL;
-}
-
 const JacobianBasis* MQuadrangle::getJacobianFuncSpace(int order) const
 {
   if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index 07f171f142e71ae57504955e1cf6c0190647962e..738ae4097c7ad4f6eec7cd09f3d98791da77019a 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -117,7 +117,6 @@ class MQuadrangle : public MElement {
   virtual const char *getStringForBDF() const { return "CQUAD4"; }
   virtual const char *getStringForDIFF() const { return "ElmB4n2D"; }
   virtual const char *getStringForINP() const { return "CPS4"/*"C2D4"*/; }
-  virtual const nodalBasis* getFunctionSpace(int o=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
   virtual void getNode(int num, double &u, double &v, double &w) const
   {
diff --git a/Geo/MSubElement.cpp b/Geo/MSubElement.cpp
index 3a4e24db8dc0c1abcd1541d1ce98d6697e9f6492..ae58a393d8ecd4785ff69da95161a1a5ab411d53 100644
--- a/Geo/MSubElement.cpp
+++ b/Geo/MSubElement.cpp
@@ -17,9 +17,9 @@ MSubTetrahedron::~MSubTetrahedron()
   if(_base) delete _base;
 }
 
-const nodalBasis* MSubTetrahedron::getFunctionSpace(int order) const
+const nodalBasis* MSubTetrahedron::getFunctionSpace(int order, bool serendip) const
 {
-  if(_orig) return _orig->getFunctionSpace(order);
+  if(_orig) return _orig->getFunctionSpace(order, serendip);
   return 0;
 }
 
@@ -209,9 +209,9 @@ MSubTriangle::~MSubTriangle()
   if(_base) delete _base;
 }
 
-const nodalBasis* MSubTriangle::getFunctionSpace(int order) const
+const nodalBasis* MSubTriangle::getFunctionSpace(int order, bool serendip) const
 {
-  if(_orig) return _orig->getFunctionSpace(order);
+  if(_orig) return _orig->getFunctionSpace(order, serendip);
   return 0;
 }
 
@@ -447,9 +447,9 @@ MSubLine::~MSubLine()
   if(_base) delete _base;
 }
 
-const nodalBasis* MSubLine::getFunctionSpace(int order) const
+const nodalBasis* MSubLine::getFunctionSpace(int order, bool serendip) const
 {
-  if(_orig) return _orig->getFunctionSpace(order);
+  if(_orig) return _orig->getFunctionSpace(order, serendip);
   return 0;
 }
 
@@ -676,9 +676,9 @@ MSubPoint::~MSubPoint()
   if(_base) delete _base;
 }
 
-const nodalBasis* MSubPoint::getFunctionSpace(int order) const
+const nodalBasis* MSubPoint::getFunctionSpace(int order, bool serendip) const
 {
-  if(_orig) return _orig->getFunctionSpace(order);
+  if(_orig) return _orig->getFunctionSpace(order, serendip);
   return 0;
 }
 
diff --git a/Geo/MSubElement.h b/Geo/MSubElement.h
index 72404d034e8126adf129c40947ba2c89abb9e535..9415ead679716660ab1d0bf9db955717a69e9034 100644
--- a/Geo/MSubElement.h
+++ b/Geo/MSubElement.h
@@ -41,7 +41,7 @@ class MSubTetrahedron : public MTetrahedron
   ~MSubTetrahedron();
 
   virtual int getTypeForMSH() const { return MSH_TET_SUB; }
-  virtual const nodalBasis* getFunctionSpace(int order=-1) const;
+  virtual const nodalBasis* getFunctionSpace(int order=-1, bool serendip=false) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
   // the parametric coordinates are the coordinates in the local parent element
   virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1) const;
@@ -99,7 +99,7 @@ class MSubTriangle : public MTriangle
   ~MSubTriangle();
 
   virtual int getTypeForMSH() const { return MSH_TRI_SUB; }
-  virtual const nodalBasis* getFunctionSpace(int order=-1) const;
+  virtual const nodalBasis* getFunctionSpace(int order=-1, bool serendip=false) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
   // the parametric coordinates are the coordinates in the local parent element
   virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1) const;
@@ -157,7 +157,7 @@ class MSubLine : public MLine
   ~MSubLine();
 
   virtual int getTypeForMSH() const { return MSH_LIN_SUB; }
-  virtual const nodalBasis* getFunctionSpace(int order=-1) const;
+  virtual const nodalBasis* getFunctionSpace(int order=-1, bool serendip=false) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
   // the parametric coordinates are the coordinates in the local parent element
   virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1) const;
@@ -215,7 +215,7 @@ class MSubPoint : public MPoint
   ~MSubPoint();
 
   virtual int getTypeForMSH() const { return MSH_PNT_SUB; }
-  virtual const nodalBasis* getFunctionSpace(int order=-1) const;
+  virtual const nodalBasis* getFunctionSpace(int order=-1, bool serendip=false) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
   // the parametric coordinates are the coordinates in the local parent element
   virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1) const;
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index 6a7ce299315a4c70ece950d5308dbc4e0d82194c..43e5587f546a76db4baf1b673e7f30ce6c19b9ae 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -100,14 +100,6 @@ void MTetrahedron::xyz2uvw(double xyz[3], double uvw[3]) const
   sys3x3(mat, b, uvw, &det);
 }
 
-const nodalBasis* MTetrahedron::getFunctionSpace(int order) const
-{
-  if (order == -1) return BasisFactory::getNodalBasis(getTypeForMSH());
-
-  int tag = ElementType::getTag(TYPE_TET, order);
-  return tag ? BasisFactory::getNodalBasis(tag) : NULL;
-}
-
 const JacobianBasis* MTetrahedron::getJacobianFuncSpace(int order) const
 {
   if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index 245fdc09a87b4b11026cc9150be86d0882972feb..721b88017b29a1dbbe092a8248d11776ea4efc07 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -129,7 +129,6 @@ class MTetrahedron : public MElement {
   virtual double getCircumRadius();
   virtual double etaShapeMeasure();
   virtual void xyz2uvw(double xyz[3], double uvw[3]) const;
-  virtual const nodalBasis* getFunctionSpace(int o=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
   virtual void getNode(int num, double &u, double &v, double &w) const
   {
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index 28ab2572276c4b8fd752493021fae3a5c0fb0319..587d7a1a520771ad273ac22aecfd4795d94e1169 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -116,14 +116,6 @@ void MTriangle::xyz2uvw(double xyz[3], double uvw[3]) const
   uvw[2] = 0.;
 }
 
-const nodalBasis* MTriangle::getFunctionSpace(int order) const
-{
-  if (order == -1) return BasisFactory::getNodalBasis(getTypeForMSH());
-
-  int tag = ElementType::getTag(TYPE_TRI, order);
-  return tag ? BasisFactory::getNodalBasis(tag) : NULL;
-}
-
 const JacobianBasis* MTriangle::getJacobianFuncSpace(int order) const
 {
   if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index dbda43a1f23410ace1b8bd70e0b743a1269810a8..0d648ac219b632cacc04c7aa8cc13411978c1ee1 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -124,7 +124,6 @@ class MTriangle : public MElement {
   {
     MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
   }
-  virtual const nodalBasis* getFunctionSpace(int o=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
   virtual void getNode(int num, double &u, double &v, double &w) const
   {