diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 6e7ffa83fd629995dda0323811748dc7a42de53f..3fda9d52e4d72be2aa4dca8998a1eae1e0ec4ee3 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -46,7 +46,7 @@ GFace::~GFace()
 }
 
 void GFace::delFreeEdge(GEdge *e)
-{ 
+{
   // delete the edge from the edge list and the orientation list
   std::list<GEdge*>::iterator ite = l_edges.begin();
   std::list<int>::iterator itd = l_dirs.begin();
@@ -57,12 +57,12 @@ void GFace::delFreeEdge(GEdge *e)
       if(itd != l_dirs.end()) l_dirs.erase(itd);
       break;
     }
-    ite++; 
+    ite++;
     if(itd != l_dirs.end()) itd++;
   }
 
   // delete the edge from the edge loops
-  for(std::list<GEdgeLoop>::iterator it = edgeLoops.begin(); 
+  for(std::list<GEdgeLoop>::iterator it = edgeLoops.begin();
       it != edgeLoops.end(); it++){
     for(GEdgeLoop::iter it2 = it->begin(); it2 != it->end(); it2++){
       if(e == it2->ge){
@@ -90,8 +90,8 @@ void GFace::deleteMesh()
 }
 
 unsigned int GFace::getNumMeshElements()
-{ 
-  return triangles.size() + quadrangles.size() + polygons.size(); 
+{
+  return triangles.size() + quadrangles.size() + polygons.size();
 }
 
 void GFace::getNumMeshElements(unsigned *const c) const
@@ -118,7 +118,7 @@ MElement *const *GFace::getStartElementType(int type) const
 }
 
 MElement *GFace::getMeshElement(unsigned int index)
-{ 
+{
   if(index < triangles.size())
     return triangles[index];
   else if(index < triangles.size() + quadrangles.size())
@@ -180,7 +180,7 @@ SOrientedBoundingBox GFace::getOBB()
         vertices.push_back(pt1);
         vertices.push_back(pt2);
       }
-    } 
+    }
     else if(buildSTLTriangulation()) {
       for (unsigned int i = 0; i < stl_vertices.size(); i++){
         GPoint p = point(stl_vertices[i]);
@@ -293,11 +293,11 @@ void GFace::writeGEO(FILE *fp)
     }
   }
 
-  for(std::list<GEdge*>::iterator it = embedded_edges.begin(); 
+  for(std::list<GEdge*>::iterator it = embedded_edges.begin();
       it != embedded_edges.end(); it++)
     fprintf(fp, "Line {%d} In Surface {%d};\n", (*it)->tag(), tag());
 
-  for(std::list<GVertex*>::iterator it = embedded_vertices.begin(); 
+  for(std::list<GVertex*>::iterator it = embedded_vertices.begin();
       it != embedded_vertices.end(); it++)
     fprintf(fp, "Point {%d} In Surface {%d};\n", (*it)->tag(), tag());
 
@@ -661,7 +661,7 @@ double GFace::curvatureDiv(const SPoint2 &param) const
 
   double ddu = dot(dndu, du);
   double ddv = dot(dndv, dv);
-  
+
   return (fabs(ddu) + fabs(ddv)) / detJ;
 }
 
@@ -687,7 +687,7 @@ double GFace::curvatures(const SPoint2 &param, SVector3 *dirMax, SVector3 *dirMi
     *curvMin = 0.;
     return 0.;
   }
-  
+
   if(geomType() == Sphere){
     *dirMax = D1.first();
     *dirMin = D1.second();
@@ -704,7 +704,7 @@ double GFace::curvatures(const SPoint2 &param, SVector3 *dirMax, SVector3 *dirMi
   *curvMin = fabs(eigVal[0]);
   *dirMax = eigVec[1] * D1.first() + eigVec[3] * D1.second();
   *dirMin = eigVec[0] * D1.first() + eigVec[2] * D1.second();
-  
+
   return *curvMax;
 }
 
@@ -756,7 +756,7 @@ void GFace::getMetricEigenVectors(const SPoint2 &param,
   N(0, 1) = inv_form1[0][0] * form2[0][1] + inv_form1[0][1] * form2[1][1];
   N(1, 0) = inv_form1[1][0] * form2[0][0] + inv_form1[1][1] * form2[1][0];
   N(1, 1) = inv_form1[1][0] * form2[0][1] + inv_form1[1][1] * form2[1][1];
-  
+
   // eigen values and vectors of N
   fullMatrix<double> vl(2, 2), vr(2, 2);
   fullVector<double> dr(2), di(2);
@@ -1013,7 +1013,7 @@ bool GFace::buildSTLTriangulation(bool force)
     }
     return true;
   }
-  
+
   // build STL for general surfaces here
 
   return false;
@@ -1096,19 +1096,19 @@ double GFace::length(const SPoint2 &pt1, const SPoint2 &pt2, int nbQuadPoints)
   return L;
 }
 
-int GFace::poincareMesh()  
+int GFace::poincareMesh()
 {
   std::set<MEdge, Less_Edge> es;
   std::set<MVertex*> vs;
-  for(unsigned int i = 0; i < getNumMeshElements(); i++){ 
+  for(unsigned int i = 0; i < getNumMeshElements(); i++){
     MElement *e = getMeshElement(i);
     for(int j = 0; j < e->getNumVertices(); j++) vs.insert(e->getVertex(j));
     for(int j = 0; j < e->getNumEdges(); j++) es.insert(e->getEdge(j));
   }
-  return vs.size() - es.size() + getNumMeshElements();  
+  return vs.size() - es.size() + getNumMeshElements();
 }
 
-int GFace::genusGeom()  
+int GFace::genusGeom()
 {
   int nSeams = 0;
   std::set<GEdge*> single_seams;
@@ -1131,7 +1131,7 @@ bool GFace::fillPointCloud(double maxDist, std::vector<SPoint3> *points,
     Msg::Error("No STL triangulation available to fill point cloud");
     return false;
   }
-  
+
   if(!points) return false;
 
   for(unsigned int i = 0; i < stl_triangles.size(); i += 3){
@@ -1149,7 +1149,7 @@ bool GFace::fillPointCloud(double maxDist, std::vector<SPoint3> *points,
       for(double v = 0.; v < 1 - u; v += 1. / N){
         SPoint2 p = p0 * (1. - u - v) + p1 * u + p2 * v;
         GPoint gp(point(p));
-        points->push_back(SPoint3(gp.x(), gp.y(), gp.z())); 
+        points->push_back(SPoint3(gp.x(), gp.y(), gp.z()));
         if(normals) normals->push_back(normal(p));
       }
     }
@@ -1174,7 +1174,7 @@ void GFace::replaceEdges (std::list<GEdge*> &new_edges)
   std::list<int> newdirs;
   for ( ; it != l_edges.end(); ++it, ++it2, ++it3){
     (*it)->delFace(this);
-    (*it2)->addFace(this);        
+    (*it2)->addFace(this);
     if ((*it2)->getBeginVertex() == (*it)->getBeginVertex())
       newdirs.push_back(*it3);
     else
@@ -1216,8 +1216,8 @@ void GFace::addLayersOfQuads(int nLayers, GVertex *gv, double hmin, double ratio
     bool found = false;
     // look if this edge loop has the GVertex as an endpoint
     for (GEdgeLoop::iter it2 = it->begin(); it2 != it->end(); ++it2){
-      if (it2->ge->getBeginVertex() == gv || it2->ge->getEndVertex() == gv) 
-	found = true;	
+      if (it2->ge->getBeginVertex() == gv || it2->ge->getEndVertex() == gv)
+	found = true;
     }
     // we found an edge loop with the GVertex that was specified
     if (found){
diff --git a/Geo/SPoint3.h b/Geo/SPoint3.h
index 4f4ad376acc7d7567e350662557c214ea848c205..c4cf69bc28dd20f52d2219e1b8345c6c82d48f2d 100644
--- a/Geo/SPoint3.h
+++ b/Geo/SPoint3.h
@@ -7,7 +7,6 @@
 #define _SPOINT3_H_
 
 #include <math.h>
-
 // A point in 3-space
 class SPoint3 {
  protected:
@@ -17,6 +16,7 @@ class SPoint3 {
   SPoint3(double x, double y, double z) { P[0] = x; P[1] = y; P[2] = z; }
   SPoint3(const double *p) { P[0] = p[0]; P[1] = p[1]; P[2] = p[2]; }
   SPoint3(const SPoint3 &pt) { P[0] = pt.P[0]; P[1] = pt.P[1]; P[2] = pt.P[2]; }
+  SPoint3(const SPoint3 &pt,const SPoint3 &dir,const double dist_) {P[0]=pt.P[0]; P[1]=pt.P[1]; P[2]=pt.P[2]; SPoint3 a(dir); a*=dist_; P[0]+=a[0]; P[1]+=a[1]; P[2]+=a[2];}
   virtual ~SPoint3() {}
   void setPosition(double xx, double yy, double zz);
   void getPosition(double *xx, double *yy, double *zz) const;
@@ -95,7 +95,7 @@ inline double SPoint3::operator[](int i) const
 { return P[i]; }
 
 inline double SPoint3::distance(const SPoint3 &p)const
-{ 
+{
   double x = P[0] - p.P[0], y = P[1] - p.P[1], z = P[2] - p.P[2];
   return sqrt(x * x + y * y + z * z);
 }
diff --git a/Geo/SVector3.h b/Geo/SVector3.h
index 842229e455665bcc244fcc59e49ed87c9967d8c4..16b72faf3fa6c3a435e16c0885dccab16affd67b 100644
--- a/Geo/SVector3.h
+++ b/Geo/SVector3.h
@@ -68,6 +68,9 @@ class SVector3 {
   operator double *() { return P; }
   void print(std::string name="") const
   { printf("Vector \'%s\':  %f  %f  %f\n",name.c_str(),P[0],P[1],P[2]); }
+
+  // Needed to allow the initialization of a SPoint3 from a SPoint3, a distance and a direction
+  SPoint3 point() const{return P;}
 };
 
 inline double dot(const SVector3 &a, const SVector3 &b)