diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.cpp b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
index 4c7dd27eb206824360caa0e5a425d7b26a93df68..aeeb08df770a0dd17f275208731729e160799084 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
@@ -132,7 +132,6 @@ void OptHOM::printProgress(const alglib::real_1d_array &x, double Obj)
   if (iter % progressInterv == 0) {
     double maxD, avgD, minJ, maxJ;
     getDistances(maxD, avgD, minJ, maxJ);
-
     printf("--> Iteration %3d --- OBJ %12.5E (relative decrease = %12.5E) -- minJ = %12.5E  maxJ = %12.5E Max D = %12.5E Avg D = %12.5E\n", iter, Obj, Obj/initObj, minJ, maxJ, maxD, avgD);
   }
 
@@ -151,15 +150,8 @@ void OptHOM::OptimPass(alglib::real_1d_array &x, const alglib::real_1d_array &in
 {
 
   static const double EPSG = 0.;
-  static const double EPSF = 1.e-8;
-//  static const double EPSF = 0.;
+  static const double EPSF = 0.;
   static const double EPSX = 0.;
-//  const double EPSX = x.length()*1.e-4/sqrt(invLengthScaleSq);
-//  std::cout << "DEBUG: EPSX = " << EPSX << ", EPSX/x.length() = " << EPSX/x.length() << std::endl;
-
-//  double iGONorm = 0;
-//  for (int i=0; i<initGradObj.length(); i++) iGONorm += initGradObj[i]*initGradObj[i];
-//  const double EPSG = 1.e-2*sqrt(iGONorm);
 
   std::cout << "--- Optimization pass with jacBar = " << jacBar << ", lambda = " << lambda << ", lambda2 = " << lambda2 << std::endl;
 
@@ -187,10 +179,10 @@ void OptHOM::OptimPass(alglib::real_1d_array &x, const alglib::real_1d_array &in
   switch(int(rep.terminationtype)) {
 //  case -2: std::cout << ", because rounding errors prevented further improvement"; break;
 //  case -1: std::cout << ", because incorrect parameters were specified"; break;
-//  case 1: std::cout << ", because relative function improvement is no more than EpsF"; break;
-//  case 2: std::cout << ", because relative step is no more than EpsX"; break;
-//  case 4: std::cout << ", because gradient norm is no more than EpsG"; break;
-//  case 5: std::cout << ", because the maximum number of steps was taken"; break;
+  case 1: std::cout << ", because relative function improvement is no more than EpsF"; break;
+  case 2: std::cout << ", because relative step is no more than EpsX"; break;
+  case 4: std::cout << ", because gradient norm is no more than EpsG"; break;
+  case 5: std::cout << ", because the maximum number of steps was taken"; break;
 //  case 7: std::cout << ", because stopping conditions are too stringent, further improvement is impossible"; break;
   default: std::cout << " with code " << int(rep.terminationtype); break;
   }
@@ -227,7 +219,6 @@ int OptHOM::optimize(double weightFixed, double weightFree, double barrier_, int
   const double jacBarStart = (minJ > 0.) ? 0.9*minJ : 1.1*minJ;
   jacBar = jacBarStart;
   setBarrierTerm(jacBarStart);
-  std::cout << "DEBUG: jacBarStart = " << jacBarStart << std::endl;
 
   // Calculate initial objective function value and gradient
   initObj = 0.;
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index a4249443864900a6c8daf6a7384bac04613b7754..58885217fc0fb10b973263764a7948c329cb7904 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -4,11 +4,14 @@
 #include "MTetrahedron.h"
 #include "ParamCoord.h"
 #include "OptHomMesh.h"
-#include "JacobianBasis.h"
+
+
 
 std::map<int, std::vector<double> > Mesh::_jacBez;
 
-static std::vector<double> _computeJB(const polynomialBasis *lagrange, const bezierBasis *bezier)
+
+
+std::vector<double> Mesh::computeJB(const polynomialBasis *lagrange, const bezierBasis *bezier)
 {
   int nbNodes = lagrange->points.size1();
   
@@ -36,7 +39,7 @@ static std::vector<double> _computeJB(const polynomialBasis *lagrange, const bez
         for (int j = 0; j < nbNodes; j++) {
           double Jij = dPsi(i, 0) * dPsi(j, 1) - dPsi(i, 1) * dPsi(j,0);
           for (int l = 0; l < bezier->points.size1(); l++) {
-            JB[(l * nbNodes + i) * nbNodes + j] += bezier->matrixLag2Bez(l, k) * Jij;
+            JB[indJB2DBase(nbNodes,l,i,j)] += bezier->matrixLag2Bez(l, k) * Jij;
           }
         }
       }
@@ -50,7 +53,7 @@ static std::vector<double> _computeJB(const polynomialBasis *lagrange, const bez
               + (dPsi(j, 2) * dPsi(m, 0) - dPsi(j, 0) * dPsi(m, 2)) * dPsi(i, 1)
               + (dPsi(j, 0) * dPsi(m, 1) - dPsi(j, 1) * dPsi(m, 0)) * dPsi(i, 2);
             for (int l = 0; l < bezier->points.size1(); l++) {
-              JB[(l * nbNodes + (i * nbNodes + j)) * nbNodes + m] += bezier->matrixLag2Bez(l, k) * Jijm;
+              JB[indJB3DBase(nbNodes,l,i,j,m)] += bezier->matrixLag2Bez(l, k) * Jijm;
             }
           }
         }
@@ -60,6 +63,8 @@ static std::vector<double> _computeJB(const polynomialBasis *lagrange, const bez
   return JB;
 }
 
+
+
 Mesh::Mesh(GEntity *ge, std::set<MVertex*> &toFix, int method) :
     _ge(ge)
 {
@@ -100,7 +105,7 @@ Mesh::Mesh(GEntity *ge, std::set<MVertex*> &toFix, int method) :
     const polynomialBasis *lagrange = el->getFunctionSpace();
     const bezierBasis *bezier = JacobianBasis::find(lagrange->type)->bezier;
     if (_jacBez.find(lagrange->type) == _jacBez.end()) {
-      _jacBez[lagrange->type] = _computeJB(lagrange, bezier);
+      _jacBez[lagrange->type] = computeJB(lagrange, bezier);
     }
     _nBezEl[iEl] = bezier->points.size1();
     _nNodEl[iEl] = lagrange->points.size1();
@@ -152,6 +157,8 @@ Mesh::Mesh(GEntity *ge, std::set<MVertex*> &toFix, int method) :
   }
 }
 
+
+
 SVector3 Mesh::getNormalEl(int iEl)
 {
 
@@ -324,7 +331,7 @@ void Mesh::scaledJac(int iEl, std::vector<double> &sJ)
         }
         sJ[l] *= n.z();
       }
- }
+  }
   else {
     for (int l = 0; l < _nBezEl[iEl]; l++) {
       sJ[l] = 0.;
@@ -421,8 +428,8 @@ void Mesh::gradScaledJac(int iEl, std::vector<double> &gSJ)
         _pc->gXyz2gUvw(_freeVert[iFVi],_uvw[iFVi],gXyzV,gUvwV);
         for (int l = 0; l < _nBezEl[iEl]; l++) {
           gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][0];
-          if (_nPCFV[iFVi] >= 2) gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][1];
-          if (_nPCFV[iFVi] == 3) gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][2];
+          if (_nPCFV[iFVi] >= 2) gSJ[indGSJ(iEl,l,iPC+1)] = gUvwV[l][1];
+          if (_nPCFV[iFVi] == 3) gSJ[indGSJ(iEl,l,iPC+2)] = gUvwV[l][2];
         }
         iPC += _nPCFV[iFVi];
       }
