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

Modified CMakeFiles and disabled features to fix compilation after r19760

parent 68f652a3
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Numeric/CMakeLists.txt
+ 2
− 2
View file @ 3445a130
......@@ -37,8 +37,8 @@ set(SRC
decasteljau.cpp
mathEvaluator.cpp
Iso.cpp
approximationError.cpp
ConjugateGradients.cpp
# approximationError.cpp
# ConjugateGradients.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
......
contrib/HighOrderMeshOptimizer/OptHOM.cpp
+ 118
− 118
View file @ 3445a130
......@@ -36,7 +36,7 @@ static int NEVAL = 0;
#include "OptHomRun.h"
#include "GmshMessage.h"
#include "GmshConfig.h"
#include "ConjugateGradients.h"
//#include "ConjugateGradients.h"
#include "MLine.h"
#include "MTriangle.h"
#include "GModel.h"
......@@ -142,21 +142,21 @@ bool OptHOM::addJacObjGrad(double &Obj, std::vector<double> &gradObj)
}
bool OptHOM::addApproximationErrorObjGrad(double factor, double &Obj, alglib::real_1d_array &gradObj, simpleFunction<double>& fct)
{
std::vector<double> gradF;
for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
double f;
mesh.approximationErrorAndGradients(iEl, f, gradF, 1.e-6, fct);
Obj += f * factor;
for (size_t i = 0; i < mesh.nPCEl(iEl); ++i){
gradObj[mesh.indPCEl(iEl, i)] += gradF[i] * factor;
}
}
// printf("DIST = %12.5E\n",DISTANCE);
return true;
}
//bool OptHOM::addApproximationErrorObjGrad(double factor, double &Obj, alglib::real_1d_array &gradObj, simpleFunction<double>& fct)
//{
// std::vector<double> gradF;
// for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
// double f;
// mesh.approximationErrorAndGradients(iEl, f, gradF, 1.e-6, fct);
// Obj += f * factor;
// for (size_t i = 0; i < mesh.nPCEl(iEl); ++i){
// gradObj[mesh.indPCEl(iEl, i)] += gradF[i] * factor;
// }
// }
// // printf("DIST = %12.5E\n",DISTANCE);
// return true;
//
//}
static void computeGradSFAtNodes (MElement *e, std::vector<std::vector<SVector3> > &gsf)
{
......@@ -804,12 +804,12 @@ void OptHOM::calcScale(alglib::real_1d_array &scale)
}
}
void OptHOM::OptimPass(std::vector<double> &x, int itMax)
{
Msg::Info("--- In-house Optimization pass with initial jac. range (%g, %g), jacBar = %g",
minJac, maxJac, jacBar);
GradientDescent (evalObjGradFunc, x, this);
}
//void OptHOM::OptimPass(std::vector<double> &x, int itMax)
//{
// Msg::Info("--- In-house Optimization pass with initial jac. range (%g, %g), jacBar = %g",
// minJac, maxJac, jacBar);
// GradientDescent (evalObjGradFunc, x, this);
//}
void OptHOM::OptimPass(alglib::real_1d_array &x, int itMax)
......@@ -880,7 +880,7 @@ void OptHOM::OptimPass(alglib::real_1d_array &x, int itMax)
}
int OptHOM::optimize_inhouse(double weightFixed, double weightFree, double weightCAD, double b_min,
int OptHOM::optimize(double weightFixed, double weightFree, double weightCAD, double b_min,
double b_max, bool optimizeMetricMin, int pInt,
int itMax, int optPassMax, int optCAD, double distanceMax, double tolerance)
{
......@@ -981,100 +981,100 @@ int OptHOM::optimize_inhouse(double weightFixed, double weightFree, double weigh
}
int OptHOM::optimize(double weightFixed, double weightFree, double weightCAD, double b_min,
double b_max, bool optimizeMetricMin, int pInt,
int itMax, int optPassMax, int optCAD, double distanceMax, double tolerance)
{
barrier_min = b_min;
barrier_max = b_max;
distance_max = distanceMax;
progressInterv = pInt;
// powM = 4;
// powP = 3;
_optimizeMetricMin = optimizeMetricMin;
_optimizeCAD = optCAD;
// Set weights & length scale for non-dimensionalization
lambda = weightFixed;
lambda2 = weightFree;
lambda3 = weightCAD;
geomTol = tolerance;
std::vector<double> dSq(mesh.nEl());
mesh.distSqToStraight(dSq);
const double maxDSq = *max_element(dSq.begin(),dSq.end());
if (maxDSq < 1.e-10) { // Length scale for non-dim. distance
std::vector<double> sSq(mesh.nEl());
mesh.elSizeSq(sSq);
const double maxSSq = *max_element(sSq.begin(),sSq.end());
invLengthScaleSq = 1./maxSSq;
}
else invLengthScaleSq = 1./maxDSq;
// Set initial guess
std::vector<double> x(mesh.nPC());
mesh.getUvw(x.data());
// Calculate initial performance
recalcJacDist();
initMaxDist = maxDist;
initAvgDist = avgDist;
const double jacBarStart = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
jacBar = jacBarStart;
_optimizeBarrierMax = false;
// Calculate initial objective function value and gradient
initObj = 0.;
std::vector<double>gradObj(mesh.nPC());
for (int i = 0; i < mesh.nPC(); i++) gradObj[i] = 0.;
evalObjGrad(x, initObj, true, gradObj);
Msg::Info("Start optimizing %d elements (%d vertices, %d free vertices, %d variables) "
"with min barrier %g and max. barrier %g", mesh.nEl(), mesh.nVert(),
mesh.nFV(), mesh.nPC(), barrier_min, barrier_max);
int ITER = 0;
bool minJacOK = true;
while (minJac < barrier_min || (maxDistCAD > distance_max && _optimizeCAD)) {
const double startMinJac = minJac;
OptimPass(x, itMax);
recalcJacDist();
jacBar = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
if (_optimizeCAD) jacBar = std::min(jacBar,barrier_min);
if (ITER ++ > optPassMax) {
minJacOK = (minJac > barrier_min && (maxDistCAD < distance_max || !_optimizeCAD));
break;
}
if (fabs((minJac-startMinJac)/startMinJac) < 0.01) {
Msg::Info("Stagnation in minJac detected, stopping optimization");
minJacOK = false;
break;
}
}
ITER = 0;
if (minJacOK && (!_optimizeMetricMin)) {
_optimizeBarrierMax = true;
jacBar = 1.1 * maxJac;
while (maxJac > barrier_max ) {
const double startMaxJac = maxJac;
OptimPass(x, itMax);
recalcJacDist();
jacBar = 1.1 * maxJac;
if (ITER ++ > optPassMax) break;
if (fabs((maxJac-startMaxJac)/startMaxJac) < 0.01) {
Msg::Info("Stagnation in maxJac detected, stopping optimization");
break;
}
}
}
Msg::Info("Optimization done Range (%g,%g)", minJac, maxJac);
if (minJac > barrier_min && maxJac < barrier_max) return 1;
if (minJac > 0.0) return 0;
return -1;
}
//int OptHOM::optimize_inhouse(double weightFixed, double weightFree, double weightCAD, double b_min,
// double b_max, bool optimizeMetricMin, int pInt,
// int itMax, int optPassMax, int optCAD, double distanceMax, double tolerance)
//{
// barrier_min = b_min;
// barrier_max = b_max;
// distance_max = distanceMax;
// progressInterv = pInt;
//// powM = 4;
//// powP = 3;
//
// _optimizeMetricMin = optimizeMetricMin;
// _optimizeCAD = optCAD;
// // Set weights & length scale for non-dimensionalization
// lambda = weightFixed;
// lambda2 = weightFree;
// lambda3 = weightCAD;
// geomTol = tolerance;
// std::vector<double> dSq(mesh.nEl());
// mesh.distSqToStraight(dSq);
// const double maxDSq = *max_element(dSq.begin(),dSq.end());
// if (maxDSq < 1.e-10) { // Length scale for non-dim. distance
// std::vector<double> sSq(mesh.nEl());
// mesh.elSizeSq(sSq);
// const double maxSSq = *max_element(sSq.begin(),sSq.end());
// invLengthScaleSq = 1./maxSSq;
// }
// else invLengthScaleSq = 1./maxDSq;
//
// // Set initial guess
// std::vector<double> x(mesh.nPC());
// mesh.getUvw(x.data());
//
// // Calculate initial performance
// recalcJacDist();
// initMaxDist = maxDist;
// initAvgDist = avgDist;
//
// const double jacBarStart = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
// jacBar = jacBarStart;
//
// _optimizeBarrierMax = false;
// // Calculate initial objective function value and gradient
// initObj = 0.;
// std::vector<double>gradObj(mesh.nPC());
// for (int i = 0; i < mesh.nPC(); i++) gradObj[i] = 0.;
// evalObjGrad(x, initObj, true, gradObj);
//
// Msg::Info("Start optimizing %d elements (%d vertices, %d free vertices, %d variables) "
// "with min barrier %g and max. barrier %g", mesh.nEl(), mesh.nVert(),
// mesh.nFV(), mesh.nPC(), barrier_min, barrier_max);
//
// int ITER = 0;
// bool minJacOK = true;
//
// while (minJac < barrier_min || (maxDistCAD > distance_max && _optimizeCAD)) {
// const double startMinJac = minJac;
// OptimPass(x, itMax);
// recalcJacDist();
// jacBar = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
// if (_optimizeCAD) jacBar = std::min(jacBar,barrier_min);
// if (ITER ++ > optPassMax) {
// minJacOK = (minJac > barrier_min && (maxDistCAD < distance_max || !_optimizeCAD));
// break;
// }
// if (fabs((minJac-startMinJac)/startMinJac) < 0.01) {
// Msg::Info("Stagnation in minJac detected, stopping optimization");
// minJacOK = false;
// break;
// }
// }
//
// ITER = 0;
// if (minJacOK && (!_optimizeMetricMin)) {
// _optimizeBarrierMax = true;
// jacBar = 1.1 * maxJac;
// while (maxJac > barrier_max ) {
// const double startMaxJac = maxJac;
// OptimPass(x, itMax);
// recalcJacDist();
// jacBar = 1.1 * maxJac;
// if (ITER ++ > optPassMax) break;
// if (fabs((maxJac-startMaxJac)/startMaxJac) < 0.01) {
// Msg::Info("Stagnation in maxJac detected, stopping optimization");
// break;
// }
// }
// }
//
// Msg::Info("Optimization done Range (%g,%g)", minJac, maxJac);
//
// if (minJac > barrier_min && maxJac < barrier_max) return 1;
// if (minJac > 0.0) return 0;
// return -1;
//}
#endif
contrib/HighOrderMeshOptimizer/OptHOM.h
+ 5
− 5
View file @ 3445a130
......@@ -53,9 +53,9 @@ public:
// mesh is invalid : some jacobians cannot be made positive
int optimize(double lambda, double lambda2, double lambda3, double barrier_min, double barrier_max,
bool optimizeMetricMin, int pInt, int itMax, int optPassMax, int optimizeCAD, double optCADDistMax, double tolerance);
int optimize_inhouse(double weightFixed, double weightFree, double weightCAD, double b_min,
double b_max, bool optimizeMetricMin, int pInt,
int itMax, int optPassMax, int optCAD, double distanceMax, double tolerance);
// int optimize_inhouse(double weightFixed, double weightFree, double weightCAD, double b_min,
// double b_max, bool optimizeMetricMin, int pInt,
// int itMax, int optPassMax, int optCAD, double distanceMax, double tolerance);
void recalcJacDist();
inline void getJacDist(double &minJ, double &maxJ, double &maxD, double &avgD);
void updateMesh(const alglib::real_1d_array &x);
......@@ -77,7 +77,7 @@ public:
// true : fixed barrier min + moving barrier max
bool _optimizeCAD; // false : do not minimize the distance between mesh and CAD
// true : minimize the distance between mesh and CAD
bool addApproximationErrorObjGrad(double Fact, double &Obj, alglib::real_1d_array &gradObj, simpleFunction<double>& fct);
// bool addApproximationErrorObjGrad(double Fact, double &Obj, alglib::real_1d_array &gradObj, simpleFunction<double>& fct);
bool addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj);
bool addJacObjGrad(double &Obj, std::vector<double> &);
bool addBndObjGrad (double Fact, double &Obj, alglib::real_1d_array &gradObj);
......@@ -90,7 +90,7 @@ public:
alglib::real_1d_array &gradObj);
void calcScale(alglib::real_1d_array &scale);
void OptimPass(alglib::real_1d_array &x, int itMax);
void OptimPass(std::vector<double> &x, int itMax);
// void OptimPass(std::vector<double> &x, int itMax);
};
void OptHOM::getJacDist(double &minJ, double &maxJ, double &maxD, double &avgD)
......
contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+ 62
− 62
View file @ 3445a130
......@@ -335,68 +335,68 @@ void Mesh::metricMinAndGradients(int iEl, std::vector<double> &lambda,
}
}
void Mesh::approximationErrorAndGradients(int iEl, double &f, std::vector<double> &gradF, double eps,
simpleFunction<double> &fct)
{
std::vector<SPoint3> _xyz_temp;
for (int iV = 0; iV < nVert(); iV++){
_xyz_temp.push_back(SPoint3( _vert[iV]->x(), _vert[iV]->y(), _vert[iV]->z()));
_vert[iV]->setXYZ(_xyz[iV].x(),_xyz[iV].y(),_xyz[iV].z());
}
MElement *element = _el[iEl];
f = approximationError (fct, element);
// FIME
// if (iEl < 1)printf("approx error elem %d = %g\n",iEl,f);
int currentId = 0;
// compute the size of the gradient
// depends on how many dofs exist per vertex (0,1,2 or 3)
for (size_t i = 0; i < element->getNumVertices(); ++i) {
if (_el2FV[iEl][i] >= 0) {// some free coordinates
currentId += _nPCFV[_el2FV[iEl][i]];
}
}
gradF.clear();
gradF.resize(currentId, 0.);
currentId = 0;
for (size_t i = 0; i < element->getNumVertices(); ++i) {
if (_el2FV[iEl][i] >= 0) {// some free coordinates
MVertex *v = element->getVertex(i);
// vertex classified on a model edge
if (_nPCFV[_el2FV[iEl][i]] == 1){
double t = _uvw[_el2FV[iEl][i]].x();
GEdge *ge = (GEdge*)v->onWhat();
SPoint3 p (v->x(),v->y(),v->z());
GPoint d = ge->point(t+eps);
v->setXYZ(d.x(),d.y(),d.z());
double f_d = approximationError (fct, element);
gradF[currentId++] = (f_d-f)/eps;
if (iEl < 1)printf("df = %g\n",(f_d-f)/eps);
v->setXYZ(p.x(),p.y(),p.z());
}
else if (_nPCFV[_el2FV[iEl][i]] == 2){
double uu = _uvw[_el2FV[iEl][i]].x();
double vv = _uvw[_el2FV[iEl][i]].y();
GFace *gf = (GFace*)v->onWhat();
SPoint3 p (v->x(),v->y(),v->z());
GPoint d = gf->point(uu+eps,vv);
v->setXYZ(d.x(),d.y(),d.z());
double f_u = approximationError (fct, element);
gradF[currentId++] = (f_u-f)/eps;
d = gf->point(uu,vv+eps);
v->setXYZ(d.x(),d.y(),d.z());
double f_v = approximationError (fct, element);
gradF[currentId++] = (f_v-f)/eps;
v->setXYZ(p.x(),p.y(),p.z());
// if (iEl < 1)printf("df = %g %g\n",(f_u-f)/eps,(f_v-f)/eps);
}
}
}
for (int iV = 0; iV < nVert(); iV++)
_vert[iV]->setXYZ(_xyz_temp[iV].x(),_xyz_temp[iV].y(),_xyz_temp[iV].z());
}
//void Mesh::approximationErrorAndGradients(int iEl, double &f, std::vector<double> &gradF, double eps,
// simpleFunction<double> &fct)
//{
// std::vector<SPoint3> _xyz_temp;
// for (int iV = 0; iV < nVert(); iV++){
// _xyz_temp.push_back(SPoint3( _vert[iV]->x(), _vert[iV]->y(), _vert[iV]->z()));
// _vert[iV]->setXYZ(_xyz[iV].x(),_xyz[iV].y(),_xyz[iV].z());
// }
//
// MElement *element = _el[iEl];
//
// f = approximationError (fct, element);
// // FIME
// // if (iEl < 1)printf("approx error elem %d = %g\n",iEl,f);
// int currentId = 0;
// // compute the size of the gradient
// // depends on how many dofs exist per vertex (0,1,2 or 3)
// for (size_t i = 0; i < element->getNumVertices(); ++i) {
// if (_el2FV[iEl][i] >= 0) {// some free coordinates
// currentId += _nPCFV[_el2FV[iEl][i]];
// }
// }
// gradF.clear();
// gradF.resize(currentId, 0.);
// currentId = 0;
// for (size_t i = 0; i < element->getNumVertices(); ++i) {
// if (_el2FV[iEl][i] >= 0) {// some free coordinates
// MVertex *v = element->getVertex(i);
// // vertex classified on a model edge
// if (_nPCFV[_el2FV[iEl][i]] == 1){
// double t = _uvw[_el2FV[iEl][i]].x();
// GEdge *ge = (GEdge*)v->onWhat();
// SPoint3 p (v->x(),v->y(),v->z());
// GPoint d = ge->point(t+eps);
// v->setXYZ(d.x(),d.y(),d.z());
// double f_d = approximationError (fct, element);
// gradF[currentId++] = (f_d-f)/eps;
// if (iEl < 1)printf("df = %g\n",(f_d-f)/eps);
// v->setXYZ(p.x(),p.y(),p.z());
// }
// else if (_nPCFV[_el2FV[iEl][i]] == 2){
// double uu = _uvw[_el2FV[iEl][i]].x();
// double vv = _uvw[_el2FV[iEl][i]].y();
// GFace *gf = (GFace*)v->onWhat();
// SPoint3 p (v->x(),v->y(),v->z());
// GPoint d = gf->point(uu+eps,vv);
// v->setXYZ(d.x(),d.y(),d.z());
// double f_u = approximationError (fct, element);
// gradF[currentId++] = (f_u-f)/eps;
// d = gf->point(uu,vv+eps);
// v->setXYZ(d.x(),d.y(),d.z());
// double f_v = approximationError (fct, element);
// gradF[currentId++] = (f_v-f)/eps;
// v->setXYZ(p.x(),p.y(),p.z());
// // if (iEl < 1)printf("df = %g %g\n",(f_u-f)/eps,(f_v-f)/eps);
// }
// }
// }
// for (int iV = 0; iV < nVert(); iV++)
// _vert[iV]->setXYZ(_xyz_temp[iV].x(),_xyz_temp[iV].y(),_xyz_temp[iV].z());
//
//}
bool Mesh::bndDistAndGradients(int iEl, double &f , std::vector<double> &gradF, double eps)
......
contrib/HighOrderMeshOptimizer/OptHomMesh.h
+ 6
− 6
View file @ 3445a130
......@@ -38,7 +38,7 @@
#include "ParamCoord.h"
#include "polynomialBasis.h"
#include "simpleFunction.h"
#include "approximationError.h"
//#include "approximationError.h"
class Mesh
{
......@@ -58,7 +58,7 @@ public:
inline const int &nBezEl(int iEl) { return _nBezEl[iEl]; }
int getFreeVertexStartIndex(MVertex* vert);
void approximationErrorAndGradients(int iEl, double &f, std::vector<double> &gradF, double eps, simpleFunction<double> &fct);
// void approximationErrorAndGradients(int iEl, double &f, std::vector<double> &gradF, double eps, simpleFunction<double> &fct);
void metricMinAndGradients(int iEl, std::vector<double> &sJ, std::vector<double> &gSJ);
void scaledJacAndGradients(int iEl, std::vector<double> &sJ, std::vector<double> &gSJ);
bool bndDistAndGradients(int iEl, double &f , std::vector<double> &gradF, double eps);
......@@ -142,10 +142,10 @@ private:
return indJB3DBase(_nNodEl[iEl],l,i,j,m);
}
public:
double approximationErr(int iEl, simpleFunction<double> &f)
{
return approximationError (f, _el[iEl]);
}
// double approximationErr(int iEl, simpleFunction<double> &f)
// {
// return approximationError (f, _el[iEl]);
// }
};
......
contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+ 4
− 4
View file @ 3445a130
......@@ -453,11 +453,11 @@ static void optimizeConnectedBlobs(const vertElVecMap &vertex2elements,
if (temp.mesh.nPC() == 0)
Msg::Info("Blob %i has no degree of freedom, skipping", i+1);
else
success = temp.optimize_inhouse(p.weightFixed, p.weightFree, p.optCADWeight, p.BARRIER_MIN,
success = temp.optimize(p.weightFixed, p.weightFree, p.optCADWeight, p.BARRIER_MIN,
p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax, p.optCAD, p.optCADDistMax, p.discrTolerance);
if (success >= 0 && p.BARRIER_MIN_METRIC > 0) {
Msg::Info("Jacobian optimization succeed, starting svd optimization");
success = temp.optimize_inhouse(p.weightFixed, p.weightFree, p.optCADWeight, p.BARRIER_MIN_METRIC, p.BARRIER_MAX,
success = temp.optimize(p.weightFixed, p.weightFree, p.optCADWeight, p.BARRIER_MIN_METRIC, p.BARRIER_MAX,
true, samples, p.itMax, p.optPassMax, p.optCAD, p.optCADDistMax,p.discrTolerance);
}
double minJac, maxJac, distMaxBND, distAvgBND;
......@@ -539,11 +539,11 @@ static void optimizeOneByOne
std::ostringstream ossI1;
ossI1 << "initial_blob-" << iBadEl << ".msh";
opt->mesh.writeMSH(ossI1.str().c_str());
success = opt->optimize_inhouse(p.weightFixed, p.weightFree, p.optCADWeight, p.BARRIER_MIN,
success = opt->optimize(p.weightFixed, p.weightFree, p.optCADWeight, p.BARRIER_MIN,
p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax, p.optCAD, p.optCADDistMax,p.discrTolerance);
if (success >= 0 && p.BARRIER_MIN_METRIC > 0) {
Msg::Info("Jacobian optimization succeed, starting svd optimization");
success = opt->optimize_inhouse(p.weightFixed, p.weightFree, p.optCADWeight, p.BARRIER_MIN_METRIC,
success = opt->optimize(p.weightFixed, p.weightFree, p.optCADWeight, p.BARRIER_MIN_METRIC,
p.BARRIER_MAX, true, samples, p.itMax, p.optPassMax, p.optCAD, p.optCADDistMax,p.discrTolerance);
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment