Added possibility of criterion on geometric entity for patch creation in mesh...

Added possibility of criterion on geometric entity for patch creation in mesh optimizer. Added possibility to exclude BL elements in mesh quality optimizer.

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;

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

contrib/MeshOptimizer/MeshOptPatch.cpp

@@ -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);
@@ -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 && (geoNorm.normSq() == 0.)) {
-//    SPoint2 param = ((GFace*)ge)->parFromPoint(_el[iEl]->barycenter(true),false);
-//    geoNorm = ((GFace*)ge)->normal(param);
-//  }
+    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);
+  }
 
   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);
   }

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);
 

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(),
@@ -429,10 +446,12 @@ void meshOptimizer(GModel *gm, MeshOptParameters &par)
     Msg::Info("Computing connectivity and bad elements for entity %d...",
               entity->tag());
     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);
     }
   }
 

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;

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),