diff --git a/Geo/MElement.h b/Geo/MElement.h
index b63efac170707c95aa9cbd9f665c658789c43f69..533c877a49979e5e99f475ba71a9dff1ec1a9601 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -161,6 +161,10 @@ class MElement
   virtual int getNumChildren() const { return 0; }
   virtual MElement *getChild(int i) const { return NULL; }
   virtual bool ownsParent() const { return false; }
+  // get base element in case of MSubElement
+  virtual const MElement *getBaseElement() const { return this; }
+  virtual MElement *getBaseElement() { return this; }
+
   // get and set domain for borders
   virtual MElement *getDomain(int i) const { return NULL; }
   virtual void setDomain (MElement *e, int i) { }
@@ -239,10 +243,10 @@ class MElement
                                      int order=-1);
   // return the Jacobian of the element evaluated at point (u,v,w) in
   // parametric coordinates
-  double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
+  virtual double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
   // To be compatible with _vgrads of functionSpace without having to put under fullMatrix form
-  double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]);
-  double getJacobian(double u, double v, double w, double jac[3][3]);
+  virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]);
+  virtual double getJacobian(double u, double v, double w, double jac[3][3]);
   inline double getJacobian(double u, double v, double w, fullMatrix<double> &j){
     double JAC[3][3];
     const double detJ = getJacobian (u,v,w,JAC);
@@ -253,7 +257,7 @@ class MElement
     }
     return detJ;
   }
-  double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
+  virtual double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
   double getJacobianDeterminant(double u, double v, double w)
   {
     double jac[3][3]; return getJacobian(u, v, w, jac);
diff --git a/Geo/MSubElement.cpp b/Geo/MSubElement.cpp
index 93534c4d0e1db6e11d619a5a30ad4186f4df2b7d..11f82cd0091e5504b10eae5cfc544b19da96aeea 100644
--- a/Geo/MSubElement.cpp
+++ b/Geo/MSubElement.cpp
@@ -7,12 +7,107 @@
 //   Frederic Duboeuf
 
 #include "MSubElement.h"
+#include "Numeric.h"
 
 // MSubTetrahedron
 
 MSubTetrahedron::~MSubTetrahedron()
 {
-  if(_intpt) delete [] _intpt;
+  if(_pts) delete [] _pts;
+  if(_base) delete _base;
+}
+
+const nodalBasis* MSubTetrahedron::getFunctionSpace(int order) const
+{
+  if(_orig) return _orig->getFunctionSpace(order);
+  return 0;
+}
+
+const JacobianBasis* MSubTetrahedron::getJacobianFuncSpace(int order) const
+{
+  if(_orig) return _orig->getJacobianFuncSpace(order);
+  return 0;
+}
+
+void MSubTetrahedron::getShapeFunctions(double u, double v, double w, double s[], int order)
+{
+  if(_orig) _orig->getShapeFunctions(u, v, w, s, order);
+}
+
+void MSubTetrahedron::getGradShapeFunctions(double u, double v, double w, double s[][3], int order)
+{
+  if(_orig) _orig->getGradShapeFunctions(u, v, w, s, order);
+}
+
+void MSubTetrahedron::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order)
+{
+  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);
+}
+
+double MSubTetrahedron::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
+{
+  if(_orig) return _orig->getJacobian(gsf, jac);
+  return 0;
+}
+double MSubTetrahedron::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])
+{
+  if(_orig) return _orig->getJacobian(gsf, jac);
+  return 0;
+}
+double MSubTetrahedron::getJacobian(double u, double v, double w, double jac[3][3])
+{
+  if(_orig) return _orig->getJacobian(u, v, w, jac);
+  return 0;
+}
+double MSubTetrahedron::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
+{
+  if(_orig) return _orig->getPrimaryJacobian(u, v, w, jac);
+  return 0;
+}
+
+int MSubTetrahedron::getNumShapeFunctions()
+{
+  if(_orig) return _orig->getNumShapeFunctions();
+  return 0;
+}
+
+int MSubTetrahedron::getNumPrimaryShapeFunctions()
+{
+  if(_orig) return _orig->getNumPrimaryShapeFunctions();
+  return 0;
+}
+
+MVertex* MSubTetrahedron::getShapeFunctionNode(int i)
+{
+  if(_orig) return _orig->getShapeFunctionNode(i);
+  return 0;
+}
+
+void MSubTetrahedron::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w)
+{
+  if(!_orig) return;
+  SPoint3 p;
+  _orig->pnt(u, v, w, p);
+  double xyz[3] = {p.x(), p.y(), p.z()};
+  double uvwE[3];
+  getBaseElement()->xyz2uvw(xyz, uvwE);
+  u = uvwE[0]; v = uvwE[1]; w = uvwE[2];
+}
+
+void MSubTetrahedron::movePointFromElementSpaceToParentSpace(double &u, double &v, double &w)
+{
+  if(!_orig) return;
+  SPoint3 p;
+  getBaseElement()->pnt(u, v, w, p);
+  double xyz[3] = {p.x(), p.y(), p.z()};
+  double uvwP[3];
+  _orig->xyz2uvw(xyz, uvwP);
+  u = uvwP[0]; v = uvwP[1]; w = uvwP[2];
 }
 
 bool MSubTetrahedron::isInside(double u, double v, double w)
