Improved OptHOM: cleaned up, added strategy for 3D and small features.

Fltk/highOrderToolsWindow.cpp

@@ -96,20 +96,33 @@ static void chooseopti_cb(Fl_Widget *w, void *data)
 
   if (elastic == 1){
     o->choice[0]->deactivate();
+    o->choice[3]->deactivate();
     for (int i=3;i<=6;i++)
       o->value[i]->deactivate();
+    for (int i=9;i<=11;i++) o->value[i]->deactivate();
     //   o->push[1]->deactivate();
   }
   else {
     o->push[0]->activate();
     o->choice[0]->activate();
+    o->choice[3]->activate();
     for (int i=3;i<=6;i++)
       o->value[i]->activate();
+    for (int i=9;i<=11;i++) o->value[i]->activate();
     //    o->push[1]->activate();
   }
 #endif
 }
 
+static void chooseopti_strategy(Fl_Widget *w, void *data)
+{
+#if defined(HAVE_OPTHOM)
+  highOrderToolsWindow *o = FlGui::instance()->highordertools;
+  if (o->choice[3]->value() == 0) for (int i=9;i<=11;i++) o->value[i]->deactivate();
+  else for (int i=9;i<=11;i++) o->value[i]->activate();
+#endif
+}
+
 static void highordertools_runelas_cb(Fl_Widget *w, void *data)
 {
 #if defined(HAVE_OPTHOM)
@@ -133,8 +146,12 @@ static void highordertools_runelas_cb(Fl_Widget *w, void *data)
     p.optPassMax = (int) o->value[4]->value();
     p.weightFixed =  o->value[5]->value();
     p.weightFree =  o->value[6]->value();
-    p.DistanceFactor =  o->value[7]->value();
-    p.method =  o->CAD ? (int)o->choice[0]->value() : 2;
+    p.distanceFactor =  o->value[7]->value();
+    p.fixBndNodes = (!o->CAD) || (o->choice[0]->value() == 0);
+    p.strategy = o->choice[3]->value();
+    p.maxAdaptBlob = o->value[9]->value();
+    p.adaptBlobLayerFact = (int) o->value[10]->value();
+    p.adaptBlobDistFact = o->value[11]->value();
     HighOrderMeshOptimizer (GModel::current(),p);
     printf("CPU TIME = %4f seconds\n",p.CPU);
   }
@@ -149,7 +166,7 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
   FL_NORMAL_SIZE -= deltaFontSize;
 
   int width = 3 * IW + 4 * WB;
-  int height = 25 * BH;
+  int height = 28 * BH;
 
   win = new paletteWindow
     (width, height, CTX::instance()->nonModalWindows ? true : false, "High Order Tools");
@@ -315,6 +332,48 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
   value[4]->align(FL_ALIGN_RIGHT);
   value[4]->value(50);
 
+  static Fl_Menu_Item menu_strategy[] = {
+    {"Connected blobs", 0, 0, 0},
+    {"Adaptive one-by-one", 0, 0, 0},
+    {0}
+  };
+
+  y += BH;
+  choice[3] = new Fl_Choice
+    (x,y, IW, BH, "Strategy");
+  choice[3]->menu(menu_strategy);
+  choice[3]->align(FL_ALIGN_RIGHT);
+  choice[3]->callback(chooseopti_strategy);
+
+  y += BH;
+  value[9] = new Fl_Value_Input
+    (x,y, IW, BH, "Max. number of blob adaptation iter.");
+  value[9]->minimum(1);
+  value[9]->maximum(100);
+  value[9]->step(1);
+  value[9]->align(FL_ALIGN_RIGHT);
+  value[9]->value(2);
+  value[9]->deactivate();
+
+  y += BH;
+  value[10] = new Fl_Value_Input
+    (x,y, IW, BH, "Num. layer adaptation factor");
+  value[10]->align(FL_ALIGN_RIGHT);
+  value[10]->minimum(1);
+  value[10]->maximum(100);
+  value[10]->step(1);
+  value[10]->value(2);
+  value[10]->deactivate();
+
+  y += BH;
+  value[11] = new Fl_Value_Input
+    (x,y, IW, BH, "Distance adaptation factor");
+  value[11]->align(FL_ALIGN_RIGHT);
+  value[11]->minimum(1.);
+  value[11]->maximum(100.);
+  value[11]->value(2.);
+  value[11]->deactivate();
+
   y += 1.5*BH;
   push[1] = new Fl_Button
     (x,y, IW, BH, "Apply");

