From 833ee6bdb414cea384e5634f9ffc809b4510509a Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Fri, 20 Apr 2012 14:27:19 +0000
Subject: [PATCH] Fixed scaling problem by preconditioning with scales based on
parametrization gradients
---
contrib/HighOrderMeshOptimizer/OptHOM.cpp | 83 ++++++++++++++-----
contrib/HighOrderMeshOptimizer/OptHOM.h | 4 +-
contrib/HighOrderMeshOptimizer/OptHomMesh.cpp | 19 +++++
contrib/HighOrderMeshOptimizer/OptHomMesh.h | 2 +
4 files changed, 83 insertions(+), 25 deletions(-)
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.cpp b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
index 355a0c443f..18fc68ee70 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
@@ -149,45 +149,82 @@ void printProgressFunc(const alglib::real_1d_array &x, double Obj, void *HOInst)
+void OptHOM::calcScale(alglib::real_1d_array &scale)
+{
+
+ scale.setlength(mesh.nPC());
+
+ // Calculate scale
+ for (int iFV = 0; iFV < mesh.nFV(); iFV++) {
+ std::vector<double> scaleFV(mesh.nPCFV(iFV),1.);
+ mesh.pcScale(iFV,scaleFV);
+ for (int iPC = 0; iPC < mesh.nPCFV(iFV); iPC++) scale[mesh.indPCFV(iFV,iPC)] = scaleFV[iPC];
+ }
+
+ // Normalize scale vector (otherwise ALGLIB routines may fail)
+ double scaleNormSq = 0.;
+ for (int i = 0; i < mesh.nPC(); i++) scaleNormSq += scale[i]*scale[i];
+ const double scaleNorm = sqrt(scaleNormSq);
+ for (int i = 0; i < mesh.nPC(); i++) scale[i] /= scaleNorm;
+
+}
+
+
+
+
void OptHOM::OptimPass(alglib::real_1d_array &x, const alglib::real_1d_array &initGradObj, int itMax)
{
static const double EPSG = 0.;
static const double EPSF = 0.;
static const double EPSX = 0.;
+ static int OPTMETHOD = 1;
std::cout << "--- Optimization pass with jacBar = " << jacBar << ", lambda = " << lambda << ", lambda2 = " << lambda2 << std::endl;
-// alglib::minlbfgsstate state;
-// alglib::minlbfgsreport rep;
- alglib::mincgstate state;
- alglib::mincgreport rep;
-
-// minlbfgscreate(3, x, state);
-// minlbfgssetcond(state, EPSG, EPSF, EPSX, itMax);
-// minlbfgssetxrep(state, true);
- mincgcreate(x, state);
- mincgsetcond(state, EPSG, EPSF, EPSX, itMax);
- mincgsetxrep(state, true);
-
iter = 0;
-// alglib::minlbfgsoptimize(state, evalObjGradFunc, printProgressFunc, this);
- alglib::mincgoptimize(state, evalObjGradFunc, printProgressFunc, this);
-
-// minlbfgsresults(state, x, rep);
- mincgresults(state, x, rep);
+ int iterationscount = 0, nfev = 0, terminationtype = -1;
+ if (OPTMETHOD == 1) {
+ alglib::mincgstate state;
+ alglib::mincgreport rep;
+ mincgcreate(x, state);
+ alglib::real_1d_array scale;
+ calcScale(scale);
+ mincgsetscale(state,scale);
+ mincgsetprecscale(state);
+ mincgsetcond(state, EPSG, EPSF, EPSX, itMax);
+ mincgsetxrep(state, true);
+ alglib::mincgoptimize(state, evalObjGradFunc, printProgressFunc, this);
+ mincgresults(state, x, rep);
+ iterationscount = rep.iterationscount;
+ nfev = rep.nfev;
+ terminationtype = rep.terminationtype;
+ }
+ else {
+ alglib::minlbfgsstate state;
+ alglib::minlbfgsreport rep;
+ minlbfgscreate(3, x, state);
+ alglib::real_1d_array scale;
+ calcScale(scale);
+ minlbfgssetscale(state,scale);
+ minlbfgssetprecscale(state);
+ minlbfgssetcond(state, EPSG, EPSF, EPSX, itMax);
+ minlbfgssetxrep(state, true);
+ alglib::minlbfgsoptimize(state, evalObjGradFunc, printProgressFunc, this);
+ minlbfgsresults(state, x, rep);
+ iterationscount = rep.iterationscount;
+ nfev = rep.nfev;
+ terminationtype = rep.terminationtype;
+ }
- std::cout << "Optimization finalized after " << rep.iterationscount << " iterations (" << rep.nfev << " functions evaluations)";
- 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;
+ std::cout << "Optimization finalized after " << iterationscount << " iterations (" << nfev << " functions evaluations)";
+ switch(int(terminationtype)) {
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;
+ default: std::cout << " with code " << int(terminationtype); break;
}
std::cout << "." << std::endl;
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.h b/contrib/HighOrderMeshOptimizer/OptHOM.h
index e16620958d..29880cfb16 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.h
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.h
@@ -39,11 +39,12 @@ private:
double initObj, initMaxDist, initAvgDist; // Values for reporting
double minJac, maxJac, maxDist, avgDist; // Values for reporting
- inline double setBarrierTerm(double jacBarrier) { bTerm = jacBarrier/(1.-jacBarrier); }
+ inline void setBarrierTerm(double jacBarrier) { bTerm = jacBarrier/(1.-jacBarrier); }
inline double compute_f(double v);
inline double compute_f1(double v);
bool addJacObjGrad(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);
};
@@ -83,7 +84,6 @@ inline double OptHOM::compute_f1(double v)
void OptHOM::getJacDist(double &minJ, double &maxJ, double &maxD, double &avgD)
{
minJ = minJac; maxJ = maxJac; maxD = maxDist; avgD = avgDist;
- std::cout << "getJacDist: minJac = " << minJac << ", maxJac = " << maxJac << ", maxDist = " << maxDist << ", avgDist = " << avgDist << std::endl;
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index 58885217fc..3b03ad6ecf 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -440,6 +440,25 @@ void Mesh::gradScaledJac(int iEl, std::vector<double> &gSJ)
+void Mesh::pcScale(int iFV, std::vector<double> &scale)
+{
+
+ // Calc. derivative of x, y & z w.r.t. parametric coordinates
+ const SPoint3 dX(1.,0.,0.), dY(0.,1.,0.), dZ(0.,0.,1.);
+ SPoint3 gX, gY, gZ;
+ _pc->gXyz2gUvw(_freeVert[iFV],_uvw[iFV],dX,gX);
+ _pc->gXyz2gUvw(_freeVert[iFV],_uvw[iFV],dY,gY);
+ _pc->gXyz2gUvw(_freeVert[iFV],_uvw[iFV],dZ,gZ);
+
+ // Scale = inverse norm. of vector (dx/du, dy/du, dz/du)
+ scale[0] = 1./sqrt(gX[0]*gX[0]+gY[0]*gY[0]+gZ[0]*gZ[0]);
+ if (_nPCFV[iFV] >= 2) scale[1] = 1./sqrt(gX[1]*gX[1]+gY[1]*gY[1]+gZ[1]*gZ[1]);
+ if (_nPCFV[iFV] == 3) scale[2] = 1./sqrt(gX[2]*gX[2]+gY[2]*gY[2]+gZ[2]*gZ[2]);
+
+}
+
+
+
void Mesh::writeMSH(const char *filename)
{
FILE *f = fopen(filename, "w");
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.h b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
index 20329df038..2e2741fbb4 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
@@ -37,6 +37,8 @@ public:
inline double distSq(int iFV);
inline void gradDistSq(int iV, std::vector<double> &gDSq);
+ void pcScale(int iFV, std::vector<double> &scale); // Calc. scale of parametric coordinates for preconditioning
+
void getUvw(double *it);
void updateMesh(const double *it);
void distSqToStraight(std::vector<double> &dSq);
--
GitLab