@@ -40,44 +135,57 @@ bool MSubTetrahedron::isInside(double u, double v, double w)
 
 void MSubTetrahedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 {
-  *npts=0;
-  if(_intpt) delete [] _intpt;
-  if(!_orig) return;
-  _intpt = new IntPt[getNGQTetPts(pOrder)];
-  int nptsi;
-  IntPt *ptsi;
+  if(_pts)
+  {
+    if(pOrder==_pOrder)
+    {
+      *npts = _npts;
+      *pts = _pts;
+      return;
+    }
+    else
+      delete [] _pts;
+  }
 
-  // work in the parametric space of the parent element
-  // (i) get the parametric coordinates of MSubElement vertices in this space
-  double v_uvw[4][3];
-  for(int i=0; i<4; ++i)
+  _pOrder = pOrder;
+
+  if(!_orig)
   {
-    MVertex *vi = getVertex(i);
-    double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
-    _orig->xyz2uvw(v_xyz, v_uvw[i]);
+    getBaseElement()->getIntegrationPoints(pOrder, &_npts, &_pts);
+    *npts = _npts;
+    *pts = _pts;
+    return;
   }
-  // (ii) build a MElement on the vertices with these coordinates in this space
-  MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
-  MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
-  MVertex v2(v_uvw[2][0], v_uvw[2][1], v_uvw[2][2]);
-  MVertex v3(v_uvw[3][0], v_uvw[3][1], v_uvw[3][2]);
-  MTetrahedron tt(&v0, &v1, &v2, &v3);
-  // (iii) get the integration points of the new element in its parametric space 
-  tt.getIntegrationPoints(pOrder, &nptsi, &ptsi);
-  // (iv) get the coordinates of these points in the parametric space of parent element
-  for(int ip=0; ip<nptsi; ++ip)
+
+  // work in the parametric space of the parent element
+  *npts=0;
+  _pts = new IntPt[getNGQTetPts(pOrder)];
+
+  // (i) get the integration points of the base element in its parametric space
+  IntPt *ptsb;
+  getBaseElement()->getIntegrationPoints(pOrder, &_npts, &ptsb);
+
+  // (ii) get the coordinates of these points in the parametric space of parent element
+  double u,v,w;
+  double jac[3][3];
+  double baseJac, origJac;
+  for(int i=0; i<_npts; ++i)
   {
-    const double u = ptsi[ip].pt[0];
-    const double v = ptsi[ip].pt[1];
-    const double w = ptsi[ip].pt[2];
-    SPoint3 p; tt.pnt(u, v, w, p);
-    _intpt[*npts + ip].pt[0] = p.x();
-    _intpt[*npts + ip].pt[1] = p.y();
-    _intpt[*npts + ip].pt[2] = p.z();
-    _intpt[*npts + ip].weight = ptsi[ip].weight;
+    u = ptsb[i].pt[0];
+    v = ptsb[i].pt[1];
+    w = ptsb[i].pt[2];
+    baseJac = getBaseElement()->getJacobian(u, v, w, jac);
+
+    movePointFromElementSpaceToParentSpace(u, v, w);
+    origJac = _orig->getJacobian(u, v, w, jac);
+
+    _pts[*npts + i].pt[0] = u;
+    _pts[*npts + i].pt[1] = v;
+    _pts[*npts + i].pt[2] = w;
+    _pts[*npts + i].weight = ptsb[i].weight * baseJac/origJac;
   }
-  *npts = nptsi;
-  *pts = _intpt;
+  *npts = _npts;
+  *pts = _pts;
 }
 
 