Geo/MPyramid.h

@@ -221,18 +221,17 @@ class MPyramidN : public MPyramid {
  protected:
   std::vector<MVertex*> _vs;
   const char _order;
-  double _disto;
  public:
   MPyramidN(MVertex* v0, MVertex* v1, MVertex* v2, MVertex* v3, MVertex* v4,
       const std::vector<MVertex*> &v, char order, int num=0, int part=0)
-    : MPyramid(v0, v1, v2, v3, v4, num, part), _vs(v), _order(order), _disto(-1e22)
+    : MPyramid(v0, v1, v2, v3, v4, num, part), _vs(v), _order(order)
   {
     for (unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
     getFunctionSpace(order);
   }
 
   MPyramidN(const std::vector<MVertex*> &v, char order, int num=0, int part=0)
-    : MPyramid(v[0], v[1], v[2], v[3], v[4], num, part), _order(order), _disto(-1e22)
+    : MPyramid(v[0], v[1], v[2], v[3], v[4], num, part), _order(order)
   {
     for (unsigned int i = 5; i < v.size(); i++ ) _vs.push_back(v[i]);
     for (unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);

contrib/HighOrderMeshOptimizer/OptHOM.cpp

@@ -42,11 +42,11 @@ static inline double compute_f1(double v, double barrier)
 
 
 // Constructor
-OptHOM::OptHOM(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> & toFix, int method) :
-       mesh(ge, els, toFix, method)
+OptHOM::OptHOM(const std::set<MElement*> &els, std::set<MVertex*> & toFix, bool fixBndNodes) :
+       mesh(els, toFix, fixBndNodes)
 {
   _optimizeMetricMin = false;
-};
+}
 
 
 
@@ -241,6 +241,8 @@ void OptHOM::OptimPass(alglib::real_1d_array &x, const alglib::real_1d_array &in
   static int OPTMETHOD = 1;
 
   Msg::Debug("--- Optimization pass with jacBar = %12.5E",jacBar);
+  std::cout << "--- Optimization pass with initial jac. range ("
+            << minJac << "," << maxJac << "), jacBar = " << jacBar << "\n";
 
   iter = 0;
 
@@ -321,7 +323,6 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double
 
   const double jacBarStart = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
   jacBar = jacBarStart;
-  setBarrierTerm(jacBarStart);
 
   _optimizeBarrierMax = false;
   // Calculate initial objective function value and gradient
@@ -337,24 +338,37 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double
             << " and max. barrier = " << barrier_max << std::endl;
 
   int ITER = 0;
+  bool minJacOK = true;
   while (minJac < barrier_min) {
+    const double startMinJac = minJac;
     OptimPass(x, gradObj, itMax);
     recalcJacDist();
     jacBar = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
-    setBarrierTerm(jacBar);
-    if (ITER ++ > optPassMax) break;
+    if (ITER ++ > optPassMax) {
+      minJacOK = (minJac > barrier_min);
+      break;
+    }
+    if (fabs((minJac-startMinJac)/startMinJac) < 0.01) {
+      std::cout << "Stagnation in minJac detected, stopping optimization\n";
+      minJacOK = false;
+      break;
+    }
   }
 
-  if (!_optimizeMetricMin) {
+  ITER = 0;
+  if (minJacOK && (!_optimizeMetricMin)) {
     _optimizeBarrierMax = true;
     jacBar =  1.1 * maxJac;
-    setBarrierTerm(jacBar);
     while (maxJac > barrier_max ) {
+      const double startMaxJac = maxJac;
       OptimPass(x, gradObj, itMax);
       recalcJacDist();
       jacBar =  1.1 * maxJac;
-      setBarrierTerm(jacBar);
       if (ITER ++ > optPassMax) break;
+      if (fabs((maxJac-startMaxJac)/startMaxJac) < 0.01) {
+        std::cout << "Stagnation in maxJac detected, stopping optimization\n";
+        break;
+      }
     }
   }
 

contrib/HighOrderMeshOptimizer/OptHOM.h

@@ -21,7 +21,7 @@ public:
 
   Mesh mesh;
 
-  OptHOM(GEntity *gf, const std::set<MElement*> &els, std::set<MVertex*> & toFix, int method);
+  OptHOM(const std::set<MElement*> &els, std::set<MVertex*> & toFix, bool fixBndNodes);
   // returns 1 if the mesh has been optimized with success i.e. all jacobians are in the range
   // returns 0 if the mesh is valid (all jacobians positive, JMIN > 0) but JMIN < barrier_min || JMAX > barrier_max
   // returns -1 if the mesh is invalid : some jacobians cannot be made positive
@@ -36,7 +36,6 @@ public:
 
 private:
 
-//  double lambda, lambda2, powM, powP, invLengthScaleSq;
   double lambda, lambda2, jacBar, invLengthScaleSq;
   int iter, progressInterv;            // Current iteration, interval of iterations for reporting
   bool _optimizeMetricMin;
@@ -45,7 +44,6 @@ private:
 
   bool _optimizeBarrierMax; // false : only moving barrier min; true : fixed barrier min + moving barrier max
 
-  inline void setBarrierTerm(double jacBarrier) {jacBar = jacBarrier;}
   bool addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj);
   bool addMetricMinObjGrad(double &Obj, alglib::real_1d_array &gradObj);
   bool addDistObjGrad(double Fact, double Fact2, double &Obj, alglib::real_1d_array &gradObj);

contrib/HighOrderMeshOptimizer/OptHomMesh.cpp

@@ -8,32 +8,21 @@
 
 
 
-Mesh::Mesh(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> &toFix, int method) :
-    _ge(ge)
+Mesh::Mesh(const std::set<MElement*> &els, std::set<MVertex*> &toFix, bool fixBndNodes)
 {
-  _dim = _ge->dim();
 
-  if (method & METHOD_PHYSCOORD) {
+  _dim = (*els.begin())->getDim();
+
+  if (fixBndNodes) {
     if (_dim == 2) _pc = new ParamCoordPhys2D;
     else _pc = new ParamCoordPhys3D;
-    Msg::Debug("METHOD: Using physical coordinates");
-  }
-  else if (method & METHOD_SURFCOORD)  {
-    if (_dim == 2) {
-      _pc = new ParamCoordSurf(_ge);
-      Msg::Debug("METHOD: Using surface parametric coordinates");
-    }
-    else Msg::Error("ERROR: Surface parametric coordinates only for 2D optimization");
+    Msg::Debug("METHOD: Fixing boundary nodes and using physical coordinates");
   }
   else {
     _pc = new ParamCoordParent;
-    Msg::Debug("METHOD: Using parent parametric coordinates");
+    Msg::Debug("METHOD: Freeing boundary nodes and using parent parametric coordinates");
   }
 
-  if (method & METHOD_RELAXBND)Msg::Debug("METHOD: Relaxing boundary vertices");
-  else if (method & METHOD_FIXBND) Msg::Debug("METHOD: Fixing all boundary vertices");
-  else Msg::Debug("METHOD: Fixing vertices on geometric points and \"toFix\" boundary");
-
   // Initialize elements, vertices, free vertices and element->vertices connectivity
   const int nElements = els.size();
   _nPC = 0;
@@ -55,8 +44,7 @@ Mesh::Mesh(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> &toFi
       _el2V[iEl].push_back(iV);
       const int nPCV = _pc->nCoord(vert);
       bool isFV = false;
-      if (method & METHOD_RELAXBND) isFV = true;
-      else if (method & METHOD_FIXBND) isFV = (vert->onWhat()->dim() == _dim) && (toFix.find(vert) == toFix.end());
+      if (fixBndNodes) isFV = (vert->onWhat()->dim() == _dim) && (toFix.find(vert) == toFix.end());
       else isFV = (vert->onWhat()->dim() >= 1) && (toFix.find(vert) == toFix.end());
       if (isFV) {
         int iFV = addFreeVert(vert,iV,nPCV,toFix);
@@ -87,14 +75,7 @@ Mesh::Mesh(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> &toFi
     double dumJac[3][3];
     for (int iEl = 0; iEl < nEl(); iEl++) _invStraightJac[iEl] = 1. / _el[iEl]->getPrimaryJacobian(0.,0.,0.,dumJac);
   }
-  if ((_dim == 2) && (method & METHOD_PROJJAC)) {
-    projJac = true;
-    Msg::Debug("METHOD: Using projected Jacobians");
-  }
-  else {
-    projJac = false;
-    Msg::Debug("METHOD: Using usual Jacobians");
-  }
+
 }
 
 
@@ -371,7 +352,7 @@ void Mesh::writeMSH(const char *filename)
   fprintf(f, "%d\n", nEl());
   for (int iEl = 0; iEl < nEl(); iEl++) {
 //    MElement *MEl = _el[iEl];
-    fprintf(f, "%d %d 2 0 %d", iEl+1, _el[iEl]->getTypeForMSH(), _ge->tag());
+    fprintf(f, "%d %d 2 0 0", iEl+1, _el[iEl]->getTypeForMSH());
     for (size_t iVEl = 0; iVEl < _el2V[iEl].size(); iVEl++) fprintf(f, " %d", _el2V[iEl][iVEl] + 1);
     fprintf(f, "\n");
   }

contrib/HighOrderMeshOptimizer/OptHomMesh.h

@@ -17,7 +17,7 @@ class Mesh
 
 public:
 
-  Mesh(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> & toFix, int method);
+  Mesh(const std::set<MElement*> &els, std::set<MVertex*> & toFix, bool fixBndNodes);
 
   inline const int &nPC() { return _nPC; }
   inline int nVert() { return _vert.size(); }
@@ -46,9 +46,6 @@ public:
   void updateGEntityPositions();
   void writeMSH(const char *filename);
 
-  enum { METHOD_RELAXBND = 1 << 0, METHOD_FIXBND = 1 << 1, METHOD_PHYSCOORD = 1 << 2,
-    METHOD_SURFCOORD = 1 << 3, METHOD_PROJJAC = 1 << 4 };
-
 //  inline double xyzDBG(int iV, int iCoord) { return _xyz[iV][iCoord]; }
 //  inline double ixyzDBG(int iV, int iCoord) { return _ixyz[iV][iCoord]; }
 //  inline double sxyzDBG(int iV, int iCoord) { return _sxyz[iV][iCoord]; }
@@ -58,7 +55,6 @@ public:
 
 private:
 
-  GEntity *_ge;
   int _dim;
   int _nPC;                                       // Total nb. of parametric coordinates
 
@@ -77,7 +73,6 @@ private:
   std::vector<std::vector<int> > _indPCEl;        // Index of parametric coord. for an el.
 
   ParamCoord *_pc;
-  bool projJac;                                   // Using "projected" Jacobians or not
 
   int addVert(MVertex* vert);
   int addFreeVert(MVertex* vert, const int iV, const int nPCV, std::set<MVertex*> &toFix);

contrib/HighOrderMeshOptimizer/OptHomRun.h

@@ -16,9 +16,14 @@ struct OptHomParameters {
   int optPassMax ; // max number of optimization passes
   double TMAX ; // max CPU time allowed
   bool onlyVisible ; // apply optimization to visible entities ONLY
-  double DistanceFactor; // filter elements such that no elements further away than
+  double distanceFactor; // filter elements such that no elements further away than
                          // DistanceFactor times the max distance to straight sided version of an element are optimized
-  int method ;  // how jacobians are computed and if points can move on boundaries
+  bool fixBndNodes;  // how jacobians are computed and if points can move on boundaries
+  int strategy;       // 0 = connected blobs, 1 = adaptive one-by-one
+  int maxAdaptBlob;   // Max. nb. of blob adaptation interations
+  int adaptBlobLayerFact;   // Growth factor in number of layers for blob adaptation
+  double adaptBlobDistFact; // Growth factor in distance factor for blob adaptation
+
   // OUTPUT ------>
   int SUCCESS ; // 0 --> success , 1 --> Not converged
   double minJac, maxJac; // after optimization, range of jacobians
@@ -30,7 +35,8 @@ struct OptHomParameters {
   OptHomParameters ()
   // default values
   : BARRIER_MIN_METRIC(-1.), BARRIER_MIN(0.1), BARRIER_MAX(2.0), weightFixed(1.e6), weightFree (1.e2),
-    nbLayers (6) , dim(3) , itMax(10000), onlyVisible(true), DistanceFactor(12), method(1)
+    nbLayers (6) , dim(3) , itMax(300), onlyVisible(true), distanceFactor(12), fixBndNodes(false),
+    strategy(0), maxAdaptBlob(3), adaptBlobLayerFact(2.), adaptBlobDistFact(2.)
   {
   }
 };

contrib/HighOrderMeshOptimizer/ParamCoord.cpp

@@ -6,53 +6,6 @@
 
 
 
-ParamCoordSurf::ParamCoordSurf(GEntity *ge)
-{
-  if ((ge->dim() == 2) && (ge->geomType() != GEntity::DiscreteSurface)) _gf = static_cast<GFace*>(ge);
-  else std::cout << "ERROR: Using surface parametric coordinates with non-surface geometric entity" << std::endl;
-}
-
-
-
-SPoint3 ParamCoordSurf::getUvw(MVertex* vert)
-{
-  SPoint2 p;
-  reparamMeshVertexOnFace(vert,_gf,p);
-  return SPoint3(p[0],p[1],0.);
-}
-
-
-
-SPoint3 ParamCoordSurf::uvw2Xyz(MVertex* vert, const SPoint3 &uvw)
-{
-  GPoint gp = _gf->point(uvw[0],uvw[1]);
-  return SPoint3(gp.x(),gp.y(),gp.z());
-}
-
-
-
-void ParamCoordSurf::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw)
-{
-  Pair<SVector3,SVector3> der = _gf->firstDer(SPoint2(uvw[0],uvw[1]));
-  gUvw[0] = gXyz.x() * der.first().x() + gXyz.y() * der.first().y() + gXyz.z() * der.first().z();
-  gUvw[1] = gXyz.x() * der.second().x() + gXyz.y() * der.second().y() + gXyz.z() * der.second().z();
-}
-
-
-
-void ParamCoordSurf::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw)
-{
-  Pair<SVector3,SVector3> der = _gf->firstDer(SPoint2(uvw[0],uvw[1]));
-  std::vector<SPoint3>::iterator itUvw=gUvw.begin();
-  for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin(); itXyz != gXyz.end(); itXyz++) {
-    (*itUvw)[0] = itXyz->x() * der.first().x() + itXyz->y() * der.first().y() + itXyz->z() * der.first().z();
-    (*itUvw)[1] = itXyz->x() * der.second().x() + itXyz->y() * der.second().y() + itXyz->z() * der.second().z();
-    itUvw++;
-  }
-}
-
-
-
 SPoint3 ParamCoordParent::getUvw(MVertex* vert)
 {
 

contrib/HighOrderMeshOptimizer/ParamCoord.h

@@ -64,26 +64,6 @@ public:
 
 
 
-class ParamCoordSurf : public ParamCoord
-{
-
-public:
-
-  ParamCoordSurf(GEntity *ge);
-  int nCoord(MVertex* vert) { return 2; }
-  virtual SPoint3 getUvw(MVertex* vert);
-  virtual SPoint3 uvw2Xyz(MVertex* vert, const SPoint3 &uvw);
-  virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw);
-  virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw);
-
-private:
-
-  GFace *_gf;
-
-};
-
-
-
 class ParamCoordParent : public ParamCoord
 {