diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index 97be674aee50d1832034b349a2b9f80f608f25d7..d62dd21e44800d3c96331f6cf62b6ec9eb255447 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -314,6 +314,46 @@ gLevelset::gLevelset(const gLevelset &lv)
   tag_ = lv.tag_;
 }
 
+gLevelsetSphere::gLevelsetSphere (const double &x, const double &y, const double &z, const double &R, int tag=1)
+  : gLevelsetPrimitive(tag), xc(x), yc(y), zc(z), r(R)
+{
+_hasDerivatives = true;
+}
+
+void gLevelsetSphere::gradient (double x, double y, double z,
+           double & dfdx, double & dfdy, double & dfdz) const
+{
+
+  const double xx = x-xc, yy = y-yc, zz = z-zc, dist = sqrt(xx*xx+yy*yy+zz*zz);
+
+  dfdx = xx/dist;
+  dfdy = yy/dist;
+  dfdz = zz/dist;
+
+}
+
+void gLevelsetSphere::hessian (double x, double y, double z,
+          double & dfdxx, double & dfdxy, double & dfdxz,
+          double & dfdyx, double & dfdyy, double & dfdyz,
+          double & dfdzx, double & dfdzy, double & dfdzz) const
+{
+
+  const double xx = x-xc, yy = y-yc, zz = z-zc;
+  const double distSq = xx*xx+yy*yy, fact = 1./(distSq*sqrt(distSq));
+
+  dfdxx = (zz*zz+yy*yy)*fact;
+  dfdxy = -xx*yy*fact;
+  dfdxz = -xx*zz*fact;
+  dfdyx = dfdxy;
+  dfdyy = (xx*xx+zz*zz)*fact;
+  dfdyz = -yy*zz*fact;
+  dfdzx = dfdxz;
+  dfdzy = dfdyz;
+  dfdzz = (xx*xx+yy*yy)*fact;
+
+}
+
+
 gLevelsetPlane::gLevelsetPlane(const std::vector<double>  &pt,
                                const std::vector<double> &norm, int tag)
   : gLevelsetPrimitive(tag)
diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h
index d40cb9c239c4fc382e8b9832adcf9a1b70ce8b64..d6716c50b151afbb80c3813578312754f3d09b1b 100644
--- a/Geo/gmshLevelset.h
+++ b/Geo/gmshLevelset.h
@@ -128,12 +128,18 @@ class gLevelsetSphere : public gLevelsetPrimitive
 protected:
   double xc, yc, zc, r;
 public:
-  gLevelsetSphere (const double &x, const double &y, const double &z, const double &R, int tag=1) : gLevelsetPrimitive(tag), xc(x), yc(y), zc(z), r(R) {}
+  gLevelsetSphere (const double &x, const double &y, const double &z, const double &R, int tag);
   virtual double operator () (const double x, const double y, const double z) const{
     if(r >= 0.)
       return sqrt((xc - x) * (xc - x) + (yc - y) * (yc - y) + (zc - z) * (zc - z)) - r;
     return (- r - sqrt((xc - x) * (xc - x) + (yc - y) * (yc - y) + (zc - z) * (zc - z)));
   }
+  void gradient (double x, double y, double z,
+    double & dfdx, double & dfdy, double & dfdz) const;
+  void hessian (double x, double y, double z,
+    double & dfdxx, double & dfdxy, double & dfdxz,
+    double & dfdyx, double & dfdyy, double & dfdyz,
+    double & dfdzx, double & dfdzy, double & dfdzz  ) const;
   int type() const {return SPHERE;}
 };