@@ -431,6 +438,8 @@ void Mesh::gradScaledJac(int iEl, std::vector<double> &gSJ)
 
 }
 
+
+
 void Mesh::writeMSH(const char *filename)
 {
   FILE *f = fopen(filename, "w");
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.h b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
index 6a9067d98ed956fe63e3d053f3265f1fd870abf1..20329df0381cca3a8f4e2a4c95914dae3873345f 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
@@ -7,6 +7,10 @@
 #include "GFace.h"
 #include "MVertex.h"
 #include "ParamCoord.h"
+#include "polynomialBasis.h"
+#include "JacobianBasis.h"
+
+
 
 class Mesh
 {
@@ -78,8 +82,11 @@ private:
   int addVert(MVertex* vert);
   int addFreeVert(MVertex* vert, const int iV, const int nPCV, std::set<MVertex*> &toFix);
   SVector3 getNormalEl(int iEl);
-  inline int indJB2D(int iEl, int l, int i, int j) { return (l*_nNodEl[iEl]+i)*_nNodEl[iEl]+j; }
-  inline int indJB3D(int iEl, int l, int i, int j, int m) { return ((l*_nNodEl[iEl]+i)*_nNodEl[iEl]+j)*_nNodEl[iEl]+m; }
+  static std::vector<double> computeJB(const polynomialBasis *lagrange, const bezierBasis *bezier);
+  static inline int indJB2DBase(int nNod, int l, int i, int j) { return (l*nNod+i)*nNod+j; }
+  inline int indJB2D(int iEl, int l, int i, int j) { return indJB2DBase(_nNodEl[iEl],l,i,j); }
+  static inline int indJB3DBase(int nNod, int l, int i, int j, int m) { return ((l*nNod+i)*nNod+j)*nNod+m; }
+  inline int indJB3D(int iEl, int l, int i, int j, int m) { return indJB3DBase(_nNodEl[iEl],l,i,j,m); }
 
 };
 
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index ef79978e23cc2666a5e340e6766937a04919a367..444544e99e38c8379218f43908792ee93e86f601 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -129,7 +129,7 @@ std::set<MVertex*> filter3D(GRegion *gr, int nbLayers, double _qual)
       for (int j = 0; j < t->getNumVertices(); j++)
         touched.insert(t->getVertex(j));
     }
-    if (ts.size() == 1) break;
+//    if (ts.size() == 1) break;
   }
 
   //  printf("FILTER3D : %d bad tets found among %6d\n", ts.size(), gr->tetrahedra.size());