@@ -85,7 +193,148 @@ void MSubTetrahedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 
 MSubTriangle::~MSubTriangle()
 {
-  if(_intpt) delete [] _intpt;
+  if(_pts) delete [] _pts;
+  if(_base) delete _base;
+}
+
+const nodalBasis* MSubTriangle::getFunctionSpace(int order) const
+{
+  if(_orig) return _orig->getFunctionSpace(order);
+  return 0;
+}
+
+const JacobianBasis* MSubTriangle::getJacobianFuncSpace(int order) const
+{
+  if(_orig) return _orig->getJacobianFuncSpace(order);
+  return 0;
+}
+
+void MSubTriangle::getShapeFunctions(double u, double v, double w, double s[], int order)
+{
+  if(_orig) _orig->getShapeFunctions(u, v, w, s, order);
+}
+
+void MSubTriangle::getGradShapeFunctions(double u, double v, double w, double s[][3], int order)
+{
+  if(!_orig)
+    return;
+
+  if (_orig->getDim()==getDim())
+    return _orig->getGradShapeFunctions(u, v, w, s, order);
+
+  int nsf = getNumShapeFunctions();
+  double gradsuvw[nsf][3];
+  _orig->getGradShapeFunctions(u, v, w, gradsuvw, order);
+
+  // work in the parametric space of the parent element
+  double jac[3][3];
+  double invjac[3][3];
+  _orig->getJacobian(u, v, w, jac);
+  inv3x3(jac, invjac);
+
+  MEdge edge[2];
+  edge[0]=getBaseElement()->getEdge(0);
+  edge[1]=getBaseElement()->getEdge(1);
+  SVector3 tang[2];
+  tang[0]=edge[0].tangent();
+  tang[1]=edge[1].tangent();
+  SVector3 vect = crossprod(tang[0],tang[1]);
+  tang[1] = crossprod(vect,tang[0]);
+
+  double gradxyz[3];
+  double projgradxyz[3];
+  for (int i=0; i<nsf; ++i)
+  {
+  // (i) get the cartesian coordinates of the gradient
+    gradxyz[0] = invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2];
+    gradxyz[1] = invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2];
+    gradxyz[2] = invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2];
+
+    // (ii) projection of the gradient on edges in the cartesian space
+    SVector3 grad(&gradxyz[0]);
+    double prodscal[2];
+    prodscal[0]=dot(tang[0],grad);
+    prodscal[1]=dot(tang[1],grad);
+    projgradxyz[0] = prodscal[0]*tang[0].x() + prodscal[1]*tang[1].x();
+    projgradxyz[1] = prodscal[0]*tang[0].y() + prodscal[1]*tang[1].y();
+    projgradxyz[2] = prodscal[0]*tang[0].z() + prodscal[1]*tang[1].z();
+
+    // (iii) get the parametric coordinates of the projection in the parametric space of the parent element
+    s[i][0] = jac[0][0] * projgradxyz[0] + jac[0][1] * projgradxyz[1] + jac[0][2] * projgradxyz[2];
+    s[i][1] = jac[1][0] * projgradxyz[0] + jac[1][1] * projgradxyz[1] + jac[1][2] * projgradxyz[2];
+    s[i][2] = jac[2][0] * projgradxyz[0] + jac[2][1] * projgradxyz[1] + jac[2][2] * projgradxyz[2];
+  }
+}
+
+void MSubTriangle::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order)
+{
+  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);
+}
+
+double MSubTriangle::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
+{
+  if(_orig) return _orig->getJacobian(gsf, jac);
+  return 0;
+}
+double MSubTriangle::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])
+{
+  if(_orig) return _orig->getJacobian(gsf, jac);
+  return 0;
+}
+double MSubTriangle::getJacobian(double u, double v, double w, double jac[3][3])
+{
+  if(_orig) return _orig->getJacobian(u, v, w, jac);
+  return 0;
+}
+double MSubTriangle::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
+{
+  if(_orig) return _orig->getPrimaryJacobian(u, v, w, jac);
+  return 0;
+}
+
+int MSubTriangle::getNumShapeFunctions()
+{
+  if(_orig) return _orig->getNumShapeFunctions();
+  return 0;
+}
+
+int MSubTriangle::getNumPrimaryShapeFunctions()
+{
+  if(_orig) return _orig->getNumPrimaryShapeFunctions();
+  return 0;
+}
+
+MVertex* MSubTriangle::getShapeFunctionNode(int i)
+{
+  if(_orig) return _orig->getShapeFunctionNode(i);
+  return 0;
+}
+
+void MSubTriangle::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w)
+{
+  if(!_orig) return;
+  SPoint3 p;
+  _orig->pnt(u, v, w, p);
+  double xyz[3] = {p.x(), p.y(), p.z()};
+  double uvwE[3];
+  getBaseElement()->xyz2uvw(xyz, uvwE);
+  u = uvwE[0]; v = uvwE[1]; w = uvwE[2];
+}
+
+void MSubTriangle::movePointFromElementSpaceToParentSpace(double &u, double &v, double &w)
+{
+  if(!_orig) return;
+  SPoint3 p;
+  getBaseElement()->pnt(u, v, w, p);
+  double xyz[3] = {p.x(), p.y(), p.z()};
+  double uvwP[3];
+  _orig->xyz2uvw(xyz, uvwP);
+  u = uvwP[0]; v = uvwP[1]; w = uvwP[2];
 }
 
 bool MSubTriangle::isInside(double u, double v, double w)
