diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 854874116e948944535f8ab94f462342d4fbb2af..40c2ca585a3e5f005f2b0fc760bf8d0efc37ec90 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1334,7 +1334,8 @@ void deMeshGFace::operator() (GFace *gf)
   gf->meshStatistics.nbTriangle = gf->meshStatistics.nbEdge = 0;
 }
 
-int debugSurface = -100;
+//For Debugging, change value from -1 to -100;
+int debugSurface = -1; 
 
 void meshGFace::operator() (GFace *gf)
 {
diff --git a/Plugin/Distance.cpp b/Plugin/Distance.cpp
index 80b0df84e3e5ffd0ba0bceb39e359e5ae561aa2e..51a4108dd9f53e40e185fb328fedd67b6d1ebda6 100644
--- a/Plugin/Distance.cpp
+++ b/Plugin/Distance.cpp
@@ -53,7 +53,7 @@ std::string GMSH_DistancePlugin::getHelp() const
     
     "Define the elementary entities to which the distance is computed. If Point=0, Line=0, and Surface=0, then the distance is computed to all the boundaries of the mesh (edges in 2D and faces in 3D)\n\n"
 
-  "Computation<0. computes the geometrical euclidian distance (warning: different than the geodesic distance), and  Computation=a>0.0 solves a PDE on the mesh with the diffusion constant mu = bbox/a, with bbox being the max size of the bounding box of the mesh (see paper Legrand 2006) \n\n"
+  "Computation<0. computes the geometrical euclidian distance (warning: different than the geodesic distance), and  Computation=a>0.0 solves a PDE on the mesh with the diffusion constant mu = a*bbox, with bbox being the max size of the bounding box of the mesh (see paper Legrand 2006) \n\n"
 
   "Plugin(Distance) creates a new distance view and also saves the view in the fileName.pos file.";
 }
@@ -258,10 +258,10 @@ PView *GMSH_DistancePlugin::execute(PView *v)
     }
 
     double L = norm(SVector3(bbox.max(), bbox.min())); 
-    double mu = L/type;
+    double mu = type*L;
 
-    simpleFunction<double> DIFF(mu * mu), MONE(1.0);
-    distanceTerm distance(GModel::current(), 1, &DIFF, &MONE);
+    simpleFunction<double> DIFF(mu*mu), ONE(1.0);
+    distanceTerm distance(GModel::current(), 1, &DIFF, &ONE);
   
     for(unsigned int ii = 0; ii < entities.size(); ii++){
       if(entities[ii]->dim() == maxDim) {
@@ -331,8 +331,54 @@ PView *GMSH_DistancePlugin::execute(PView *v)
 #endif
   }
 
-  //FILE *fp = fopen(fileName.c_str(), "rb");
-  //data->readPOS(fp, 1.0, false);
+  //compute also orthogonal vector to distance field
+  //------------------------------------------------
+// #if defined(HAVE_SOLVER)
+  
+// #ifdef HAVE_TAUCS
+//     linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
+// #else
+//     linearSystemCSRGmm<double> *lsys = new linearSystemCSRGmm<double>;
+//     lsys->setNoisy(1);
+//     lsys->setGmres(1);
+//     lsys->setPrec(5.e-8);
+// #endif
+
+//     dofManager<double> myAssembler(lsys);
+
+//     for(unsigned int ii = 0; ii < entities.size(); ii++){
+//       if(entities[ii]->dim() == maxDim) {
+// 	GEntity *ge = entities[ii];
+// 	for(unsigned int i = 0; i < ge->getNumMeshElements(); ++i){
+// 	  MElement *t = ge->getMeshElement(i);
+// 	  for(int k = 0; k < t->getNumVertices(); k++){
+// 	    myAssembler.numberVertex(t->getVertex(k), 0, 1);
+// 	  }
+// 	}    
+//       }  
+//     }
+
+//     simpleFunction<double> ONE(1.0);
+//     simpleFunction<double> MONE(-1.0 );
+//     laplaceTerm laplace(model(), 1, &ONE);
+//     crossConfTerm cross12(model(), 1, 1, &ONE);
+
+//     for(unsigned int ii = 0; ii < entities.size(); ii++){
+//       if(entities[ii]->dim() == maxDim) {
+// 	GEntity *ge = entities[ii];
+// 	for(unsigned int i = 0; i < ge->getNumMeshElements(); ++i){
+// 	  SElement se(ge->getMeshElement(i));
+// 	  laplace.addToMatrix(myAssembler, &se);
+// 	}
+// 	//groupOfElements g((GFace*)ge); //WARNING GFACE HE
+// 	//distance.addToRightHandSide(myAssembler, g);
+//       }
+//     }
+
+// #endif
+
+    //-------------------------------------------------
+
 
   data->Time.push_back(0);
   data->setFileName(fileName.c_str());
diff --git a/Solver/distanceTerm.h b/Solver/distanceTerm.h
index 3ab7f6b0dbf9ce9484a63cbea99300cacea3017d..dcea0b969ba50d4cd31d503aa1334a51f1d7f3c7 100644
--- a/Solver/distanceTerm.h
+++ b/Solver/distanceTerm.h
@@ -23,9 +23,7 @@ class distanceTerm : public helmholtzTerm<double> {
     double jac[3][3];
     double ff[256];
     e->getIntegrationPoints(integrationOrder, &npts, &GP);
-  
-    m.scale(0.);
-    
+    m.scale(0.); 
     for (int i = 0; i < npts; i++){
       const double u = GP[i].pt[0];
       const double v = GP[i].pt[1];
diff --git a/Solver/helmholtzTerm.h b/Solver/helmholtzTerm.h
index 64b6d1890b7764dc9402d7c7b17f9b741f1f0e6a..2af2587d623e7fcef8c1294d572577bcc0e9a9e4 100644
--- a/Solver/helmholtzTerm.h
+++ b/Solver/helmholtzTerm.h
@@ -15,7 +15,7 @@
 #include "fullMatrix.h"
 #include "Numeric.h"
 
-// \nabla \cdot k \nabla U - a U 
+// \nabla \cdot k \nabla U - a U = 0 
 template<class scalar>
 class helmholtzTerm : public femTerm<scalar> {
  protected: