diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.cpp b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
index 0197263f11b33680aa8e0d3802c9196916890cb0..6c9582c26a0a79f130a77fdb8e296fb74f9d7d41 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
@@ -25,6 +25,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);
@@ -35,6 +38,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]);
}
}
@@ -48,13 +53,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;
@@ -70,8 +84,12 @@ 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_min) && (maxJac < barrier_max))
+ printf("INFO: reached Jacobian requirements, setting null gradient\n");
+ else {
+ addJacObjGrad(Obj, gradObj);
+ addDistObjGrad(lambda, lambda2, Obj, gradObj);
+ }
}
@@ -79,38 +97,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_min && maxJac < HO.barrier_max) {
- 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);
@@ -130,11 +140,8 @@ 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);
- Msg::Debug("--> Iteration %3d --- OBJ %12.5E (relative decrease = %12.5E) -- minJ = %12.5E maxJ = %12.5E Max D = %12.5E Avg D = %12.5E", iter, Obj, Obj/initObj, minJ, maxJ, maxD, avgD);
+ 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);
+ Msg::Debug("--> Iteration %3d --- OBJ %12.5E (relative decrease = %12.5E) -- minJ = %12.5E maxJ = %12.5E Max D = %12.5E Avg D = %12.5E", iter, Obj, Obj/initObj, minJac, maxJac, maxDist, avgDist);
}
}
@@ -148,60 +155,89 @@ 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 = 1.e-8;
-// static const double EPSF = 0.;
+ static const double EPSF = 0.;
static const double EPSX = 0.;
-// const double EPSX = x.length()*1.e-4/sqrt(invLengthScaleSq);
-// std::cout << "DEBUG: EPSX = " << EPSX << ", EPSX/x.length() = " << EPSX/x.length() << std::endl;
-
-// double iGONorm = 0;
-// for (int i=0; i<initGradObj.length(); i++) iGONorm += initGradObj[i]*initGradObj[i];
-// const double EPSG = 1.e-2*sqrt(iGONorm);
+ static int OPTMETHOD = 1;
Msg::Debug("--- Optimization pass with jacBar = %12.5E",jacBar);
-// 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);
-
- // 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;
-// 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;
-// }
-// std::cout << "." << std::endl;
+ 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 " << 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;
+ default: std::cout << " with code " << int(terminationtype); break;
+ }
+ std::cout << "." << std::endl;
}
-int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double b_max, int pInt, int itMax, double &minJ, double &maxJ)
+int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double b_max, int pInt, int itMax)
{
barrier_min = b_min;
barrier_max = b_max;
@@ -223,14 +259,13 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double
mesh.getUvw(x.getcontent());
// Calculate initial performance
- // double minJ, maxJ;
- double initMaxD, initAvgD;
- 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);
- // std::cout << "DEBUG: jacBarStart = " << jacBarStart << std::endl;
// Calculate initial objective function value and gradient
initObj = 0.;
@@ -240,14 +275,13 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double
evalObjGrad(x, initObj, gradObj);
- // 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;
int ITER = 0;
- while (minJ < barrier_min) {
+ while (minJac < barrier_min) {
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);
if (ITER ++ > 5) break;
}
@@ -257,12 +291,11 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double
// OptimPass(x, gradObj, itMax);
// }
- double finalObj = 0., finalMaxD, finalAvgD;
+ double finalObj = 0.;
evalObjGrad(x, finalObj, gradObj);
- getDistances(finalMaxD, finalAvgD, minJ, maxJ);
- Msg::Info("Optimization finished : Avg distance to bnd = %12.5E MinJac %12.5E MaxJac %12.5E",finalAvgD,minJ,maxJ);
+ Msg::Info("Optimization finished : Avg distance to bnd = %12.5E MinJac %12.5E MaxJac %12.5E",avgDist,minJac,maxJac);
- return (minJ > barrier_min && maxJ < barrier_max );
+ return (minJac > barrier_min && maxJac < barrier_max );
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.h b/contrib/HighOrderMeshOptimizer/OptHOM.h
index b05b5d8fa679ccfc4a40cbc4547b4848af8844f2..f2b18530e2ff576fb92cc7e2ac9055b41790e98a 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.h
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.h
@@ -22,12 +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_min, double barrier_max,int pInt, int itMax, double &minJ, double &maxJ); // optimize one list of elements
- // OptHOM(GEntity *gf, const std::set<MElement*> &els, std::set<MVertex*> & toFix, int method);
- // 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
+ int optimize(double lambda, double lambda2, double barrier_min, double barrier_max, 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);
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);
@@ -84,10 +81,10 @@ 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;
-//}
+void OptHOM::getJacDist(double &minJ, double &maxJ, double &maxD, double &avgD)
+{
+ minJ = minJac; maxJ = maxJac; maxD = maxDist; avgD = avgDist;
+}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index 7ea06cfa7eba5a7f57e9a81c89dc4b94ec703c4f..6f1f46a5b70748101a584d65c4c044e9cf64012a 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -5,11 +5,14 @@
#include "MTetrahedron.h"
#include "ParamCoord.h"
#include "OptHomMesh.h"
-#include "JacobianBasis.h"
+
+
std::map<int, std::vector<double> > Mesh::_jacBez;
-static std::vector<double> _computeJB(const polynomialBasis *lagrange, const bezierBasis *bezier)
+
+
+std::vector<double> Mesh::computeJB(const polynomialBasis *lagrange, const bezierBasis *bezier)
{
int nbNodes = lagrange->points.size1();
@@ -37,7 +40,7 @@ static std::vector<double> _computeJB(const polynomialBasis *lagrange, const bez
for (int j = 0; j < nbNodes; j++) {
double Jij = dPsi(i, 0) * dPsi(j, 1) - dPsi(i, 1) * dPsi(j,0);
for (int l = 0; l < bezier->points.size1(); l++) {
- JB[(l * nbNodes + i) * nbNodes + j] += bezier->matrixLag2Bez(l, k) * Jij;
+ JB[indJB2DBase(nbNodes,l,i,j)] += bezier->matrixLag2Bez(l, k) * Jij;
}
}
}
@@ -51,7 +54,7 @@ static std::vector<double> _computeJB(const polynomialBasis *lagrange, const bez
+ (dPsi(j, 2) * dPsi(m, 0) - dPsi(j, 0) * dPsi(m, 2)) * dPsi(i, 1)
+ (dPsi(j, 0) * dPsi(m, 1) - dPsi(j, 1) * dPsi(m, 0)) * dPsi(i, 2);
for (int l = 0; l < bezier->points.size1(); l++) {
- JB[(l * nbNodes + (i * nbNodes + j)) * nbNodes + m] += bezier->matrixLag2Bez(l, k) * Jijm;
+ JB[indJB3DBase(nbNodes,l,i,j,m)] += bezier->matrixLag2Bez(l, k) * Jijm;
}
}
}
@@ -61,6 +64,8 @@ static std::vector<double> _computeJB(const polynomialBasis *lagrange, const bez
return JB;
}
+
+
Mesh::Mesh(GEntity *ge, std::set<MVertex*> &toFix, int method) :
_ge(ge)
{
@@ -101,7 +106,7 @@ Mesh::Mesh(GEntity *ge, std::set<MVertex*> &toFix, int method) :
const polynomialBasis *lagrange = el->getFunctionSpace();
const bezierBasis *bezier = JacobianBasis::find(lagrange->type)->bezier;
if (_jacBez.find(lagrange->type) == _jacBez.end()) {
- _jacBez[lagrange->type] = _computeJB(lagrange, bezier);
+ _jacBez[lagrange->type] = computeJB(lagrange, bezier);
}
_nBezEl[iEl] = bezier->points.size1();
_nNodEl[iEl] = lagrange->points.size1();
@@ -153,6 +158,8 @@ Mesh::Mesh(GEntity *ge, std::set<MVertex*> &toFix, int method) :
}
}
+
+
SVector3 Mesh::getNormalEl(int iEl)
{
@@ -421,8 +428,8 @@ void Mesh::gradScaledJac(int iEl, std::vector<double> &gSJ)
_pc->gXyz2gUvw(_freeVert[iFVi],_uvw[iFVi],gXyzV,gUvwV);
for (int l = 0; l < _nBezEl[iEl]; l++) {
gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][0];
- if (_nPCFV[iFVi] >= 2) gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][1];
- if (_nPCFV[iFVi] == 3) gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][2];
+ if (_nPCFV[iFVi] >= 2) gSJ[indGSJ(iEl,l,iPC+1)] = gUvwV[l][1];
+ if (_nPCFV[iFVi] == 3) gSJ[indGSJ(iEl,l,iPC+2)] = gUvwV[l][2];
}
iPC += _nPCFV[iFVi];
}
@@ -431,6 +438,27 @@ 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/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 88acf1dd7400ed2578c9d5d0400233264a759bd0..022e1b00d7d008d4ad4c37386a23771ed241cc24 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -465,7 +465,8 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
OptHOM *temp = new OptHOM(*itf, toFix, method);
double distMaxBND, distAvgBND, minJac, maxJac;
- temp->getDistances(distMaxBND, distAvgBND, minJac, maxJac);
+ temp->recalcJacDist();
+ temp->getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
if (minJac < 1.e2)Msg::Info("Optimizing a blob of %4d elements minJ %12.5E -- maxJ %12.5E",(*itf)->getNumMeshElements(), minJac, maxJac);
(*itf)->triangles = tris;
(*itf)->quadrangles = quas;
@@ -474,7 +475,8 @@ 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 = temp->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, samples, p.itMax, p.minJac, p.maxJac);
+ p.SUCCESS = temp->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, samples, p.itMax);
+ temp->getJacDist(p.minJac, p.maxJac, distMaxBND, distAvgBND);
temp->mesh.updateGEntityPositions();
if (!p.SUCCESS) break;
ITER ++;
@@ -499,13 +501,13 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
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->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 (minJac > p.BARRIER_MIN && maxJac < p.BARRIER_MAX) continue;
- p.SUCCESS = temp->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, samples, p.itMax, p.minJac, p.maxJac);
+ p.SUCCESS = temp->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, 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 76970347a1ccf8316076fd96407bbe81c8f9ad41..83e3107bde73ba9f7b943d046433f0396c6fbbce 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.h
@@ -9,7 +9,6 @@ struct OptHomParameters {
double BARRIER_MAX ; // maximum scaled jcaobian
double weightFixed ; // weight of the energy for fixed nodes
double weightFree ; // weight of the energy for free nodes
- double STOP;
int nbLayers ; // number of layers taken around a bad element
int dim ; // which dimension to optimize
int itMax ; // max number of iterations in the optimization process
@@ -26,7 +25,7 @@ struct OptHomParameters {
OptHomParameters ()
// default values
- : STOP (0.01) , BARRIER_MIN (0.1), BARRIER_MAX (2.0) , weightFixed (1.e6), weightFree (1.e2),
+ : 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), filter(1)
{
}