diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 2d54c47d4f4af6e5b5f2beb36576fa21de212417..ca89017d999ce28202f23ecb2b1998a87e3f3dc7 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -369,6 +369,24 @@ double MElement::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
   return _computeDeterminantAndRegularize(this, jac);
 }
 
+double MElement::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])
+{
+  jac[0][0] = jac[0][1] = jac[0][2] = 0.;
+  jac[1][0] = jac[1][1] = jac[1][2] = 0.;
+  jac[2][0] = jac[2][1] = jac[2][2] = 0.;
+
+  for (int i = 0; i < getNumShapeFunctions(); i++) {
+    const MVertex *v = getShapeFunctionNode(i);
+    for (int j = 0; j < 3; j++) {
+      double mult = gsf[i][j];
+      jac[j][0] += v->x() * mult;
+      jac[j][1] += v->y() * mult;
+      jac[j][2] += v->z() * mult;
+    }
+  }
+  return _computeDeterminantAndRegularize(this, jac);
+}
+
 double MElement::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
 {
   jac[0][0] = jac[0][1] = jac[0][2] = 0.;
@@ -404,6 +422,18 @@ void MElement::pnt(double u, double v, double w, SPoint3 &p)
   p = SPoint3(x, y, z);
 }
 
+void MElement::pnt(const std::vector<double> &sf, SPoint3 &p)
+{
+  double x = 0., y = 0., z = 0.;
+  for (int j = 0; j < getNumShapeFunctions(); j++) {
+    const MVertex *v = getShapeFunctionNode(j);
+    x += sf[j] * v->x();
+    y += sf[j] * v->y();
+    z += sf[j] * v->z();
+  }
+  p = SPoint3(x, y, z);
+}
+
 void MElement::primaryPnt(double u, double v, double w, SPoint3 &p)
 {
   double x = 0., y = 0., z = 0.;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 5d048996537be150eb00753162d7334b6eda4f79..f7b824e3b91215a9189b950680fa18b86629e4f2 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -243,6 +243,8 @@ class MElement
   // 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]);
+  // 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]);
   inline double getJacobian(double u, double v, double w, fullMatrix<double> &j){
     double JAC[3][3];
@@ -266,6 +268,7 @@ class MElement
   // get the point in cartesian coordinates corresponding to the point
   // (u,v,w) in parametric coordinates
   virtual void pnt(double u, double v, double w, SPoint3 &p);
+  virtual void pnt(const std::vector<double> &sf,SPoint3 &p); // To be compatible with functionSpace without changing form
   virtual void primaryPnt(double u, double v, double w, SPoint3 &p);
 
   // invert the parametrisation