diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 57f09d9e2385214e9d0b02a01195fb442a237b5b..718f08216d19d8cd4df9ad3e994e6fd9640e2e55 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -701,10 +701,11 @@ struct HOPatchDefParameters : public MeshOptPatchDef
 {
   HOPatchDefParameters(const OptHomParameters &p);
   virtual ~HOPatchDefParameters() {}
-  virtual double elBadness(MElement *el) const;
+  virtual double elBadness(MElement *el, GEntity* gEnt) const;
   virtual double maxDistance(MElement *el) const;
   virtual int inPatch(const SPoint3 &badBary,
-                      double limDist, MElement *el) const;
+                      double limDist, MElement *el,
+                      GEntity* gEnt) const;
 private:
   double jacMin, jacMax;
   double distanceFactor;
@@ -736,7 +737,7 @@ HOPatchDefParameters::HOPatchDefParameters(const OptHomParameters &p)
 }
 
 
-double HOPatchDefParameters::elBadness(MElement *el) const
+double HOPatchDefParameters::elBadness(MElement *el, GEntity* gEnt) const
 {
   double jmin, jmax;
   el->scaledJacRange(jmin, jmax);
@@ -756,7 +757,8 @@ double HOPatchDefParameters::maxDistance(MElement *el) const
 
 
 int HOPatchDefParameters::inPatch(const SPoint3 &badBary,
-                                  double limDist, MElement *el) const
+                                  double limDist, MElement *el,
+                                  GEntity* gEnt) const
 {
   return testElInDist(badBary, limDist, el) ? 1 : 0;
 }
@@ -770,6 +772,7 @@ void HighOrderMeshOptimizerNew(GModel *gm, OptHomParameters &p)
   par.dim = p.dim;
   par.onlyVisible = p.onlyVisible;
   par.fixBndNodes = p.fixBndNodes;
+  par.useGeom = false;
   HOPatchDefParameters patchDef(p);
   par.patchDef = &patchDef;
   par.optDisplay = 30;
diff --git a/contrib/MeshOptimizer/MeshOpt.h b/contrib/MeshOptimizer/MeshOpt.h
index f0ba9f8416ed60269431140afa4f034034721ef2..844d0ca61e54b4e43438277f9a6dde991a6bf380 100644
--- a/contrib/MeshOptimizer/MeshOpt.h
+++ b/contrib/MeshOptimizer/MeshOpt.h
@@ -47,7 +47,7 @@ class MeshOpt
 {
 public:
   Patch patch;
-  MeshOpt(const std::map<MElement*,GEntity*> &element2entity,
+  MeshOpt(const std::map<MElement*, GEntity*> &element2entity,
           const std::set<MElement*> &els, std::set<MVertex*> &toFix,
           const MeshOptParameters &par);
   ~MeshOpt();
diff --git a/contrib/MeshOptimizer/MeshOptCommon.h b/contrib/MeshOptimizer/MeshOptCommon.h
index 189c7689d5030a4d4c009f794c3ebe208a2ebfce..dbb81abf1918c777a480453cee5e9a1bea6c8779 100644
--- a/contrib/MeshOptimizer/MeshOptCommon.h
+++ b/contrib/MeshOptimizer/MeshOptCommon.h
@@ -34,6 +34,7 @@
 
 
 class MElement;
+class GEntity;
 class SPoint3;
 class ObjContrib;
 
@@ -51,11 +52,13 @@ public:
     bool weakMerge;                                     // If connected strategy: weak or strong merging of patches
   };
   virtual ~MeshOptPatchDef() {}
-  virtual double elBadness(MElement *el) const = 0;     // Determine "badness" of a given element (for patch creation)
+  virtual double elBadness(MElement *el,                // Determine "badness" of a given element (for patch creation)
+                           GEntity* gEnt) const = 0;
   virtual double maxDistance(MElement *el) const = 0;   // Compute max. distance to a given bad element for elements in patch
   virtual int inPatch(const SPoint3 &badBary,           // Determine whether a given element should be included in the patch around a...
                       double limDist,                   // ... given bad element barycenter, with a limit distance if needed. Output: ...
-                      MElement *el) const = 0;          // ... -1 = excluded, 0 = included only up to minLayers, 1 = included up to maxLayers
+                      MElement *el,                     // ... -1 = excluded, 0 = included only up to minLayers, 1 = included up to maxLayers
+                      GEntity* gEnt) const = 0;
 protected:
   bool testElInDist(const SPoint3 &P, double limDist,   // Test whether an element is within a certain distance from a point
                     MElement *el) const;
