From d992e38e9767f1f0ce1c4a19e4697cd129bcd9bd Mon Sep 17 00:00:00 2001
From: Sebastien Blaise <sebastien.blaise@uclouvain.be>
Date: Thu, 22 Oct 2015 07:57:34 +0000
Subject: [PATCH] meshOptimize: fixes + added ncurses interface
---
contrib/MeshOptimizer/CMakeLists.txt | 20 ++
contrib/MeshOptimizer/MeshOpt.cpp | 112 ++++++++-
contrib/MeshOptimizer/MeshOpt.h | 8 +-
contrib/MeshOptimizer/MeshOptCommon.cpp | 215 +++++++++++++++-
contrib/MeshOptimizer/MeshOptCommon.h | 32 ++-
.../MeshOptObjectiveFunction.cpp | 24 +-
.../MeshOptimizer/MeshOptObjectiveFunction.h | 2 +
contrib/MeshOptimizer/MeshOptimizer.cpp | 237 ++++++++++++++++--
.../MeshOptimizer/MeshOptimizerConfig.h.in | 6 +
.../MeshOptimizer/doc/mesh_optimization.tex | 2 +-
10 files changed, 618 insertions(+), 40 deletions(-)
create mode 100644 contrib/MeshOptimizer/MeshOptimizerConfig.h.in
diff --git a/contrib/MeshOptimizer/CMakeLists.txt b/contrib/MeshOptimizer/CMakeLists.txt
index f96d3f6053..b30ab35afa 100644
--- a/contrib/MeshOptimizer/CMakeLists.txt
+++ b/contrib/MeshOptimizer/CMakeLists.txt
@@ -14,5 +14,25 @@ set(SRC
MeshOptVertexCoord.cpp
)
+
+
+opt(NCURSES "Enable ncurses console tools" OFF)
+
+if(ENABLE_NCURSES)
+ find_library(NCURSES_LIB ncurses)
+ find_path(NCURSES_INC "ncurses.h" PATH_SUFFIXES src include)
+endif(ENABLE_NCURSES)
+if(NCURSES_LIB AND NCURSES_INC)
+ set_config_option(HAVE_NCURSES "NCURSES")
+ list(APPEND EXTERNAL_LIBRARIES ${NCURSES_LIB})
+ set(EXTERNAL_LIBRARIES ${EXTERNAL_LIBRARIES} PARENT_SCOPE)
+ list(APPEND EXTERNAL_INCLUDES ${NCURSES_INC})
+endif(NCURSES_LIB AND NCURSES_INC)
+
+mark_as_advanced(NCURSES_INC NCURSES_LIB)
+
file(GLOB_RECURSE HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.hpp)
append_gmsh_src(contrib/MeshOptimizer "${SRC};${HDR}")
+
+configure_file(${CMAKE_CURRENT_SOURCE_DIR}/MeshOptimizerConfig.h.in
+ ${CMAKE_CURRENT_BINARY_DIR}/MeshOptimizerConfig.h)
diff --git a/contrib/MeshOptimizer/MeshOpt.cpp b/contrib/MeshOptimizer/MeshOpt.cpp
index a35c3125ee..ad8e9c4665 100644
--- a/contrib/MeshOptimizer/MeshOpt.cpp
+++ b/contrib/MeshOptimizer/MeshOpt.cpp
@@ -36,6 +36,8 @@
#include "MeshOptCommon.h"
#include "MeshOpt.h"
+#include <unistd.h>
+
#if defined(HAVE_BFGS)
#include "ap.h"
@@ -47,6 +49,7 @@
namespace {
+
void evalObjGradFunc(const alglib::real_1d_array &x, double &Obj,
alglib::real_1d_array &gradObj, void *MOInst)
@@ -68,7 +71,7 @@ MeshOpt::MeshOpt(const std::map<MElement*,GEntity*> &element2entity,
const std::map<MElement*, GEntity*> &bndEl2Ent,
const std::set<MElement*> &els, std::set<MVertex*> &toFix,
const std::set<MElement*> &bndEls, const MeshOptParameters &par) :
- patch(element2entity, bndEl2Ent, els, toFix, bndEls, par.fixBndNodes), _verbose(0),
+ patch(element2entity, bndEl2Ent, els, toFix, bndEls, par.fixBndNodes), _verbose(0), _nCurses(1),
_iPass(0), _iter(0), _intervDisplay(0), _initObj(0)
{
_allObjFunc.resize(par.pass.size());
@@ -78,6 +81,33 @@ MeshOpt::MeshOpt(const std::map<MElement*,GEntity*> &element2entity,
_allObjFunc[iPass][iC] = par.pass[iPass].contrib[iC]->copy();
}
_objFunc = &_allObjFunc[0];
+ if (par.nCurses){
+ int nbRow, nbCol;
+ mvgetScreenSize(nbRow, nbCol);
+ int minCol = 0;
+ for (int i=0; i < _objFunc->names().size(); i++){
+ if (i > 0)
+ minCol += 5;
+ minCol += 34;
+ minCol += _objFunc->names()[i].size();
+ }
+ minCol = std::max(minCol, 96);
+ int minRow = std::max(34, 34+(int)_objFunc->names().size()-5);
+ if (nbRow < minRow || nbCol < minCol){
+ for (int i=0; i<nbRow; i++){
+ mvfillRow(i,' ');
+ }
+ int firstRow = (nbRow-5)/2;
+ mvprintCenter(firstRow,"Given the number of contributions to the objective function, ncurses");
+ mvprintCenter(firstRow+1,"enhanced interface requires a terminal window of minimal size");
+ mvprintCenter(firstRow+2,"%4ix%4i characters. Yours is %4ix%4i, try resizing the window", minRow, minCol, nbRow, nbCol);
+ mvprintCenter(firstRow+3,"or disabling the ncurses interface.");
+ mvpause();
+ mvterminate();
+ Msg::SetCallback(NULL);
+ Msg::Exit(EXIT_FAILURE);
+ }
+ }
}
@@ -85,6 +115,14 @@ MeshOpt::~MeshOpt()
{
for (int iPass=0; iPass<_allObjFunc.size(); iPass++)
for (int iC=0; iC<_allObjFunc[iPass].size(); iC++) delete _allObjFunc[iPass][iC];
+ while (_optHistory.size() > 0){
+ delete[] _optHistory.back();
+ _optHistory.pop_back();
+ }
+ while (_iterHistory.size() > 0){
+ delete[] _iterHistory.back();
+ _iterHistory.pop_back();
+ }
}
@@ -101,10 +139,26 @@ void MeshOpt::evalObjGrad(const alglib::real_1d_array &x, double &obj,
}
+ObjectiveFunction *MeshOpt::objFunction()
+{
+ return _objFunc;
+}
+
void MeshOpt::printProgress(const alglib::real_1d_array &x, double Obj)
{
_iter++;
+ if (_nCurses){
+ mvprintCenter(21, "Iteration %3d --- OBJ %12.5E (relative decrease = %12.5E)", _iter, Obj, Obj/_initObj);
+ if (_iterHistory.size() < 5){
+ _iterHistory.push_back(new char[1000]);
+ } else {
+ _iterHistory.push_back(_iterHistory.front());
+ _iterHistory.pop_front();
+ }
+ sprintf(_iterHistory.back(), _objFunc->minMaxStr().c_str());
+ mvprintList(22, 5, _iterHistory, 1);
+ }
if ((_verbose > 2) && (_iter % _intervDisplay == 0))
Msg::Info(("--> Iteration %3d --- OBJ %12.5E (relative decrease = %12.5E)" +
_objFunc->minMaxStr()).c_str(), _iter, Obj, Obj/_initObj);
@@ -132,7 +186,7 @@ void MeshOpt::updateResults()
void MeshOpt::runOptim(alglib::real_1d_array &x,
- const alglib::real_1d_array &initGradObj, int itMax)
+ const alglib::real_1d_array &initGradObj, int itMax, int iBar)
{
static const double EPSG = 0.;
static const double EPSF = 0.;
@@ -161,7 +215,34 @@ void MeshOpt::runOptim(alglib::real_1d_array &x,
iterationscount = rep.iterationscount;
nfev = rep.nfev;
terminationtype = rep.terminationtype;
-
+ if (_nCurses){
+ if (_optHistory.size() < 8){
+ _optHistory.push_front(new char[1000]);
+ } else {
+ _optHistory.push_front(_optHistory.back());
+ _optHistory.pop_back();
+ }
+ switch(int(terminationtype)) {
+ case 1: sprintf(_optHistory.front(),"Optimization run %3d (%3d iterations, %3d function evaluations): rel function improvement <= EpsF", iBar, iterationscount, nfev); break;
+ case 2: sprintf(_optHistory.front(),"Optimization run %3d (%3d iterations, %3d function evaluations): rel step <= EpsX ", iBar, iterationscount, nfev); break;
+ case 4: sprintf(_optHistory.front(),"Optimization run %3d (%3d iterations, %3d function evaluations): gradient norm <= EpsG ", iBar, iterationscount, nfev); break;
+ case 5: sprintf(_optHistory.front(),"Optimization run %3d (%3d iterations, %3d function evaluations): max number of steps taken ", iBar, iterationscount, nfev); break;
+ default: sprintf(_optHistory.front(),"Optimization run %3d (%3d iterations, %3d function evaluations): code %d ", iBar, iterationscount, nfev, int(terminationtype)); break;
+ }
+ if (_optHistory.size() < 8){
+ _optHistory.push_front(new char[1000]);
+ } else {
+ _optHistory.push_front(_optHistory.back());
+ _optHistory.pop_back();
+ }
+ sprintf(_optHistory.front(), _iterHistory.back());
+ mvprintList(19, -8, _optHistory, 2);
+ while (_iterHistory.size() > 0){
+ delete[] _iterHistory.back();
+ _iterHistory.pop_back();
+ }
+ mvprintList(22, 5, _iterHistory, 1);
+ }
if (_verbose > 2) {
Msg::Info("Optimization finalized after %d iterations (%d function evaluations),",
iterationscount, nfev);
@@ -180,16 +261,19 @@ int MeshOpt::optimize(const MeshOptParameters &par)
{
_intervDisplay = par.displayInterv;
_verbose = par.verbose;
-
+ _nCurses = par.nCurses;
// Set initial guess & result
int result = 1;
alglib::real_1d_array x;
x.setlength(patch.nPC());
patch.getUvw(x.getcontent());
-
+ if (_nCurses){
+ mvprintCenter(11, "%7i elements, %7i vertices, %7i free vertices, %7i variables",
+ patch.nEl(), patch.nVert(), patch.nFV(), patch.nPC());
+ }
if (_verbose > 2)
- Msg::Info("Patch composed of %i elements, %i vertices, %i free vertices, %i variables",
- patch.nEl(), patch.nVert(), patch.nFV(), patch.nPC());
+ Msg::Info("Patch composed of %i elements, %i vertices, %i free vertices, %i variables",
+ patch.nEl(), patch.nVert(), patch.nFV(), patch.nPC());
// Loop on passes
for (_iPass=0; _iPass<par.pass.size(); _iPass++) {
@@ -215,9 +299,21 @@ int MeshOpt::optimize(const MeshOptParameters &par)
// Loop for update of objective function parameters (barrier movement)
bool targetReached = _objFunc->targetReached();
for (int iBar=0; (iBar<par.pass[_iPass].maxParamUpdates) && (!targetReached); iBar++) {
+ if (_nCurses){
+ mvbold(true);
+ mvprintCenter(20, " OPTIMIZATION RUN %3d ", iBar);
+ mvbold(false);
+ if (iBar == 0){
+ while (_optHistory.size() > 0){
+ delete[] _optHistory.back();
+ _optHistory.pop_back();
+ }
+ }
+ mvprintList(19, -8, _optHistory, 2);
+ }
if (_verbose > 2) Msg::Info("--- Optimization run %d", iBar);
_objFunc->updateParameters();
- runOptim(x, gradObj, par.pass[_iPass].maxOptIter);
+ runOptim(x, gradObj, par.pass[_iPass].maxOptIter, iBar);
_objFunc->updateMinMax();
targetReached = _objFunc->targetReached();
if (_objFunc->stagnated()) {
diff --git a/contrib/MeshOptimizer/MeshOpt.h b/contrib/MeshOptimizer/MeshOpt.h
index 9c7a21c808..d132bb9410 100644
--- a/contrib/MeshOptimizer/MeshOpt.h
+++ b/contrib/MeshOptimizer/MeshOpt.h
@@ -30,6 +30,8 @@
#ifndef _MESHOPT_H_
#define _MESHOPT_H_
+#include <iostream>
+#include <fstream>
#include <string>
#include <math.h>
#include "MeshOptObjectiveFunction.h"
@@ -58,9 +60,11 @@ public:
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);
-
+ ObjectiveFunction *objFunction();
private:
int _verbose;
+ bool _nCurses;
+ std::list<char*> _iterHistory, _optHistory;
int _iPass;
std::vector<ObjectiveFunction> _allObjFunc; // Contributions to objective function for current pass
ObjectiveFunction *_objFunc; // Contributions to objective function for current pass
@@ -68,7 +72,7 @@ public:
double _initObj; // Values for reporting
void calcScale(alglib::real_1d_array &scale);
void runOptim(alglib::real_1d_array &x,
- const alglib::real_1d_array &initGradObj, int itMax);
+ const alglib::real_1d_array &initGradObj, int itMax, int iBar);
};
diff --git a/contrib/MeshOptimizer/MeshOptCommon.cpp b/contrib/MeshOptimizer/MeshOptCommon.cpp
index 4f82cba30d..62ff7b0177 100644
--- a/contrib/MeshOptimizer/MeshOptCommon.cpp
+++ b/contrib/MeshOptimizer/MeshOptCommon.cpp
@@ -1,4 +1,4 @@
-// Copyright (C) 2013 ULg-UCL
+ // Copyright (C) 2013 ULg-UCL
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@@ -27,12 +27,223 @@
//
// Contributors: Thomas Toulorge, Jonathan Lambrechts
+#include <iostream>
+#include <fstream>
+#include <list>
+#include <string.h>
#include "MElement.h"
#include "MeshOptCommon.h"
+#ifdef HAVE_NCURSES
-namespace {
+#include "ncurses.h"
+void catchPause(){
+ timeout(0);
+ if (getch() == ' '){
+ mvpause();
+ }
+};
+
+void mvinit(){
+ initscr();
+ start_color();
+ noecho();
+ curs_set(FALSE);
+ timeout(0);
+ init_pair(0, COLOR_WHITE, COLOR_BLACK);
+ init_pair(1, COLOR_BLACK, COLOR_WHITE);
+ init_pair(2, COLOR_BLACK, COLOR_YELLOW);
+ init_pair(3, COLOR_RED, COLOR_WHITE);
+ init_pair(4, COLOR_GREEN, COLOR_WHITE);
+ init_pair(5, COLOR_BLUE, COLOR_WHITE);
+ init_pair(6, COLOR_RED, COLOR_BLACK);
+ init_pair(7, COLOR_GREEN, COLOR_BLACK);
+ init_pair(8, COLOR_BLUE, COLOR_BLACK);
+}
+void mvterminate(){
+ endwin();
+}
+
+void mvgetScreenSize(int &nbRow, int &nbCol){
+ getmaxyx(stdscr,nbRow,nbCol);
+}
+
+void mvbold(bool on){
+ if (on)
+ attron(A_BOLD);
+ else
+ attroff(A_BOLD);
+}
+
+void mvcolor(int colorScheme, bool on){
+ if (on)
+ attron(COLOR_PAIR(colorScheme));
+ else
+ attroff(COLOR_PAIR(colorScheme));
+}
+void mvpause(){
+ attron(COLOR_PAIR(1));
+ attron(A_BOLD);
+ mvprintCenter(-1, " PAUSED (PRESS SPACE TO CONTINUE) ");
+ attroff(COLOR_PAIR(1));
+ attroff(A_BOLD);
+ timeout(-1);
+ while (getch() != ' '){}
+ mvfillRow(-1);
+}
+
+void mvprintCenter(int row, const char* fmt, ...){
+ catchPause();
+ char str[1000];
+ va_list args;
+ va_start(args, fmt);
+ vsnprintf(str, sizeof(str), fmt, args);
+ va_end(args);
+ int nbRow, nbCol;
+ getmaxyx(stdscr,nbRow,nbCol);
+ if (row<0)
+ row = nbRow+row;
+ mvprintw(row, (nbCol-strlen(str))/2, str, args);
+ refresh();
+}
+
+void mvprintLeft(int row, const char* fmt, ...){
+ catchPause();
+ char str[1000];
+ va_list args;
+ va_start(args, fmt);
+ vsnprintf(str, sizeof(str), fmt, args);
+ va_end(args);
+ int nbRow, nbCol;
+ getmaxyx(stdscr,nbRow,nbCol);
+ if (row<0)
+ row = nbRow+row;
+ mvprintw(row, 0, str, args);
+ refresh();
+}
+
+void mvprintRight(int row, const char* fmt, ...){
+ catchPause();
+ char str[1000];
+ va_list args;
+ va_start(args, fmt);
+ vsnprintf(str, sizeof(str), fmt, args);
+ va_end(args);
+ int nbRow, nbCol;
+ getmaxyx(stdscr,nbRow,nbCol);
+ if (row<0)
+ row = nbRow+row;
+ mvprintw(row, nbCol-strlen(str), str, args);
+ refresh();
+}
+
+void mvprintXY(int row, int col, const char* fmt, ...){
+ catchPause();
+ char str[1000];
+ va_list args;
+ va_start(args, fmt);
+ vsnprintf(str, sizeof(str), fmt, args);
+ va_end(args);
+ int nbRow, nbCol;
+ getmaxyx(stdscr,nbRow,nbCol);
+ if (row<0)
+ row = nbRow+row;
+ if (col<0)
+ col = nbCol+col;
+ mvprintw(row, col, str, args);
+ refresh();
+}
+
+void mvprintList(int row, int maxSize, std::list<char*> listStr, int colorScheme){
+ int nbRow, nbCol;
+ getmaxyx(stdscr,nbRow,nbCol);
+ if (row<0)
+ row = nbRow+row;
+ int i=0;
+ for (std::list<char*>::iterator it=listStr.begin(); it != listStr.end(); it++){
+ if (i >= abs(maxSize)) break;
+ if (colorScheme==1){
+ if (*it == listStr.back())
+ attron(COLOR_PAIR(2));
+ else
+ attron(COLOR_PAIR(1));
+ }
+ if (colorScheme==2){
+ if (i%2==0)
+ attron(COLOR_PAIR(1));
+ }
+ mvprintLeft(row + maxSize/abs(maxSize)*i, *it);
+ if (colorScheme==1){
+ if (*it == listStr.back())
+ attroff(COLOR_PAIR(2));
+ else
+ attroff(COLOR_PAIR(1));
+ }
+ if (colorScheme==2){
+ if (i%2==0)
+ attroff(COLOR_PAIR(1));
+ }
+ i++;
+ }
+ while (i < abs(maxSize)) {
+ mvfillRow(row + maxSize/abs(maxSize)*i++);
+ }
+}
+
+void mvfillRow(int row, char fillWith){
+ int nbRow, nbCol;
+ getmaxyx(stdscr,nbRow,nbCol);
+ if (row<0)
+ row = nbRow+row;
+ char toFill[1] = {fillWith};
+ for (int k = 0; k < nbCol; k++)
+ mvprintXY(row, k, toFill);
+}
+
+#else
+
+void catchPause(){Msg::Fatal("Gmsh should be configured with ncurses to use enhanced interface");}
+void mvinit(){Msg::Fatal("Gmsh should be configured with ncurses to use enhanced interface");}
+void mvterminate(){Msg::Fatal("Gmsh should be configured with ncurses to use enhanced interface");}
+void mvgetScreenSize(int &nbRow, int &nbCol){Msg::Fatal("Gmsh should be configured with ncurses to use enhanced interface");}
+void mvbold(bool on){Msg::Fatal("Gmsh should be configured with ncurses to use enhanced interface");}
+void mvcolor(int colorScheme, bool on){Msg::Fatal("Gmsh should be configured with ncurses to use enhanced interface");}
+void mvpause(){Msg::Fatal("Gmsh should be configured with ncurses to use enhanced interface");}
+void mvprintCenter(int row, const char* fmt, ...){Msg::Fatal("Gmsh should be configured with ncurses to use enhanced interface");}
+void mvprintLeft(int row, const char* fmt, ...){Msg::Fatal("Gmsh should be configured with ncurses to use enhanced interface");}
+void mvprintRight(int row, const char* fmt, ...){Msg::Fatal("Gmsh should be configured with ncurses to use enhanced interface");}
+void mvprintXY(int row, int col, const char* fmt, ...){Msg::Fatal("Gmsh should be configured with ncurses to use enhanced interface");}
+void mvprintList(int row, int maxSize, std::list<char*> listStr, int colorScheme){Msg::Fatal("Gmsh should be configured with ncurses to use enhanced interface");}
+void mvfillRow(int row, char fillWith){Msg::Fatal("Gmsh should be configured with ncurses to use enhanced interface");}
+
+#endif
+
+redirectMessage::redirectMessage(std::string logFileName, bool console){
+ std::ofstream logFile;
+ _logFileName = logFileName;
+ _console = console;
+ if(logFileName.compare("") != 0){
+ logFile.open(_logFileName.c_str());
+ logFile.close();
+ }
+}
+
+void redirectMessage::operator()(std::string level, std::string message){
+ std::ofstream logFile;
+ if(_logFileName.compare("") != 0){
+ logFile.open(_logFileName.c_str(), std::ios::app);
+ logFile << level << " : "<< message << std::endl;
+ logFile.close();
+ }
+ if (_console){
+ fprintf(stdout, "%s : %s\n", level.c_str(), message.c_str());
+ fflush(stdout);
+ }
+}
+
+
+namespace {
// Test intersection between sphere and segment
bool testSegSphereIntersect(SPoint3 A, SPoint3 B, const SPoint3& P, const double rr)
diff --git a/contrib/MeshOptimizer/MeshOptCommon.h b/contrib/MeshOptimizer/MeshOptCommon.h
index eddbcc042f..232fd359c2 100644
--- a/contrib/MeshOptimizer/MeshOptCommon.h
+++ b/contrib/MeshOptimizer/MeshOptCommon.h
@@ -30,14 +30,42 @@
#ifndef _MESHOPTCOMMON_H_
#define _MESHOPTCOMMON_H_
-#include <vector>
+#include <vector>
+#include "GmshMessage.h"
+#include "MeshOptimizerConfig.h"
class MElement;
class GEntity;
class SPoint3;
class ObjContrib;
+//ncurses shortcuts
+void mvinit();
+void mvterminate();
+void mvpause();
+void mvgetScreenSize(int &nbRow, int &nbCol);
+void mvprintCenter(int row, const char* fmt, ...);
+void mvprintLeft(int row, const char* fmt, ...);
+void mvprintRight(int row, const char* fmt, ...);
+void mvprintXY(int row, int col, const char* fmt, ...);
+//color scheme: 0=default, 1=last in yellow back, others in white back, 2=even numbers in white back
+void mvprintList(int row, int maxSize, std::list<char*> listStr, int colorScheme=0);
+void mvfillRow(int row, char fillWith=' ');
+void mvbold(bool on);
+void mvcolor(int colorScheme, bool on);
+
+
+
+class redirectMessage : public GmshMessage
+{
+ private:
+ std::string _logFileName;
+ bool _console;
+ public:
+ virtual void operator()(std::string level, std::string message);
+ redirectMessage(std::string logFileName , bool console);
+};
class MeshOptPatchDef {
public:
@@ -85,6 +113,8 @@ struct MeshOptParameters { // Parameters controllin
std::vector<MeshOptPass> pass;
int displayInterv; // Sampling rate in opt. iterations for display
int verbose; // Level of information displayed and written to disk
+ int nCurses; // Enhanced text output (not affected by verbose)
+ std::string logFileName; // External log file (affected by verbose)
int success; // Success flag: -1 = fail, 0 = partial fail (target not reached), 1 = success
double CPU; // Time for optimization
};
diff --git a/contrib/MeshOptimizer/MeshOptObjectiveFunction.cpp b/contrib/MeshOptimizer/MeshOptObjectiveFunction.cpp
index 45e339924f..a8dc3acd00 100644
--- a/contrib/MeshOptimizer/MeshOptObjectiveFunction.cpp
+++ b/contrib/MeshOptimizer/MeshOptObjectiveFunction.cpp
@@ -1,6 +1,8 @@
// TODO: Copyright
#include <sstream>
+#include <iostream>
+#include <iomanip>
#include "MeshOptObjContrib.h"
#include "MeshOptObjectiveFunction.h"
@@ -26,13 +28,31 @@ std::string ObjectiveFunction::minMaxStr()
std::string str;
for (std::vector<ObjContrib*>::iterator it=begin(); it!=end(); it++) {
std::ostringstream oss;
- oss << " -- Min. " + (*it)->getMeasureName() + " = " << (*it)->getMin();
- oss << " -- Max. " + (*it)->getMeasureName() + " = " << (*it)->getMax();
+ if (it != begin())
+ oss << " | ";
+ oss << std::scientific << std::setw(13) << (*it)->getMin() << " <= " << (*it)->getMeasureName() << " <= " << std::setw(13) << (*it)->getMax();
str += oss.str();
}
return str;
}
+std::vector<std::pair<double,double> > ObjectiveFunction::minMax() {
+ std::vector<std::pair<double,double> > range;
+ for (std::vector<ObjContrib*>::iterator it=begin(); it!=end(); it++) {
+ std::pair<double,double> oneRange = std::make_pair((*it)->getMin(),(*it)->getMax());
+ range.push_back(oneRange);
+ }
+ return range;
+}
+
+std::vector<std::string> ObjectiveFunction::names(){
+ std::vector<std::string > namesStr;
+ for (std::vector<ObjContrib*>::iterator it=begin(); it!=end(); it++) {
+ std::string name = (*it)->getMeasureName();
+ namesStr.push_back(name);
+ }
+ return namesStr;
+}
void ObjectiveFunction::updateMinMax()
{
diff --git a/contrib/MeshOptimizer/MeshOptObjectiveFunction.h b/contrib/MeshOptimizer/MeshOptObjectiveFunction.h
index 29d87ec3a7..f44c9ab553 100644
--- a/contrib/MeshOptimizer/MeshOptObjectiveFunction.h
+++ b/contrib/MeshOptimizer/MeshOptObjectiveFunction.h
@@ -17,6 +17,8 @@ public:
void initialize(Patch *mesh);
std::string contribNames();
std::string minMaxStr();
+ std::vector< std::pair<double,double> > minMax();
+ std::vector<std::string> names();
void updateMinMax();
void updateParameters();
void updateResults();
diff --git a/contrib/MeshOptimizer/MeshOptimizer.cpp b/contrib/MeshOptimizer/MeshOptimizer.cpp
index b8445e941e..924081d243 100644
--- a/contrib/MeshOptimizer/MeshOptimizer.cpp
+++ b/contrib/MeshOptimizer/MeshOptimizer.cpp
@@ -30,8 +30,10 @@
#include <stdio.h>
#include <sstream>
#include <iterator>
+#include <limits>
#include <string.h>
#include <stack>
+#include "Context.h"
#include "GmshConfig.h"
#include "Gmsh.h"
#include "GModel.h"
@@ -274,6 +276,9 @@ std::vector<elSetVertSetPair> getDisjointPatches(const vertElVecMap &vertex2elem
const elSet &badElements,
const MeshOptParameters &par)
{
+ if (par.nCurses) {
+ mvprintCenter(16, "Starting patch generation from %5i bad elements...", badElements.size());
+ }
Msg::Info("Starting patch generation from %i bad elements...", badElements.size());
elElSetMap element2elements; // Element to element connectivity, built progressively
@@ -358,6 +363,63 @@ void getAdjacentBndElts(const elElMap &el2BndEl, const elEntMap &bndEl2Ent,
}
}
+void displayMinMaxVal(int nbPatchSuccess[3], std::vector<std::string > &objFunctionNames, std::vector<std::pair<double,double> > newObjFunctionRange){
+ for (int i = 0; i < objFunctionNames.size(); i++){
+ if (nbPatchSuccess[0] > 0)
+ mvcolor(6, true);
+ else if (nbPatchSuccess[1] > 0)
+ mvcolor(8, true);
+ else
+ mvcolor(7, true);
+ mvprintRight(28+i, "%s on optimized patches: min %+.6e max %+.6e",objFunctionNames[i].c_str(), newObjFunctionRange[i].first, newObjFunctionRange[i].second);
+ if (nbPatchSuccess[0] > 0)
+ mvcolor(6, false);
+ else if (nbPatchSuccess[1] > 0)
+ mvcolor(8, false);
+ else
+ mvcolor(7, false);
+ }
+}
+
+void displayResultTable(int nbPatchSuccess[3], int nbPatch){
+ mvcolor(1, true);
+ mvprintLeft(28," TOTAL PATCHES : %4i ",nbPatch);
+ mvcolor(1, false);
+ mvcolor(3, true);
+ mvprintLeft(29," Fail : %4i ",nbPatchSuccess[0]);
+ mvcolor(3, false);
+ mvcolor(5, true);
+ mvprintLeft(30," Partial fail : %4i ",nbPatchSuccess[1]);
+ mvcolor(5, false);
+ mvcolor(4, true);
+ mvprintLeft(31," Success : %4i ",nbPatchSuccess[2]);
+ mvcolor(4, false);
+ mvcolor(1, true);
+ mvprintLeft(32," REMAINING : %4i ",nbPatch-nbPatchSuccess[0]-nbPatchSuccess[1]-nbPatchSuccess[2]);
+ mvcolor(1, false);
+}
+
+void updateDisplayPatchHistory(std::list<char*> &_patchHistory, const std::string &objFunctionStr, int iPatch, int iAdapt){
+ if (_patchHistory.size() < 8){
+ _patchHistory.push_front(new char[1000]);
+ } else {
+ _patchHistory.push_front(_patchHistory.back());
+ _patchHistory.pop_back();
+ }
+ if (iAdapt < 0)
+ sprintf(_patchHistory.front(),"Patch %i", iPatch);
+ else
+ sprintf(_patchHistory.front(),"Patch %i - Adaptation step %i", iPatch, iAdapt);
+ if (_patchHistory.size() < 8){
+ _patchHistory.push_front(new char[1000]);
+ } else {
+ _patchHistory.push_front(_patchHistory.back());
+ _patchHistory.pop_back();
+ }
+ sprintf(_patchHistory.front(), objFunctionStr.c_str());
+ mvprintList(9, -8, _patchHistory, 2);
+}
+
void optimizeDisjointPatches(const vertElVecMap &vertex2elements,
const elEntMap &element2entity,
@@ -365,6 +427,12 @@ void optimizeDisjointPatches(const vertElVecMap &vertex2elements,
const elEntMap &bndEl2Ent,
elSet &badasses, MeshOptParameters &par)
{
+ int nbPatchSuccess[3] = {0, 0, 0}; //0: fail, 1: partial fail, 2: success
+ std::list<char*> _patchHistory;
+ std::vector<std::pair<double,double> > newObjFunctionRange;
+ std::vector<std::string > objFunctionNames;
+
+
par.success = 1;
const elEntMap &e2ePatch = par.useGeomForPatches ? element2entity : elEntMap();
@@ -381,10 +449,16 @@ void optimizeDisjointPatches(const vertElVecMap &vertex2elements,
for (int iPatch = 0; iPatch < toOptimize.size(); ++iPatch)
getAdjacentBndElts(el2BndEl, bndEl2Ent, toOptimize[iPatch].first, bndElts[iPatch], par);
}
-
+ if (par.nCurses)
+ displayResultTable(nbPatchSuccess, toOptimize.size());
for (int iPatch = 0; iPatch < toOptimize.size(); ++iPatch) {
// Initialize optimization and output if asked
+ if (par.nCurses){
+ mvbold(true);
+ mvprintCenter(10, " PATCH %5i ", iPatch);
+ mvbold(false);
+ }
if (par.verbose > 1)
Msg::Info("Optimizing patch %i/%i composed of %i elements, "
"%i boundary elements", iPatch, toOptimize.size()-1,
@@ -412,10 +486,30 @@ void optimizeDisjointPatches(const vertElVecMap &vertex2elements,
// Evaluate mesh and update it if (partial) success
opt.updateResults();
+ if (newObjFunctionRange.size() == 0){
+ newObjFunctionRange = opt.objFunction()->minMax();
+ objFunctionNames = opt.objFunction()->names();
+ } else {
+ for (int i = 0; i < newObjFunctionRange.size(); i++){
+ newObjFunctionRange[i].first = std::min(newObjFunctionRange[i].first, opt.objFunction()->minMax()[i].first);
+ newObjFunctionRange[i].second = std::max(newObjFunctionRange[i].second, opt.objFunction()->minMax()[i].second);
+ }
+ }
if (success >= 0) opt.patch.updateGEntityPositions();
//#pragma omp critical
par.success = std::min(par.success, success);
+
+ nbPatchSuccess[success+1]++;
+ if (par.nCurses){
+ displayMinMaxVal(nbPatchSuccess, objFunctionNames, newObjFunctionRange);
+ displayResultTable(nbPatchSuccess, toOptimize.size());
+ updateDisplayPatchHistory(_patchHistory, opt.objFunction()->minMaxStr(), iPatch, -1);
+ }
+ }
+ while (_patchHistory.size() > 0){
+ delete[] _patchHistory.back();
+ _patchHistory.pop_back();
}
}
@@ -443,28 +537,36 @@ MElement *getWorstElement(elSet &badElts,
return worstEl;
}
-
+
void optimizeOneByOne(const vertElVecMap &vertex2elements,
const elEntMap &element2entity,
const elElMap &el2BndEl,
const elEntMap &bndEl2Ent,
elSet badElts, MeshOptParameters &par)
{
+ int nbPatchSuccess[3] = {0, 0, 0}; //0: fail, 1: partial fail, 2: success
+ std::list<char*> _patchHistory;
+ std::vector<std::pair<double,double> > newObjFunctionRange;
+ std::vector<std::string > objFunctionNames;
+
par.success = 1;
const elEntMap &e2ePatch = par.useGeomForPatches ? element2entity : elEntMap();
const elEntMap &e2eOpt = par.useGeomForOpt ? element2entity : elEntMap();
const int initNumBadElts = badElts.size();
+ if (par.nCurses)
+ displayResultTable(nbPatchSuccess, initNumBadElts);
if (par.verbose > 0) Msg::Info("%d bad elements, starting to iterate...", initNumBadElts);
elElSetMap element2elements; // Element to element connectivity, built progressively
-
+
// Loop over bad elements
for (int iBadEl=0; iBadEl<initNumBadElts; iBadEl++) {
+ int success;
if (badElts.empty()) break;
-
+
// Create patch around worst element and remove it from badElts
MElement *worstEl = getWorstElement(badElts, e2ePatch, par);
badElts.erase(worstEl);
@@ -472,8 +574,7 @@ void optimizeOneByOne(const vertElVecMap &vertex2elements,
// Initialize patch size to be adapted
int maxLayers = par.patchDef->maxLayers;
double distanceFactor = 1.;
- int success;
-
+
// Patch adaptation loop
for (int iAdapt=0; iAdapt<par.patchDef->maxPatchAdapt; iAdapt++) {
@@ -495,6 +596,11 @@ void optimizeOneByOne(const vertElVecMap &vertex2elements,
}
// Initialize optimization and output if asked
+ if (par.nCurses){
+ mvbold(true);
+ mvprintCenter(10, " PATCH %5i - ADAPTATION STEP %i ", iBadEl, iAdapt);
+ mvbold(false);
+ }
if (par.verbose > 1)
Msg::Info("Optimizing patch %i (max. %i remaining) composed of %4d elements",
iBadEl, badElts.size(), toOptimize.size());
@@ -520,9 +626,21 @@ void optimizeOneByOne(const vertElVecMap &vertex2elements,
opt.patch.writeMSH(ossI2.str().c_str());
}
+ if (par.nCurses)
+ updateDisplayPatchHistory(_patchHistory, opt.objFunction()->minMaxStr(), iBadEl, iAdapt);
+
// If (partial) success, update mesh and break adaptation loop, otherwise adapt
if ((success > 0) || (iAdapt == par.patchDef->maxPatchAdapt-1)) {
opt.updateResults();
+ if (newObjFunctionRange.size() == 0){
+ newObjFunctionRange = opt.objFunction()->minMax();
+ objFunctionNames = opt.objFunction()->names();
+ } else {
+ for (int i = 0; i < newObjFunctionRange.size(); i++){
+ newObjFunctionRange[i].first = std::min(newObjFunctionRange[i].first, opt.objFunction()->minMax()[i].first);
+ newObjFunctionRange[i].second = std::max(newObjFunctionRange[i].second, opt.objFunction()->minMax()[i].second);
+ }
+ }
if (success >= 0) {
opt.patch.updateGEntityPositions();
break;
@@ -535,20 +653,31 @@ void optimizeOneByOne(const vertElVecMap &vertex2elements,
}
}
+
} // End of adaptation loop
-
+
+ nbPatchSuccess[success+1]++;
+ if (par.nCurses){
+ displayMinMaxVal(nbPatchSuccess, objFunctionNames, newObjFunctionRange);
+ displayResultTable(nbPatchSuccess, initNumBadElts);
+ }
if (par.verbose > 1)
switch (success) {
- case 1: Msg::Info("Patch %i succeeded", iBadEl); break;
- case 0:
- Msg::Info("Patch %i partially failed (measure "
- "above critical value but below target)", iBadEl);
- break;
- case -1: Msg::Info("Patch %i failed", iBadEl); break;
+ case 1: Msg::Info("Patch %i succeeded", iBadEl); break;
+ case 0:
+ Msg::Info("Patch %i partially failed (measure "
+ "above critical value but below target)", iBadEl);
+ break;
+ case -1: Msg::Info("Patch %i failed", iBadEl); break;
}
-
+
par.success = std::min(par.success, success);
}
+
+ while (_patchHistory.size() > 0){
+ delete[] _patchHistory.back();
+ _patchHistory.pop_back();
+ }
}
@@ -561,8 +690,23 @@ void optimizeOneByOne(const vertElVecMap &vertex2elements,
void meshOptimizer(GModel *gm, MeshOptParameters &par)
{
#if defined(HAVE_BFGS)
+ if (par.nCurses)
+ mvinit();
+ redirectMessage _logWriter(par.logFileName, !par.nCurses);
+ if (par.logFileName.compare("") != 0 || par.nCurses)
+ Msg::SetCallback(&_logWriter);
+
double startTime = Cpu();
+ if (par.nCurses) {
+ mvbold(true);
+ mvprintCenter(0, "OPTIMIZING MESH");
+ mvfillRow(1,'-');
+ mvfillRow(10,'-');
+ mvfillRow(20,'-');
+ mvfillRow(27,'-');
+ mvbold(false);
+ }
if (par.verbose > 0) Msg::StatusBar(true, "Optimizing mesh...");
std::vector<GEntity*> entities;
@@ -576,6 +720,9 @@ void meshOptimizer(GModel *gm, MeshOptParameters &par)
GEntity* &entity = entities[iEnt];
if (entity->dim() != par.dim ||
(par.onlyVisible && !entity->getVisibility())) continue;
+ if (par.nCurses) {
+ mvprintCenter(15, "Computing connectivity and bad elements for entity %3d...", entity->tag());
+ }
Msg::Info("Computing connectivity and bad elements for entity %d...",
entity->tag());
calcVertex2Elements(par.dim, entity, vertex2elements);
@@ -603,21 +750,63 @@ void meshOptimizer(GModel *gm, MeshOptParameters &par)
else if (par.patchDef->strategy == MeshOptPatchDef::STRAT_ONEBYONE)
optimizeOneByOne(vertex2elements, element2entity,
el2BndEl, bndEl2Ent, badElts, par);
- else
- Msg::Error("Unknown strategy %d for mesh optimization", par.patchDef->strategy);
+ else {
+ if (par.nCurses){
+ mvcolor(2,true);
+ mvbold(true);
+ mvprintCenter(-2, " ERROR: Unknown strategy %d for mesh optimization ", par.patchDef->strategy);
+ mvcolor(2,false);
+ mvbold(false);
+ }
+ else
+ Msg::Error("Unknown strategy %d for mesh optimization", par.patchDef->strategy);
+ }
if (par.verbose > 0) {
- if (par.success == 1)
- Msg::Info("Optimization succeeded");
- else if (par.success == 0)
- Msg::Warning("Optimization partially failed (all measures above critical "
- "value, but some below target)");
- else if (par.success == -1)
- Msg::Error("Optimization failed (some measures below critical value)");
+ if (par.success == 1){
+ if (par.nCurses){
+ mvcolor(4,true);
+ mvbold(true);
+ mvprintCenter(-2, " Optimization succeeded ");
+ mvcolor(4,false);
+ mvbold(false);
+ }
+ else
+ Msg::Info("Optimization succeeded");
+ }
+ else if (par.success == 0){
+ if (par.nCurses){
+ mvcolor(5,true);
+ mvbold(true);
+ mvprintCenter(-2, " Optimization partially failed (all measures above critical value, but some below target) ");
+ mvcolor(5,false);
+ mvbold(false);
+ }
+ else
+ Msg::Warning("Optimization partially failed (all measures above critical "
+ "value, but some below target)");
+ }
+ else if (par.success == -1){
+ if (par.nCurses){
+ mvcolor(3,true);
+ mvbold(true);
+ mvprintCenter(-2, "Optimization Failed");
+ mvcolor(3,false);
+ mvbold(false);
+ }
+ else
+ Msg::Error("Optimization failed (some measures below critical value)");
+ }
par.CPU = Cpu()-startTime;
Msg::StatusBar(true, "Done optimizing mesh (%g s)", par.CPU);
}
-
+ if (par.nCurses){
+ mvpause();
+ mvterminate();
+ }
+ if (par.logFileName.compare("") != 0 || par.nCurses)
+ Msg::SetCallback(NULL);
+
#else
Msg::Error("Mesh optimizer requires BFGS");
#endif
diff --git a/contrib/MeshOptimizer/MeshOptimizerConfig.h.in b/contrib/MeshOptimizer/MeshOptimizerConfig.h.in
new file mode 100644
index 0000000000..84172db05e
--- /dev/null
+++ b/contrib/MeshOptimizer/MeshOptimizerConfig.h.in
@@ -0,0 +1,6 @@
+#ifndef _MESH_OPTIMIZER_CONFIG_H_
+#define _MESH_OPTIMIZER_H_
+
+#cmakedefine HAVE_NCURSES
+
+#endif
diff --git a/contrib/MeshOptimizer/doc/mesh_optimization.tex b/contrib/MeshOptimizer/doc/mesh_optimization.tex
index 87b6038294..d813f2cf41 100644
--- a/contrib/MeshOptimizer/doc/mesh_optimization.tex
+++ b/contrib/MeshOptimizer/doc/mesh_optimization.tex
@@ -27,7 +27,7 @@
\usepackage{parskip}
%\usepackage{alltt}
\usepackage{verbatim}
-\usepackage[latin1]{inputenc}
+%\usepackage[latin1]{inputenc}
\usepackage{float}
\usepackage{url}
%\usepackage[francais]{babel}
--
GitLab