@@ -112,51 +361,195 @@ bool MSubTriangle::isInside(double u, double v, double w)
 
 void MSubTriangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 {
-  *npts=0;
-  if(_intpt) delete [] _intpt;
-  if(!_orig) return;
-  _intpt = new IntPt[getNGQTPts(pOrder)];
-  int nptsi;
-  IntPt *ptsi;
+  if(_pts)
+  {
+    if(pOrder==_pOrder)
+    {
+      *npts = _npts;
+      *pts = _pts;
+      return;
+    }
+    else
+      delete [] _pts;
+  }
 
-  // work in the parametric space of the parent element
-  // (i) get the parametric coordinates of MSubElement vertices in this space
-  double v_uvw[3][3];
-  for(int i=0; i<3; ++i)
+  _pOrder = pOrder;
+
+  if(!_orig)
   {
-    MVertex *vi = getVertex(i);
-    double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
-    _orig->xyz2uvw(v_xyz, v_uvw[i]);
+    getBaseElement()->getIntegrationPoints(pOrder, &_npts, &_pts);
+    *npts = _npts;
+    *pts = _pts;
+    return;
   }
-  // (ii) build a MElement on the vertices with these coordinates in this space
-  MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
-  MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
-  MVertex v2(v_uvw[2][0], v_uvw[2][1], v_uvw[2][2]);
-  MTriangle t(&v0, &v1, &v2);
-  // (iii) get the integration points of the new element in its parametric space 
-  t.getIntegrationPoints(pOrder, &nptsi, &ptsi);
-  // (iv) get the coordinates of these points in the parametric space of parent element
-  for(int ip=0; ip<nptsi; ++ip)
+
+  // work in the parametric space of the parent element
+  *npts=0;
+  _pts = new IntPt[getNGQTetPts(pOrder)];
+
+  // (i) get the integration points of the base element in its parametric space
+  IntPt *ptsb;
+  getBaseElement()->getIntegrationPoints(pOrder, &_npts, &ptsb);
+
+  // (ii) get the coordinates of these points in the parametric space of parent element
+  double u,v,w;
+  double jac[3][3];
+  double baseJac, origJac;
+  for(int i=0; i<_npts; ++i)
   {
-    const double u = ptsi[ip].pt[0];
-    const double v = ptsi[ip].pt[1];
-    const double w = ptsi[ip].pt[2];
-    SPoint3 p; t.pnt(u, v, w, p);
-    _intpt[*npts + ip].pt[0] = p.x();
-    _intpt[*npts + ip].pt[1] = p.y();
-    _intpt[*npts + ip].pt[2] = p.z();
-    _intpt[*npts + ip].weight = ptsi[ip].weight;
+    u = ptsb[i].pt[0];
+    v = ptsb[i].pt[1];
+    w = ptsb[i].pt[2];
+    baseJac = getBaseElement()->getJacobian(u, v, w, jac);
+
+    movePointFromElementSpaceToParentSpace(u, v, w);
+    origJac = _orig->getJacobian(u, v, w, jac);
+
+    _pts[*npts + i].pt[0] = u;
+    _pts[*npts + i].pt[1] = v;
+    _pts[*npts + i].pt[2] = w;
+    _pts[*npts + i].weight = ptsb[i].weight * baseJac/origJac;
   }
-  *npts = nptsi;
-  *pts = _intpt;
+  *npts = _npts;
+  *pts = _pts;
 }
 
-
 // MSubLine
 
 MSubLine::~MSubLine()
 {
-  if(_intpt) delete [] _intpt;
+  if(_pts) delete [] _pts;
+  if(_base) delete _base;
+}
+
+const nodalBasis* MSubLine::getFunctionSpace(int order) const
+{
+  if(_orig) return _orig->getFunctionSpace(order);
+  return 0;
+}
+
+const JacobianBasis* MSubLine::getJacobianFuncSpace(int order) const
+{
+  if(_orig) return _orig->getJacobianFuncSpace(order);
+  return 0;
+}
+
+void MSubLine::getShapeFunctions(double u, double v, double w, double s[], int order)
+{
+  if(_orig) _orig->getShapeFunctions(u, v, w, s, order);
+}
+
+void MSubLine::getGradShapeFunctions(double u, double v, double w, double s[][3], int order)
+{
+  if(!_orig)
+    return;
+
+  if (_orig->getDim()==getDim())
+    return _orig->getGradShapeFunctions(u, v, w, s, order);
+
+  int nsf = _orig->getNumShapeFunctions();
+  double gradsuvw[nsf][3];
+  _orig->getGradShapeFunctions(u, v, w, gradsuvw, order);
+
+  double jac[3][3];
+  double invjac[3][3];
+  _orig->getJacobian(u, v, w, jac);
+  inv3x3(jac, invjac);
+  MEdge edge = getBaseElement()->getEdge(0);
+  SVector3 tang = edge.tangent();
+
+  double gradxyz[3];
+  double projgradxyz[3];
+  for (int i=0; i<nsf; ++i)
+  {
+    // (i) get the cartesian coordinates of the gradient
+    gradxyz[0] = invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2];
+    gradxyz[1] = invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2];
+    gradxyz[2] = invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2];
+
+    // (ii) projection of the gradient on edges in the cartesian space
+    SVector3 grad(&gradxyz[0]);
+    double prodscal = dot(tang,grad);
+    projgradxyz[0] = prodscal * tang.x();
+    projgradxyz[1] = prodscal * tang.y();
+    projgradxyz[2] = prodscal * tang.z();
+
+    // (iii) get the parametric coordinates of the projection in the parametric space of the parent element
+    s[i][0] = jac[0][0] * projgradxyz[0] + jac[0][1] * projgradxyz[1] + jac[0][2] * projgradxyz[2];
+    s[i][1] = jac[1][0] * projgradxyz[0] + jac[1][1] * projgradxyz[1] + jac[1][2] * projgradxyz[2];
+    s[i][2] = jac[2][0] * projgradxyz[0] + jac[2][1] * projgradxyz[1] + jac[2][2] * projgradxyz[2];
+  }
+}
+
+void MSubLine::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order)
+{
+  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);
+}
+
+double MSubLine::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
+{
+  if(_orig) return _orig->getJacobian(gsf, jac);
+  return 0;
+}
+double MSubLine::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])
+{
+  if(_orig) return _orig->getJacobian(gsf, jac);
+  return 0;
+}
+double MSubLine::getJacobian(double u, double v, double w, double jac[3][3])
+{
+  if(_orig) return _orig->getJacobian(u, v, w, jac);
+  return 0;
+}
+double MSubLine::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
+{
+  if(_orig) return _orig->getPrimaryJacobian(u, v, w, jac);
+  return 0;
+}
+
+int MSubLine::getNumShapeFunctions()
+{
+  if(_orig) return _orig->getNumShapeFunctions();
+  return 0;
+}
+
+int MSubLine::getNumPrimaryShapeFunctions()
+{
+  if(_orig) return _orig->getNumPrimaryShapeFunctions();
+  return 0;
+}
+
+MVertex* MSubLine::getShapeFunctionNode(int i)
+{
+  if(_orig) return _orig->getShapeFunctionNode(i);
+  return 0;
+}
+
+void MSubLine::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w)
+{
+  if(!_orig) return;
+  SPoint3 p;
+  _orig->pnt(u, v, w, p);
+  double xyz[3] = {p.x(), p.y(), p.z()};
+  double uvwE[3];
+  getBaseElement()->xyz2uvw(xyz, uvwE);
+  u = uvwE[0]; v = uvwE[1]; w = uvwE[2];
+}
+
+void MSubLine::movePointFromElementSpaceToParentSpace(double &u, double &v, double &w)
+{
+  if(!_orig) return;
+  SPoint3 p;
+  getBaseElement()->pnt(u, v, w, p);
+  double xyz[3] = {p.x(), p.y(), p.z()};
+  double uvwP[3];
+  _orig->xyz2uvw(xyz, uvwP);
+  u = uvwP[0]; v = uvwP[1]; w = uvwP[2];
 }
 
 bool MSubLine::isInside(double u, double v, double w)