@@ -74,6 +77,7 @@ struct MeshOptParameters {                              // Parameters controllin
   int dim ;                                             // Which dimension to optimize
   bool onlyVisible ;                                    // Apply optimization to visible entities ONLY
   bool fixBndNodes;                                     // If points can move on boundaries
+  bool useGeom;                                         // Compute and use info from geometric (CAD model) entities where helpful
   MeshOptPatchDef *patchDef;
   std::vector<MeshOptPass> pass;
   int optDisplay;                                       // Sampling rate in opt. iterations for display
diff --git a/contrib/MeshOptimizer/MeshOptPatch.cpp b/contrib/MeshOptimizer/MeshOptPatch.cpp
index 6b4e866c696187b1e1671c11bdc2c90613576b01..ce0aa46d77da05e3777efc2602f953f4f7a3c96a 100644
--- a/contrib/MeshOptimizer/MeshOptPatch.cpp
+++ b/contrib/MeshOptimizer/MeshOptPatch.cpp
@@ -40,7 +40,7 @@
 #include "MeshOptPatch.h"
 
 
-Patch::Patch(const std::map<MElement*,GEntity*> &element2entity,
+Patch::Patch(const std::map<MElement*, GEntity*> &element2entity,
            const std::set<MElement*> &els, std::set<MVertex*> &toFix,
            bool fixBndNodes) :
   _typeLengthScale(LS_NONE), _invLengthScaleSq(0.)
@@ -57,11 +57,16 @@ Patch::Patch(const std::map<MElement*,GEntity*> &element2entity,
   _el2V.resize(nElements);
   _nNodEl.resize(nElements);
   _indPCEl.resize(nElements);
+  if (!element2entity.empty()) _gEnt.resize(nElements);
   int iEl = 0;
   bool nonGeoMove = false;
   for(std::set<MElement*>::const_iterator it = els.begin();
       it != els.end(); ++it, ++iEl) {
     _el[iEl] = *it;
+    if (!element2entity.empty()) {
+      std::map<MElement*,GEntity*>::const_iterator itEl2Ent = element2entity.find(*it);
+      _gEnt[iEl] = (itEl2Ent == element2entity.end()) ? 0 : itEl2Ent->second;
+    }
     _nNodEl[iEl] = _el[iEl]->getNumVertices();
     for (int iVEl = 0; iVEl < _nNodEl[iEl]; iVEl++) {
       MVertex *vert = _el[iEl]->getVertex(iVEl);
@@ -88,7 +93,7 @@ Patch::Patch(const std::map<MElement*,GEntity*> &element2entity,
   }
 
   if (nonGeoMove) Msg::Warning("Some vertices will be moved along local lines "
-                            "or planes, they may not remain on the exact geometry");
+                               "or planes, they may not remain on the exact geometry");
 
   // Initial coordinates
   _ixyz.resize(nVert());
@@ -223,21 +228,21 @@ void Patch::initScaledNodeDispSq(LengthScaling scaling)
     switch(scaling) {
       case LS_MAXNODEDIST : {
         for (int iEl = 0; iEl < nEl(); iEl++) {
-          const double d = el(iEl)->maxDistToStraight(), dd = d*d;
+          const double d = _el[iEl]->maxDistToStraight(), dd = d*d;
           if (dd > maxDSq) maxDSq = dd;
         }
         break;
       }
       case LS_MAXOUTERRADIUS : {
         for (int iEl = 0; iEl < nEl(); iEl++) {
-          const double d = el(iEl)->getOuterRadius(), dd = d*d;
+          const double d = _el[iEl]->getOuterRadius(), dd = d*d;
           if (dd > maxDSq) maxDSq = dd;
         }
         break;
       }
       case LS_MINEDGELENGTH : {
         for (int iEl = 0; iEl < nEl(); iEl++) {
-          const double d = el(iEl)->minEdge(), dd = d*d;
+          const double d = _el[iEl]->minEdge(), dd = d*d;
           if (dd > maxDSq) maxDSq = dd;
         }
         break;
@@ -282,7 +287,6 @@ void Patch::initScaledJac()
   // Jacobians of 3D elements
   if ((_dim == 2) && _JacNormEl.empty()) {
     _JacNormEl.resize(nEl());
-//    for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(element2entity,iEl);
     for (int iEl = 0; iEl < nEl(); iEl++)
       calcNormalEl2D(iEl, NS_INVNORM, _JacNormEl[iEl], false);
   }
@@ -306,34 +310,31 @@ void Patch::initMetricMin()
 }
 
 
-// TODO: Re-introduce normal to geometry?
-//void Mesh::calcScaledNormalEl2D(const std::map<MElement*,GEntity*> &element2entity, int iEl)
 void Patch::calcNormalEl2D(int iEl, NormalScaling scaling,
                            fullMatrix<double> &elNorm, bool ideal)
 {
   const JacobianBasis *jac = _el[iEl]->getJacobianFuncSpace();
 
   fullMatrix<double> primNodesXYZ(jac->getNumPrimMapNodes(),3);
-//  SVector3 geoNorm(0.,0.,0.);
-//  std::map<MElement*,GEntity*>::const_iterator itEl2ent = element2entity.find(_el[iEl]);
-//  GEntity *ge = (itEl2ent == element2entity.end()) ? 0 : itEl2ent->second;
-//  const bool hasGeoNorm = ge && (ge->dim() == 2) && ge->haveParametrization();
+  SVector3 geoNorm(0.,0.,0.);
+  GEntity *ge = (_gEnt.empty()) ? 0 : _gEnt[iEl];
+  const bool hasGeoNorm = ge && (ge->dim() == 2) && ge->haveParametrization();
   for (int i=0; i<jac->getNumPrimMapNodes(); i++) {
     const int &iV = _el2V[iEl][i];
     primNodesXYZ(i, 0) = _xyz[iV].x();
     primNodesXYZ(i, 1) = _xyz[iV].y();
     primNodesXYZ(i, 2) = _xyz[iV].z();
-//    if (hasGeoNorm && (_vert[iV]->onWhat() == ge)) {
-//      double u, v;
-//      _vert[iV]->getParameter(0,u);
-//      _vert[iV]->getParameter(1,v);
-//      geoNorm += ((GFace*)ge)->normal(SPoint2(u,v));
-//    }
+    if (hasGeoNorm && (_vert[iV]->onWhat() == ge)) {
+      double u, v;
+      _vert[iV]->getParameter(0,u);
+      _vert[iV]->getParameter(1,v);
+      geoNorm += ((GFace*)ge)->normal(SPoint2(u,v));
+    }
+  }
+  if (hasGeoNorm && (geoNorm.normSq() == 0.)) {
+    SPoint2 param = ((GFace*)ge)->parFromPoint(_el[iEl]->barycenter(true),false);
+    geoNorm = ((GFace*)ge)->normal(param);
   }
-//  if (hasGeoNorm && (geoNorm.normSq() == 0.)) {
-//    SPoint2 param = ((GFace*)ge)->parFromPoint(_el[iEl]->barycenter(true),false);
-//    geoNorm = ((GFace*)ge)->normal(param);
-//  }
 
   elNorm.resize(1, 3);
   const double norm = jac->getPrimNormal2D(primNodesXYZ, elNorm, ideal);
@@ -349,10 +350,10 @@ void Patch::calcNormalEl2D(int iEl, NormalScaling scaling,
       factor = sqrt(norm);
       break;
   }
-//  if (hasGeoNorm) {
-//    const double scal = geoNorm(0)*elNorm(0,0)+geoNorm(1)*elNorm(0,1)+geoNorm(2)*elNorm(0,2);
-//    if (scal < 0.) factor = -factor;
-//  }
+  if (hasGeoNorm) {
+    const double scal = geoNorm(0)*elNorm(0,0)+geoNorm(1)*elNorm(0,1)+geoNorm(2)*elNorm(0,2);
+    if (scal < 0.) factor = -factor;
+  }
   elNorm.scale(factor);   // Re-scaling normal here is faster than an extra scaling operation on the Jacobian
 }
 
@@ -537,7 +538,6 @@ void Patch::initIdealJac()
   // Jacobians of 3D elements
   if ((_dim == 2) && _IJacNormEl.empty()) {
     _IJacNormEl.resize(nEl());
-//    for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(element2entity,iEl);
     for (int iEl = 0; iEl < nEl(); iEl++)
       calcNormalEl2D(iEl, NS_INVNORM, _IJacNormEl[iEl], true);
   }
@@ -575,18 +575,6 @@ void Patch::idealJacAndGradients(int iEl, std::vector<double> &iJ, std::vector<d
   // regularization normals in 2D)
   jacBasis->getSignedIdealJacAndGradients(nodesXYZ,_IJacNormEl[iEl],JDJ);
   if (_dim == 3) JDJ.scale(_invIJac[iEl]);
-//  if (_el[iEl]->getNum() == 90370) {
-//    std::cout << "DBGTT: bad el.: " << _el[iEl]->getNum() << "\n";
-//    for (int i = 0; i < numMapNodes; i++)
-//      std::cout << "DBGTT: {x,y,z}" << i << " = (" << nodesXYZ(i,0)
-//                << ", " << nodesXYZ(i,1) << ", " << nodesXYZ(i,2) << ")\n";
-//   for (int l = 0; l < numJacNodes; l++) {
-//      for (int i = 0; i < numMapNodes; i++)
-//        std::cout << "DBGTT: dJ" << l << "d{x,y,z}" << i << " = (" << JDJ(l, i)
-//                  << ", " << JDJ(l, i+numMapNodes) << ", " << JDJ(l, i+2*numMapNodes)<< ")\n";
-//      std::cout << "DBGTT: J" << l << " = " << JDJ(l, 3*numMapNodes)<< "\n";
-//    }
-//  }
 
   // Transform Jacobian and gradients from Lagrangian to Bezier basis
   jacBasis->lag2Bez(JDJ,BDB);
@@ -629,7 +617,6 @@ void Patch::initInvCondNum()
   // Set normals to 2D elements
   if ((_dim == 2) && _condNormEl.empty()) {
     _condNormEl.resize(nEl());
-//    for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(element2entity,iEl);
     for (int iEl = 0; iEl < nEl(); iEl++)
       calcNormalEl2D(iEl, NS_UNIT, _condNormEl[iEl], true);
   }
diff --git a/contrib/MeshOptimizer/MeshOptPatch.h b/contrib/MeshOptimizer/MeshOptPatch.h
index 7e18e8409ba209a420ed5f8c3ae118e48a56165e..6dea34f62de5bda5af56c0c04d65ae19d327f315 100644
--- a/contrib/MeshOptimizer/MeshOptPatch.h
+++ b/contrib/MeshOptimizer/MeshOptPatch.h
@@ -63,7 +63,6 @@ public:
   inline const int &indPCEl(int iEl, int iPC) { return _indPCEl[iEl][iPC]; }
   inline const int &el2V(int iEl, int i) { return _el2V[iEl][i]; }
   inline const int &el2FV(int iEl, int i) { return _el2FV[iEl][i]; }
-  inline MElement *el(int iEl) { return _el[iEl]; }
   inline const int &fv2V(int iFV) { return _fv2V[iFV]; }
   inline const SPoint3 &xyz(int iV) { return _xyz[iV]; }
   inline const SPoint3 &ixyz(int iV) { return _ixyz[iV]; }
@@ -109,6 +108,7 @@ private:
   int _dim;
   int _nPC;                                         // Total nb. of parametric coordinates
   std::vector<MElement*> _el;                       // List of elements
+  std::vector<GEntity*> _gEnt;                      // Geometric entity corresponding to each element
   std::vector<MVertex*> _vert, _freeVert;           // List of vert., free vert.
   std::vector<int> _fv2V;                           // Index of free vert. -> index of vert.
   std::vector<bool> _forced;                        // Is vertex forced?
@@ -149,7 +149,6 @@ private:
   std::vector<int> _nBezEl;                                     // Number of Bezier poly. for an el.
   std::vector<fullMatrix<double> > _JacNormEl;                  // Normals to 2D elements for Jacobian regularization and scaling
   std::vector<double> _invStraightJac;                          // Initial Jacobians for 3D elements
-//  void calcScaledNormalEl2D(const std::map<MElement*,GEntity*> &element2entity, int iEl);
   void calcNormalEl2D(int iEl, NormalScaling scaling,
                       fullMatrix<double> &elNorm, bool ideal);
 
diff --git a/contrib/MeshOptimizer/MeshOptimizer.cpp b/contrib/MeshOptimizer/MeshOptimizer.cpp
index 9694635bb647b9a2051b4df2ada3eb33c43f36d0..5ddd0848f88e5c60b4f32cb018c06f56af5e1125 100644
--- a/contrib/MeshOptimizer/MeshOptimizer.cpp
+++ b/contrib/MeshOptimizer/MeshOptimizer.cpp
@@ -105,7 +105,8 @@ void getElementNeighbours(MElement *el, const vertElVecMap &v2e,
 elSet getSurroundingPatch(MElement *el, const MeshOptPatchDef *patchDef,
                           double limDist, int maxLayers,
                           const vertElVecMap &vertex2elements,
-                          elElSetMap &element2elements)
+                          elElSetMap &element2elements,
+                          const elEntMap &element2entity)
 {
   const SPoint3 pnt = el->barycenter(true);
 
@@ -121,7 +122,12 @@ elSet getSurroundingPatch(MElement *el, const MeshOptPatchDef *patchDef,
       for (elSetIter itN = neighbours.begin(); itN != neighbours.end(); ++itN) {      // Loop over neighbours
         if ((lastLayer.find(*itN) == lastLayer.end()) &&                              // If neighbour already in last layer...
             (excluded.find(*itN) == excluded.end())) {                                // ... or marked as excluded, skip
-          const int elIn = patchDef->inPatch(pnt, limDist, *itN);                     // Test if element in patch according to user-defined criteria
+          GEntity *gEnt = 0;
+          if (!element2entity.empty()) {
+            elEntMap::const_iterator itEl2Ent = element2entity.find(el);
+            if (itEl2Ent != element2entity.end()) gEnt = itEl2Ent->second;
+          }
+          const int elIn = patchDef->inPatch(pnt, limDist, *itN, gEnt);               // Test if element in patch according to user-defined criteria
           if ((elIn > 0) || ((d < patchDef->minLayers) && (elIn == 0))) {
             if (patch.insert(*itN).second) currentLayer.insert(*itN);                 // If element in, insert in patch and in current layer...
           }
@@ -173,6 +179,7 @@ void calcElement2Entity(GEntity *entity, elEntMap &element2entity)
 
 
 std::vector<elSetVertSetPair> getConnectedPatches(const vertElVecMap &vertex2elements,
+                                                  const elEntMap &element2entity,
                                                   const elSet &badElements,
                                                   const MeshOptParameters &par)
 {
@@ -188,7 +195,8 @@ std::vector<elSetVertSetPair> getConnectedPatches(const vertElVecMap &vertex2ele
     const double limDist = par.patchDef->maxDistance(*it);
     primPatches.push_back(getSurroundingPatch(*it, par.patchDef, limDist,
                                               par.patchDef->maxLayers,
-                                              vertex2elements, element2elements));
+                                              vertex2elements, element2elements,
+                                              element2entity));
   }
 
   // Compute patch connectivity
@@ -252,6 +260,7 @@ void optimizeConnectedPatches(const vertElVecMap &vertex2elements,
   par.success = 1;
 
   std::vector<elSetVertSetPair> toOptimize = getConnectedPatches(vertex2elements,
+                                                                 element2entity,
                                                                  badasses, par);
 
   for (int iPatch = 0; iPatch < toOptimize.size(); ++iPatch) {
@@ -290,13 +299,20 @@ void optimizeConnectedPatches(const vertElVecMap &vertex2elements,
 }
 
 
-MElement *getWorstElement(elSet &badElts, const MeshOptParameters &par)
+MElement *getWorstElement(elSet &badElts,
+                          const elEntMap &element2entity,
+                          const MeshOptParameters &par)
 {
   double worst = 1.e300;
   MElement *worstEl = 0;
 
   for (elSetIter it=badElts.begin(); it!=badElts.end(); it++) {
-    const double val = par.patchDef->elBadness(*it);
+    GEntity *gEnt = 0;
+    if (!element2entity.empty()) {
+      elEntMap::const_iterator itEl2Ent = element2entity.find(*it);
+      if (itEl2Ent != element2entity.end()) gEnt = itEl2Ent->second;
+    }
+    const double val = par.patchDef->elBadness(*it, gEnt);
     if (val < worst) {
       worst = val;
       worstEl = *it;
@@ -324,7 +340,7 @@ void optimizeOneByOne(const vertElVecMap &vertex2elements,
     if (badElts.empty()) break;
 
     // Create patch around worst element and remove it from badElts
-    MElement *worstEl = getWorstElement(badElts, par);
+    MElement *worstEl = getWorstElement(badElts, element2entity, par);
     badElts.erase(worstEl);
 
     // Initialize patch size to be adapted
@@ -338,7 +354,8 @@ void optimizeOneByOne(const vertElVecMap &vertex2elements,
       // Set up patch
       const double limDist = par.patchDef->maxDistance(worstEl);
       elSet toOptimizePrim = getSurroundingPatch(worstEl, par.patchDef, limDist,
-                                                 maxLayers, vertex2elements, element2elements);
+                                                 maxLayers, vertex2elements,
+                                                 element2elements, element2entity);
       vertSet toFix = getAllBndVertices(toOptimizePrim, vertex2elements);
       elSet toOptimize;
       std::set_difference(toOptimizePrim.begin(),toOptimizePrim.end(),
@@ -428,11 +445,13 @@ void meshOptimizer(GModel *gm, MeshOptParameters &par)
         (par.onlyVisible && !entity->getVisibility())) continue;
     Msg::Info("Computing connectivity and bad elements for entity %d...",
               entity->tag());
-    calcVertex2Elements(par.dim,entity,vertex2elements);
+    calcVertex2Elements(par.dim, entity, vertex2elements);
+    if (par.useGeom) calcElement2Entity(entity, element2entity);
     for (int iEl = 0; iEl < entity->getNumMeshElements();iEl++) {                               // Detect bad elements
       double jmin, jmax;
       MElement *el = entity->getMeshElement(iEl);
-      if ((el->getDim() == par.dim) && (par.patchDef->elBadness(el) < 0.)) badElts.insert(el);
+      if ((el->getDim() == par.dim) && (par.patchDef->elBadness(el, entity) < 0.))
+        badElts.insert(el);
     }
   }
 
diff --git a/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp b/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp
index e5c75535ac03327b3c9c457429c2c542ae18c18e..f0ff5c58e276cc64e7421302315d6d885ff0627b 100644
--- a/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp
+++ b/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp
@@ -2,6 +2,8 @@
 
 //#include "GModel.h"
 #include "GEntity.h"
+#include "GFace.h"
+#include "GRegion.h"
 #include "MElement.h"
 #include "MTriangle.h"
 #include "MQuadrangle.h"
@@ -20,12 +22,12 @@ struct QualPatchDefParameters : public MeshOptPatchDef
 {
   QualPatchDefParameters(const MeshQualOptParameters &p);
   virtual ~QualPatchDefParameters() {}
-  virtual double elBadness(MElement *el) const;
+  virtual double elBadness(MElement *el, GEntity* gEnt) const;
   virtual double maxDistance(MElement *el) const;
-  virtual int inPatch(const SPoint3 &badBary,
-                      double limDist, MElement *el) const;
+  virtual int inPatch(const SPoint3 &badBary, double limDist,
+                      MElement *el, GEntity* gEnt) const;
 private:
-  bool _excludeQuad, _excludeHex, _excludePrism;
+  bool _excludeQuad, _excludeHex, _excludePrism, _excludeBL;
   double _idealJacMin, _invCondNumMin;
   double _distanceFactor;
 };
@@ -36,6 +38,7 @@ QualPatchDefParameters::QualPatchDefParameters(const MeshQualOptParameters &p)
   _excludeQuad = p.excludeQuad;
   _excludeHex = p.excludeHex;
   _excludePrism = p.excludePrism;
+  _excludeBL = p.excludeBL;
   _idealJacMin = p.minTargetIdealJac;
   _invCondNumMin = p.minTargetInvCondNum;
   strategy = (p.strategy == 1) ? MeshOptParameters::STRAT_ONEBYONE :
@@ -53,12 +56,23 @@ QualPatchDefParameters::QualPatchDefParameters(const MeshQualOptParameters &p)
 }
 
 
-double QualPatchDefParameters::elBadness(MElement *el) const
+double QualPatchDefParameters::elBadness(MElement *el, GEntity* gEnt) const
 {
   const int typ = el->getType();
   if (_excludeQuad && (typ == TYPE_QUA)) return 1.;
   if (_excludeHex && (typ == TYPE_HEX)) return 1.;
   if (_excludePrism && (typ == TYPE_PRI)) return 1.;
+  if (_excludeBL) {
+    BoundaryLayerColumns *blc = 0;
+    if (gEnt->dim() == 2)
+      blc = static_cast<GFace*>(gEnt)->getColumns();
+    else if (gEnt->dim() == 3)
+      blc = static_cast<GRegion*>(gEnt)->getColumns();
+    if (blc) {
+      std::map<MElement*, MElement*>::iterator itBLEl = blc->_toFirst.find(el);
+      if (itBLEl != blc->_toFirst.end()) return 1.;
+    }
+  }
 //  double jMin, jMax;
 //  el->idealJacRange(jMin, jMax);
 //  return jMin-_idealJacMin;
@@ -74,13 +88,24 @@ double QualPatchDefParameters::maxDistance(MElement *el) const
 }
 
 
-int QualPatchDefParameters::inPatch(const SPoint3 &badBary,
-                                    double limDist, MElement *el) const
+int QualPatchDefParameters::inPatch(const SPoint3 &badBary, double limDist,
+                                    MElement *el, GEntity* gEnt) const
 {
   const int typ = el->getType();
   if (_excludeQuad && (typ == TYPE_QUA)) return -1;
   if (_excludeHex && (typ == TYPE_HEX)) return -1;
   if (_excludePrism && (typ == TYPE_PRI)) return -1;
+  if (_excludeBL) {
+    BoundaryLayerColumns *blc = 0;
+    if (gEnt->dim() == 2)
+      blc = static_cast<GFace*>(gEnt)->getColumns();
+    else if (gEnt->dim() == 3)
+      blc = static_cast<GRegion*>(gEnt)->getColumns();
+    if (blc) {
+      std::map<MElement*, MElement*>::iterator itBLEl = blc->_toFirst.find(el);
+      if (itBLEl != blc->_toFirst.end()) return -1;
+    }
+  }
   return testElInDist(badBary, limDist, el) ? 1 : 0;
 }
 
@@ -93,6 +118,7 @@ void MeshQualityOptimizer(GModel *gm, MeshQualOptParameters &p)
   par.dim = p.dim;
   par.onlyVisible = p.onlyVisible;
   par.fixBndNodes = p.fixBndNodes;
+  par.useGeom = p.excludeBL;
   QualPatchDefParameters patchDef(p);
   par.patchDef = &patchDef;
   par.optDisplay = 20;
diff --git a/contrib/MeshQualityOptimizer/MeshQualityOptimizer.h b/contrib/MeshQualityOptimizer/MeshQualityOptimizer.h
index 8353638e4782aba17eb133d6289c39d7a317cc70..d722f456e84a20949e7408fe6004a505f839bf87 100644
--- a/contrib/MeshQualityOptimizer/MeshQualityOptimizer.h
+++ b/contrib/MeshQualityOptimizer/MeshQualityOptimizer.h
@@ -34,7 +34,7 @@ class GModel;
 
 struct MeshQualOptParameters {
   bool onlyValidity;
-  bool excludeQuad, excludeHex, excludePrism;
+  bool excludeQuad, excludeHex, excludePrism, excludeBL;
   double minTargetIdealJac;
   double minTargetInvCondNum;
   double weightFixed;                                                 // weight of the energy for fixed nodes
@@ -58,7 +58,7 @@ struct MeshQualOptParameters {
 
   MeshQualOptParameters ()
     : onlyValidity(false), excludeQuad(false),
-      excludeHex(false), excludePrism(false),
+      excludeHex(false), excludePrism(false), excludeBL(false),
       minTargetIdealJac(0.1), minTargetInvCondNum(0.1), weightFixed(1000.),
       weightFree (1.), nbLayers (6) , dim(3) , itMax(300), optPassMax(50),
       onlyVisible(true), distanceFactor(12), fixBndNodes(false), strategy(0),