diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.cpp b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
index aeeb08df770a0dd17f275208731729e160799084..355a0c443f6728ce802756bca5fb881f02959c54 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
@@ -24,6 +24,9 @@ OptHOM::OptHOM(GEntity *ge, std::set<MVertex*> & toFix, int method) :
bool OptHOM::addJacObjGrad(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
mesh.scaledJac(iEl,sJ);
@@ -34,6 +37,8 @@ bool OptHOM::addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj)
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]);
}
}
@@ -47,13 +52,22 @@ bool OptHOM::addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj)
bool OptHOM::addDistObjGrad(double Fact, double Fact2, double &Obj, alglib::real_1d_array &gradObj)
{
+ maxDist = 0;
+ avgDist = 0;
+ int nbBnd = 0;
+
for (int iFV = 0; iFV < mesh.nFV(); iFV++) {
const double Factor = invLengthScaleSq*(mesh.forced(iFV) ? Fact : Fact2);
- Obj += Factor * mesh.distSq(iFV);
+ const double dSq = mesh.distSq(iFV), dist = sqrt(dSq);
+ Obj += Factor * dSq;
std::vector<double> gDSq(mesh.nPCFV(iFV));
mesh.gradDistSq(iFV,gDSq);
for (int iPC = 0; iPC < mesh.nPCFV(iFV); iPC++) gradObj[mesh.indPCFV(iFV,iPC)] += Factor*gDSq[iPC];
+ maxDist = std::max(maxDist, dist);
+ avgDist += dist;
+ nbBnd++;
}
+ if (nbBnd != 0) avgDist /= nbBnd;
return true;
@@ -69,8 +83,11 @@ void OptHOM::evalObjGrad(const alglib::real_1d_array &x, double &Obj, alglib::re
Obj = 0.;
for (int i = 0; i < gradObj.length(); i++) gradObj[i] = 0.;
- addJacObjGrad(Obj, gradObj);
- addDistObjGrad(lambda, lambda2, Obj, gradObj);
+ if (minJac > barrier) printf("INFO: reached minimum Jacobian requirement, setting null gradient\n");
+ else {
+ addJacObjGrad(Obj, gradObj);
+ addDistObjGrad(lambda, lambda2, Obj, gradObj);
+ }
}
@@ -78,39 +95,30 @@ void OptHOM::evalObjGrad(const alglib::real_1d_array &x, double &Obj, alglib::re
void evalObjGradFunc(const alglib::real_1d_array &x, double &Obj, alglib::real_1d_array &gradObj, void *HOInst)
{
- OptHOM &HO = *static_cast<OptHOM*> (HOInst);
- HO.evalObjGrad(x, Obj, gradObj);
- double distMaxBnd, distAvgBnd, minJac, maxJac;
- HO.getDistances(distMaxBnd, distAvgBnd, minJac, maxJac);
- if (minJac > HO.barrier) {
- printf("STOP\n");
- for (int i = 0; i < gradObj.length(); ++i) {
- gradObj[i] = 0;
- }
- }
+ (static_cast<OptHOM*>(HOInst))->evalObjGrad(x, Obj, gradObj);
}
-void OptHOM::getDistances(double &distMaxBND, double &distAvgBND, double &minJac, double &maxJac)
+void OptHOM::recalcJacDist()
{
- distMaxBND = 0;
- distAvgBND = 0;
+ maxDist = 0;
+ avgDist = 0;
int nbBnd = 0;
for (int iFV = 0; iFV < mesh.nFV(); iFV++) {
if (mesh.forced(iFV)) {
double dSq = mesh.distSq(iFV);
- distMaxBND = std::max(distMaxBND, sqrt(dSq));
- distAvgBND += sqrt(dSq);
+ maxDist = std::max(maxDist, sqrt(dSq));
+ avgDist += sqrt(dSq);
nbBnd++;
}
}
- if (nbBnd != 0) distAvgBND /= nbBnd;
+ if (nbBnd != 0) avgDist /= nbBnd;
- minJac = 1000;
- maxJac = -1000;
+ minJac = 1.e300;
+ maxJac = -1.e300;
for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
std::vector<double> sJ(mesh.nBezEl(iEl)); // Scaled Jacobians
mesh.scaledJac(iEl,sJ);
@@ -128,12 +136,7 @@ void OptHOM::printProgress(const alglib::real_1d_array &x, double Obj)
{
iter++;
-
- if (iter % progressInterv == 0) {
- double maxD, avgD, minJ, maxJ;
- getDistances(maxD, avgD, minJ, maxJ);
- printf("--> Iteration %3d --- OBJ %12.5E (relative decrease = %12.5E) -- minJ = %12.5E maxJ = %12.5E Max D = %12.5E Avg D = %12.5E\n", iter, Obj, Obj/initObj, minJ, maxJ, maxD, avgD);
- }
+ if (iter % progressInterv == 0) printf("--> Iteration %3d --- OBJ %12.5E (relative decrease = %12.5E) -- minJ = %12.5E maxJ = %12.5E Max D = %12.5E Avg D = %12.5E\n", iter, Obj, Obj/initObj, minJac, maxJac, maxDist, avgDist);
}
@@ -192,7 +195,7 @@ void OptHOM::OptimPass(alglib::real_1d_array &x, const alglib::real_1d_array &in
-int OptHOM::optimize(double weightFixed, double weightFree, double barrier_, int pInt, int itMax, double &minJ, double &maxJ)
+int OptHOM::optimize(double weightFixed, double weightFree, double barrier_, int pInt, int itMax)
{
barrier = barrier_;
progressInterv = pInt;
@@ -213,10 +216,11 @@ int OptHOM::optimize(double weightFixed, double weightFree, double barrier_, int
mesh.getUvw(x.getcontent());
// Calculate initial performance
- // double minJ, maxJ;
- getDistances(initMaxD, initAvgD, minJ, maxJ);
+ recalcJacDist();
+ initMaxDist = maxDist;
+ initAvgDist = avgDist;
- const double jacBarStart = (minJ > 0.) ? 0.9*minJ : 1.1*minJ;
+ const double jacBarStart = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
jacBar = jacBarStart;
setBarrierTerm(jacBarStart);
@@ -227,15 +231,14 @@ int OptHOM::optimize(double weightFixed, double weightFree, double barrier_, int
for (int i = 0; i < mesh.nPC(); i++) gradObj[i] = 0.;
evalObjGrad(x, initObj, gradObj);
- std::cout << "Initial mesh: Obj = " << initObj << ", minJ = " << minJ << ", maxJ = " << maxJ << ", maxD = " << initMaxD << ", avgD = " << initAvgD << std::endl;
+ std::cout << "Initial mesh: Obj = " << initObj << ", minJ = " << minJac << ", maxJ = " << maxJac << ", maxD = " << initMaxDist << ", avgD = " << initAvgDist << std::endl;
- std::cout << "Start optimizing with barrier = " << barrier << std::endl;
+ std::cout << "Start optimizing " << mesh.nEl() << " elements (" << mesh.nVert() << " vertices, "
+ << mesh.nFV() << " free vertices, " << mesh.nPC() << " variables) with barrier = " << barrier << std::endl;
- while (minJ < barrier) {
+ while (minJac < barrier) {
OptimPass(x, gradObj, itMax);
- double dumMaxD, dumAvgD;
- getDistances(dumMaxD, dumAvgD, minJ, maxJ);
- jacBar = (minJ > 0.) ? 0.9*minJ : 1.1*minJ;
+ jacBar = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
setBarrierTerm(jacBar);
}
@@ -244,10 +247,9 @@ int OptHOM::optimize(double weightFixed, double weightFree, double barrier_, int
// OptimPass(x, gradObj, itMax);
// }
- double finalObj = 0., finalMaxD, finalAvgD;
+ double finalObj = 0.;
evalObjGrad(x, finalObj, gradObj);
- getDistances(finalMaxD, finalAvgD, minJ, maxJ);
- std::cout << "Final mesh: Obj = " << finalObj << ", maxD = " << finalMaxD << ", avgD = " << finalAvgD << ", minJ = " << minJ << ", maxJ = " << maxJ << std::endl;
+ std::cout << "Final mesh: Obj = " << finalObj << ", maxD = " << maxDist << ", avgD = " << avgDist << ", minJ = " << minJac << ", maxJ = " << maxJac << std::endl;
return 0;
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.h b/contrib/HighOrderMeshOptimizer/OptHOM.h
index 23c7efaab920da8ba8b0ca1ad69adafafdbd5f75..e16620958dcb83f84cd92f83163e7b490ee507da 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.h
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.h
@@ -22,8 +22,9 @@ public:
Mesh mesh;
OptHOM(GEntity *gf, std::set<MVertex*> & toFix, int method);
- void getDistances(double &distMaxBND, double &distAvgBND, double &minJac, double &maxJac);
- int optimize(double lambda, double lambda2, double barrier, int pInt, int itMax, double &minJ, double &maxJ); // optimize one list of elements
+ void recalcJacDist();
+ inline void getJacDist(double &minJ, double &maxJ, double &maxD, double &avgD);
+ int optimize(double lambda, double lambda2, double barrier, int pInt, int itMax); // optimize one list of elements
void updateMesh(const alglib::real_1d_array &x);
void evalObjGrad(const alglib::real_1d_array &x, double &Obj, alglib::real_1d_array &gradObj);
void printProgress(const alglib::real_1d_array &x, double Obj);
@@ -35,7 +36,8 @@ private:
// double lambda, lambda2, powM, powP, invLengthScaleSq;
double lambda, lambda2, jacBar, bTerm, invLengthScaleSq;
int iter, progressInterv; // Current iteration, interval of iterations for reporting
- double initObj, initMaxD, initAvgD; // Values for reporting
+ 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 double compute_f(double v);
@@ -78,6 +80,14 @@ 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;
+}
+
+
+
#endif
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 444544e99e38c8379218f43908792ee93e86f601..1bac92a11bb759c05a7cacef9bdd1e1ca7087e12 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -171,7 +171,7 @@ std::set<MVertex*> filter3D(GRegion *gr, int nbLayers, double _qual)
void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
{
- int samples = 30;
+ int samples = 20;
// int method = Mesh::METHOD_RELAXBND | Mesh::METHOD_PHYSCOORD | Mesh::METHOD_PROJJAC;
// int method = Mesh::METHOD_FIXBND | Mesh::METHOD_PHYSCOORD | Mesh::METHOD_PROJJAC;
@@ -196,13 +196,15 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
(*itf)->quadrangles = quas;
double distMaxBND, distAvgBND, minJac, maxJac;
- temp->getDistances(distMaxBND, distAvgBND, minJac, maxJac);
+ temp->recalcJacDist();
+ temp->getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
// std::ostringstream ossI;
// ossI << "initial_" << (*itf)->tag() << ".msh";
// temp->mesh.writeMSH(ossI.str().c_str());
// printf("--> Mesh Loaded for GFace %d : minJ %12.5E -- maxJ %12.5E\n", (*itf)->tag(), minJac, maxJac);
if (distMaxBND < 1.e-8 * SIZE && minJac > p.BARRIER) continue;
- p.SUCCESS = temp->optimize(p.weightFixed, p.weightFree, p.BARRIER, samples, p.itMax, p.minJac, p.maxJac);
+ p.SUCCESS = temp->optimize(p.weightFixed, p.weightFree, p.BARRIER, samples, p.itMax);
+ temp->getJacDist(p.minJac, p.maxJac, distMaxBND, distAvgBND);
temp->mesh.updateGEntityPositions();
// std::ostringstream ossF;
// ossF << "final_" << (*itf)->tag() << ".msh";
@@ -217,11 +219,13 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
filter3D(*itr, p.nbLayers, p.BARRIER);
OptHOM *temp = new OptHOM(*itr, toFix, method);
double distMaxBND, distAvgBND, minJac, maxJac;
- temp->getDistances(distMaxBND, distAvgBND, minJac, maxJac);
+ temp->recalcJacDist();
+ temp->getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
// temp->mesh.writeMSH("initial.msh");
// printf("--> Mesh Loaded for GRegion %d : minJ %12.5E -- maxJ %12.5E\n", (*itr)->tag(), minJac, maxJac);
if (distMaxBND < 1.e-8 * SIZE && minJac > p.BARRIER) continue;
- p.SUCCESS = temp->optimize(p.weightFixed, p.weightFree, p.BARRIER, samples, p.itMax, p.minJac, p.maxJac);
+ p.SUCCESS = temp->optimize(p.weightFixed, p.weightFree, p.BARRIER, samples, p.itMax);
+ temp->getJacDist(p.minJac, p.maxJac, distMaxBND, distAvgBND);
temp->mesh.updateGEntityPositions();
(*itr)->tetrahedra = tets;
// temp->mesh.writeMSH("final.msh");
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.h b/contrib/HighOrderMeshOptimizer/OptHomRun.h
index 10ad8a3f86390d0327571c97eedfde2ee8b4d399..c77229145d4ff2eea93f2a4343af823c279cb2c6 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.h
@@ -5,7 +5,6 @@ class GModel;
struct OptHomParameters {
// INPUT ------>
- double STOP ; // stop criterion
double BARRIER ; // minimum scaled jcaobian
double weightFixed ; // weight of the energy for fixed nodes
double weightFree ; // weight of the energy for free nodes
@@ -21,8 +20,8 @@ struct OptHomParameters {
OptHomParameters ()
// default values
- : STOP (0.01) , BARRIER (0.1) , weightFixed (10.), weightFree (1.e-3),
- nbLayers (6) , dim(3) , itMax(10000), onlyVisible(true)
+ : BARRIER (0.1) , weightFixed (10.), weightFree (1.e-3),
+ nbLayers (6) , dim(3) , itMax(1000), onlyVisible(true)
{
}
};