@@ -182,42 +575,57 @@ bool MSubLine::isInside(double u, double v, double w)
 
 void MSubLine::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 {
-  *npts=0;
-  if(_intpt) delete [] _intpt;
-  if(!_orig) return;
-  _intpt = new IntPt[getNGQLPts(pOrder)];
-  int nptsi;
-  IntPt *ptsi;
+  if(_pts)
+  {
+    if(pOrder==_pOrder)
+    {
+      *npts = _npts;
+      *pts = _pts;
+      return;
+    }
+    else
+      delete [] _pts;
+  }
 
-  // work in the parametric space of the parent element
-  // (i) get the parametric coordinates of MSubElement vertices in this space
-  double v_uvw[2][3];
-  for(int i=0; i<2; ++i)
+  _pOrder = pOrder;
+
+  if(!_orig)
   {
-    MVertex *vi = getVertex(i);
-    double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
-    _orig->xyz2uvw(v_xyz, v_uvw[i]);
+    getBaseElement()->getIntegrationPoints(pOrder, &_npts, &_pts);
+    *npts = _npts;
+    *pts = _pts;
+    return;
   }
-  // (ii) build a MElement on the vertices with these coordinates in this space
-  MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
-  MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
-  MLine l(&v0, &v1);
-  // (iii) get the integration points of the new element in its parametric space 
-  l.getIntegrationPoints(pOrder, &nptsi, &ptsi);
-  // (iv) get the coordinates of these points in the parametric space of parent element
-  for(int ip=0; ip<nptsi; ++ip)
+
+  // work in the parametric space of the parent element
+  *npts=0;
+  _pts = new IntPt[getNGQTetPts(pOrder)];
+
+  // (i) get the integration points of the base element in its parametric space
+  IntPt *ptsb;
+  getBaseElement()->getIntegrationPoints(pOrder, &_npts, &ptsb);
+
+  // (ii) get the coordinates of these points in the parametric space of parent element
+  double u,v,w;
+  double jac[3][3];
+  double baseJac, origJac;
+  for(int i=0; i<_npts; ++i)
   {
-    const double u = ptsi[ip].pt[0];
-    const double v = ptsi[ip].pt[1];
-    const double w = ptsi[ip].pt[2];
-    SPoint3 p; l.pnt(u, v, w, p);
-    _intpt[*npts + ip].pt[0] = p.x();
-    _intpt[*npts + ip].pt[1] = p.y();
-    _intpt[*npts + ip].pt[2] = p.z();
-    _intpt[*npts + ip].weight = ptsi[ip].weight;
+    u = ptsb[i].pt[0];
+    v = ptsb[i].pt[1];
+    w = ptsb[i].pt[2];
+    baseJac = getBaseElement()->getJacobian(u, v, w, jac);
+
+    movePointFromElementSpaceToParentSpace(u, v, w);
+    origJac = _orig->getJacobian(u, v, w, jac);
+
+    _pts[*npts + i].pt[0] = u;
+    _pts[*npts + i].pt[1] = v;
+    _pts[*npts + i].pt[2] = w;
+    _pts[*npts + i].weight = ptsb[i].weight * baseJac/origJac;
   }
-  *npts = nptsi;
-  *pts = _intpt;
+  *npts = _npts;
+  *pts = _pts;
 }
 
 
