diff --git a/Fltk/graphicWindow.cpp b/Fltk/graphicWindow.cpp
index d6ada68fc30df33ef9891b5f92345d577e16f36a..03bbbe73dba2cf2a18f33736a1a253ecdb90778f 100644
--- a/Fltk/graphicWindow.cpp
+++ b/Fltk/graphicWindow.cpp
@@ -3289,7 +3289,7 @@ static menuItem static_modules[] = {
(Fl_Callback *)mesh_degree_cb, (void*)2},
{"0Modules/Mesh/Set order 3",
(Fl_Callback *)mesh_degree_cb, (void*)3},
- {"0Modules/Mesh/Optimize high order",
+ {"0Modules/Mesh/High order tools",
(Fl_Callback *)highordertools_cb},
{"0Modules/Mesh/Inspect",
(Fl_Callback *)mesh_inspect_cb} ,
diff --git a/Fltk/highOrderToolsWindow.cpp b/Fltk/highOrderToolsWindow.cpp
index 566a5b5e3616f69b9e6218aa460e7fbe320bb530..cc556d3d120dfc465907d84f4619b5b994405c8a 100644
--- a/Fltk/highOrderToolsWindow.cpp
+++ b/Fltk/highOrderToolsWindow.cpp
@@ -4,11 +4,6 @@
// bugs and problems to the public mailing list <gmsh@geuz.org>.
#include "GmshConfig.h"
-#if !defined(HAVE_NO_STDINT_H)
-#include <stdint.h>
-#elif defined(HAVE_NO_INTPTR_T)
-typedef unsigned long intptr_t;
-#endif
#include <string>
#include <sstream>
#include <map>
@@ -20,19 +15,13 @@ typedef unsigned long intptr_t;
#include <FL/Fl_Value_Input.H>
#include "GmshConfig.h"
#include "FlGui.h"
+#include "graphicWindow.h"
#include "drawContext.h"
#include "highOrderToolsWindow.h"
#include "paletteWindow.h"
-#include "contextWindow.h"
-#include "graphicWindow.h"
-#include "GmshDefines.h"
#include "GmshMessage.h"
#include "GModel.h"
#include "MElement.h"
-#include "PView.h"
-#include "PViewData.h"
-#include "GeoStringInterface.h"
-#include "Options.h"
#include "Context.h"
#include "HighOrder.h"
@@ -40,10 +29,6 @@ typedef unsigned long intptr_t;
#include "OptHomRun.h"
#endif
-#if defined(HAVE_PARSER)
-#include "Parser.h"
-#endif
-
static void change_completeness_cb(Fl_Widget *w, void *data)
{
highOrderToolsWindow *o = FlGui::instance()->highordertools;
@@ -62,7 +47,6 @@ static void change_completeness_cb(Fl_Widget *w, void *data)
static void highordertools_runp_cb(Fl_Widget *w, void *data)
{
-#if defined(HAVE_OPTHOM)
highOrderToolsWindow *o = FlGui::instance()->highordertools;
int order = (int)o->value[0]->value();
@@ -79,18 +63,16 @@ static void highordertools_runp_cb(Fl_Widget *w, void *data)
computeDistanceFromMeshToGeometry (GModel::current(), dist);
for (std::map<GEntity*, double> ::iterator it = dist.d2.begin();
it !=dist.d2.end();++it){
- printf ("GEntity %d of dim %d : dist %12.5E\n",
- it->first->tag(), it->first->dim(), it->second);
+ printf("GEntity %d of dim %d : dist %12.5E\n",
+ it->first->tag(), it->first->dim(), it->second);
}
CTX::instance()->mesh.changed |= (ENT_LINE | ENT_SURFACE | ENT_VOLUME);
drawContext::global()->draw();
-#endif
}
static void chooseopti_cb(Fl_Widget *w, void *data)
{
-#if defined(HAVE_OPTHOM)
highOrderToolsWindow *o = FlGui::instance()->highordertools;
int elastic = o->choice[2]->value();
@@ -111,37 +93,39 @@ static void chooseopti_cb(Fl_Widget *w, void *data)
for (int i=9;i<=11;i++) o->value[i]->activate();
// o->push[1]->activate();
}
-#endif
}
static void chooseopti_strategy(Fl_Widget *w, void *data)
{
-#if defined(HAVE_OPTHOM)
highOrderToolsWindow *o = FlGui::instance()->highordertools;
if (o->choice[3]->value() == 0) for (int i=9;i<=11;i++) o->value[i]->deactivate();
else for (int i=9;i<=11;i++) o->value[i]->activate();
-#endif
}
-static void highordertools_runelas_cb(Fl_Widget *w, void *data)
+static void highordertools_runopti_cb(Fl_Widget *w, void *data)
{
-#if defined(HAVE_OPTHOM)
highOrderToolsWindow *o = FlGui::instance()->highordertools;
+ if(o->butt[3]->value())
+ FlGui::instance()->graph[0]->showMessages();
+
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();
double threshold_max = o->value[8]->value();
- if(elastic)ElasticAnalogy(GModel::current(), threshold_min, onlyVisible);
+ if(elastic){
+ ElasticAnalogy(GModel::current(), threshold_min, onlyVisible);
+ }
else {
+#if defined(HAVE_OPTHOM)
OptHomParameters p;
p.nbLayers = nbLayers;
p.BARRIER_MIN = threshold_min;
p.BARRIER_MAX = threshold_max;
p.onlyVisible = onlyVisible;
- p.dim = GModel::current()->getDim();//o->meshOrder;
+ p.dim = GModel::current()->getDim();
p.itMax = (int) o->value[3]->value();
p.optPassMax = (int) o->value[4]->value();
p.weightFixed = o->value[5]->value();
@@ -153,12 +137,11 @@ static void highordertools_runelas_cb(Fl_Widget *w, void *data)
p.adaptBlobLayerFact = (int) o->value[10]->value();
p.adaptBlobDistFact = o->value[11]->value();
HighOrderMeshOptimizer (GModel::current(),p);
- printf("CPU TIME = %4f seconds\n",p.CPU);
+#endif
}
CTX::instance()->mesh.changed |= (ENT_LINE | ENT_SURFACE | ENT_VOLUME);
drawContext::global()->draw();
-#endif
}
highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
@@ -166,24 +149,31 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
FL_NORMAL_SIZE -= deltaFontSize;
int width = 3 * IW + 4 * WB;
- int height = 28 * BH;
+ int height = 24 * BH;
win = new paletteWindow
- (width, height, CTX::instance()->nonModalWindows ? true : false, "High Order Tools");
+ (width, height, CTX::instance()->nonModalWindows ? true : false, "High order tools");
win->box(GMSH_WINDOW_BOX);
int y = WB;
int x = 2 * WB;
+ box[0] = new Fl_Box(x, y, width - 4 * WB, BH);
+ box[0]->align(FL_ALIGN_LEFT|FL_ALIGN_INSIDE);
+
+ y += BH;
+
butt[1] = new Fl_Check_Button
- (x,y, 1.5*IW-WB, BH, "Visible entities only");
+ (x, y, width - 4 * WB, BH, "Only apply high order tools to visible entities");
butt[1]->type(FL_TOGGLE_BUTTON);
butt[1]->value(1);
- output[0] = new Fl_Output
- (width/2,y, IW, BH, "CAD");
- output[0]->align(FL_ALIGN_RIGHT);
- output[0]->value("Available");
+ y += BH;
+
+ butt[3] = new Fl_Check_Button
+ (x, y, width - 4 * WB, BH, "Show detailed log messages");
+ butt[3]->type(FL_TOGGLE_BUTTON);
+ butt[3]->value(1);
{
y += BH / 2;
@@ -195,14 +185,14 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
{
y += BH;
Fl_Box *b = new Fl_Box
- (x - WB, y, width, BH, "1. Generation of high order nodes");
+ (x - WB, y, width, BH, "1. Generation of high order vertices");
b->align(FL_ALIGN_LEFT | FL_ALIGN_INSIDE);
}
y += BH;
value[0] = new Fl_Value_Input
- (x,y, IW, BH, "Polynomial order");
+ (x, y, IW, BH, "Polynomial order");
value[0]->minimum(1);
value[0]->maximum(10);
value[0]->step(1);
@@ -212,7 +202,7 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
y += BH;
butt[0] = new Fl_Check_Button
- (x,y, 1.5*IW-WB, BH, "Use incomplete elements");
+ (x, y, width - 4 * WB, BH, "Use incomplete elements");
butt[0]->type(FL_TOGGLE_BUTTON);
butt[0]->value(!complete);
butt[0]->callback(change_completeness_cb);
@@ -220,12 +210,14 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
y += BH;
butt[2] = new Fl_Check_Button
- (x,y, 1.5*IW-WB, BH, "Generate curvilinear elements");
+ (x, y, width - 4 * WB, BH, "Generate curvilinear elements");
butt[2]->type(FL_TOGGLE_BUTTON);
butt[2]->value(1);
+ y += BH;
+
push[0] = new Fl_Button
- (width - BB - 2 * WB, y, BB, BH, "Apply");
+ (width - BB - 2 * WB, y, BB, BH, "Generate");
push[0]->callback(highordertools_runp_cb);
{
@@ -238,14 +230,14 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
{
y += BH;
Fl_Box *b = new Fl_Box
- (x - WB, y, width, BH, "2. Optimization");
+ (x - WB, y, width, BH, "2. Optimization of high order elements");
b->align(FL_ALIGN_LEFT | FL_ALIGN_INSIDE);
}
y += BH;
value[1] = new Fl_Value_Input
- (x,y, IW/2.0, BH);
+ (x, y, IW/2.0, BH);
value[1]->minimum(0);
value[1]->maximum(1);
value[1]->step(.01);
@@ -253,7 +245,7 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
value[1]->value(0.1);
value[8] = new Fl_Value_Input
- (x+IW/2.0,y, IW/2.0, BH, "Jacobian range");
+ (x+IW/2.0,y, IW/2.0, BH, "Target Jacobian range");
value[8]->minimum(1);
value[8]->maximum(10);
value[8]->step(.01);
@@ -262,7 +254,7 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
y += BH;
value[2] = new Fl_Value_Input
- (x,y, IW, BH, "Number of layers");
+ (x, y, IW, BH, "Number of layers");
value[2]->minimum(1);
value[2]->maximum(20);
value[2]->step(1);
@@ -271,7 +263,7 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
y += BH;
value[7] = new Fl_Value_Input
- (x,y, IW, BH, "Distance factor");
+ (x, y, IW, BH, "Distance factor");
value[7]->minimum(1);
value[7]->maximum(20000);
value[7]->step(1);
@@ -286,7 +278,7 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
y += BH;
choice[2] = new Fl_Choice
- (x,y, IW, BH, "Algorithm");
+ (x, y, IW, BH, "Algorithm");
choice[2]->align(FL_ALIGN_RIGHT);
choice[2]->menu(menu_method);
choice[2]->callback(chooseopti_cb);
@@ -299,24 +291,23 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
y += BH;
choice[0] = new Fl_Choice
- (x,y, IW, BH, "Boundary nodes");
+ (x, y, IW, BH, "Boundary vertices");
choice[0]->menu(menu_objf);
choice[0]->align(FL_ALIGN_RIGHT);
y += BH;
value[5] = new Fl_Value_Input
- (x,y, IW, BH, "W fixed");
+ (x, y, IW/2, BH);
value[5]->align(FL_ALIGN_RIGHT);
value[5]->value(1.e+5);
-
value[6] = new Fl_Value_Input
- (x+IW*1.5,y, IW, BH, "W free");
+ (x+IW/2,y, IW/2, BH, "W fixed / W free");
value[6]->align(FL_ALIGN_RIGHT);
value[6]->value(1.e+2);
y += BH;
value[3] = new Fl_Value_Input
- (x,y, IW, BH, "Max. number of iterations");
+ (x, y, IW, BH, "Maximum number of iterations");
value[3]->minimum(1);
value[3]->maximum(10000);
value[3]->step(10);
@@ -325,7 +316,7 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
y += BH;
value[4] = new Fl_Value_Input
- (x,y, IW, BH, "Max. number of optimization passes");
+ (x, y, IW, BH, "Max. number of optimization passes");
value[4]->minimum(1);
value[4]->maximum(100);
value[4]->step(1);
@@ -340,14 +331,14 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
y += BH;
choice[3] = new Fl_Choice
- (x,y, IW, BH, "Strategy");
+ (x, y, IW, BH, "Strategy");
choice[3]->menu(menu_strategy);
choice[3]->align(FL_ALIGN_RIGHT);
choice[3]->callback(chooseopti_strategy);
y += BH;
value[9] = new Fl_Value_Input
- (x,y, IW, BH, "Max. number of blob adaptation iter.");
+ (x, y, IW, BH, "Max. number of blob adaptation iter.");
value[9]->minimum(1);
value[9]->maximum(100);
value[9]->step(1);
@@ -357,7 +348,7 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
y += BH;
value[10] = new Fl_Value_Input
- (x,y, IW, BH, "Num. layer adaptation factor");
+ (x, y, IW, BH, "Num. layer adaptation factor");
value[10]->align(FL_ALIGN_RIGHT);
value[10]->minimum(1);
value[10]->maximum(100);
@@ -367,7 +358,7 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
y += BH;
value[11] = new Fl_Value_Input
- (x,y, IW, BH, "Distance adaptation factor");
+ (x, y, IW, BH, "Distance adaptation factor");
value[11]->align(FL_ALIGN_RIGHT);
value[11]->minimum(1.);
value[11]->maximum(100.);
@@ -376,16 +367,8 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
y += 1.5*BH;
push[1] = new Fl_Button
- (x,y, IW, BH, "Apply");
- push[1]->callback(highordertools_runelas_cb);
-
- 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);
+ (width - BB - 2 * WB, y, BB, BH, "Optimize");
+ push[1]->callback(highordertools_runopti_cb);
// win->resizable(o);
win->position(CTX::instance()->hotPosition[0], CTX::instance()->hotPosition[1]);
@@ -403,11 +386,11 @@ void highOrderToolsWindow::show(bool redrawOnly)
value[0]->value(meshOrder);
butt[0]->value(!complete);
if (CAD) {
- output[0]->value("Available");
+ box[0]->label("CAD model is available");
choice[0]->value(1);
}
else {
- output[0]->value("Not Available");
+ box[0]->label("CAD model is not available");
choice[0]->deactivate();
}
win->show();
diff --git a/Fltk/highOrderToolsWindow.h b/Fltk/highOrderToolsWindow.h
index eedf7d24d9653d081b3249f59381432370423ce0..21e4f24594e7f16d3130995e94e4ab8e520478fa 100644
--- a/Fltk/highOrderToolsWindow.h
+++ b/Fltk/highOrderToolsWindow.h
@@ -7,6 +7,7 @@
#define _HIGHORDERTOOLS_WINDOW_H_
#include <FL/Fl_Window.H>
+#include <FL/Fl_Box.H>
#include <FL/Fl_Choice.H>
#include <FL/Fl_Browser.H>
#include <FL/Fl_Multi_Browser.H>
@@ -20,12 +21,11 @@ class highOrderToolsWindow{
bool CAD, complete;
int meshOrder;
Fl_Window *win;
+ Fl_Box *box[20];
Fl_Check_Button *butt[20];
Fl_Value_Input *value[20];
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/contrib/HighOrderMeshOptimizer/OptHOM.cpp b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
index e722f24a545dc4a9073f5c68bf826e5bc84417ef..25b72d1e25d4d81f322ae553114545bb726965c6 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
@@ -6,7 +6,7 @@
#include "GmshMessage.h"
#include "GmshConfig.h"
-#ifdef HAVE_BFGS
+#if defined(HAVE_BFGS)
#include "ap.h"
#include "alglibinternal.h"
@@ -22,9 +22,9 @@ static inline double compute_f(double v, double barrier)
return l * l + m * m;
}
else return 1.e300;
-// if (v < 1.) return pow(1.-v,powM);
-// if (v < 1.) return exp((long double)pow(1.-v,3));
-// else return pow(v-1.,powP);
+ // if (v < 1.) return pow(1.-v,powM);
+ // if (v < 1.) return exp((long double)pow(1.-v,3));
+ // else return pow(v-1.,powP);
}
static inline double compute_f1(double v, double barrier)
@@ -33,35 +33,32 @@ static inline double compute_f1(double v, double barrier)
return 2 * (v - 1) + 2 * log((v - barrier) / (1 - barrier)) / (v - barrier);
}
else return barrier < 1 ? -1.e300 : 1e300;
-// if (v < 1.) return -powM*pow(1.-v,powM-1.);
-// if (v < 1.) return -3.*pow(1.-v,2)*exp((long double)pow(1.-v,3));
-// else return powP*pow(v-1.,powP-1.);
+ // if (v < 1.) return -powM*pow(1.-v,powM-1.);
+ // if (v < 1.) return -3.*pow(1.-v,2)*exp((long double)pow(1.-v,3));
+ // else return powP*pow(v-1.,powP-1.);
}
-
-
-
-// Constructor
-OptHOM::OptHOM(const std::set<MElement*> &els, std::set<MVertex*> & toFix, bool fixBndNodes) :
- mesh(els, toFix, fixBndNodes)
+OptHOM::OptHOM(const std::set<MElement*> &els, std::set<MVertex*> & toFix,
+ bool fixBndNodes) :
+ mesh(els, toFix, fixBndNodes)
{
_optimizeMetricMin = false;
}
-
-
-// Contribution of the element Jacobians to the objective function value and gradients (2D version)
+// Contribution of the element Jacobians to the objective function value and
+// gradients (2D version)
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
- std::vector<double> gSJ(mesh.nBezEl(iEl)*mesh.nPCEl(iEl)); // Gradients of scaled Jacobians
+ // Scaled Jacobians
+ std::vector<double> sJ(mesh.nBezEl(iEl));
+ // Gradients of scaled Jacobians
+ std::vector<double> gSJ(mesh.nBezEl(iEl)*mesh.nPCEl(iEl));
mesh.scaledJacAndGradients (iEl,sJ,gSJ);
-
+
for (int l = 0; l < mesh.nBezEl(iEl); l++) {
double f1 = compute_f1(sJ[l], jacBar);
Obj += compute_f(sJ[l], jacBar);
@@ -77,7 +74,6 @@ bool OptHOM::addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj)
}
return true;
-
}
bool OptHOM::addMetricMinObjGrad(double &Obj, alglib::real_1d_array &gradObj)
@@ -87,10 +83,12 @@ bool OptHOM::addMetricMinObjGrad(double &Obj, alglib::real_1d_array &gradObj)
maxJac = -1.e300;
for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
- std::vector<double> sJ(mesh.nBezEl(iEl)); // Scaled Jacobians
- std::vector<double> gSJ(mesh.nBezEl(iEl)*mesh.nPCEl(iEl)); // Gradients of scaled Jacobians
+ // Scaled Jacobians
+ std::vector<double> sJ(mesh.nBezEl(iEl));
+ // Gradients of scaled Jacobians
+ std::vector<double> gSJ(mesh.nBezEl(iEl)*mesh.nPCEl(iEl));
mesh.metricMinAndGradients (iEl,sJ,gSJ);
-
+
for (int l = 0; l < mesh.nBezEl(iEl); l++) {
Obj += compute_f(sJ[l], jacBar);
const double f1 = compute_f1(sJ[l], jacBar);
@@ -103,12 +101,11 @@ bool OptHOM::addMetricMinObjGrad(double &Obj, alglib::real_1d_array &gradObj)
return true;
}
-
-
-// Contribution of the vertex distance to the objective function value and gradients (2D version)
-bool OptHOM::addDistObjGrad(double Fact, double Fact2, double &Obj, alglib::real_1d_array &gradObj)
+// Contribution of the vertex distance to the objective function value and
+// gradients (2D version)
+bool OptHOM::addDistObjGrad(double Fact, double Fact2, double &Obj,
+ alglib::real_1d_array &gradObj)
{
-
maxDist = 0;
avgDist = 0;
int nbBnd = 0;
@@ -119,7 +116,8 @@ bool OptHOM::addDistObjGrad(double Fact, double Fact2, double &Obj, alglib::real
Obj += Factor * dSq;
std::vector<double> gDSq(mesh.nPCFV(iFV));
mesh.gradDistSq(iFV,gDSq);
- for (int iPC = 0; iPC < mesh.nPCFV(iFV); iPC++) gradObj[mesh.indPCFV(iFV,iPC)] += Factor*gDSq[iPC];
+ 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++;
@@ -127,14 +125,11 @@ bool OptHOM::addDistObjGrad(double Fact, double Fact2, double &Obj, alglib::real
if (nbBnd != 0) avgDist /= nbBnd;
return true;
-
}
-
-
-void OptHOM::evalObjGrad(const alglib::real_1d_array &x, double &Obj, alglib::real_1d_array &gradObj)
+void OptHOM::evalObjGrad(const alglib::real_1d_array &x, double &Obj,
+ alglib::real_1d_array &gradObj)
{
-
mesh.updateMesh(x.getcontent());
Obj = 0.;
@@ -145,25 +140,21 @@ void OptHOM::evalObjGrad(const alglib::real_1d_array &x, double &Obj, alglib::re
if(_optimizeMetricMin)
addMetricMinObjGrad(Obj, gradObj);
if ((minJac > barrier_min) && (maxJac < barrier_max || !_optimizeBarrierMax)) {
- printf("INFO: reached %s (%g %g) requirements, setting null gradient\n", _optimizeMetricMin ? "svd" : "jacobian", minJac, maxJac);
+ Msg::Info("Reached %s (%g %g) requirements, setting null gradient",
+ _optimizeMetricMin ? "svd" : "jacobian", minJac, maxJac);
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)
+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);
}
-
-
void OptHOM::recalcJacDist()
- {
-
+{
maxDist = 0;
avgDist = 0;
int nbBnd = 0;
@@ -181,8 +172,10 @@ void OptHOM::recalcJacDist()
minJac = 1.e300;
maxJac = -1.e300;
for (int iEl = 0; iEl < mesh.nEl(); iEl++) {
- std::vector<double> sJ(mesh.nBezEl(iEl)); // Scaled Jacobians
- std::vector<double> dumGSJ(mesh.nBezEl(iEl)*mesh.nPCEl(iEl)); // (Dummy) gradients of scaled Jacobians
+ // Scaled Jacobians
+ std::vector<double> sJ(mesh.nBezEl(iEl));
+ // (Dummy) gradients of scaled Jacobians
+ std::vector<double> dumGSJ(mesh.nBezEl(iEl)*mesh.nPCEl(iEl));
mesh.scaledJacAndGradients (iEl,sJ,dumGSJ);
if(_optimizeMetricMin)
mesh.metricMinAndGradients (iEl,sJ,dumGSJ);
@@ -193,56 +186,45 @@ void OptHOM::recalcJacDist()
}
}
-
-
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);
+ Msg::Info("--> 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);
}
-
}
-
-
void printProgressFunc(const alglib::real_1d_array &x, double Obj, void *HOInst)
{
((OptHOM*)HOInst)->printProgress(x,Obj);
}
-
-
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];
+ for (int iPC = 0; iPC < mesh.nPCFV(iFV); iPC++)
+ scale[mesh.indPCFV(iFV,iPC)] = scaleFV[iPC];
}
-
}
-
-
-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 EPSF = 0.;
static const double EPSX = 0.;
static int OPTMETHOD = 1;
- Msg::Debug("--- Optimization pass with jacBar = %12.5E",jacBar);
- std::cout << "--- Optimization pass with initial jac. range ("
- << minJac << "," << maxJac << "), jacBar = " << jacBar << "\n";
+ Msg::Info("--- Optimization pass with initial jac. range (%g, %g), jacBar = %g",
+ minJac, maxJac, jacBar);
iter = 0;
@@ -279,23 +261,21 @@ void OptHOM::OptimPass(alglib::real_1d_array &x, const alglib::real_1d_array &in
terminationtype = rep.terminationtype;
}
- std::cout << "Optimization finalized after " << iterationscount << " iterations (" << nfev << " functions evaluations)";
+ Msg::Info("Optimization finalized after %d iterations (%d function evaluations),",
+ iterationscount, nfev);
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;
+ case 1: Msg::Info("because relative function improvement is no more than EpsF"); break;
+ case 2: Msg::Info("because relative step is no more than EpsX"); break;
+ case 4: Msg::Info("because gradient norm is no more than EpsG"); break;
+ case 5: Msg::Info("because the maximum number of steps was taken"); break;
+ default: Msg::Info("with code %d", int(terminationtype)); break;
}
- std::cout << "." << std::endl;
-
}
-
-
-int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double b_max, bool optimizeMetricMin, int pInt, int itMax, int optPassMax)
+int OptHOM::optimize(double weightFixed, double weightFree, double b_min,
+ double b_max, bool optimizeMetricMin, int pInt,
+ int itMax, int optPassMax)
{
-
barrier_min = b_min;
barrier_max = b_max;
progressInterv = pInt;
@@ -332,10 +312,9 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double
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;
+ 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;
@@ -349,7 +328,7 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double
break;
}
if (fabs((minJac-startMinJac)/startMinJac) < 0.01) {
- std::cout << "Stagnation in minJac detected, stopping optimization\n";
+ Msg::Info("Stagnation in minJac detected, stopping optimization");
minJacOK = false;
break;
}
@@ -366,20 +345,17 @@ int OptHOM::optimize(double weightFixed, double weightFree, double b_min, double
jacBar = 1.1 * maxJac;
if (ITER ++ > optPassMax) break;
if (fabs((maxJac-startMaxJac)/startMaxJac) < 0.01) {
- std::cout << "Stagnation in maxJac detected, stopping optimization\n";
+ Msg::Info("Stagnation in maxJac detected, stopping optimization");
break;
}
}
}
- OptHomMessage("Optimization done Range (%g,%g)",minJac,maxJac);
+ 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 53e60216ec7084353b8fc2d0cb86f5412f90bde2..e7bb8d4aca286f23587f52ba9a010097fc07fdb2 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.h
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.h
@@ -3,66 +3,57 @@
#include <string>
#include <math.h>
+#include "GmshConfig.h"
-
-
-#ifdef HAVE_BFGS
+#if defined(HAVE_BFGS)
#include "ap.h"
-
-
class Mesh;
-
class OptHOM
{
public:
-
Mesh mesh;
-
- OptHOM(const std::set<MElement*> &els, std::set<MVertex*> & toFix, bool fixBndNodes);
- // 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, bool optimizeMetricMin, int pInt, int itMax, int optPassMax); // optimize one list of elements
+ OptHOM(const std::set<MElement*> &els, std::set<MVertex*> & toFix,
+ bool fixBndNodes);
+ // 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,
+ bool optimizeMetricMin, int pInt, int itMax, int optPassMax);
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 evalObjGrad(const alglib::real_1d_array &x, double &Obj,
+ alglib::real_1d_array &gradObj);
void printProgress(const alglib::real_1d_array &x, double Obj);
double barrier_min, barrier_max;
-private:
-
+ private:
double lambda, lambda2, jacBar, invLengthScaleSq;
- int iter, progressInterv; // Current iteration, interval of iterations for reporting
+ int iter, progressInterv; // Current iteration, interval of iterations for reporting
bool _optimizeMetricMin;
- double initObj, initMaxDist, initAvgDist; // Values for reporting
- double minJac, maxJac, maxDist, avgDist; // Values for reporting
-
- bool _optimizeBarrierMax; // false : only moving barrier min; true : fixed barrier min + moving barrier max
-
+ double initObj, initMaxDist, initAvgDist; // Values for reporting
+ double minJac, maxJac, maxDist, avgDist; // Values for reporting
+ bool _optimizeBarrierMax; // false : only moving barrier min;
+ // true : fixed barrier min + moving barrier max
bool addJacObjGrad(double &Obj, alglib::real_1d_array &gradObj);
bool addMetricMinObjGrad(double &Obj, alglib::real_1d_array &gradObj);
- bool addDistObjGrad(double Fact, double Fact2, double &Obj, alglib::real_1d_array &gradObj);
+ bool addDistObjGrad(double Fact, double Fact2, double &Obj,
+ alglib::real_1d_array &gradObj);
void calcScale(alglib::real_1d_array &scale);
- void OptimPass(alglib::real_1d_array &x, const alglib::real_1d_array &initGradObj, int itMax);
+ void OptimPass(alglib::real_1d_array &x, const alglib::real_1d_array &initGradObj,
+ int itMax);
};
-
-
-
void OptHOM::getJacDist(double &minJ, double &maxJ, double &maxD, double &avgD)
{
minJ = minJac; maxJ = maxJac; maxD = maxDist; avgD = avgDist;
}
-
-
#endif
-
-
#endif
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index ca60318901b226100e15407f34c12c66b81a7318..72f04ce243f5faec7cafc2d9a6142d994c82d741 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -6,8 +6,6 @@
#include "ParamCoord.h"
#include "OptHomMesh.h"
-
-
Mesh::Mesh(const std::set<MElement*> &els, std::set<MVertex*> &toFix, bool fixBndNodes)
{
@@ -23,7 +21,8 @@ Mesh::Mesh(const std::set<MElement*> &els, std::set<MVertex*> &toFix, bool fixBn
Msg::Debug("METHOD: Freeing boundary nodes and using parent parametric coordinates");
}
- // Initialize elements, vertices, free vertices and element->vertices connectivity
+ // Initialize elements, vertices, free vertices and element->vertices
+ // connectivity
const int nElements = els.size();
_nPC = 0;
_el.resize(nElements);
@@ -33,7 +32,8 @@ Mesh::Mesh(const std::set<MElement*> &els, std::set<MVertex*> &toFix, bool fixBn
_nNodEl.resize(nElements);
_indPCEl.resize(nElements);
int iEl = 0;
- for(std::set<MElement*>::const_iterator it = els.begin(); it != els.end(); ++it, ++iEl) {
+ for(std::set<MElement*>::const_iterator it = els.begin();
+ it != els.end(); ++it, ++iEl) {
_el[iEl] = *it;
const JacobianBasis *jac = _el[iEl]->getJacobianFuncSpace();
_nBezEl[iEl] = jac->getNumJacNodes();
@@ -44,12 +44,15 @@ Mesh::Mesh(const std::set<MElement*> &els, std::set<MVertex*> &toFix, bool fixBn
_el2V[iEl].push_back(iV);
const int nPCV = _pc->nCoord(vert);
bool isFV = false;
- if (fixBndNodes) isFV = (vert->onWhat()->dim() == _dim) && (toFix.find(vert) == toFix.end());
- else isFV = (vert->onWhat()->dim() >= 1) && (toFix.find(vert) == toFix.end());
+ if (fixBndNodes)
+ isFV = (vert->onWhat()->dim() == _dim) && (toFix.find(vert) == toFix.end());
+ else
+ isFV = (vert->onWhat()->dim() >= 1) && (toFix.find(vert) == toFix.end());
if (isFV) {
int iFV = addFreeVert(vert,iV,nPCV,toFix);
_el2FV[iEl].push_back(iFV);
- for (int i=_startPCFV[iFV]; i<_startPCFV[iFV]+nPCV; i++) _indPCEl[iEl].push_back(i);
+ for (int i=_startPCFV[iFV]; i<_startPCFV[iFV]+nPCV; i++)
+ _indPCEl[iEl].push_back(i);
}
else _el2FV[iEl].push_back(-1);
}
@@ -65,7 +68,8 @@ Mesh::Mesh(const std::set<MElement*> &els, std::set<MVertex*> &toFix, bool fixBn
_xyz = _ixyz;
_uvw = _iuvw;
- // Set normals to 2D elements (with magnitude of inverse Jacobian) or initial Jacobians of 3D elements
+ // Set normals to 2D elements (with magnitude of inverse Jacobian) or initial
+ // Jacobians of 3D elements
if (_dim == 2) {
_scaledNormEl.resize(nEl());
for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(iEl);
@@ -73,13 +77,12 @@ Mesh::Mesh(const std::set<MElement*> &els, std::set<MVertex*> &toFix, bool fixBn
else {
_invStraightJac.resize(nEl(),1.);
double dumJac[3][3];
- for (int iEl = 0; iEl < nEl(); iEl++) _invStraightJac[iEl] = 1. / _el[iEl]->getPrimaryJacobian(0.,0.,0.,dumJac);
+ for (int iEl = 0; iEl < nEl(); iEl++)
+ _invStraightJac[iEl] = 1. / _el[iEl]->getPrimaryJacobian(0.,0.,0.,dumJac);
}
}
-
-
void Mesh::calcScaledNormalEl2D(int iEl)
{
@@ -102,11 +105,8 @@ void Mesh::calcScaledNormalEl2D(int iEl)
}
-
-
int Mesh::addVert(MVertex* vert)
{
-
std::vector<MVertex*>::iterator itVert = find(_vert.begin(),_vert.end(),vert);
if (itVert == _vert.end()) {
_vert.push_back(vert);
@@ -116,12 +116,11 @@ int Mesh::addVert(MVertex* vert)
}
-
-
-int Mesh::addFreeVert(MVertex* vert, const int iV, const int nPCV, std::set<MVertex*> &toFix)
+int Mesh::addFreeVert(MVertex* vert, const int iV, const int nPCV,
+ std::set<MVertex*> &toFix)
{
-
- std::vector<MVertex*>::iterator itVert = find(_freeVert.begin(),_freeVert.end(),vert);
+ std::vector<MVertex*>::iterator itVert = find(_freeVert.begin(),
+ _freeVert.end(),vert);
if (itVert == _freeVert.end()) {
const int iStart = (_startPCFV.size() == 0)? 0 : _startPCFV.back()+_nPCFV.back();
const bool forcedV = (vert->onWhat()->dim() < 2) || (toFix.find(vert) != toFix.end());
@@ -137,12 +136,8 @@ int Mesh::addFreeVert(MVertex* vert, const int iV, const int nPCV, std::set<MVer
}
-
-
void Mesh::getUvw(double *it)
{
-
-// std::vector<double>::iterator it = uvw.begin();
for (int iFV = 0; iFV < nFV(); iFV++) {
SPoint3 &uvwV = _uvw[iFV];
*it = uvwV[0]; it++;
@@ -152,12 +147,8 @@ void Mesh::getUvw(double *it)
}
-
-
void Mesh::updateMesh(const double *it)
{
-
-// std::vector<double>::const_iterator it = uvw.begin();
for (int iFV = 0; iFV < nFV(); iFV++) {
int iV = _fv2V[iFV];
SPoint3 &uvwV = _uvw[iFV];
@@ -169,7 +160,6 @@ void Mesh::updateMesh(const double *it)
}
-
void Mesh::distSqToStraight(std::vector<double> &dSq)
{
std::vector<SPoint3> sxyz(nVert());
@@ -185,7 +175,8 @@ void Mesh::distSqToStraight(std::vector<double> &dSq)
int dim = lagrange->points.size2();
for (int i = nV1; i < nV; ++i) {
double f[256];
- lagrange1->f(lagrange->points(i, 0), dim > 1 ? lagrange->points(i, 1) : 0., dim > 2 ? lagrange->points(i, 2) : 0., f);
+ lagrange1->f(lagrange->points(i, 0), dim > 1 ? lagrange->points(i, 1) : 0.,
+ dim > 2 ? lagrange->points(i, 2) : 0., f);
for (int j = 0; j < nV1; ++j)
sxyz[_el2V[iEl][i]] += sxyz[_el2V[iEl][j]] * f[j];
}
@@ -197,28 +188,25 @@ void Mesh::distSqToStraight(std::vector<double> &dSq)
}
}
-
-
void Mesh::updateGEntityPositions()
{
for (int iV = 0; iV < nVert(); iV++)
- _vert[iV]->setXYZ(_xyz[iV].x(),_xyz[iV].y(),_xyz[iV].z());
+ _vert[iV]->setXYZ(_xyz[iV].x(),_xyz[iV].y(),_xyz[iV].z());
for (int iFV = 0; iFV < nFV(); iFV++)
_pc->exportParamCoord(_freeVert[iFV], _uvw[iFV]);
}
-
-
-void Mesh::metricMinAndGradients(int iEl, std::vector<double> &lambda , std::vector<double> &gradLambda)
+void Mesh::metricMinAndGradients(int iEl, std::vector<double> &lambda,
+ std::vector<double> &gradLambda)
{
-
const JacobianBasis *jacBasis = _el[iEl]->getJacobianFuncSpace();
const int &numJacNodes = jacBasis->getNumJacNodes();
- const int &numMapNodes = jacBasis->getNumMapNodes(), &numPrimMapNodes = jacBasis->getNumPrimMapNodes();
+ const int &numMapNodes = jacBasis->getNumMapNodes();
+ const int &numPrimMapNodes = jacBasis->getNumPrimMapNodes();
fullVector<double> lambdaJ(numJacNodes), lambdaB(numJacNodes);
fullMatrix<double> gradLambdaJ(numJacNodes, 2 * numMapNodes);
fullMatrix<double> gradLambdaB(numJacNodes, 2 * numMapNodes);
-
+
// Coordinates of nodes
fullMatrix<double> nodesXYZ(numMapNodes,3), nodesXYZStraight(numPrimMapNodes,3);
for (int i = 0; i < numMapNodes; i++) {
@@ -234,7 +222,7 @@ void Mesh::metricMinAndGradients(int iEl, std::vector<double> &lambda , std::vec
}
jacBasis->getMetricMinAndGradients(nodesXYZ,nodesXYZStraight,lambdaJ,gradLambdaJ);
-
+
//l2b.mult(lambdaJ, lambdaB);
//l2b.mult(gradLambdaJ, gradLambdaB);
lambdaB = lambdaJ;
@@ -250,7 +238,8 @@ void Mesh::metricMinAndGradients(int iEl, std::vector<double> &lambda , std::vec
int &iFVi = _el2FV[iEl][i];
if (iFVi >= 0) {
for (int l = 0; l < numJacNodes; l++) {
- gXyzV [l] = SPoint3(gradLambdaB(l,i+0*numMapNodes),gradLambdaB(l,i+1*numMapNodes),/*BDB(l,i+2*nbNod)*/ 0.);
+ gXyzV [l] = SPoint3(gradLambdaB(l,i+0*numMapNodes),
+ gradLambdaB(l,i+1*numMapNodes),/*BDB(l,i+2*nbNod)*/ 0.);
}
_pc->gXyz2gUvw(_freeVert[iFVi],_uvw[iFVi],gXyzV,gUvwV);
for (int l = 0; l < numJacNodes; l++) {
@@ -263,16 +252,15 @@ void Mesh::metricMinAndGradients(int iEl, std::vector<double> &lambda , std::vec
}
}
-
-
-void Mesh::scaledJacAndGradients(int iEl, std::vector<double> &sJ , std::vector<double> &gSJ) {
-
+void Mesh::scaledJacAndGradients(int iEl, std::vector<double> &sJ,
+ std::vector<double> &gSJ)
+{
const JacobianBasis *jacBasis = _el[iEl]->getJacobianFuncSpace();
const int &numJacNodes = jacBasis->getNumJacNodes();
const int &numMapNodes = jacBasis->getNumMapNodes();
fullMatrix<double> JDJ (numJacNodes,3*numMapNodes+1);
fullMatrix<double> BDB (numJacNodes,3*numMapNodes+1);
-
+
// Coordinates of nodes
fullMatrix<double> nodesXYZ(numMapNodes,3), normals(_dim,3);
for (int i = 0; i < numMapNodes; i++) {
@@ -282,10 +270,11 @@ void Mesh::scaledJacAndGradients(int iEl, std::vector<double> &sJ , std::vector<
nodesXYZ(i,2) = _xyz[iVi].z();
}
- // Calculate Jacobian and gradients, scale if 3D (already scaled by regularization normals in 2D)
+ // Calculate Jacobian and gradients, scale if 3D (already scaled by
+ // regularization normals in 2D)
jacBasis->getSignedJacAndGradients(nodesXYZ,_scaledNormEl[iEl],JDJ);
if (_dim == 3) JDJ.scale(_invStraightJac[iEl]);
-
+
// Transform Jacobian and gradients from Lagrangian to Bezier basis
jacBasis->lag2Bez(JDJ,BDB);
@@ -300,7 +289,8 @@ void Mesh::scaledJacAndGradients(int iEl, std::vector<double> &sJ , std::vector<
int &iFVi = _el2FV[iEl][i];
if (iFVi >= 0) {
for (int l = 0; l < numJacNodes; l++)
- gXyzV [l] = SPoint3(BDB(l,i+0*numMapNodes),BDB(l,i+1*numMapNodes),BDB(l,i+2*numMapNodes));
+ gXyzV [l] = SPoint3(BDB(l,i+0*numMapNodes), BDB(l,i+1*numMapNodes),
+ BDB(l,i+2*numMapNodes));
_pc->gXyz2gUvw(_freeVert[iFVi],_uvw[iFVi],gXyzV,gUvwV);
for (int l = 0; l < numJacNodes; l++) {
gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][0];
@@ -308,16 +298,13 @@ void Mesh::scaledJacAndGradients(int iEl, std::vector<double> &sJ , std::vector<
if (_nPCFV[iFVi] == 3) gSJ[indGSJ(iEl,l,iPC+2)] = gUvwV[l][2];
}
iPC += _nPCFV[iFVi];
- }
+ }
}
}
-
-
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;
@@ -329,11 +316,8 @@ void Mesh::pcScale(int iFV, std::vector<double> &scale)
scale[0] = 1./sqrt(gX[0]*gX[0]+gY[0]*gY[0]+gZ[0]*gZ[0]);
if (_nPCFV[iFV] >= 2) scale[1] = 1./sqrt(gX[1]*gX[1]+gY[1]*gY[1]+gZ[1]*gZ[1]);
if (_nPCFV[iFV] == 3) scale[2] = 1./sqrt(gX[2]*gX[2]+gY[2]*gY[2]+gZ[2]*gZ[2]);
-
}
-
-
void Mesh::writeMSH(const char *filename)
{
FILE *f = fopen(filename, "w");
@@ -351,13 +335,12 @@ void Mesh::writeMSH(const char *filename)
fprintf(f, "$Elements\n");
fprintf(f, "%d\n", nEl());
for (int iEl = 0; iEl < nEl(); iEl++) {
-// MElement *MEl = _el[iEl];
fprintf(f, "%d %d 2 0 0", iEl+1, _el[iEl]->getTypeForMSH());
- for (size_t iVEl = 0; iVEl < _el2V[iEl].size(); iVEl++) fprintf(f, " %d", _el2V[iEl][iVEl] + 1);
+ for (size_t iVEl = 0; iVEl < _el2V[iEl].size(); iVEl++)
+ fprintf(f, " %d", _el2V[iEl][iVEl] + 1);
fprintf(f, "\n");
}
fprintf(f, "$EndElements\n");
fclose(f);
-
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 78fd0f844d993929187aaee2ca4467f9a1f11827..7ee49c94b42f9240acc9cc120d8d5f864d385233 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -1,4 +1,3 @@
-//#include <fenv.h>
#include <stdio.h>
#include <sstream>
#include <iterator>
@@ -23,32 +22,7 @@
#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);
-#else
- fprintf(stdout,"%s\n", str);
-#endif
-}
-
-double distMaxStraight (MElement *el)
+double distMaxStraight(MElement *el)
{
const polynomialBasis *lagrange = (polynomialBasis*)el->getFunctionSpace();
const polynomialBasis *lagrange1 = (polynomialBasis*)el->getFunctionSpace(1);
@@ -76,7 +50,7 @@ double distMaxStraight (MElement *el)
return maxdx;
}
-void exportMeshToDassault (GModel *gm, const std::string &fn, int dim)
+void exportMeshToDassault(GModel *gm, const std::string &fn, int dim)
{
FILE *f = fopen(fn.c_str(),"w");
@@ -224,10 +198,10 @@ static std::vector<std::pair<std::set<MElement*>, std::set<MVertex*> > > getConn
bool weakMerge = false)
{
- OptHomMessage("Starting blob generation from %i bad elements...",badElements.size());
+ Msg::Info("Starting blob generation from %i bad elements...", badElements.size());
// Contruct primary blobs
- std::cout << "Constructing " << badElements.size() << " primary blobs...\n";
+ Msg::Info("Constructing %i primary blobs", badElements.size());
std::vector<std::set<MElement*> > primBlobs;
primBlobs.reserve(badElements.size());
for (std::set<MElement*>::const_iterator it = badElements.begin();
@@ -235,7 +209,7 @@ static std::vector<std::pair<std::set<MElement*>, std::set<MVertex*> > > getConn
primBlobs.push_back(getSurroundingBlob(*it, depth, vertex2elements, distFactor));
// Compute blob connectivity
- std::cout << "Computing blob connectivity...\n";
+ Msg::Info("Computing blob connectivity...");
std::map<MElement*, std::set<int> > tags;
std::vector<std::set<int> > blobConnect(primBlobs.size());
for (int iB = 0; iB < primBlobs.size(); ++iB) {
@@ -253,7 +227,7 @@ static std::vector<std::pair<std::set<MElement*>, std::set<MVertex*> > > getConn
}
// Identify groups of connected blobs
- std::cout << "Identifying groups of primary blobs...\n";
+ Msg::Info("Identifying groups of primary blobs...");
std::list<std::set<int> > groups;
std::vector<bool> todoPB(primBlobs.size(), true);
for (int iB = 0; iB < primBlobs.size(); ++iB)
@@ -264,7 +238,7 @@ static std::vector<std::pair<std::set<MElement*>, std::set<MVertex*> > > getConn
}
// Merge primary blobs according to groups
- std::cout << "Merging primary blobs into " << groups.size() << " blobs...\n";
+ Msg::Info("Merging primary blobs into %i blobs...", groups.size());
std::list<std::set<MElement*> > blobs;
for (std::list<std::set<int> >::iterator itG = groups.begin();
itG != groups.end(); ++itG) {
@@ -276,14 +250,14 @@ static std::vector<std::pair<std::set<MElement*>, std::set<MVertex*> > > getConn
}
// Store and compute blob boundaries
- std::cout << "Computing boundaries for " << blobs.size() << " blobs...\n";
+ Msg::Info("Computing boundaries for %i blobs...", blobs.size());
std::vector<std::pair<std::set<MElement *>, std::set<MVertex*> > > result;
for (std::list<std::set<MElement*> >::iterator itB = blobs.begin();
itB != blobs.end(); ++itB)
result.push_back(std::pair<std::set<MElement*>, std::set<MVertex*> >
(*itB, getAllBndVertices(*itB, vertex2elements)));
- OptHomMessage("Generated %i blobs",blobs.size());
+ Msg::Info("Generated %i blobs", blobs.size());
return result;
}
@@ -302,7 +276,8 @@ static void optimizeConnectedBlobs
//#pragma omp parallel for schedule(dynamic, 1)
for (int i = 0; i < toOptimize.size(); ++i) {
- OptHomMessage("Optimizing a blob %i/%i composed of %4d elements", i+1, toOptimize.size(), toOptimize[i].first.size());
+ Msg::Info("Optimizing a blob %i/%i composed of %4d elements", i+1,
+ toOptimize.size(), toOptimize[i].first.size());
fflush(stdout);
OptHOM temp(toOptimize[i].first, toOptimize[i].second, p.fixBndNodes);
//std::ostringstream ossI1;
@@ -311,7 +286,7 @@ static void optimizeConnectedBlobs
int success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX,
false, samples, p.itMax, p.optPassMax);
if (success >= 0 && p.BARRIER_MIN_METRIC > 0) {
- OptHomMessage("Jacobian optimization succeed, starting svd optimization");
+ Msg::Info("Jacobian optimization succeed, starting svd optimization");
success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC, p.BARRIER_MAX,
true, samples, p.itMax, p.optPassMax);
}
@@ -460,7 +435,7 @@ static void optimizeOneByOne
p.maxJac = -1.e100;
const int initNumBadElts = badElts.size();
- std::cout << "DBGTT: " << initNumBadElts << " badasses, starting to iterate...\n";
+ Msg::Info("%d badasses, starting to iterate...", initNumBadElts);
for (int iBadEl=0; iBadEl<initNumBadElts; iBadEl++) {
@@ -488,8 +463,8 @@ static void optimizeOneByOne
std::set_difference(toOptimizePrim.begin(),toOptimizePrim.end(),
badElts.begin(),badElts.end(), // Do not optimize badElts
std::inserter(toOptimize1, toOptimize1.end()));
- OptHomMessage("Optimizing primary blob %i (max. %i remaining) composed of"
- " %4d elements", iBadEl, badElts.size(), toOptimize1.size());
+ Msg::Info("Optimizing primary blob %i (max. %i remaining) composed of"
+ " %4d elements", iBadEl, badElts.size(), toOptimize1.size());
fflush(stdout);
opt = new OptHOM(toOptimize1, toFix1, p.fixBndNodes);
//std::ostringstream ossI1;
@@ -510,9 +485,9 @@ static void optimizeOneByOne
std::set_difference(toOptimizePrim.begin(),toOptimizePrim.end(),
badElts.begin(),badElts.end(),
std::inserter(toOptimize2, toOptimize2.end()));
- OptHomMessage("Optimizing corrective blob %i (max. %i remaining) "
- "composed of %4d elements", iBadEl, badElts.size(),
- toOptimize2.size());
+ Msg::Info("Optimizing corrective blob %i (max. %i remaining) "
+ "composed of %4d elements", iBadEl, badElts.size(),
+ toOptimize2.size());
fflush(stdout);
opt = new OptHOM(toOptimize2, toFix2, p.fixBndNodes);
//std::ostringstream ossI1;
@@ -521,14 +496,14 @@ static void optimizeOneByOne
success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN,
p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax);
if (success >= 0 && p.BARRIER_MIN_METRIC > 0) {
- OptHomMessage("Jacobian optimization succeed, starting svd optimization");
+ Msg::Info("Jacobian optimization succeed, starting svd optimization");
success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC,
p.BARRIER_MAX, true, samples, p.itMax, p.optPassMax);
}
}
else {
- OptHomMessage("Primary blob %i did not break new elements, "
- "no correction needed", iBadEl);
+ Msg::Info("Primary blob %i did not break new elements, "
+ "no correction needed", iBadEl);
fflush(stdout);
}
}
@@ -549,8 +524,8 @@ static void optimizeOneByOne
if (success <= 0) {
distanceFactor *= p.adaptBlobDistFact;
nbLayers *= p.adaptBlobLayerFact;
- OptHomMessage("Blob %i failed (adapt #%i), adapting with increased size",
- iBadEl, iterBlob);
+ Msg::Info("Blob %i failed (adapt #%i), adapting with increased size",
+ iBadEl, iterBlob);
if (iterBlob == p.maxAdaptBlob-1) {
std::ostringstream ossI2;
ossI2 << "final_" << iBadEl << ".msh";
@@ -566,189 +541,7 @@ static void optimizeOneByOne
}
}
-//static void optimizeOneByOneTest(const std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
-// std::set<MElement*> badasses, OptHomParameters &p, int method, int samples)
-//{
-//
-// p.SUCCESS = 1;
-// p.minJac = 1.e100;
-// p.maxJac = -1.e100;
-//
-// const int initNumBadAsses = badasses.size();
-// std::cout << "DBGTT: " << initNumBadAsses << " badasses, starting to iterate...\n";
-//
-// for (int iter=0; iter<initNumBadAsses; iter++) {
-//
-// if (badasses.empty()) break;
-//
-// // Create blob around worst element
-// MElement *worstEl = getWorstElement(badasses,p);
-//
-// const int maxRepeatBlob = 1;
-// const int repeatBlobLayerFact = 2;
-// const int repeatBlobDistFact = 2.;
-//
-// int nbLayers = p.nbLayers, minNbLayers = 1;
-// double distanceFactor = p.distanceFactor;
-//
-// int success;
-//
-// for (int iterBlob=0; iterBlob<maxRepeatBlob; iterBlob++) {
-//
-// // TODO: First opt. with only 1st-order bnd. vertices fixed,
-// // then add external layer and opt. with all bnd. vert. fixed
-//// std::set<MElement*> toOptimizePrim = getSurroundingBlob(worstEl, nbLayers, vertex2elements, distanceFactor, 0);
-//// std::set<MVertex*> toFix = getBndVertices3D(toOptimizePrim, vertex2elements);
-// std::set<MElement*> toOptimizePrim = getSurroundingBlob3D(worstEl, nbLayers, vertex2elements, distanceFactor);
-// std::set<MVertex*> toFix = getAllBndVertices(toOptimizePrim, vertex2elements);
-// badasses.erase(worstEl);
-// std::set<MElement*> toOptimize;
-// std::set_difference(toOptimizePrim.begin(),toOptimizePrim.end(),badasses.begin(),badasses.end(),
-// std::inserter(toOptimize, toOptimize.end()));
-//
-// // Optimize blob
-// OptHomMessage("Optimizing blob %i (max. %i remaining) composed of %4d elements", iter,
-// badasses.size(), toOptimize.size());
-// fflush(stdout);
-// OptHOM temp(toOptimize, toFix, method);
-////std::cout << "DBGTT: toOptimize:\n";
-////for (std::set<MElement*>::iterator it=toOptimize.begin(); it!=toOptimize.end(); it++ ) {
-//// SPoint3 bar = (*it)->barycenter();
-//// std::cout << "DBGTT: -> (" << bar.x() << "," << bar.y() << "," << bar.z() << ")\n";
-//// std::cout << "DBGTT: -> (" << (*it)->getVertex(0)->onWhat()->dim();
-//// for (int j=1; j<(*it)->getNumVertices(); j++) std::cout << "," << (*it)->getVertex(j)->onWhat()->dim();
-//// std::cout << ")\n";
-////}
-////std::cout << "DBGTT: toFix:\n";
-////for (std::set<MVertex*>::iterator it=toFix.begin(); it!=toFix.end(); it++ )
-//// std::cout << "DBGTT: -> (" << (*it)->x() << "," << (*it)->y() << "," << (*it)->z() << ")\n";
-//std::ostringstream ossI1;
-//ossI1 << "initial_ITER_" << iter << ".msh";
-//temp.mesh.writeMSH(ossI1.str().c_str());
-// success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX,
-// false, samples, p.itMax, p.optPassMax);
-// if (success >= 0 && p.BARRIER_MIN_METRIC > 0) {
-// OptHomMessage("Jacobian optimization succeed, starting svd optimization");
-// success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC, p.BARRIER_MAX,
-// true, samples, p.itMax, p.optPassMax);
-// }
-//
-// // Measure min and max Jac, update mesh
-// double minJac, maxJac, distMaxBND, distAvgBND;
-// temp.recalcJacDist();
-// temp.getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
-// p.minJac = std::min(p.minJac,minJac);
-// p.maxJac = std::max(p.maxJac,maxJac);
-// temp.mesh.updateGEntityPositions();
-//
-////std::ostringstream ossI2;
-////ossI2 << "final_ITER_" << iter << ".msh";
-////temp.mesh.writeMSH(ossI2.str().c_str());
-// if (success <= 0) {
-// distanceFactor *= repeatBlobDistFact;
-// nbLayers *= repeatBlobLayerFact;
-// OptHomMessage("Blob %i failed (repeat %i), repeating with increased size", iter, iterBlob);
-// if (iterBlob == maxRepeatBlob-1) {
-// std::ostringstream ossI2;
-// ossI2 << "final_ITER_" << iter << ".msh";
-// temp.mesh.writeMSH(ossI2.str().c_str());
-// }
-// }
-// else break;
-//
-// }
-//
-// //#pragma omp critical
-// p.SUCCESS = std::min(p.SUCCESS, success);
-//
-// }
-//
-//}
-//
-//
-//
-//static void optimizePropagate(const std::map<MVertex*, std::vector<MElement*> > &vertex2elements,
-// std::set<MElement*> badasses, OptHomParameters &p, int method, int samples)
-//{
-//
-// p.SUCCESS = 1;
-// p.minJac = 1.e100;
-// p.maxJac = -1.e100;
-//
-// const int initNumBadAsses = badasses.size();
-// std::cout << "DBGTT: " << initNumBadAsses << " badasses, starting to iterate...\n";
-//
-// std::set<MVertex*> toFix;
-// std::set<MElement*> done;
-//
-// while (!badasses.empty()) {
-//
-// OptHomMessage("Optimizing first of %i bad elements, %i already done", badasses.size(), done.size());
-// fflush(stdout);
-////char dum;
-////std::cin >> dum;
-//
-// MElement *el = *badasses.begin();
-//
-// // Optimize element
-// std::set<MElement*> blob;
-// blob.insert(el);
-// OptHOM temp(blob, toFix, method);
-////std::ostringstream ossI1;
-////ossI1 << "initial.msh";
-////temp.mesh.writeMSH(ossI1.str().c_str());
-// int success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX,
-// false, samples, p.itMax, p.optPassMax);
-// if (success <= 0) {
-// OptHomMessage("Error, element optimization failed");
-// fflush(stdout);
-// std::ostringstream ossI1;
-// ossI1 << "failed_" << el->getNum() << ".msh";
-// temp.mesh.writeMSH(ossI1.str().c_str());
-// }
-////std::ostringstream ossI2;
-////ossI2 << "final.msh";
-////temp.mesh.writeMSH(ossI2.str().c_str());
-//
-// // Move element from list of bad elements to liste of elements done
-// badasses.erase(el);
-// done.insert(el);
-//std::cout << "DBGTT: Removing el. " << el->getNum() << "\n";
-//
-// // Measure success, min and max Jac and update mesh
-// double minJac, maxJac, distMaxBND, distAvgBND;
-// temp.recalcJacDist();
-// temp.getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
-// p.minJac = std::min(p.minJac,minJac);
-// p.maxJac = std::max(p.maxJac,maxJac);
-// p.SUCCESS = std::min(p.SUCCESS, success);
-// temp.mesh.updateGEntityPositions();
-//
-// // Fix elements nodes
-// for (int i = 0; i < el->getNumVertices(); ++i) toFix.insert(el->getVertex(i));
-//
-// // Add elements that have been broken by optimization to list of bad elements
-// std::set<MElement*> layer;
-// for (int i = 0; i < el->getNumPrimaryVertices(); ++i) {
-// const std::vector<MElement*> &neighbours = vertex2elements.find(el->getVertex(i))->second;
-// for (size_t k = 0; k < neighbours.size(); ++k) layer.insert(neighbours[k]);
-// }
-// for (std::set<MElement*>::const_iterator it = layer.begin(); it != layer.end(); ++it)
-// if (done.find(*it) == done.end()) { // Add only if not done
-// double jmin, jmax;
-// (*it)->scaledJacRange(jmin, jmax);
-// if (p.BARRIER_MIN_METRIC > 0) jmax = jmin;
-// if (jmin < p.BARRIER_MIN || jmax > p.BARRIER_MAX) {
-//std::cout << "DBGTT: Adding el. " << (*it)->getNum() << "\n";
-// badasses.insert(*it);
-// }
-// }
-//
-// }
-//
-//}
-
-void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
+void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p)
{
double t1 = Cpu();
@@ -756,7 +549,7 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
double tf1 = Cpu();
- OptHomMessage("High order mesh optimizer starts");
+ Msg::StatusBar(true, "Optimizing high order mesh...");
std::vector<GEntity*> entities;
gm->getEntities(entities);
@@ -765,11 +558,9 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
for (int iEnt = 0; iEnt < entities.size(); ++iEnt) {
GEntity* &entity = entities[iEnt];
if (entity->dim() != p.dim || (p.onlyVisible && !entity->getVisibility())) continue;
- std::cout << "Computing connectivity and bad elements for entity "
- << entity->tag() << " ...\n";
+ Msg::Info("Computing connectivity and bad elements for entity %d...",
+ entity->tag());
calcVertex2Elements(p.dim,entity,vertex2elements); // Compute vert. -> elt. connectivity
- //std::cout << " -> " << entity->getNumMeshElements()
- // << " elements, vertex2elements.size() = " << vertex2elements.size() << "\n";
for (int iEl = 0; iEl < entity->getNumMeshElements();iEl++) { // Detect bad elements
double jmin, jmax;
MElement *el = entity->getMeshElement(iEl);
@@ -788,15 +579,15 @@ void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p)
else
optimizeOneByOne(vertex2elements, badasses, p, samples);
- double DTF = Cpu()-tf1;
if (p.SUCCESS == 1)
- OptHomMessage("Optimization succeeded (CPU %g sec)", DTF);
+ Msg::Info("Optimization succeeded");
else if (p.SUCCESS == 0)
- OptHomMessage("Warning : All jacobians positive but not all in the range"
- " (CPU %g sec)", DTF);
+ Msg::Warning("All jacobians positive but not all in the range");
else if (p.SUCCESS == -1)
- OptHomMessage("Error : Still negative jacobians (CPU %g sec)", DTF);
+ Msg::Error("Still negative jacobians");
double t2 = Cpu();
p.CPU = t2-t1;
+
+ Msg::StatusBar(true, "Done optimizing high order mesh (%g s)", p.CPU);
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.h b/contrib/HighOrderMeshOptimizer/OptHomRun.h
index 2382923312a25b0ddce40c373af01dc368c93057..4d9c207e9ebf4d73aee3c71c9092897db169fc96 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.h
@@ -16,12 +16,13 @@ struct OptHomParameters {
int optPassMax ; // max number of optimization passes
double TMAX ; // max CPU time allowed
bool onlyVisible ; // apply optimization to visible entities ONLY
- double distanceFactor; // filter elements such that no elements further away than
- // DistanceFactor times the max distance to straight sided version of an element are optimized
+ double distanceFactor; // filter elements such that no elements further away
+ // than DistanceFactor times the max distance to
+ // straight sided version of an element are optimized
bool fixBndNodes; // how jacobians are computed and if points can move on boundaries
- int strategy; // 0 = connected blobs, 1 = adaptive one-by-one
- int maxAdaptBlob; // Max. nb. of blob adaptation interations
- int adaptBlobLayerFact; // Growth factor in number of layers for blob adaptation
+ int strategy; // 0 = connected blobs, 1 = adaptive one-by-one
+ int maxAdaptBlob; // Max. nb. of blob adaptation interations
+ int adaptBlobLayerFact; // Growth factor in number of layers for blob adaptation
double adaptBlobDistFact; // Growth factor in distance factor for blob adaptation
// OUTPUT ------>
@@ -29,20 +30,15 @@ struct OptHomParameters {
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
- : BARRIER_MIN_METRIC(-1.), BARRIER_MIN(0.1), BARRIER_MAX(2.0), weightFixed(1.e6), weightFree (1.e2),
- nbLayers (6) , dim(3) , itMax(300), onlyVisible(true), distanceFactor(12), fixBndNodes(false),
- strategy(0), maxAdaptBlob(3), adaptBlobLayerFact(2.), adaptBlobDistFact(2.)
+ : BARRIER_MIN_METRIC(-1.), BARRIER_MIN(0.1), BARRIER_MAX(2.0), weightFixed(1.e6),
+ weightFree (1.e2), nbLayers (6) , dim(3) , itMax(300), onlyVisible(true),
+ distanceFactor(12), fixBndNodes(false), strategy(0), maxAdaptBlob(3),
+ adaptBlobLayerFact(2.), adaptBlobDistFact(2.)
{
}
};
-void HighOrderMeshOptimizer (GModel *gm, OptHomParameters &p);
-void OptHomMessage (const char *s, ...);
-
+void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p);
#endif
diff --git a/contrib/HighOrderMeshOptimizer/ParamCoord.cpp b/contrib/HighOrderMeshOptimizer/ParamCoord.cpp
index b8f2518e194f6418097daa71a99ccd85b03540aa..6dac43cf0037e30b0043599b5ce8217f3e5ea7cc 100644
--- a/contrib/HighOrderMeshOptimizer/ParamCoord.cpp
+++ b/contrib/HighOrderMeshOptimizer/ParamCoord.cpp
@@ -3,15 +3,14 @@
#include "GEdge.h"
#include "MVertex.h"
#include "ParamCoord.h"
-
-
+#include "GmshMessage.h"
SPoint3 ParamCoordParent::getUvw(MVertex* vert)
{
-
GEntity *ge = vert->onWhat();
- if ((ge->geomType() == GEntity::DiscreteCurve) || (ge->geomType() == GEntity::DiscreteSurface))
- std::cout << "ERROR: using parent coordinates on discrete curve or surface" << std::endl;
+ if ((ge->geomType() == GEntity::DiscreteCurve) ||
+ (ge->geomType() == GEntity::DiscreteSurface))
+ Msg::Error("Using parent coordinates on discrete curve or surface");
switch (ge->dim()) {
case 1: {
@@ -31,14 +30,10 @@ SPoint3 ParamCoordParent::getUvw(MVertex* vert)
break;
}
}
-
}
-
-
SPoint3 ParamCoordParent::uvw2Xyz(MVertex* vert, const SPoint3 &uvw)
{
-
GEntity *ge = vert->onWhat();
switch (ge->dim()) {
@@ -57,14 +52,11 @@ SPoint3 ParamCoordParent::uvw2Xyz(MVertex* vert, const SPoint3 &uvw)
break;
}
}
-
}
-
-
-void ParamCoordParent::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw)
+void ParamCoordParent::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw,
+ const SPoint3 &gXyz, SPoint3 &gUvw)
{
-
GEntity *ge = vert->onWhat();
switch (ge->dim()) {
@@ -75,8 +67,10 @@ void ParamCoordParent::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const SPoint
}
case 2: {
Pair<SVector3,SVector3> der = static_cast<GFace*>(ge)->firstDer(SPoint2(uvw[0],uvw[1]));
- gUvw[0] = gXyz.x() * der.first().x() + gXyz.y() * der.first().y() + gXyz.z() * der.first().z();
- gUvw[1] = gXyz.x() * der.second().x() + gXyz.y() * der.second().y() + gXyz.z() * der.second().z();
+ gUvw[0] = gXyz.x() * der.first().x() + gXyz.y() * der.first().y() +
+ gXyz.z() * der.first().z();
+ gUvw[1] = gXyz.x() * der.second().x() + gXyz.y() * der.second().y() +
+ gXyz.z() * der.second().z();
break;
}
case 3: {
@@ -87,9 +81,9 @@ void ParamCoordParent::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const SPoint
}
-
-
-void ParamCoordParent::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw)
+void ParamCoordParent::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw,
+ const std::vector<SPoint3> &gXyz,
+ std::vector<SPoint3> &gUvw)
{
GEntity *ge = vert->onWhat();
@@ -97,26 +91,32 @@ void ParamCoordParent::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const std::v
switch (ge->dim()) {
case 1: {
SVector3 der = static_cast<GEdge*>(ge)->firstDer(uvw[0]);
- std::vector<SPoint3>::iterator itUvw=gUvw.begin();
- for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin(); itXyz != gXyz.end(); itXyz++) {
+ std::vector<SPoint3>::iterator itUvw = gUvw.begin();
+ for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin();
+ itXyz != gXyz.end(); itXyz++) {
(*itUvw)[0] = itXyz->x() * der.x() + itXyz->y() * der.y() + itXyz->z() * der.z();
itUvw++;
}
break;
}
case 2: {
- Pair<SVector3,SVector3> der = static_cast<GFace*>(ge)->firstDer(SPoint2(uvw[0],uvw[1]));
+ Pair<SVector3,SVector3> der = static_cast<GFace*>(ge)->firstDer
+ (SPoint2(uvw[0],uvw[1]));
std::vector<SPoint3>::iterator itUvw=gUvw.begin();
- for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin(); itXyz != gXyz.end(); itXyz++) {
- (*itUvw)[0] = itXyz->x() * der.first().x() + itXyz->y() * der.first().y() + itXyz->z() * der.first().z();
- (*itUvw)[1] = itXyz->x() * der.second().x() + itXyz->y() * der.second().y() + itXyz->z() * der.second().z();
+ for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin();
+ itXyz != gXyz.end(); itXyz++) {
+ (*itUvw)[0] = itXyz->x() * der.first().x() + itXyz->y() * der.first().y() +
+ itXyz->z() * der.first().z();
+ (*itUvw)[1] = itXyz->x() * der.second().x() + itXyz->y() * der.second().y() +
+ itXyz->z() * der.second().z();
itUvw++;
}
break;
}
case 3: {
std::vector<SPoint3>::iterator itUvw=gUvw.begin();
- for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin(); itXyz != gXyz.end(); itXyz++) {
+ for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin();
+ itXyz != gXyz.end(); itXyz++) {
*itUvw = *itXyz;
itUvw++;
}
@@ -126,8 +126,6 @@ void ParamCoordParent::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const std::v
}
-
-
void ParamCoordParent::exportParamCoord(MVertex *v, const SPoint3 &uvw)
{
for (int d = 0; d < v->onWhat()->dim(); ++d) {