OptHOM : add metric based optimization (only for planar 2D without

boundary layers)

contrib/HighOrderMeshOptimizer/OptHOM.cpp

@@ -20,6 +20,7 @@
 OptHOM::OptHOM(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> & toFix, int method) :
        mesh(ge, els, toFix, method)
 {
+  _optimizeMetricMin = false;
 };
 
 
@@ -50,6 +51,29 @@ bool OptHOM::addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj)
 
 }
 
+bool OptHOM::addMetricMinObjGrad(double &Obj, alglib::real_1d_array &gradObj)
+{
+
+  minJac = 1.e300;
+  maxJac = -1.e300;
+
+  for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
+    std::vector<double> sJ(mesh.nBezEl(iEl));                   // Scaled Jacobians
+    std::vector<double> gSJ(mesh.nBezEl(iEl)*mesh.nPCEl(iEl));  // Gradients of scaled Jacobians
+    mesh.metricMinAndGradients (iEl,sJ,gSJ);
+    
+    for (int l = 0; l < mesh.nBezEl(iEl); l++) {
+      Obj += compute_f(sJ[l]);
+      const double f1 = compute_f1(sJ[l]);
+      for (int iPC = 0; iPC < mesh.nPCEl(iEl); iPC++)
+        gradObj[mesh.indPCEl(iEl,iPC)] += f1*gSJ[mesh.indGSJ(iEl,l,iPC)];
+      minJac = std::min(minJac, sJ[l]);
+      maxJac = std::max(maxJac, sJ[l]);
+    }
+  }
+  return true;
+}
+
 
 
 // Contribution of the vertex distance to the objective function value and gradients (2D version)
@@ -89,9 +113,12 @@ void OptHOM::evalObjGrad(const alglib::real_1d_array &x, double &Obj, alglib::re
 
   addJacObjGrad(Obj, gradObj);
   addDistObjGrad(lambda, lambda2, Obj, gradObj);
+  if(_optimizeMetricMin)
+    addMetricMinObjGrad(Obj, gradObj);
+  
 
   if ((minJac > barrier_min) && (maxJac < barrier_max)) {
-    printf("INFO: reached Jacobian requirements, setting null gradient\n");
+    printf("INFO: reached %s (%g %g) requirements, setting null gradient\n", _optimizeMetricMin ? "svd" : "jacobian", minJac, maxJac);
     Obj = 0.;
     for (int i = 0; i < gradObj.length(); i++) gradObj[i] = 0.;
   }
@@ -130,12 +157,13 @@ void OptHOM::recalcJacDist()
     std::vector<double> sJ(mesh.nBezEl(iEl));                       // Scaled Jacobians
     std::vector<double> dumGSJ(mesh.nBezEl(iEl)*mesh.nPCEl(iEl));   // (Dummy) gradients of scaled Jacobians
     mesh.scaledJacAndGradients (iEl,sJ,dumGSJ);
+    if(_optimizeMetricMin)
+      mesh.metricMinAndGradients (iEl,sJ,dumGSJ);
     for (int l = 0; l < mesh.nBezEl(iEl); l++) {
       minJac = std::min(minJac, sJ[l]);
       maxJac = std::max(maxJac, sJ[l]);
     }
   }
-
 }
 
 
@@ -236,7 +264,7 @@ void OptHOM::OptimPass(alglib::real_1d_array &x, const alglib::real_1d_array &in
 
 
 
-int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double b_max, int pInt, int itMax)
+int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double b_max, bool optimizeMetricMin, int pInt, int itMax)
 {
 
   barrier_min = b_min;
@@ -245,6 +273,7 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double
 //  powM = 4;
 //  powP = 3;
 
+  _optimizeMetricMin = optimizeMetricMin;
   // Set weights & length scale for non-dimensionalization
   lambda = weightFixed;
   lambda2 = weightFree;

contrib/HighOrderMeshOptimizer/OptHOM.h

@@ -25,7 +25,7 @@ public:
   // 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
-  int optimize(double lambda, double lambda2, double barrier_min, double barrier_max, int pInt, int itMax);  // optimize one list of elements
+  int optimize(double lambda, double lambda2, double barrier_min, double barrier_max, bool optimizeMetricMin, int pInt, int itMax);  // optimize one list of elements
   void recalcJacDist();
   inline void getJacDist(double &minJ, double &maxJ, double &maxD, double &avgD);
   void updateMesh(const alglib::real_1d_array &x);
@@ -39,6 +39,7 @@ private:
 //  double lambda, lambda2, powM, powP, invLengthScaleSq;
   double lambda, lambda2, jacBar, invLengthScaleSq;
   int iter, progressInterv;            // Current iteration, interval of iterations for reporting
+  bool _optimizeMetricMin;
   double initObj, initMaxDist, initAvgDist;  // Values for reporting
   double minJac, maxJac, maxDist, avgDist;  // Values for reporting
 
@@ -46,6 +47,7 @@ private:
   inline double compute_f(double v);
   inline double compute_f1(double v);
   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);
   void calcScale(alglib::real_1d_array &scale);
   void OptimPass(alglib::real_1d_array &x, const alglib::real_1d_array &initGradObj, int itMax);