@@ -225,9 +633,102 @@ void MSubLine::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 
 MSubPoint::~MSubPoint()
 {
-  if(_intpt) delete [] _intpt;
+  if(_pts) delete [] _pts;
+  if(_base) delete _base;
+}
+
+const nodalBasis* MSubPoint::getFunctionSpace(int order) const
+{
+  if(_orig) return _orig->getFunctionSpace(order);
+  return 0;
+}
+
+const JacobianBasis* MSubPoint::getJacobianFuncSpace(int order) const
+{
+  if(_orig) return _orig->getJacobianFuncSpace(order);
+  return 0;
+}
+
+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);
+}
+
+double MSubPoint::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
+{
+  if(_orig) return _orig->getJacobian(gsf, jac);
+  return 0;
+}
+double MSubPoint::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])
+{
+  if(_orig) return _orig->getJacobian(gsf, jac);
+  return 0;
+}
+double MSubPoint::getJacobian(double u, double v, double w, double jac[3][3])
+{
+  if(_orig) return _orig->getJacobian(u, v, w, jac);
+  return 0;
+}
+double MSubPoint::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
+{
+  if(_orig) return _orig->getPrimaryJacobian(u, v, w, jac);
+  return 0;
+}
+
+int MSubPoint::getNumShapeFunctions()
+{
+  if(_orig) return _orig->getNumShapeFunctions();
+  return 0;
+}
+
+int MSubPoint::getNumPrimaryShapeFunctions()
+{
+  if(_orig) return _orig->getNumPrimaryShapeFunctions();
+  return 0;
+}
+
+MVertex* MSubPoint::getShapeFunctionNode(int i)
+{
+  if(_orig) return _orig->getShapeFunctionNode(i);
+  return 0;
+}
+
+void MSubPoint::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w)
+{
+  if(!_orig) return;
+  SPoint3 p;
+  _orig->pnt(u, v, w, p);
+  double xyz[3] = {p.x(), p.y(), p.z()};
+  double uvwE[3];
+  getBaseElement()->xyz2uvw(xyz, uvwE);
+  u = uvwE[0]; v = uvwE[1]; w = uvwE[2];
+}
+
+void MSubPoint::movePointFromElementSpaceToParentSpace(double &u, double &v, double &w)
+{
+  if(!_orig) return;
+  SPoint3 p;
+  getBaseElement()->pnt(u, v, w, p);
+  double xyz[3] = {p.x(), p.y(), p.z()};
+  double uvwP[3];
+  _orig->xyz2uvw(xyz, uvwP);
+  u = uvwP[0]; v = uvwP[1]; w = uvwP[2];
+}
 
 bool MSubPoint::isInside(double u, double v, double w)
 {
@@ -248,20 +749,23 @@ bool MSubPoint::isInside(double u, double v, double w)
 void MSubPoint::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 {
   // invariable regardless of the order
-  if(!_intpt)
+  if(!_pts)
   {
     if(!_orig) return;
-    _intpt = new IntPt[1];
+
+    _pts = new IntPt[1];
     // work in the parametric space of the parent element
     MVertex *vi = getVertex(0);
     double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
     double v_uvw[3];
     _orig->xyz2uvw(v_xyz, v_uvw);
-    _intpt[0].pt[0] = v_uvw[0];
-    _intpt[0].pt[1] = v_uvw[1];
-    _intpt[0].pt[2] = v_uvw[2];
-    _intpt[0].weight = 1;
+    double jac[3][3];
+    double origJac = _orig->getJacobian(v_uvw[0], v_uvw[1], v_uvw[2], jac);
+    _pts[0].pt[0] = v_uvw[0];
+    _pts[0].pt[1] = v_uvw[1];
+    _pts[0].pt[2] = v_uvw[2];
+    _pts[0].weight = 1./origJac;
   }
   *npts = 1;
-  *pts = _intpt;
+  *pts = _pts;
 }
diff --git a/Geo/MSubElement.h b/Geo/MSubElement.h
index 37d6c56fa0d87bc3dd0b5b4d1bcd000783e2c6ba..53f191b8fe02af25d1277bd8e1d5abcb544432e9 100644
--- a/Geo/MSubElement.h
+++ b/Geo/MSubElement.h
@@ -23,22 +23,43 @@ class MSubTetrahedron : public MTetrahedron
   bool _owner;
   MElement* _orig;
   std::vector<MElement*> _parents;
-  IntPt *_intpt;
+  mutable MElement* _base;
+
+  int _pOrder;
+  int _npts;
+  IntPt *_pts;
+
  public:
-  MSubTetrahedron(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, int num=0,
-                  int part=0, bool owner=false, MElement* orig=NULL)
-    : MTetrahedron(v0, v1, v2, v3, num, part), _owner(owner), _orig(orig), _intpt(0) {}
-  MSubTetrahedron(std::vector<MVertex*> v, int num=0, int part=0, bool owner=false,
-                  MElement* orig=NULL)
-    : MTetrahedron(v, num, part), _owner(owner), _orig(orig), _intpt(0) {}
+  MSubTetrahedron(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, int num=0, int part=0,
+                  bool owner=false, MElement* orig=NULL)
+    : MTetrahedron(v0, v1, v2, v3, num, part), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
+  MSubTetrahedron(std::vector<MVertex*> v, int num=0, int part=0,
+                  bool owner=false, MElement* orig=NULL)
+    : MTetrahedron(v, num, part), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
   MSubTetrahedron(const MTetrahedron &tet, bool owner=false, MElement* orig=NULL)
-    : MTetrahedron(tet), _owner(owner), _orig(orig), _intpt(0) {}
-
+    : MTetrahedron(tet), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
   ~MSubTetrahedron();
+
   virtual int getTypeForMSH() const { return MSH_TET_SUB; }
+  virtual const nodalBasis* getFunctionSpace(int order=-1) 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);
+  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 double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
+  virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]);
+  virtual double getJacobian(double u, double v, double w, double jac[3][3]);
+  virtual double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
+  virtual int getNumShapeFunctions();
+  virtual int getNumPrimaryShapeFunctions();
+  virtual MVertex* getShapeFunctionNode(int i);
+  virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w);
+  virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w);
   virtual bool isInside(double u, double v, double w);
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
+
   virtual MElement *getParent() const { return _orig; }
   virtual void setParent(MElement *p, bool owner = false) { _orig = p; _owner = owner; }
   virtual bool ownsParent() const { return _owner; }
