Skip to content
Snippets Groups Projects
Commit d992e38e authored by Sebastien Blaise's avatar Sebastien Blaise
Browse files

meshOptimize: fixes + added ncurses interface

parent 939906a5
No related branches found
No related tags found
No related merge requests found
......@@ -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)
......@@ -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()) {
......
......@@ -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);
};
......
// 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)
......
......@@ -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
};
......
// 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()
{
......
......@@ -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();
......
......@@ -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
......
#ifndef _MESH_OPTIMIZER_CONFIG_H_
#define _MESH_OPTIMIZER_H_
#cmakedefine HAVE_NCURSES
#endif
......@@ -27,7 +27,7 @@
\usepackage{parskip}
%\usepackage{alltt}
\usepackage{verbatim}
\usepackage[latin1]{inputenc}
%\usepackage[latin1]{inputenc}
\usepackage{float}
\usepackage{url}
%\usepackage[francais]{babel}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment