diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.cpp b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
index 3d5fc413fd3fec4f94bf857d6a37f2d7a7017a45..a9a0954ccae4aacbde25428a856a07236075d834 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.cpp
+++ b/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;
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.h b/contrib/HighOrderMeshOptimizer/OptHOM.h
index a0865c7ef244c697d8baf313f34e997a0ecb258d..431c7413eb837cfe501d35f9f613a67bf7bef498 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.h
+++ b/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);
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index 23c635348d9d9542958af57c10aa9c9f737757ec..20f7035bc0555ea2e53bdb976e5a9a07c286bb6c 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/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
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.h b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
index bac927e1a3139010324bfe48c3296cd07359f980..61c89a72ca24d57379c0b3082dbc5362ebed01b0 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.h
+++ b/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; }
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index cff74b8f5fa422e27eda6ac7cb881b4c4e54c43e..d6ac67769efda27e1d60f2e3940ee56689843865 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/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();
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.h b/contrib/HighOrderMeshOptimizer/OptHomRun.h
index 5b6cebb9460ca1438c8d79559503dae4e7b14426..66fa29a22e3378158726f199383b005cbe6739d3 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.h
+++ b/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)
{
}