@@ -47,6 +68,9 @@ class MSubTetrahedron : public MTetrahedron
   {
     _parents = parents; _orig = _parents[0]; _owner = owner;
   }
+
+  virtual const MElement *getBaseElement() const {if(!_base) _base=new MTetrahedron(*this); return _base;}
+  virtual MElement *getBaseElement() {if(!_base) _base=new MTetrahedron(*this); return _base;}
 };
 
 // A sub triangle, contained in another element
@@ -56,21 +80,43 @@ class MSubTriangle : public MTriangle
   bool _owner;
   MElement* _orig;
   std::vector<MElement*> _parents;
-  IntPt *_intpt;
+  mutable MElement* _base;
+
+  int _pOrder;
+  int _npts;
+  IntPt *_pts;
+
  public:
   MSubTriangle(MVertex *v0, MVertex *v1, MVertex *v2, int num=0, int part=0,
                bool owner=false, MElement* orig=NULL)
-    : MTriangle(v0, v1, v2, num, part), _owner(owner), _orig(orig), _intpt(0) {}
-  MSubTriangle(std::vector<MVertex*> v, int num=0, int part=0, bool owner=false,
-               MElement* orig=NULL)
-    : MTriangle(v, num, part), _owner(owner), _orig(orig), _intpt(0) {}
+    : MTriangle(v0, v1, v2, num, part), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
+  MSubTriangle(std::vector<MVertex*> v, int num=0, int part=0,
+               bool owner=false, MElement* orig=NULL)
+    : MTriangle(v, num, part), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
   MSubTriangle(const MTriangle &tri, bool owner=false, MElement* orig=NULL)
-    : MTriangle(tri), _owner(owner), _orig(orig), _intpt(0) {}
+    : MTriangle(tri), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
   ~MSubTriangle();
+
   virtual int getTypeForMSH() const { return MSH_TRI_SUB; }
+  virtual const nodalBasis* getFunctionSpace(int order=-1) 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);
+  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 double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
+  virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]);
+  virtual double getJacobian(double u, double v, double w, double jac[3][3]);
+  virtual double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
+  virtual int getNumShapeFunctions();
+  virtual int getNumPrimaryShapeFunctions();
+  virtual MVertex* getShapeFunctionNode(int i);
+  virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w);
+  virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w);
   virtual bool isInside(double u, double v, double w);
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
+
   virtual MElement *getParent() const { return _orig; }
   virtual void setParent(MElement *p, bool owner = false) { _orig = p; _owner = owner; }
   virtual bool ownsParent() const { return _owner; }
@@ -79,6 +125,9 @@ class MSubTriangle : public MTriangle
   {
     _parents = parents; _orig = _parents[0]; _owner = owner;
   }
+
+  virtual const MElement *getBaseElement() const {if(!_base) _base=new MTriangle(*this); return _base;}
+  virtual MElement *getBaseElement() {if(!_base) _base=new MTriangle(*this); return _base;}
 };
 
 // A sub line, contained in another element
