diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index 97a3c6fec00afe0e01d5242fcbc40086a0829563..e4449b7522e9d0458c1fcc04bdd4eee30d9873e0 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -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)
diff --git a/Numeric/ConjugateGradients.cpp b/Numeric/ConjugateGradients.cpp
index 9720f14e3e9eea6ff15aab0395b86c81d88b9324..666267c31e9696693be81f02fe3883654d4a8013 100644
--- a/Numeric/ConjugateGradients.cpp
+++ b/Numeric/ConjugateGradients.cpp
@@ -113,7 +113,7 @@ double GradientDescent(void (*func)(std::vector<double> &x,
std::vector<double> dir(N);
double f;
- printf("entering gradient descent (%d unknowns)...\n",N);
+// printf("entering gradient descent (%d unknowns)...\n",N);
for (int iter = 0; iter < MAXIT; iter++){
// compute gradient of func
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.cpp b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
index 0684d5fca1e4f7ac3686e79a2e8179f57a6fa735..eeecc02386cec7c00f32d8b1084fa04a23d1085b 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
@@ -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)
@@ -981,100 +981,100 @@ int OptHOM::optimize(double weightFixed, double weightFree, double weightCAD, do
}
-//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;
-//}
+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
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.h b/contrib/HighOrderMeshOptimizer/OptHOM.h
index 4cb3fc7bf91489a858cdd601e37361a03df551c0..3eaf09d89c35f325f15b60934dc4c530dcd282a5 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.h
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.h
@@ -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)
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index ef905a18a55e48e6d20c46136b2d6d2747f68e6f..dbea55ed27232527a68af49daad108a5ab78548c 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -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)
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.h b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
index 894cde318160544e2ede08e1cb3ecd4674d0bc42..696a19b0c286562d3db711d8cf5d0bf9f473b292 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
@@ -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]);
+ }
};