contrib/HighOrderMeshOptimizer/OptHomMesh.cpp

@@ -269,6 +269,85 @@ void Mesh::updateGEntityPositions()
 }
 
 
+void Mesh::metricMinAndGradients(int iEl, std::vector<double> &lambda , std::vector<double> &gradLambda)
+{
+  fullMatrix<double> &gsf = _gradShapeFunctions[_el[iEl]->getTypeForMSH()];
+  //const fullMatrix<double> &l2b = _lag2Bez[_el[iEl]->getTypeForMSH()];
+  const int nbBez = _nBezEl[iEl];
+  const int nbNod = _nNodEl[iEl];
+  fullVector<double> lambdaJ(nbBez), lambdaB(nbBez);
+  fullMatrix<double> gradLambdaJ(nbBez, 2 * nbNod);
+  fullMatrix<double> gradLambdaB(nbBez, 2 * nbNod);
+
+  // jacobian of the straight elements (only triangles for now)
+  SPoint3 *IXYZ[3] = {&_ixyz[_el2V[iEl][0]], &_ixyz[_el2V[iEl][1]], &_ixyz[_el2V[iEl][2]]};
+  double jaci[2][2] = {
+    {IXYZ[1]->x() - IXYZ[0]->x(), IXYZ[2]->x() - IXYZ[0]->x()},
+    {IXYZ[1]->y() - IXYZ[0]->y(), IXYZ[2]->y() - IXYZ[0]->y()}
+  };
+  double invJaci[2][2];
+  inv2x2(jaci, invJaci);
+  
+  for (int l = 0; l < nbBez; l++) {
+    fullMatrix<double> dPsi(gsf, l * 3, 3);
+    double jac[2][2] = {{0., 0.}, {0., 0.}};
+    for (int i = 0; i < nbNod; i++) {
+      int &iVi = _el2V[iEl][i];
+      const double dpsidx = dPsi(i, 0) * invJaci[0][0] + dPsi(i, 1) * invJaci[1][0];
+      const double dpsidy = dPsi(i, 0) * invJaci[0][1] + dPsi(i, 1) * invJaci[1][1];
+      jac[0][0] += _xyz[iVi].x() * dpsidx;
+      jac[0][1] += _xyz[iVi].x() * dpsidy;
+      jac[1][0] += _xyz[iVi].y() * dpsidx;
+      jac[1][1] += _xyz[iVi].y() * dpsidy;
+    }
+    const double dxdx = jac[0][0] * jac[0][0] + jac[0][1] * jac[0][1];
+    const double dydy = jac[1][0] * jac[1][0] + jac[1][1] * jac[1][1];
+    const double dxdy = jac[0][0] * jac[1][0] + jac[0][1] * jac[1][1];
+    const double sqr = sqrt((dxdx - dydy) * (dxdx - dydy) + 4 * dxdy * dxdy);
+    const double osqr = sqr > 1e-8 ? 1/sqr : 0;
+    lambdaJ(l) = 0.5 * (dxdx + dydy - sqr);
+    const double axx = (1 - (dxdx - dydy) * osqr) * jac[0][0] - 2 * dxdy * osqr * jac[1][0];
+    const double axy = (1 - (dxdx - dydy) * osqr) * jac[0][1] - 2 * dxdy * osqr * jac[1][1];
+    const double ayx = -2 * dxdy * osqr * jac[0][0] + (1 - (dydy - dxdx) * osqr) * jac[1][0];
+    const double ayy = -2 * dxdy * osqr * jac[0][1] + (1 - (dydy - dxdx) * osqr) * jac[1][1];
+    const double axixi   = axx * invJaci[0][0] + axy * invJaci[0][1];
+    const double aetaeta = ayx * invJaci[1][0] + ayy * invJaci[1][1];
+    const double aetaxi  = ayx * invJaci[0][0] + ayy * invJaci[0][1];
+    const double axieta  = axx * invJaci[1][0] + axy * invJaci[1][1];
+    for (int i = 0; i < nbNod; i++) {
+      gradLambdaJ(l, i + 0 * nbNod) = axixi * dPsi(i, 0) + axieta * dPsi(i, 1);
+      gradLambdaJ(l, i + 1 * nbNod) = aetaxi * dPsi(i, 0) + aetaeta * dPsi(i, 1);
+    }
+  }
+  
+  //l2b.mult(lambdaJ, lambdaB);
+  //l2b.mult(gradLambdaJ, gradLambdaB);
+  lambdaB = lambdaJ;
+  gradLambdaB = gradLambdaJ;
+
+  int iPC = 0;
+  std::vector<SPoint3> gXyzV(nbBez);
+  std::vector<SPoint3> gUvwV(nbBez);
+  for (int l = 0; l < nbBez; l++) {
+    lambda[l] = lambdaB(l);
+  }
+  for (int i = 0; i < nbNod; i++) {
+    int &iFVi = _el2FV[iEl][i];
+    if (iFVi >= 0) {
+      for (int l = 0; l < nbBez; l++) {
+        gXyzV [l] = SPoint3(gradLambdaB(l,i+0*nbNod),gradLambdaB(l,i+1*nbNod),/*BDB(l,i+2*nbNod)*/ 0.);
+      }
+      _pc->gXyz2gUvw(_freeVert[iFVi],_uvw[iFVi],gXyzV,gUvwV);
+      for (int l = 0; l < nbBez; l++) {
+        gradLambda[indGSJ(iEl,l,iPC)] = gUvwV[l][0];
+        if (_nPCFV[iFVi] >= 2) gradLambda[indGSJ(iEl,l,iPC+1)] = gUvwV[l][1];
+        if (_nPCFV[iFVi] == 3) gradLambda[indGSJ(iEl,l,iPC+2)] = gUvwV[l][2];
+      }
+      iPC += _nPCFV[iFVi];
+    }
+  }
+}
+
 /*
   A faster version that computes jacobians and their gradients
 

contrib/HighOrderMeshOptimizer/OptHomMesh.h

@@ -30,6 +30,7 @@ public:
   inline const int &indPCEl(int iEl, int iPC) { return _indPCEl[iEl][iPC]; }
   inline const int &nBezEl(int iEl) { return _nBezEl[iEl]; }
 
+  void metricMinAndGradients(int iEl, std::vector<double> &sJ, std::vector<double> &gSJ);
   void scaledJacAndGradients(int iEl, std::vector<double> &sJ, std::vector<double> &gSJ);   
   inline int indGSJ(int iEl, int l, int iPC) { return iPC*_nBezEl[iEl]+l; }
 

contrib/HighOrderMeshOptimizer/OptHomRun.cpp

@@ -457,7 +457,11 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
         OptHomMessage("Optimizing a blob %i/%i composed of %4d elements", i+1, toOptimize.size(), toOptimize[i].first.size());
         fflush(stdout);
         OptHOM temp(&entity, toOptimize[i].first, toOptimize[i].second, method);
-        int success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, samples, p.itMax);
+        int success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, false, samples, p.itMax);
+        if (success >= 0 && p.BARRIER_MIN_METRIC > 0) {
+          OptHomMessage("jacobian optimization succeed, starting svd optimization");
+          success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC, p.BARRIER_MAX, true, samples, p.itMax);
+        }
         temp.mesh.updateGEntityPositions();
         if (success <= 0) {
           std::ostringstream ossI2;
@@ -508,7 +512,7 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
           temp.mesh.writeMSH(ossI.str().c_str());
           if (minJac > p.BARRIER_MIN && maxJac < p.BARRIER_MAX) break;
 
-          p.SUCCESS = std::min(p.SUCCESS,temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, samples, p.itMax));
+          p.SUCCESS = std::min(p.SUCCESS,temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, false, samples, p.itMax));
 
           //  temp.recalcJacDist();
           //  temp.getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
@@ -556,7 +560,7 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
           ossI << "initial_" << (*itr)->tag() << "ITER_" << ITER << ".msh";
           temp.mesh.writeMSH(ossI.str().c_str());
           if (minJac > p.BARRIER_MIN  && maxJac < p.BARRIER_MAX) break;
-          p.SUCCESS = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, samples, p.itMax);
+          p.SUCCESS = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, false, samples, p.itMax);
           temp.recalcJacDist();
           temp.getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
           temp.mesh.updateGEntityPositions();

contrib/HighOrderMeshOptimizer/OptHomRun.h

@@ -5,6 +5,7 @@ class GModel;
 
 struct OptHomParameters {
   // INPUT ------> 
+  double BARRIER_MIN_METRIC ; // minimum scaled jcaobian
   double BARRIER_MIN ; // minimum scaled jcaobian
   double BARRIER_MAX ; // maximum scaled jcaobian
   double weightFixed ; // weight of the energy for fixed nodes
@@ -28,7 +29,7 @@ struct OptHomParameters {
   
   OptHomParameters () 
   // default values    
-  : BARRIER_MIN (0.1), BARRIER_MAX (2.0) , weightFixed (1.e6),  weightFree (1.e2),
+  : BARRIER_MIN (0.1), BARRIER_MAX (2.0) , BARRIER_MIN_METRIC (-1.), weightFixed (1.e6),  weightFree (1.e2),
     nbLayers (6) , dim(3) , itMax(10000), onlyVisible(true), DistanceFactor(12), method(1), filter(1)
   {
   }