@@ -88,22 +137,43 @@ class MSubLine : public MLine
   bool _owner;
   MElement* _orig;
   std::vector<MElement*> _parents;
-  IntPt *_intpt;
+  mutable MElement* _base;
+
+  int _pOrder;
+  int _npts;
+  IntPt *_pts;
+
  public:
-  MSubLine(MVertex *v0, MVertex *v1, int num=0, int part=0, bool owner=false,
-           MElement* orig=NULL)
-    : MLine(v0, v1, num, part), _owner(owner), _orig(orig), _intpt(0) {}
-  MSubLine(std::vector<MVertex*> v, int num=0, int part=0, bool owner=false,
-           MElement* orig=NULL)
-    : MLine(v, num, part), _owner(owner), _orig(orig), _intpt(0) {}
+  MSubLine(MVertex *v0, MVertex *v1, int num=0, int part=0,
+           bool owner=false, MElement* orig=NULL)
+    : MLine(v0, v1, num, part), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
+  MSubLine(std::vector<MVertex*> v, int num=0, int part=0,
+           bool owner=false, MElement* orig=NULL)
+    : MLine(v, num, part), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
   MSubLine(const MLine &lin, bool owner=false, MElement* orig=NULL)
-    : MLine(lin), _owner(owner), _orig(orig), _intpt(0) {}
-
+    : MLine(lin), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
   ~MSubLine();
+
   virtual int getTypeForMSH() const { return MSH_LIN_SUB; }
+  virtual const nodalBasis* getFunctionSpace(int order=-1) 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);
+  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 double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
+  virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]);
+  virtual double getJacobian(double u, double v, double w, double jac[3][3]);
+  virtual double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
+  virtual int getNumShapeFunctions();
+  virtual int getNumPrimaryShapeFunctions();
+  virtual MVertex* getShapeFunctionNode(int i);
+  virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w);
+  virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w);
   virtual bool isInside(double u, double v, double w);
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
+
   virtual MElement *getParent() const { return _orig; }
   virtual void setParent(MElement *p, bool owner = false) { _orig = p; _owner = owner; }
   virtual bool ownsParent() const { return _owner; }
@@ -112,6 +182,9 @@ class MSubLine : public MLine
   {
     _parents = parents; _orig = _parents[0]; _owner = owner;
   }
+
+  virtual const MElement *getBaseElement() const {if(!_base) _base=new MLine(*this); return _base;}
+  virtual MElement *getBaseElement() {if(!_base) _base=new MLine(*this); return _base;}
 };
 
 // A sub point, contained in another element
@@ -121,21 +194,43 @@ class MSubPoint : public MPoint
   bool _owner;
   MElement* _orig;
   std::vector<MElement*> _parents;
-  IntPt *_intpt;
+  mutable MElement* _base;
+
+  int _pOrder;
+  int _npts;
+  IntPt *_pts;
+
  public:
-  MSubPoint(MVertex *v0, int num=0, int part=0, bool owner=false, MElement* orig=NULL)
-    : MPoint(v0, num, part), _owner(owner), _orig(orig), _intpt(0) {}
-  MSubPoint(std::vector<MVertex*> v, int num=0, int part=0, bool owner=false,
-            MElement* orig=NULL)
-    : MPoint(v, num, part), _owner(owner), _orig(orig), _intpt(0) {}
+  MSubPoint(MVertex *v0, int num=0, int part=0,
+            bool owner=false, MElement* orig=NULL)
+    : MPoint(v0, num, part), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
+  MSubPoint(std::vector<MVertex*> v, int num=0, int part=0,
+            bool owner=false, MElement* orig=NULL)
+    : MPoint(v, num, part), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
   MSubPoint(const MPoint &pt, bool owner=false, MElement* orig=NULL)
-    : MPoint(pt), _owner(owner), _orig(orig), _intpt(0) {}
-
+    : MPoint(pt), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
   ~MSubPoint();
+
   virtual int getTypeForMSH() const { return MSH_PNT_SUB; }
+  virtual const nodalBasis* getFunctionSpace(int order=-1) 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);
+  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 double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
+  virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]);
+  virtual double getJacobian(double u, double v, double w, double jac[3][3]);
+  virtual double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
+  virtual int getNumShapeFunctions();
+  virtual int getNumPrimaryShapeFunctions();
+  virtual MVertex* getShapeFunctionNode(int i);
+  virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w);
+  virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w);
   virtual bool isInside(double u, double v, double w);
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
+  
   virtual MElement *getParent() const { return _orig; }
   virtual void setParent(MElement *p, bool owner = false) { _orig = p; _owner = owner; }
   virtual bool ownsParent() const { return _owner; }
@@ -144,6 +239,9 @@ class MSubPoint : public MPoint
   {
     _parents = parents; _orig = _parents[0]; _owner = owner;
   }
+
+  virtual const MElement *getBaseElement() const {if(!_base) _base=new MPoint(*this); return _base;}
+  virtual MElement *getBaseElement() {if(!_base) _base=new MPoint(*this); return _base;}
 };
 
 #endif