Skip to content
Snippets Groups Projects
Commit 09dd18f1 authored by Thomas Toulorge's avatar Thomas Toulorge
Browse files

Reintroduced in HighOrderMeshOptimizer bug fixes and improvements of r11889,... 

Reintroduced in HighOrderMeshOptimizer bug fixes and improvements of r11889, r11890 and r11896 that had been erased in r11918
parent 27b4b552
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
contrib/HighOrderMeshOptimizer/OptHOM.cpp
+ 112
− 79
View file @ 09dd18f1
...@@ -25,6 +25,9 @@ OptHOM::OptHOM(GEntity *ge, std::set<MVertex*> & toFix, int method) : ...@@ -25,6 +25,9 @@ OptHOM::OptHOM(GEntity *ge, std::set<MVertex*> & toFix, int method) :
bool OptHOM::addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj) bool OptHOM::addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj)
{ {
minJac = 1.e300;
maxJac = -1.e300;
for (int iEl = 0; iEl < mesh.nEl(); iEl++) { for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
std::vector<double> sJ(mesh.nBezEl(iEl)); // Scaled Jacobians std::vector<double> sJ(mesh.nBezEl(iEl)); // Scaled Jacobians
mesh.scaledJac(iEl,sJ); mesh.scaledJac(iEl,sJ);
...@@ -35,6 +38,8 @@ bool OptHOM::addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj) ...@@ -35,6 +38,8 @@ bool OptHOM::addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj)
const double f1 = compute_f1(sJ[l]); const double f1 = compute_f1(sJ[l]);
for (int iPC = 0; iPC < mesh.nPCEl(iEl); iPC++) for (int iPC = 0; iPC < mesh.nPCEl(iEl); iPC++)
gradObj[mesh.indPCEl(iEl,iPC)] += f1*gSJ[mesh.indGSJ(iEl,l,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) ...@@ -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) 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++) { for (int iFV = 0; iFV < mesh.nFV(); iFV++) {
const double Factor = invLengthScaleSq*(mesh.forced(iFV) ? Fact : Fact2); 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)); std::vector<double> gDSq(mesh.nPCFV(iFV));
mesh.gradDistSq(iFV,gDSq); mesh.gradDistSq(iFV,gDSq);
for (int iPC = 0; iPC < mesh.nPCFV(iFV); iPC++) gradObj[mesh.indPCFV(iFV,iPC)] += Factor*gDSq[iPC]; 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; return true;
...@@ -70,8 +84,12 @@ void OptHOM::evalObjGrad(const alglib::real_1d_array &x, double &Obj, alglib::re ...@@ -70,8 +84,12 @@ void OptHOM::evalObjGrad(const alglib::real_1d_array &x, double &Obj, alglib::re
Obj = 0.; Obj = 0.;
for (int i = 0; i < gradObj.length(); i++) gradObj[i] = 0.; for (int i = 0; i < gradObj.length(); i++) gradObj[i] = 0.;
if ((minJac > barrier_min) && (maxJac < barrier_max))
printf("INFO: reached Jacobian requirements, setting null gradient\n");
else {
addJacObjGrad(Obj, gradObj); addJacObjGrad(Obj, gradObj);
addDistObjGrad(lambda, lambda2, Obj, gradObj); addDistObjGrad(lambda, lambda2, Obj, gradObj);
}
} }
...@@ -79,38 +97,30 @@ void OptHOM::evalObjGrad(const alglib::real_1d_array &x, double &Obj, alglib::re ...@@ -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) void evalObjGradFunc(const alglib::real_1d_array &x, double &Obj, alglib::real_1d_array &gradObj, void *HOInst)
{ {
OptHOM &HO = *static_cast<OptHOM*> (HOInst); (static_cast<OptHOM*>(HOInst))->evalObjGrad(x, Obj, gradObj);
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;
}
}
} }
void OptHOM::getDistances(double &distMaxBND, double &distAvgBND, double &minJac, double &maxJac) void OptHOM::recalcJacDist()
{ {
distMaxBND = 0; maxDist = 0;
distAvgBND = 0; avgDist = 0;
int nbBnd = 0; int nbBnd = 0;
for (int iFV = 0; iFV < mesh.nFV(); iFV++) { for (int iFV = 0; iFV < mesh.nFV(); iFV++) {
if (mesh.forced(iFV)) { if (mesh.forced(iFV)) {
double dSq = mesh.distSq(iFV); double dSq = mesh.distSq(iFV);
distMaxBND = std::max(distMaxBND, sqrt(dSq)); maxDist = std::max(maxDist, sqrt(dSq));
distAvgBND += sqrt(dSq); avgDist += sqrt(dSq);
nbBnd++; nbBnd++;
} }
} }
if (nbBnd != 0) distAvgBND /= nbBnd; if (nbBnd != 0) avgDist /= nbBnd;
minJac = 1000; minJac = 1.e300;
maxJac = -1000; maxJac = -1.e300;
for (int iEl = 0; iEl < mesh.nEl(); iEl++) { for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
std::vector<double> sJ(mesh.nBezEl(iEl)); // Scaled Jacobians std::vector<double> sJ(mesh.nBezEl(iEl)); // Scaled Jacobians
mesh.scaledJac(iEl,sJ); mesh.scaledJac(iEl,sJ);
...@@ -130,11 +140,8 @@ void OptHOM::printProgress(const alglib::real_1d_array &x, double Obj) ...@@ -130,11 +140,8 @@ void OptHOM::printProgress(const alglib::real_1d_array &x, double Obj)
iter++; iter++;
if (iter % progressInterv == 0) { if (iter % progressInterv == 0) {
double 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, minJac, maxJac, maxDist, avgDist);
getDistances(maxD, avgD, minJ, maxJ); 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);
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);
} }
} }
...@@ -148,60 +155,89 @@ void printProgressFunc(const alglib::real_1d_array &x, double Obj, void *HOInst) ...@@ -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) void OptHOM::OptimPass(alglib::real_1d_array &x, const alglib::real_1d_array &initGradObj, int itMax)
{ {
static const double EPSG = 0.; 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.; static const double EPSX = 0.;
// const double EPSX = x.length()*1.e-4/sqrt(invLengthScaleSq); static int OPTMETHOD = 1;
// 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);
Msg::Debug("--- Optimization pass with jacBar = %12.5E",jacBar); Msg::Debug("--- Optimization pass with jacBar = %12.5E",jacBar);
// alglib::minlbfgsstate state; iter = 0;
// alglib::minlbfgsreport rep;
int iterationscount = 0, nfev = 0, terminationtype = -1;
if (OPTMETHOD == 1) {
alglib::mincgstate state; alglib::mincgstate state;
alglib::mincgreport rep; alglib::mincgreport rep;
// minlbfgscreate(3, x, state);
// minlbfgssetcond(state, EPSG, EPSF, EPSX, itMax);
// minlbfgssetxrep(state, true);
mincgcreate(x, state); mincgcreate(x, state);
alglib::real_1d_array scale;
calcScale(scale);
mincgsetscale(state,scale);
mincgsetprecscale(state);
mincgsetcond(state, EPSG, EPSF, EPSX, itMax); mincgsetcond(state, EPSG, EPSF, EPSX, itMax);
mincgsetxrep(state, true); mincgsetxrep(state, true);
iter = 0;
// alglib::minlbfgsoptimize(state, evalObjGradFunc, printProgressFunc, this);
alglib::mincgoptimize(state, evalObjGradFunc, printProgressFunc, this); alglib::mincgoptimize(state, evalObjGradFunc, printProgressFunc, this);
// minlbfgsresults(state, x, rep);
mincgresults(state, x, rep); 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)"; std::cout << "Optimization finalized after " << iterationscount << " iterations (" << nfev << " functions evaluations)";
// switch(int(rep.terminationtype)) { switch(int(terminationtype)) {
// case -2: std::cout << ", because rounding errors prevented further improvement"; break; case 1: std::cout << ", because relative function improvement is no more than EpsF"; break;
// case -1: std::cout << ", because incorrect parameters were specified"; break; case 2: std::cout << ", because relative step is no more than EpsX"; break;
// case 1: std::cout << ", because relative function improvement is no more than EpsF"; break; case 4: std::cout << ", because gradient norm is no more than EpsG"; break;
// case 2: std::cout << ", because relative step is no more than EpsX"; break; case 5: std::cout << ", because the maximum number of steps was taken"; break;
// case 4: std::cout << ", because gradient norm is no more than EpsG"; break; default: std::cout << " with code " << int(terminationtype); 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; std::cout << "." << std::endl;
// default: std::cout << " with code " << int(rep.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_min = b_min;
barrier_max = b_max; barrier_max = b_max;
...@@ -223,14 +259,13 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double ...@@ -223,14 +259,13 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double
mesh.getUvw(x.getcontent()); mesh.getUvw(x.getcontent());
// Calculate initial performance // Calculate initial performance
// double minJ, maxJ; recalcJacDist();
double initMaxD, initAvgD; initMaxDist = maxDist;
getDistances(initMaxD, initAvgD, minJ, maxJ); 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; jacBar = jacBarStart;
setBarrierTerm(jacBarStart); setBarrierTerm(jacBarStart);
// std::cout << "DEBUG: jacBarStart = " << jacBarStart << std::endl;
// Calculate initial objective function value and gradient // Calculate initial objective function value and gradient
initObj = 0.; initObj = 0.;
...@@ -240,14 +275,13 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double ...@@ -240,14 +275,13 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double
evalObjGrad(x, initObj, gradObj); 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; int ITER = 0;
while (minJ < barrier_min) { while (minJac < barrier_min) {
OptimPass(x, gradObj, itMax); OptimPass(x, gradObj, itMax);
double dumMaxD, dumAvgD; jacBar = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
getDistances(dumMaxD, dumAvgD, minJ, maxJ);
jacBar = (minJ > 0.) ? 0.9*minJ : 1.1*minJ;
setBarrierTerm(jacBar); setBarrierTerm(jacBar);
if (ITER ++ > 5) break; if (ITER ++ > 5) break;
} }
...@@ -257,12 +291,11 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double ...@@ -257,12 +291,11 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double
// OptimPass(x, gradObj, itMax); // OptimPass(x, gradObj, itMax);
// } // }
double finalObj = 0., finalMaxD, finalAvgD; double finalObj = 0.;
evalObjGrad(x, finalObj, gradObj); 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",avgDist,minJac,maxJac);
Msg::Info("Optimization finished : Avg distance to bnd = %12.5E MinJac %12.5E MaxJac %12.5E",finalAvgD,minJ,maxJ);
return (minJ > barrier_min && maxJ < barrier_max ); return (minJac > barrier_min && maxJac < barrier_max );
} }
......
contrib/HighOrderMeshOptimizer/OptHOM.h
+ 7
− 10
View file @ 09dd18f1
...@@ -22,12 +22,9 @@ public: ...@@ -22,12 +22,9 @@ public:
Mesh mesh; Mesh mesh;
OptHOM(GEntity *gf, std::set<MVertex*> & toFix, int method); 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); // optimize one list of elements
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 void recalcJacDist();
// OptHOM(GEntity *gf, const std::set<MElement*> &els, std::set<MVertex*> & toFix, int method); inline void getJacDist(double &minJ, double &maxJ, double &maxD, double &avgD);
// 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 updateMesh(const alglib::real_1d_array &x);
void evalObjGrad(const alglib::real_1d_array &x, double &Obj, alglib::real_1d_array &gradObj); 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); void printProgress(const alglib::real_1d_array &x, double Obj);
...@@ -84,10 +81,10 @@ inline double OptHOM::compute_f1(double v) ...@@ -84,10 +81,10 @@ inline double OptHOM::compute_f1(double v)
//void OptHOM::getJacDist(double &minJ, double &maxJ, double &maxD, double &avgD) void OptHOM::getJacDist(double &minJ, double &maxJ, double &maxD, double &avgD)
//{ {
// minJ = minJac; maxJ = maxJac; maxD = maxDist; avgD = avgDist; minJ = minJac; maxJ = maxJac; maxD = maxDist; avgD = avgDist;
//} }
......
contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+ 35
− 7
View file @ 09dd18f1
...@@ -5,11 +5,14 @@ ...@@ -5,11 +5,14 @@
#include "MTetrahedron.h" #include "MTetrahedron.h"
#include "ParamCoord.h" #include "ParamCoord.h"
#include "OptHomMesh.h" #include "OptHomMesh.h"
#include "JacobianBasis.h"
std::map<int, std::vector<double> > Mesh::_jacBez; 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(); int nbNodes = lagrange->points.size1();
...@@ -37,7 +40,7 @@ static std::vector<double> _computeJB(const polynomialBasis *lagrange, const bez ...@@ -37,7 +40,7 @@ static std::vector<double> _computeJB(const polynomialBasis *lagrange, const bez
for (int j = 0; j < nbNodes; j++) { for (int j = 0; j < nbNodes; j++) {
double Jij = dPsi(i, 0) * dPsi(j, 1) - dPsi(i, 1) * dPsi(j,0); double Jij = dPsi(i, 0) * dPsi(j, 1) - dPsi(i, 1) * dPsi(j,0);
for (int l = 0; l < bezier->points.size1(); l++) { 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 ...@@ -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, 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); + (dPsi(j, 0) * dPsi(m, 1) - dPsi(j, 1) * dPsi(m, 0)) * dPsi(i, 2);
for (int l = 0; l < bezier->points.size1(); l++) { 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 ...@@ -61,6 +64,8 @@ static std::vector<double> _computeJB(const polynomialBasis *lagrange, const bez
return JB; return JB;
} }
Mesh::Mesh(GEntity *ge, std::set<MVertex*> &toFix, int method) : Mesh::Mesh(GEntity *ge, std::set<MVertex*> &toFix, int method) :
_ge(ge) _ge(ge)
{ {
...@@ -101,7 +106,7 @@ Mesh::Mesh(GEntity *ge, std::set<MVertex*> &toFix, int method) : ...@@ -101,7 +106,7 @@ Mesh::Mesh(GEntity *ge, std::set<MVertex*> &toFix, int method) :
const polynomialBasis *lagrange = el->getFunctionSpace(); const polynomialBasis *lagrange = el->getFunctionSpace();
const bezierBasis *bezier = JacobianBasis::find(lagrange->type)->bezier; const bezierBasis *bezier = JacobianBasis::find(lagrange->type)->bezier;
if (_jacBez.find(lagrange->type) == _jacBez.end()) { if (_jacBez.find(lagrange->type) == _jacBez.end()) {
_jacBez[lagrange->type] = _computeJB(lagrange, bezier); _jacBez[lagrange->type] = computeJB(lagrange, bezier);
} }
_nBezEl[iEl] = bezier->points.size1(); _nBezEl[iEl] = bezier->points.size1();
_nNodEl[iEl] = lagrange->points.size1(); _nNodEl[iEl] = lagrange->points.size1();
...@@ -153,6 +158,8 @@ Mesh::Mesh(GEntity *ge, std::set<MVertex*> &toFix, int method) : ...@@ -153,6 +158,8 @@ Mesh::Mesh(GEntity *ge, std::set<MVertex*> &toFix, int method) :
} }
} }
SVector3 Mesh::getNormalEl(int iEl) SVector3 Mesh::getNormalEl(int iEl)
{ {
...@@ -421,8 +428,8 @@ void Mesh::gradScaledJac(int iEl, std::vector<double> &gSJ) ...@@ -421,8 +428,8 @@ void Mesh::gradScaledJac(int iEl, std::vector<double> &gSJ)
_pc->gXyz2gUvw(_freeVert[iFVi],_uvw[iFVi],gXyzV,gUvwV); _pc->gXyz2gUvw(_freeVert[iFVi],_uvw[iFVi],gXyzV,gUvwV);
for (int l = 0; l < _nBezEl[iEl]; l++) { for (int l = 0; l < _nBezEl[iEl]; l++) {
gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][0]; gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][0];
if (_nPCFV[iFVi] >= 2) gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][1]; if (_nPCFV[iFVi] >= 2) gSJ[indGSJ(iEl,l,iPC+1)] = gUvwV[l][1];
if (_nPCFV[iFVi] == 3) gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][2]; if (_nPCFV[iFVi] == 3) gSJ[indGSJ(iEl,l,iPC+2)] = gUvwV[l][2];
} }
iPC += _nPCFV[iFVi]; iPC += _nPCFV[iFVi];
} }
...@@ -431,6 +438,27 @@ void Mesh::gradScaledJac(int iEl, std::vector<double> &gSJ) ...@@ -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) void Mesh::writeMSH(const char *filename)
{ {
FILE *f = fopen(filename, "w"); FILE *f = fopen(filename, "w");
......
contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+ 8
− 6
View file @ 09dd18f1
...@@ -465,7 +465,8 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p) ...@@ -465,7 +465,8 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
OptHOM *temp = new OptHOM(*itf, toFix, method); OptHOM *temp = new OptHOM(*itf, toFix, method);
double distMaxBND, distAvgBND, minJac, maxJac; 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); 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)->triangles = tris;
(*itf)->quadrangles = quas; (*itf)->quadrangles = quas;
...@@ -474,7 +475,8 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p) ...@@ -474,7 +475,8 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
temp->mesh.writeMSH(ossI.str().c_str()); temp->mesh.writeMSH(ossI.str().c_str());
if (minJac > p.BARRIER_MIN && maxJac < p.BARRIER_MAX) break; 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(); temp->mesh.updateGEntityPositions();
if (!p.SUCCESS) break; if (!p.SUCCESS) break;
ITER ++; ITER ++;
...@@ -499,13 +501,13 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p) ...@@ -499,13 +501,13 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
OptHOM *temp = new OptHOM(*itr, toFix, method); OptHOM *temp = new OptHOM(*itr, toFix, method);
double distMaxBND, distAvgBND, minJac, maxJac; double distMaxBND, distAvgBND, minJac, maxJac;
temp->getDistances(distMaxBND, distAvgBND, minJac, maxJac); temp->recalcJacDist();
// temp->recalcJacDist(); temp->getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
// temp->getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
// temp->mesh.writeMSH("initial.msh"); // temp->mesh.writeMSH("initial.msh");
// printf("--> Mesh Loaded for GRegion %d : minJ %12.5E -- maxJ %12.5E\n", (*itr)->tag(), minJac, maxJac); // 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; 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(); temp->mesh.updateGEntityPositions();
(*itr)->tetrahedra = tets; (*itr)->tetrahedra = tets;
// temp->mesh.writeMSH("final.msh"); // temp->mesh.writeMSH("final.msh");
......
contrib/HighOrderMeshOptimizer/OptHomRun.h
+ 1
− 2
View file @ 09dd18f1
...@@ -9,7 +9,6 @@ struct OptHomParameters { ...@@ -9,7 +9,6 @@ struct OptHomParameters {
double BARRIER_MAX ; // maximum scaled jcaobian double BARRIER_MAX ; // maximum scaled jcaobian
double weightFixed ; // weight of the energy for fixed nodes double weightFixed ; // weight of the energy for fixed nodes
double weightFree ; // weight of the energy for free nodes double weightFree ; // weight of the energy for free nodes
double STOP;
int nbLayers ; // number of layers taken around a bad element int nbLayers ; // number of layers taken around a bad element
int dim ; // which dimension to optimize int dim ; // which dimension to optimize
int itMax ; // max number of iterations in the optimization process int itMax ; // max number of iterations in the optimization process
...@@ -26,7 +25,7 @@ struct OptHomParameters { ...@@ -26,7 +25,7 @@ struct OptHomParameters {
OptHomParameters () OptHomParameters ()
// default values // 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) nbLayers (6) , dim(3) , itMax(10000), onlyVisible(true), DistanceFactor(12), method(1), filter(1)
{ {
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment