diff --git a/Fltk/highOrderToolsWindow.cpp b/Fltk/highOrderToolsWindow.cpp
index 6b6a5bcd98a721b030446fd05b5a64d6f0d1cc81..4181ea08b5465fc9fe53d3e5ba24309eb3c31116 100644
--- a/Fltk/highOrderToolsWindow.cpp
+++ b/Fltk/highOrderToolsWindow.cpp
@@ -65,10 +65,7 @@ static void highordertools_runp_cb(Fl_Widget *w, void *data)
int order = (int)o->value[0]->value();
bool linear = !(bool)o->butt[2]->value();
bool incomplete = (bool)o->butt[0]->value();
- bool elastic = (bool)o->butt[3]->value();
- double threshold = o->value[1]->value();
bool onlyVisible = (bool)o->butt[1]->value();
- int nLayers = (int) o->value[2]->value();
if (order == 1)
SetOrder1(GModel::current());
@@ -89,85 +86,31 @@ static void highordertools_runp_cb(Fl_Widget *w, void *data)
static void chooseopti_cb(Fl_Widget *w, void *data){
highOrderToolsWindow *o = FlGui::instance()->highordertools;
- bool elastic = (bool)o->butt[3]->value();
+ int elastic = o->choice[2]->value();
- if (elastic){
- o->butt[4]->value(0);
+ if (elastic == 1){
o->choice[0]->deactivate();
+ o->choice[1]->deactivate();
for (int i=3;i<=6;i++)
o->value[i]->deactivate();
- o->push[1]->deactivate();
+ // o->push[1]->deactivate();
}
else {
- o->butt[3]->value(0);
o->push[0]->activate();
o->choice[0]->activate();
+ o->choice[1]->activate();
for (int i=3;i<=6;i++)
o->value[i]->activate();
- o->push[1]->activate();
+ // o->push[1]->activate();
}
}
-
-static void chooseopti2_cb(Fl_Widget *w, void *data){
- highOrderToolsWindow *o = FlGui::instance()->highordertools;
- bool elastic = !(bool)o->butt[4]->value();
-
- if (elastic){
- o->butt[4]->value(0);
- o->choice[0]->deactivate();
- for (int i=3;i<=6;i++)
- o->value[i]->deactivate();
- o->push[1]->deactivate();
- }
- else {
- o->butt[3]->value(0);
- o->push[0]->activate();
- o->choice[0]->activate();
- for (int i=3;i<=6;i++)
- o->value[i]->activate();
- o->push[1]->activate();
- }
-
-}
-
-static void getMeshInfo (GModel *gm,
- int &meshOrder,
- bool &complete,
- bool &CAD){
- meshOrder = -1;
- CAD = true;
- for (GModel::riter itr = gm->firstRegion(); itr != gm->lastRegion(); ++itr) {
- if ((*itr)->getNumMeshElements()){
- meshOrder = (*itr)->getMeshElement(0)->getPolynomialOrder();
- complete = (meshOrder <= 2) ? 1 : (*itr)->getMeshElement(0)->getNumVolumeVertices();
- break;
- }
- }
- for (GModel::fiter itf = gm->firstFace(); itf != gm->lastFace(); ++itf) {
- printf("%d\n",(*itf)->getNumMeshElements());
- if ((*itf)->getNumMeshElements()){
- if (meshOrder == -1) {
- meshOrder = (*itf)->getMeshElement(0)->getPolynomialOrder();
- complete = (meshOrder <= 2) ? 1 : (*itf)->getMeshElement(0)->getNumFaceVertices();
- if ((*itf)->geomType() == GEntity::DiscreteSurface)CAD = false;
- printf("%d %d %d\n",meshOrder,complete, CAD);
- break;
- }
- }
- }
-}
-
-
static void highordertools_runelas_cb(Fl_Widget *w, void *data)
{
highOrderToolsWindow *o = FlGui::instance()->highordertools;
- int order = (int)o->value[0]->value();
- bool linear = !(bool)o->butt[2]->value();
- bool incomplete = (bool)o->butt[0]->value();
- bool elastic = (bool)o->butt[3]->value();
+ bool elastic = o->choice[2]->value() == 1;
double threshold_min = o->value[1]->value();
bool onlyVisible = (bool)o->butt[1]->value();
int nbLayers = (int) o->value[2]->value();
@@ -180,7 +123,7 @@ static void highordertools_runelas_cb(Fl_Widget *w, void *data)
p.BARRIER_MIN = threshold_min;
p.BARRIER_MAX = threshold_max;
p.onlyVisible = onlyVisible;
- p.dim = GModel::current()->getNumRegions() ? 3 : 2;
+ p.dim = GModel::current()->getDim();//o->meshOrder;
p.itMax = (int) o->value[3]->value();
p.weightFixed = o->value[5]->value();
p.weightFree = o->value[6]->value();
@@ -197,6 +140,7 @@ static void highordertools_runelas_cb(Fl_Widget *w, void *data)
}
+
highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
{
FL_NORMAL_SIZE -= deltaFontSize;
@@ -342,22 +286,21 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
}
{
- y += BH;
- butt[3] = new Fl_Check_Button
- (x,y, IW, BH, "Relocate mesh vertices through Elastic Analogy");
- butt[3]->type(FL_TOGGLE_BUTTON);
- butt[3]->value(0);
- butt[3]->callback(chooseopti_cb);
+ static Fl_Menu_Item menu_method[] = {
+ {"Optimization", 0, 0, 0},
+ {"ElasticAnalogy", 0, 0, 0},
+ {0}
+ };
+
+ y += BH;
+ choice[2] = new Fl_Choice
+ (x,y, IW, BH, "Algorithm");
+ choice[2]->align(FL_ALIGN_RIGHT);
+ choice[2]->menu(menu_method);
+ choice[2]->callback(chooseopti_cb);
}
{
- y += BH;
- butt[4] = new Fl_Check_Button
- (x,y, IW, BH, "Use Optimization for relocating mesh vertices");
- butt[4]->type(FL_TOGGLE_BUTTON);
- butt[4]->value(1);
- butt[4]->callback(chooseopti2_cb);
-
static Fl_Menu_Item menu_objf[] = {
{"CAD + FIXBND", 0, 0, 0},
{"CAD + FREEBND", 0, 0, 0},
@@ -422,9 +365,14 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
}
-
- // g->end();
- // Fl_Group::current()->resizable(g);
+ y += 1.5*BH;
+ messages = new Fl_Browser
+ (x,y, width-2*x, height - y - 2*WB);
+ messages->box(FL_THIN_DOWN_BOX);
+ messages->textfont(FL_COURIER);
+ messages->textsize(FL_NORMAL_SIZE - 1);
+ messages->type(FL_MULTI_BROWSER);
+ // messages->callback(message_browser_cb, this);
}
// win->resizable(o);
@@ -435,7 +383,7 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
void highOrderToolsWindow::show(bool redrawOnly)
{
- getMeshInfo (GModel::current(),meshOrder,complete, CAD);
+ getMeshInfoForHighOrder (GModel::current(),meshOrder,complete, CAD);
if(win->shown() && redrawOnly)
win->redraw();
@@ -443,7 +391,10 @@ void highOrderToolsWindow::show(bool redrawOnly)
value[0]->value(meshOrder);
butt[0]->value(!complete);
if (CAD)output[0]->value("AVAILABLE");
- else output[0]->value("NOT AVAILABLE");
+ else {
+ output[0]->value("NOT AVAILABLE");
+ choice[0]->value(2);
+ }
win->show();
}
}
diff --git a/Fltk/highOrderToolsWindow.h b/Fltk/highOrderToolsWindow.h
index 43862d8daa1c4ff2769d8f59ff4648c8de5ca54c..6e7124d5fe5701488849db32553a9d03c9fe0ead 100644
--- a/Fltk/highOrderToolsWindow.h
+++ b/Fltk/highOrderToolsWindow.h
@@ -12,7 +12,7 @@
#include <FL/Fl_Multi_Browser.H>
#include <FL/Fl_Button.H>
#include <FL/Fl_Check_Button.H>
-#include <FL/Fl_Input.H>
+#include <FL/Fl_Value_Input.H>
#include <FL/Fl_Output.H>
#include "GmshConfig.h"
@@ -26,6 +26,7 @@ class highOrderToolsWindow{
Fl_Choice *choice[20];
Fl_Button *push[20];
Fl_Output *output[20];
+ Fl_Browser *messages;
public:
highOrderToolsWindow(int deltaFontSize=0);
void show(bool redrawOnly);
diff --git a/Fltk/menuWindow.cpp b/Fltk/menuWindow.cpp
index bd6cd213fee29ab3abf03beb4d35e998e9904009..7a8e317a28b87557c35342df0a93ea56e396e518 100644
--- a/Fltk/menuWindow.cpp
+++ b/Fltk/menuWindow.cpp
@@ -1820,7 +1820,8 @@ static void mesh_inspect_cb(Fl_Widget *w, void *data)
Msg::Direct(" Rho: %g", ele->rhoShapeMeasure());
Msg::Direct(" Gamma: %g", ele->gammaShapeMeasure());
Msg::Direct(" Eta: %g", ele->etaShapeMeasure());
- Msg::Direct(" Disto: %g", ele->distoShapeMeasure());
+ double jmin,jmax;ele->scaledJacRange(jmin,jmax);
+ Msg::Direct(" Scaled Jacobian Range : %g %g",jmin,jmax );
CTX::instance()->mesh.changed = ENT_ALL;
drawContext::global()->draw();
FlGui::instance()->showMessages();
diff --git a/Fltk/statisticsWindow.cpp b/Fltk/statisticsWindow.cpp
index 78167989a72703e3f11e2bea41aa2b2ce32bddd3..4efe5f4bc219dee81aeb89792325cc53641a6aef 100644
--- a/Fltk/statisticsWindow.cpp
+++ b/Fltk/statisticsWindow.cpp
@@ -77,8 +77,10 @@ static void statistics_histogram_cb(Fl_Widget *w, void *data)
d[e->getNum()].push_back(e->etaShapeMeasure());
else if(name == "Rho3D")
d[e->getNum()].push_back(e->rhoShapeMeasure());
- else
- d[e->getNum()].push_back(e->distoShapeMeasure());
+ else{
+ double jmin,jmax; e->scaledJacRange(jmin,jmax);
+ d[e->getNum()].push_back(std::min(jmin,1./jmax));
+ }
}
}
@@ -137,7 +139,7 @@ statisticsWindow::statisticsWindow(int deltaFontSize)
value[num] = new Fl_Output(2 * WB, 2 * WB + 15 * BH, IW, BH, "Rho");
value[num]->tooltip("~ min_edge_length / max_edge_length"); num++;
value[num] = new Fl_Output(2 * WB, 2 * WB + 16 * BH, IW, BH, "Disto");
- value[num]->tooltip("~ min (J0/J, J/J0)"); num++;
+ value[num]->tooltip("~ min (J_min/J_0, J_0/J_max)"); num++;
for(int i = 0; i < 4; i++){
int ww = 3 * FL_NORMAL_SIZE;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 65c48e5a36fa01fb9fa5718b841994193bc2a557..5d048996537be150eb00753162d7334b6eda4f79 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -181,6 +181,7 @@ class MElement
virtual double etaShapeMeasure(){ return 0.; }
virtual double distoShapeMeasure(){ return 1.0; }
virtual double angleShapeMeasure() { return 1.0; }
+ virtual void scaledJacRange(double &jmin, double &jmax);
// get the radius of the inscribed circle/sphere if it exists,
// otherwise get the minimum radius of all the circles/spheres
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index f89e9a4622cc1204938d9d619c6a0159d046d0fd..9906814b6de41b495cc15181364f97530805404e 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -171,6 +171,7 @@ static void GetQualityMeasure(std::vector<T*> &ele,
double &eta, double &etaMin, double &etaMax,
double &rho, double &rhoMin, double &rhoMax,
double &disto, double &distoMin, double &distoMax,
+ double &scaledJacMin, double &scaledJacMax,
double quality[4][100])
{
for(unsigned int i = 0; i < ele.size(); i++){
@@ -186,10 +187,14 @@ static void GetQualityMeasure(std::vector<T*> &ele,
rho += r;
rhoMin = std::min(rhoMin, r);
rhoMax = std::max(rhoMax, r);
- double d = ele[i]->distoShapeMeasure();
+ double jmin,jmax; ele[i]->scaledJacRange(jmin,jmax);
+ double d = std::min(jmin,1./jmax);
disto += d;
distoMin = std::min(distoMin, d);
distoMax = std::max(distoMax, d);
+ scaledJacMin = std::min(scaledJacMin, jmin);
+ scaledJacMax = std::max(scaledJacMax, jmax);
+
for(int j = 0; j < 100; j++){
if(g > j / 100. && g <= (j + 1) / 100.) quality[0][j]++;
if(e > j / 100. && e <= (j + 1) / 100.) quality[1][j]++;
@@ -246,15 +251,16 @@ void GetStatistics(double stat[50], double quality[4][100])
double eta=0., etaMin=1., etaMax=0.;
double rho=0., rhoMin=1., rhoMax=0.;
double disto=0., distoMin=1., distoMax=0.;
+ double jmin = 1.e22, jmax = -1.e22;
double N = 0.0;
if (m->firstRegion() == m->lastRegion()){
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
GetQualityMeasure((*it)->quadrangles, gamma, gammaMin, gammaMax,
eta, etaMin, etaMax, rho, rhoMin, rhoMax,
- disto, distoMin, distoMax, quality);
+ disto, distoMin, distoMax, jmin,jmax, quality);
GetQualityMeasure((*it)->triangles, gamma, gammaMin, gammaMax,
eta, etaMin, etaMax, rho, rhoMin, rhoMax,
- disto, distoMin, distoMax, quality);
+ disto, distoMin, distoMax, jmin,jmax, quality);
}
N = stat[7] + stat[8];
}
@@ -262,16 +268,16 @@ void GetStatistics(double stat[50], double quality[4][100])
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
GetQualityMeasure((*it)->tetrahedra, gamma, gammaMin, gammaMax,
eta, etaMin, etaMax, rho, rhoMin, rhoMax,
- disto, distoMin, distoMax, quality);
+ disto, distoMin, distoMax, jmin,jmax, quality);
GetQualityMeasure((*it)->hexahedra, gamma, gammaMin, gammaMax,
eta, etaMin, etaMax, rho, rhoMin, rhoMax,
- disto, distoMin, distoMax, quality);
+ disto, distoMin, distoMax, jmin,jmax, quality);
GetQualityMeasure((*it)->prisms, gamma, gammaMin, gammaMax,
eta, etaMin, etaMax, rho, rhoMin, rhoMax,
- disto, distoMin, distoMax,quality);
+ disto, distoMin, distoMax,jmin,jmax, quality);
GetQualityMeasure((*it)->pyramids, gamma, gammaMin, gammaMax,
eta, etaMin, etaMax, rho, rhoMin, rhoMax,
- disto, distoMin, distoMax, quality);
+ disto, distoMin, distoMax, jmin,jmax, quality);
}
N = stat[9] + stat[10] + stat[11] + stat[12];
}
@@ -284,6 +290,9 @@ void GetStatistics(double stat[50], double quality[4][100])
stat[23] = rho / N ;
stat[24] = rhoMin;
stat[25] = rhoMax;
+
+ stat[45] = jmin;
+ stat[46] = jmax;
stat[46] = disto / N ;
stat[47] = distoMin;
stat[48] = distoMax;
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 43b1f16f2dd34140a3f718cecd1752e6de866ef2..cb2d12993e12a1a261762f9026913e5637857a41 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1234,8 +1234,40 @@ void printJacobians(GModel *m, const char *nm)
fclose(f);
}
+void getMeshInfoForHighOrder (GModel *gm,
+ int &meshOrder,
+ bool &complete,
+ bool &CAD){
+ meshOrder = -1;
+ CAD = true;
+ for (GModel::riter itr = gm->firstRegion(); itr != gm->lastRegion(); ++itr) {
+ if ((*itr)->getNumMeshElements()){
+ meshOrder = (*itr)->getMeshElement(0)->getPolynomialOrder();
+ complete = (meshOrder <= 2) ? 1 : (*itr)->getMeshElement(0)->getNumVolumeVertices();
+ break;
+ }
+ }
+ for (GModel::fiter itf = gm->firstFace(); itf != gm->lastFace(); ++itf) {
+ if ((*itf)->getNumMeshElements()){
+ if (meshOrder == -1) {
+ meshOrder = (*itf)->getMeshElement(0)->getPolynomialOrder();
+ complete = (meshOrder <= 2) ? 1 : (*itf)->getMeshElement(0)->getNumFaceVertices();
+ if ((*itf)->geomType() == GEntity::DiscreteSurface)CAD = false;
+ break;
+ }
+ }
+ }
+}
+
+
+
+
void ElasticAnalogy ( GModel *m, double threshold, bool onlyVisible) {
+ bool CAD, complete;
+ int meshOrder;
+
+ getMeshInfoForHighOrder (m,meshOrder, complete, CAD);
highOrderTools hot(m);
// now we smooth mesh the internal vertices of the faces
// we do that model face by model face
@@ -1248,7 +1280,8 @@ void ElasticAnalogy ( GModel *m, double threshold, bool onlyVisible) {
std::vector<MElement*> v;
v.insert(v.begin(), (*it)->triangles.begin(), (*it)->triangles.end());
v.insert(v.end(), (*it)->quadrangles.begin(), (*it)->quadrangles.end());
- hot.applySmoothingTo(v, (*it));
+ if (CAD)hot.applySmoothingTo(v, (*it));
+ else hot.applySmoothingTo(v, 1.e32, false);
}
}
// hot.ensureMinimumDistorsion(0.1);
@@ -1432,13 +1465,13 @@ void SetHighOrderInComplete (GModel *m, bool onlyVisible){
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
if (onlyVisible && !(*it)->getVisibility())continue;
std::vector<MTriangle*> newT;
- std::vector<MQuadrangle*> newQ;
+
for (int i=0;i<(*it)->triangles.size();i++){
MTriangle *t = (*it)->triangles[i];
std::vector<MVertex*> vt;
int order = t->getPolynomialOrder();
for (int j=3;j<t->getNumVertices()-t->getNumFaceVertices();j++)vt.push_back(t->getVertex(j));
- for (int j=t->getNumVertices()-t->getNumFaceVertices();j<t->getNumVertices();j++)toDelete.insert(t->getVertex(j));
+ for (int j=t->getNumVertices()-t->getNumFaceVertices();j < t->getNumVertices();j++)toDelete.insert(t->getVertex(j));
newT.push_back(new MTriangleN(t->getVertex(0), t->getVertex(1), t->getVertex(2),
vt, order));
@@ -1446,6 +1479,7 @@ void SetHighOrderInComplete (GModel *m, bool onlyVisible){
}
(*it)->triangles = newT;
+ std::vector<MQuadrangle*> newQ;
for (int i=0;i<(*it)->quadrangles.size();i++){
MQuadrangle *t = (*it)->quadrangles[i];
std::vector<MVertex*> vt;
@@ -1456,13 +1490,18 @@ void SetHighOrderInComplete (GModel *m, bool onlyVisible){
delete t;
}
(*it)->quadrangles = newQ;
+
std::vector<MVertex*> newV;
+ int numd = 0;
for (int i=0;i<(*it)->mesh_vertices.size();++i){
if (toDelete.find((*it)->mesh_vertices[i]) == toDelete.end())
newV.push_back((*it)->mesh_vertices[i]);
- else
+ else{
delete (*it)->mesh_vertices[i];
+ numd++;
+ }
}
+ printf("%d vertices deleted\n");
(*it)->mesh_vertices = newV;
}
diff --git a/Mesh/HighOrder.h b/Mesh/HighOrder.h
index 49be09e410236f4d40eab5a08c280b340f05d402..cd19ea3e2d131cac7c05e6caf2be5f52deb2d7a9 100644
--- a/Mesh/HighOrder.h
+++ b/Mesh/HighOrder.h
@@ -36,5 +36,5 @@ struct distanceFromMeshToGeometry_t {
};
void computeDistanceFromMeshToGeometry (GModel *m, distanceFromMeshToGeometry_t &dist);
-
+void getMeshInfoForHighOrder (GModel *gm, int &meshOrder, bool &complete, bool &CAD);
#endif
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 31ee8bfcd4e1a6fc8985cce2b9e25274dc711211..d55777041c54bc925383b3a922dae157c3310da0 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -307,15 +307,11 @@ double mesh_functional_distorsion_p2_exact(MTriangle *t)
double mesh_functional_distorsion_pN(MElement *t)
{
- // printf("element disto -->\n");
const bezierBasis *jac = t->getJacobianFuncSpace()->bezier;
fullVector<double>Ji(jac->points.size1());
- // printf("%d points for bez \n",jac->points.size1());
for (int i=0;i<jac->points.size1();i++){
double u = jac->points(i,0);
double v = jac->points(i,1);
-
- // JF : bezier points are defined in the [0,1] x [0,1] quad
if (t->getType() == TYPE_QUA){
u = -1 + 2*u;
v = -1 + 2*v;
@@ -326,23 +322,31 @@ double mesh_functional_distorsion_pN(MElement *t)
fullVector<double> Bi( jac->matrixLag2Bez.size1() );
jac->matrixLag2Bez.mult(Ji,Bi);
+ return *std::min_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
+}
- // for (int i=0;i<jac->points.size1();i++){
- // printf("J(%d) = %12.5E B(%d) = %12.5E\n",i,Ji(i),i,Bi(i));
- // }
- /*
- if (t->getType() == TYPE_QUA){
- jac->matrixLag2Bez.print("lag2bez");
+void MElement::scaledJacRange(double &jmin, double &jmax)
+{
+ jmin = jmax = 1.0;
+ if (getPolynomialOrder() == 1) return;
+ const bezierBasis *jac = getJacobianFuncSpace()->bezier;
+ fullVector<double>Ji(jac->points.size1());
+ for (int i=0;i<jac->points.size1();i++){
+ double u = jac->points(i,0);
+ double v = jac->points(i,1);
+ if (getType() == TYPE_QUA){
+ u = -1 + 2*u;
+ v = -1 + 2*v;
+ }
- jac->points.print("lagrangianNodesBezierQ");
- t->getFunctionSpace(t->getPolynomialOrder())->points.print("lagrangianNodesQ");
- t->getFunctionSpace(t->getPolynomialOrder())->monomials.print("MonomialsQ");
- t->getFunctionSpace(t->getPolynomialOrder())->coefficients.print("shapeFunctionCoeffQ");
+ Ji(i) = mesh_functional_distorsion(this,u,v);
}
- */
-
- return *std::min_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
+
+ fullVector<double> Bi( jac->matrixLag2Bez.size1() );
+ jac->matrixLag2Bez.mult(Ji,Bi);
+ jmin = *std::min_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
+ jmax = *std::max_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.cpp b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
index 7c677ea5d1c38a2e409850e1b43888e5eb7ac5bf..1ea6bb6851288475d13d638d2ce02afdbd814630 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
@@ -2,6 +2,7 @@
#include <algorithm>
#include "OptHomMesh.h"
#include "OptHOM.h"
+#include "OptHOMRun.h"
#include "GmshMessage.h"
#include "GmshConfig.h"
@@ -14,7 +15,6 @@
#include "optimization.h"
-
// Constructor
OptHOM::OptHOM(GEntity *ge, std::set<MVertex*> & toFix, int method) :
mesh(ge, toFix, method)
@@ -25,9 +25,6 @@ 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);
@@ -38,8 +35,6 @@ 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]);
}
}
@@ -53,22 +48,13 @@ 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);
- const double dSq = mesh.distSq(iFV), dist = sqrt(dSq);
- Obj += Factor * dSq;
+ Obj += Factor * mesh.distSq(iFV);
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;
@@ -87,42 +73,49 @@ void OptHOM::evalObjGrad(const alglib::real_1d_array &x, double &Obj, alglib::re
addJacObjGrad(Obj, gradObj);
addDistObjGrad(lambda, lambda2, Obj, gradObj);
- if ((minJac > barrier_min) && (maxJac < barrier_max)) {
- printf("INFO: reached Jacobian requirements, setting null gradient\n");
- Obj = 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");
+ // Obj = 0.;
+ // for (int i = 0; i < gradObj.length(); i++) gradObj[i] = 0.;
+ // }
}
void evalObjGradFunc(const alglib::real_1d_array &x, double &Obj, alglib::real_1d_array &gradObj, void *HOInst)
{
- (static_cast<OptHOM*>(HOInst))->evalObjGrad(x, Obj, gradObj);
+ 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;
+ }
+ }
}
-void OptHOM::recalcJacDist()
- {
+void OptHOM::getDistances(double &distMaxBND, double &distAvgBND, double &minJac, double &maxJac)
+{
- maxDist = 0;
- avgDist = 0;
+ distMaxBND = 0;
+ distAvgBND = 0;
int nbBnd = 0;
for (int iFV = 0; iFV < mesh.nFV(); iFV++) {
if (mesh.forced(iFV)) {
double dSq = mesh.distSq(iFV);
- maxDist = std::max(maxDist, sqrt(dSq));
- avgDist += sqrt(dSq);
+ distMaxBND = std::max(distMaxBND, sqrt(dSq));
+ distAvgBND += sqrt(dSq);
nbBnd++;
}
}
- if (nbBnd != 0) avgDist /= nbBnd;
+ if (nbBnd != 0) distAvgBND /= nbBnd;
- minJac = 1.e300;
- maxJac = -1.e300;
+ minJac = 1000;
+ maxJac = -1000;
for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
std::vector<double> sJ(mesh.nBezEl(iEl)); // Scaled Jacobians
mesh.scaledJac(iEl,sJ);
@@ -142,8 +135,11 @@ void OptHOM::printProgress(const alglib::real_1d_array &x, double Obj)
iter++;
if (iter % progressInterv == 0) {
- 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);
+ 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);
}
}
@@ -157,83 +153,53 @@ 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 = 0.;
+ static const double EPSF = 1.e-8;
+// static const double EPSF = 0.;
static const double EPSX = 0.;
- static int OPTMETHOD = 1;
+// 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);
Msg::Debug("--- Optimization pass with jacBar = %12.5E",jacBar);
- iter = 0;
+// alglib::minlbfgsstate state;
+// alglib::minlbfgsreport rep;
+ alglib::mincgstate state;
+ alglib::mincgreport rep;
- 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;
- }
+// 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);
- 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;
+ 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;
}
@@ -241,6 +207,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)
{
+ double minJ,maxJ;
barrier_min = b_min;
barrier_max = b_max;
progressInterv = pInt;
@@ -254,6 +221,8 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double
mesh.distSqToStraight(dSq);
const double maxDSq = *max_element(dSq.begin(),dSq.end());
invLengthScaleSq = 1./maxDSq; // Length scale for non-dimensional distance
+ // printf("DISTANCE TO STRAIGHT = %12.5E\n",maxDSq);
+
// Set initial guess
alglib::real_1d_array x;
@@ -261,32 +230,34 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double
mesh.getUvw(x.getcontent());
// Calculate initial performance
- recalcJacDist();
- initMaxDist = maxDist;
- initAvgDist = avgDist;
+ // double minJ, maxJ;
+ double initMaxD, initAvgD;
+ getDistances(initMaxD, initAvgD, minJ, maxJ);
- const double jacBarStart = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
+ const double jacBarStart = (minJ > 0.) ? 0.9*minJ : 1.1*minJ;
jacBar = jacBarStart;
setBarrierTerm(jacBarStart);
+ // std::cout << "DEBUG: jacBarStart = " << jacBarStart << std::endl;
// Calculate initial objective function value and gradient
initObj = 0.;
alglib::real_1d_array gradObj;
gradObj.setlength(mesh.nPC());
+ // printf("mesh.npc = %d\n",mesh.nPC());
for (int i = 0; i < mesh.nPC(); i++) gradObj[i] = 0.;
evalObjGrad(x, initObj, gradObj);
- std::cout << "Start optimizing " << mesh.nEl() << " elements (" << mesh.nVert() << " vertices, "
- << mesh.nFV() << " free vertices, " << mesh.nPC() << " variables) with min. barrier = " << barrier_min
- << " and max. barrier = " << barrier_max << std::endl;
+ // std::cout << "Start optimizing with barrier = " << barrier_min << std::endl;
int ITER = 0;
- while (minJac < barrier_min) {
+ while (minJ < barrier_min || maxJ > barrier_max ) {
OptimPass(x, gradObj, itMax);
- jacBar = (minJac > 0.) ? 0.9*minJac : 1.1*minJac;
+ double dumMaxD, dumAvgD;
+ getDistances(dumMaxD, dumAvgD, minJ, maxJ);
+ jacBar = (minJ > 0.) ? 0.9*minJ : 1.1*minJ;
setBarrierTerm(jacBar);
- if (ITER ++ > 10) break;
+ if (ITER ++ > 15) break;
}
// for (int i = 0; i<3; i++) {
@@ -294,12 +265,14 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double
// OptimPass(x, gradObj, itMax);
// }
- double finalObj = 0.;
+ double finalObj = 0., finalMaxD, finalAvgD;
evalObjGrad(x, finalObj, gradObj);
- Msg::Info("Optimization finished : Avg distance to bnd = %12.5E MinJac %12.5E MaxJac %12.5E",avgDist,minJac,maxJac);
-
- return (minJac > barrier_min && maxJac < barrier_max );
+ getDistances(finalMaxD, finalAvgD, minJ, maxJ);
+ OptHomMessage("Optimization done Range (%g,%g)",minJ,maxJ);
+ if (minJ < barrier_min && maxJ > barrier_max) return 0;
+ if (minJ > 0.0) return 1;
+ return -1;
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.h b/contrib/HighOrderMeshOptimizer/OptHOM.h
index f2b18530e2ff576fb92cc7e2ac9055b41790e98a..093411aabb1624fe9c6e63454cc56b3cda36d10a 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.h
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.h
@@ -19,12 +19,17 @@ class OptHOM
public:
+
Mesh mesh;
OptHOM(GEntity *gf, std::set<MVertex*> & toFix, int method);
+ // 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
- void recalcJacDist();
- inline void getJacDist(double &minJ, double &maxJ, double &maxD, double &avgD);
+ void getDistances(double &distMaxBND, double &distAvgBND, double &minJac, double &maxJac);
+ // 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);
@@ -81,10 +86,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/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 46fffaaebe1fad33207c5afec52edadbfbe5d1a2..abb960a711888c76be0f4d16d9b24755124ec245 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -16,6 +16,33 @@
#include "OptHOM.h"
#include <stack>
+#ifdef HAVE_FLTK
+#include "highOrderToolsWindow.h"
+#include "FLGui.h"
+#endif
+
+void OptHomMessage (const char *s, ...) {
+ char str[1024];
+ va_list args;
+ va_start(args, s);
+ vsnprintf(str, sizeof(str), s, args);
+ va_end(args);
+#ifdef HAVE_FLTK
+ if(FlGui::available()){
+ FlGui::instance()->check();
+ FlGui::instance()->highordertools->messages->add(str, 0);
+ if(FlGui::instance()->highordertools->win->shown() &&
+ FlGui::instance()->highordertools->messages->h() >= 10){
+ FlGui::instance()->highordertools->messages->bottomline(FlGui::instance()->highordertools->messages->size());
+ FlGui::instance()->highordertools->messages->show();
+ }
+ }
+#else
+ fprintf(stdout,"%s\n",str);
+#endif
+}
+
+
// get all elements that are neighbors
double distMaxStraight (MElement *el){
@@ -170,21 +197,39 @@ static MElement * compare_worst (MElement *a, MElement *b){
return b;
}
-template <class T>
-MElement* getTheWorstElement (const std::vector<T*> &a) {
- T *worst = 0;
- double q = 1.e22;
- for (int i=0;i<a.size();i++){
- T *t = a[i];
- double Q = t->distoShapeMeasure();
- if (Q < q) {
- worst = t;q = Q;
+template <class ITERATOR>
+MElement* getTheWorstElementDown (const ITERATOR &beg, const ITERATOR &end, double &q) {
+ MElement *worst = 0;
+ q = 1.e22;
+ ITERATOR it = beg;
+ for (;it != end;++it){
+ MElement *t = *it;
+ double jmin,jmax;
+ t->scaledJacRange(jmin,jmax);
+ if (jmin < q) {
+ worst = t;q = jmin;
}
}
- // printf("worst = %12.5E\n",q);
return worst;
}
+template <class ITERATOR>
+MElement* getTheWorstElementUp (const ITERATOR &beg, const ITERATOR &end, double &q) {
+ MElement *worst = 0;
+ q = -1.e22;
+ ITERATOR it = beg;
+ for (;it != end;++it){
+ MElement *t = *it;
+ double jmin,jmax;
+ t->scaledJacRange(jmin,jmax);
+ if (jmax > q) {
+ worst = t;q = jmax;
+ }
+ }
+ return worst;
+}
+
+
std::set<MVertex*> filter3D(GRegion *gr, int nbLayers, double _qual)
{
std::set<MVertex*> touched;
@@ -196,7 +241,7 @@ std::set<MVertex*> filter3D(GRegion *gr, int nbLayers, double _qual)
ts.insert(t);
for(int j=0;j<t->getNumVertices();j++)touched.insert(t->getVertex(j));
}
- if (ts.size() == 1)break;
+ // if (ts.size() == 1)break;
}
printf("FILTER3D : %d bad tets found among %6d\n",ts.size(),gr->tetrahedra.size());
@@ -242,9 +287,18 @@ std::set<MVertex*> filter2D_boundaryLayer(GFace *gf,
int nbLayers,
double _qual_min,
double _qual_max,
- double F ) {
- MElement *worst = compare_worst (getTheWorstElement(gf->triangles),
- getTheWorstElement(gf->quadrangles));
+ double F ,
+ std::set<MElement*> & badasses) {
+
+ double jmin, jmax;
+ MElement *worstDown = getTheWorstElementDown(badasses.begin(),badasses.end(),jmin);
+ MElement *worstUp = getTheWorstElementUp(badasses.begin(),badasses.end(),jmax);
+ MElement *worst;
+ if (jmin < _qual_min)worst = worstDown;
+ else worst = worstUp;
+
+ // MElement *worst = compare_worst (getTheWorstElement(gf->triangles),
+ // getTheWorstElement(gf->quadrangles));
std::vector<MElement*> vworst; vworst.push_back(worst);
std::vector<MElement*> all;
@@ -284,9 +338,10 @@ std::set<MVertex*> filter3D_boundaryLayer(GRegion *gr,
double _qual_min,
double _qual_max,
double F ) {
- MElement *worst = compare_worst (getTheWorstElement(gr->tetrahedra),
- getTheWorstElement(gr->prisms));
- worst = compare_worst (worst,getTheWorstElement(gr->hexahedra));
+ double jmin,jmax;
+ MElement *worst = compare_worst (getTheWorstElementDown(gr->tetrahedra.begin(),gr->tetrahedra.end(),jmin),
+ getTheWorstElementDown(gr->prisms.begin(),gr->prisms.end(),jmin));
+ worst = compare_worst (worst,getTheWorstElementDown(gr->hexahedra.begin(),gr->hexahedra.end(),jmin));
std::vector<MElement*> vworst; vworst.push_back(worst);
std::vector<MElement*> all;
@@ -327,27 +382,32 @@ std::set<MVertex*> filter3D_boundaryLayer(GRegion *gr,
std::set<MVertex*> filter2D(GFace *gf,
int nbLayers,
double _qual_min,
- double _qual_max,
- double F = 1.e22)
+ double _qual_max)
{
std::set<MVertex*> touched;
std::set<MVertex*> not_touched;
std::set<MElement*> elements;
+ std::set<MTriangle*> tt;
+ std::set<MQuadrangle*> qq;
for (int i = 0; i < gf->triangles.size(); i++) {
MTriangle *t = gf->triangles[i];
- double Q = t->distoShapeMeasure();
- if (Q < _qual_min || Q > _qual_max) {
+ double jmin,jmax;
+ t->scaledJacRange(jmin,jmax);
+ if (jmin < _qual_min || jmax > _qual_max) {
elements.insert(t);
+ tt.insert(t);
for (int j = 0; j < t->getNumVertices(); j++)
touched.insert(t->getVertex(j));
}
}
for (int i = 0; i < gf->quadrangles.size(); i++) {
MQuadrangle *q = gf->quadrangles[i];
- double Q = q->distoShapeMeasure();
- if (Q < _qual_min || Q > _qual_max) {
+ double jmin,jmax;
+ q->scaledJacRange(jmin,jmax);
+ if (jmin < _qual_min || jmax > _qual_max) {
elements.insert(q);
+ qq.insert(q);
for (int j = 0; j < q->getNumVertices(); j++)
touched.insert(q->getVertex(j));
}
@@ -360,14 +420,16 @@ std::set<MVertex*> filter2D(GFace *gf,
for (int j = 0; j < el->getNumVertices(); ++j) {
if (touched.find(el->getVertex(j)) != touched.end()) {
elements.insert(el);
+ if (el->getType() == TYPE_TRI)tt.insert((MTriangle*)el);
+ if (el->getType() == TYPE_QUA)qq.insert((MQuadrangle*)el);
break;
}
}
}
for (std::set<MElement *>::iterator it = elements.begin(); it != elements.end(); ++it) {
- MElement &el = **it;
- for (int j = 0; j < el.getNumVertices(); ++j) {
- touched.insert(el.getVertex(j));
+ MElement *el = *it;
+ for (int j = 0; j < el->getNumVertices(); ++j) {
+ touched.insert(el->getVertex(j));
}
}
}
@@ -379,6 +441,10 @@ std::set<MVertex*> filter2D(GFace *gf,
}
}
}
+ gf->triangles.clear();
+ gf->triangles.insert(gf->triangles.begin(),tt.begin(),tt.end());
+ gf->quadrangles.clear();
+ gf->quadrangles.insert(gf->quadrangles.begin(),qq.begin(),qq.end());
return not_touched;
}
@@ -427,7 +493,7 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
clock_t t1 = clock();
- int samples = 20;
+ int samples = 100;
// int method = Mesh::METHOD_RELAXBND | Mesh::METHOD_PHYSCOORD | Mesh::METHOD_PROJJAC;
// int method = Mesh::METHOD_FIXBND | Mesh::METHOD_PHYSCOORD | Mesh::METHOD_PROJJAC;
@@ -449,47 +515,68 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
SVector3 BND(gm->bounds().max(), gm->bounds().min());
double SIZE = BND.norm();
- Msg::Info("High order mesh optimizer starts");
-
+ OptHomMessage("High order mesh optimizer starts");
+ double distMaxBND, distAvgBND, minJac, maxJac;
if (p.dim == 2) {
+ clock_t tf1 = clock();;
for (GModel::fiter itf = gm->firstFace(); itf != gm->lastFace(); ++itf) {
if (p.onlyVisible && !(*itf)->getVisibility())continue;
int ITER = 0;
- Msg::Info("Model face %4d is considered",(*itf)->tag());
- p.SUCCESS = true;
+ OptHomMessage("Optimizing Model face %4d...",(*itf)->tag());
+ p.SUCCESS = 1;
+ std::set<MElement*> badasses;
+ if (p.filter == 1){
+ for (int i=0;i<(*itf)->getNumMeshElements();i++){
+ double jmin,jmax;
+ (*itf)->getMeshElement(i)->scaledJacRange(jmin,jmax);
+ if (jmin < p.BARRIER_MIN || jmax > p.BARRIER_MAX)badasses.insert((*itf)->getMeshElement(i));
+ }
+ // printf("START WITH %d bad asses\n",badasses.size());
+ if (badasses.size() == 0)continue;
+ }
while (1){
std::vector<MTriangle*> tris = (*itf)->triangles;
std::vector<MQuadrangle*> quas = (*itf)->quadrangles;
std::set<MVertex*> toFix;
- if (p.filter == 1) toFix = filter2D_boundaryLayer(*itf, p.nbLayers, p.BARRIER_MIN, p.BARRIER_MAX, p.DistanceFactor);
- else toFix = filter2D(*itf, p.nbLayers, p.BARRIER_MIN, p.BARRIER_MAX, p.DistanceFactor);
-
+ if (p.filter == 1) toFix = filter2D_boundaryLayer(*itf, p.nbLayers, p.BARRIER_MIN, p.BARRIER_MAX, p.DistanceFactor, badasses);
+ else toFix = filter2D(*itf, p.nbLayers, p.BARRIER_MIN, p.BARRIER_MAX);
OptHOM *temp = new OptHOM(*itf, toFix, method);
- double 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;
+ // temp->recalcJacDist();
+
+ temp->getDistances(distMaxBND, distAvgBND,minJac, maxJac);
+ // if (minJac < 1.e2)OptHomMessage("Optimizing a blob of %4d elements minJ %12.5E -- maxJ %12.5E",(*itf)->getNumMeshElements(), minJac, maxJac);
std::ostringstream ossI;
ossI << "initial_" << (*itf)->tag() << "ITER_" << ITER << ".msh";
temp->mesh.writeMSH(ossI.str().c_str());
+ (*itf)->triangles = tris;
+ (*itf)->quadrangles = quas;
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);
- temp->getJacDist(p.minJac, p.maxJac, distMaxBND, distAvgBND);
+ p.SUCCESS = std::min(p.SUCCESS,temp->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, samples, p.itMax));
+
+ temp->getDistances(distMaxBND, distAvgBND,minJac, maxJac);
+ // temp->getJacDist(p.minJac, p.maxJac, distMaxBND, distAvgBND);
temp->mesh.updateGEntityPositions();
- if (!p.SUCCESS) break;
+ if (p.SUCCESS == -1) break;
ITER ++;
+ if (p.filter == 1 && ITER > badasses.size() * 2)break;
// std::ostringstream ossF;
// ossF << "final_" << (*itf)->tag() << ".msh";
// temp->mesh.writeMSH(ossF.str().c_str());
}
- Msg::Info("Model face %4d is done (%d subdomains were computed) SUCCESS %1d",(*itf)->tag(),ITER,p.SUCCESS);
- Msg::Info("----------------------------------------------------------------");
+ double DTF = (double)(clock()-tf1) / CLOCKS_PER_SEC;
+ if (p.SUCCESS == 1){
+ OptHomMessage("Optimization succeeded (CPU %g sec)",DTF);
+ }
+ else if (p.SUCCESS == 0)
+ OptHomMessage("Warning : Model face %4d has all jacobians positive but not all in the range CPU %g",(*itf)->tag(),DTF);
+ else if (p.SUCCESS == -1)
+ OptHomMessage("Error : Model face %4d has still negative jacobians",(*itf)->tag());
+
+ Msg::Info("---------------------------------------------------------------------------------------------------------------");
}
exportMeshToDassault (gm,gm->getName() + "_opti.das", 2);
}
@@ -509,8 +596,9 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
// if ((*itr)->tetrahedra.size() > 0) {
OptHOM *temp = new OptHOM(*itr, toFix, method);
double distMaxBND, distAvgBND, minJac, maxJac;
- temp->recalcJacDist();
- temp->getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
+ 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",(*itr)->getNumMeshElements(), minJac, maxJac);
(*itr)->tetrahedra = tets;
std::ostringstream ossI;
@@ -518,7 +606,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);
- temp->getJacDist(p.minJac, p.maxJac, distMaxBND, distAvgBND);
+ temp->getDistances(distMaxBND, distAvgBND,minJac, maxJac);
+ // temp->getJacDist(p.minJac, p.maxJac, distMaxBND, distAvgBND);
temp->mesh.updateGEntityPositions();
if (!p.SUCCESS) break;
// }
@@ -530,5 +619,5 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
}
}
clock_t t2 = clock();
- p.CPU = (double)(t2-t2)/CLOCKS_PER_SEC;
+ p.CPU = (double)(t2-t1)/CLOCKS_PER_SEC;
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.h b/contrib/HighOrderMeshOptimizer/OptHomRun.h
index feac05b006f3b45f0ecac6ad5d0595748b926c85..5b6cebb9460ca1438c8d79559503dae4e7b14426 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.h
@@ -22,6 +22,9 @@ struct OptHomParameters {
int SUCCESS ; // 0 --> success , 1 --> Not converged
double minJac, maxJac; // after optimization, range of jacobians
double CPU; // Time for optimization
+
+ // is called by the optimizer at every stage of the
+ // virtual void progressCallback (int percetageDone) const {}
OptHomParameters ()
// default values
@@ -32,6 +35,7 @@ struct OptHomParameters {
};
void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p);
+void OptHomMessage (const char *s, ...);
#endif