diff --git a/Mesh/meshRecombine2D.cpp b/Mesh/meshRecombine2D.cpp
index c28768c0b3d22998c42f669faeee1801fcb8dfda..36f2e7a8318efb9b5431ce442f873a3bf0469e8b 100644
--- a/Mesh/meshRecombine2D.cpp
+++ b/Mesh/meshRecombine2D.cpp
@@ -326,9 +326,9 @@ int Recombine2D::apply()
_gf->quadrangles = _quads;
_applied = true;
- _gf->model()->writeMSH("recSquare.msh");
- laplaceSmoothing(_gf,100);
- _gf->model()->writeMSH("aftSquare.msh");
+ //_gf->model()->writeMSH("recSquare.msh");
+ laplaceSmoothing(_gf,10);
+ //_gf->model()->writeMSH("aftSquare.msh");
return 1;
}
diff --git a/Mesh/meshRecombine2D.h b/Mesh/meshRecombine2D.h
index 5c624db8e584fd81aa8e4f3c8c1fd470b2bab031..eecfd57ca12b0b18859730aa7e3e86ea549f2411 100644
--- a/Mesh/meshRecombine2D.h
+++ b/Mesh/meshRecombine2D.h
@@ -22,16 +22,184 @@ class RecombTriangle;
class Recomb2D_Node;
class TrianglesUnion;
class Rec2d_vertex;
+class list_actions;
+class Recombine2D;
+//class list_actions;
+//class list_actions::iterator;
struct lessRecombTri {
bool operator()(RecombTriangle *rt1, RecombTriangle *rt2) const;
};
struct lessTriUnion {
bool operator()(TrianglesUnion *, TrianglesUnion *) const;
};
-
+struct lessTriUnionInv {
+ bool operator()(TrianglesUnion *, TrianglesUnion *) const;
+};
typedef std::set<RecombTriangle*,lessRecombTri> Set_Recomb;
typedef std::map<MElement*,std::set<RecombTriangle*> > Map_Tri_Recomb;
typedef std::map<MElement*,Recomb2D_Node*> Map_Tri_Node;
+typedef std::map<MTriangle*, std::set<TrianglesUnion*> > mapTriUnion;
+typedef std::set<TrianglesUnion*, lessTriUnionInv> setTriUnion;
+typedef std::vector<setTriUnion> vectSetTriUnion;
+
+class list_actions {
+ private :
+ vectSetTriUnion _cont;
+
+ public :
+ list_actions(int h) {_cont.reserve(h);}
+
+ class iterator {
+ private :
+ int pos;
+ list_actions *la;
+ setTriUnion::iterator it;
+
+ public :
+ void print() {
+ Msg::Info("iterator(pos %d, la %d, tu %d)", pos, la, *it);
+ }
+ iterator(list_actions *l, int i, setTriUnion::iterator t) {
+ pos = i;
+ la = l;
+ it = t;
+ }
+ iterator() {
+ pos = -1;
+ }
+ bool operator==(iterator a) const {
+ if (pos != a.pos)
+ return false;
+ if (pos < 0)
+ return true;
+ return it == a.it;
+ }
+ bool operator!=(iterator a) const {
+ if (pos != a.pos)
+ return true;
+ if (pos < 0)
+ return false;
+ return it != a.it;
+ }
+ iterator& operator=(const iterator &a) {
+ if (this != &a) {
+ pos = a.pos;
+ la = a.la;
+ it = a.it;
+ }
+ return *this;
+ }
+ TrianglesUnion* operator*() const {
+ //Msg::Info("pos %d", pos);
+ //Msg::Info(" ");
+ if (pos < 0)
+ return NULL;
+ return *it;
+ }
+ iterator& operator++() {
+ int p = pos;
+ ++it;
+ while (it == la->_end(p) && ++p != la->_size()) {
+ ++pos;
+ it = la->_begin(pos);
+ if (p != pos)
+ Msg::Error(" ");
+ }
+ return *this;
+ }
+ iterator operator++(int n) {
+ iterator t = *this;
+ int p = pos;
+ ++it;
+ while (it == la->_end(p) && ++p != la->_size()) {
+ ++pos;
+ it = la->_begin(pos);
+ }
+ return t;
+ }
+ inline int getPos() {return pos;}
+ };
+
+ int add(setTriUnion &s) {
+ setTriUnion::iterator it = s.begin();
+ while (it != s.end()) {
+ if (find(*it) != end())
+ s.erase(it++);
+ else
+ ++it;
+ }
+ _cont.push_back(s);
+ return _cont.size() - 1;
+ }
+ void remove(int a) {
+ for (; a < _cont.size(); a++)
+ _cont.pop_back();
+ }
+ iterator find(TrianglesUnion* tu) {
+ if (_cont.empty())
+ return end();
+ for (int i = 0; i < _cont.size(); ++i) {
+ setTriUnion::iterator it = _cont[i].find(tu);
+ if (it != _cont[i].end())
+ return iterator(this, i, it);
+ }
+ return end();
+ }
+ iterator begin() {
+ if (_cont.empty())
+ return iterator(this, -1, _begin());
+ return iterator(this, 0, _begin());
+ }
+ iterator end() {
+ return iterator(this, _cont.size() - 1, _end());
+ }
+ void pop_back() {
+ _cont.pop_back();
+ }
+
+ void sizes(){
+ switch (_cont.size()) {
+ case 3 :
+ Msg::Info("--Actions : %d + %d + %d", _cont[0].size(), _cont[1].size(), _cont[2].size());
+ break;
+ case 2 :
+ Msg::Info("--Actions : %d + %d", _cont[0].size(), _cont[1].size());
+ break;
+ case 1 :
+ Msg::Info("--Actions : %d", _cont[0].size());
+ break;
+ default :
+ break;
+ }
+ //Msg::Info("%d sets", _cont.size());
+ //for (int i = 0; i < _cont.size(); ++i)
+ // Msg::Info("[%d] %d", i+1, _cont[i].size());
+ }
+
+ private :
+ setTriUnion::iterator _end(int pos = -1) {
+ if (_cont.size() == 0) {
+ setTriUnion::iterator it;
+ return it; // Not clean !
+ }
+ if (pos < 0)
+ return _cont.back().end();
+ return _cont[pos].end();
+ }
+ setTriUnion::iterator _begin(int pos = -1) {
+ if (_cont.size() == 0) {
+ setTriUnion::iterator it;
+ return it; // Not clean !
+ }
+ if (pos < 0)
+ return _cont.front().begin();
+ return _cont[pos].begin();
+ }
+ inline int _size() {return _cont.size();}
+
+};
+
+
class Recombine2D {
private :
@@ -43,6 +211,7 @@ class Recombine2D {
Map_Tri_Recomb _possibleRecomb;
std::set<MElement*> _isolated;
std::vector<MQuadrangle*> _quads;
+ static int ra;
template <class E>
void _buildEdgeToElement(std::vector<E*> &, e2t_cont &);
@@ -57,15 +226,27 @@ class Recombine2D {
~Recombine2D();
int apply();
- int apply(bool);
inline double getBenef() const {return _benef;}
inline int numTriangle() const {return _isolated.size();}
private :
//std::set<TrianglesUnion*, lessTriUnion> _setQuads;
std::list<TrianglesUnion*> _setQuads;
- void _removeImpossible(std::list<TrianglesUnion*>::iterator);
+ std::list<TrianglesUnion*> _lastQuad;
+ mapTriUnion _possibleActions;
+ mapTriUnion _lessActions;
+ void _removeImpossible(TrianglesUnion*);
void _recombine(bool);
+
+ void _lookAhead(std::list<TrianglesUnion*>&, int horiz);
+ void _getIncompatibles(const TrianglesUnion*, setTriUnion&);
+ void _getNeighbors(const TrianglesUnion*, setTriUnion&);
+ void _addNeighbors(int horiz, std::vector<list_actions::iterator>&, int current, list_actions*);
+ void _removeNeighbors(int horiz, int current, list_actions*);
+ int _checkIsolatedTri(std::vector<list_actions::iterator>&, int size, std::set<MTriangle*>&);
+
+ public :
+ int apply(bool);
};
class RecombTriangle {
@@ -135,9 +316,6 @@ class Rec2d_vertex {
// GEntity : virtual Range<double> parBounds(int i)
static double **_Vvalues;
static double **_Vbenefs;
-
- void _initStaticTable();
- double _computeAngle(MVertex *, std::list<MTriangle*> &, std::set<GEdge*> &);
public :
Rec2d_vertex(MVertex *, std::list<MTriangle*> &,
@@ -146,6 +324,12 @@ class Rec2d_vertex {
inline void changeNumEl(int c) {_numEl += c;}
double getReward();
double getReward(int);
+ //int onwhat() {return _onWhat;}
+ //int numel() {return _numEl;}
+
+ private :
+ void _initStaticTable();
+ double _computeAngle(MVertex *, std::list<MTriangle*> &, std::set<GEdge*> &);
};
class TrianglesUnion {
@@ -155,24 +339,31 @@ class TrianglesUnion {
Rec2d_vertex **_vertices;
MTriangle **_triangles;
static int _RECOMPUTE; //incremented when a recombination is executed
+
public:
static int _NumEdge, _NumVert;
static double _ValEdge, _ValVert;
private:
- double _computeEdgeLength(GFace *, MVertex **, double *u, double *v,
- int numIntermedPoint= 1);
- double _computeAlignment(MEdge &, std::list<MTriangle*> &);
- void _update();
-
public :
TrianglesUnion(GFace *, std::list<MTriangle*> &, std::list<MEdge> &,
std::list<Rec2d_vertex*> &, std::map<MVertex*,double> &);
bool operator<(TrianglesUnion &other);
- inline double getEdgeValue() {return _embEdgeValue;}
- inline double getNumTriangles() {return _numTri;}
- inline MTriangle *getTriangle(int num) {return _triangles[num];}
+ void addTri(std::set<MTriangle*>&) const;
+ void removeTri(std::set<MTriangle*>&) const;
+ bool isIn(std::set<MTriangle*>&) const;
+ void addInfluence(int*, double*, std::map<Rec2d_vertex*, int>&) const;
+ void removeInfluence(int*, double*, std::map<Rec2d_vertex*, int>&) const;
+ static double computeReturn(const int*, const double*, const std::map<Rec2d_vertex*, int>&);
+ inline double getEdgeValue() const {return _embEdgeValue;}
+ inline double getVertValue() const {return _embVertValue;}
+ inline int getNumVerts() const {return _numEmbVert;}
+ inline int getNumEdges() const {return _numEmbEdge;}
+ inline int getNumVertices() const {return _numBoundVert;}
+ inline Rec2d_vertex* getVertex(int num) const {return _vertices[num];}
+ inline int getNumTriangles() const {return _numTri;}
+ inline MTriangle* getTriangle(int num) const {return _triangles[num];}
void select() {
_ValEdge -= _embEdgeValue;
_NumEdge -= _numEmbEdge;
@@ -184,7 +375,7 @@ class TrianglesUnion {
}
_RECOMPUTE++;
}
- MQuadrangle *createQuad() const;
+ MQuadrangle* createQuad() const;
void print() {
Msg::Info("Printing TU (%d,%d,%d,%d)", _numEmbVert, _numBoundVert, _numEmbEdge, _numTri);
for (int i = 0; i < _numTri; i++)
@@ -192,5 +383,11 @@ class TrianglesUnion {
}
static inline void addRec() {_RECOMPUTE++;}
inline double getReturn() {return _globValIfExecuted;}
+
+ private:
+ double _computeEdgeLength(GFace *, MVertex **, double *u, double *v,
+ int numIntermedPoint= 1);
+ double _computeAlignment(MEdge &, std::list<MTriangle*> &);
+ void _update();
};
#endif
diff --git a/Mesh/meshRecombine2D_2.cpp b/Mesh/meshRecombine2D_2.cpp
index 1543764f200ee781e8a56af7255e7ec5ff81d01c..1fcb2584beff77438ad65d1fc24eb5afd273008e 100644
--- a/Mesh/meshRecombine2D_2.cpp
+++ b/Mesh/meshRecombine2D_2.cpp
@@ -17,8 +17,8 @@
#include "Context.h"
#include "OpenFile.h"//pas propre
#include "Field.h"
-
-static int HORIZ = 20;
+#include "OS.h"
+#define HORIZ2 7
double **Rec2d_vertex::_Vvalues = NULL;
double **Rec2d_vertex::_Vbenefs = NULL;
@@ -27,6 +27,16 @@ int TrianglesUnion::_NumEdge = 0;
int TrianglesUnion::_NumVert = 0;
double TrianglesUnion::_ValEdge = .0;
double TrianglesUnion::_ValVert = .0;
+int Recombine2D::ra = 0;
+
+/*
+ a faire :
+ - modifier gain des edges
+ - fonction retournant qualite
+ - verifier la qualite
+ - action supplementaire (faire classe de base)
+ - Faire un super conteneur pour les actions
+*/
Recombine2D::Recombine2D(GFace *gf)
: _horizon(0), _gf(gf), _benef(.0), _applied(true)
@@ -69,10 +79,10 @@ Recombine2D::Recombine2D(GFace *gf)
std::map<MVertex*, Rec2d_vertex*>::iterator itV2N = mapV2N.begin();
std::list<std::map<MVertex*, std::list<MTriangle*> >::iterator> fourTri;
for (; itvert != mapVert.end(); itvert++) {
- if ((*itvert).second.size() == 4 && (*itvert).first->onWhat()->dim() == 2)
+ if (itvert->second.size() == 4 && itvert->first->onWhat()->dim() == 2)
fourTri.push_back(itvert);
- Rec2d_vertex *n = new Rec2d_vertex((*itvert).first, (*itvert).second, mapCornerVert);
- itV2N = mapV2N.insert(itV2N, std::make_pair((*itvert).first,n));
+ Rec2d_vertex *n = new Rec2d_vertex(itvert->first, itvert->second, mapCornerVert);
+ itV2N = mapV2N.insert(itV2N, std::make_pair(itvert->first,n));
TrianglesUnion::_NumVert++;
TrianglesUnion::_ValVert += n->getReward();
}
@@ -98,13 +108,13 @@ Recombine2D::Recombine2D(GFace *gf)
{
std::set<MVertex*> setVert;
std::multiset<MEdge, Less_Edge> msetEdge;
- std::list<MTriangle*>::iterator itTri = (*(*it4)).second.begin();
- for (; itTri != (*(*it4)).second.end(); itTri++) {
+ std::list<MTriangle*>::iterator itTri = (*it4)->second.begin();
+ for (; itTri != (*it4)->second.end(); itTri++) {
MVertex *vt;
for (int i = 0; i < 3; i++) {
msetEdge.insert((*itTri)->getEdge(i));
vt = (*itTri)->getVertex(i);
- if (vt != (*(*it4)).first)
+ if (vt != (*it4)->first)
setVert.insert(vt);
}
}
@@ -118,10 +128,14 @@ Recombine2D::Recombine2D(GFace *gf)
std::set<MVertex*>::iterator itsVert = setVert.begin();
for (; itsVert != setVert.end(); itsVert++)
listVert.push_back(mapV2N[*itsVert]);
- listVert.push_back(mapV2N[(*(*it4)).first]);
+ listVert.push_back(mapV2N[(*it4)->first]);
}
Msg::Info("%d",++j);
- _setQuads.push_back(new TrianglesUnion(gf, (*(*it4)).second, listEmbEdge, listVert, mapV2LcUV));
+ TrianglesUnion *tu = new TrianglesUnion(gf, (*it4)->second, listEmbEdge, listVert, mapV2LcUV);
+ std::list<MTriangle*>::iterator itTri = (*it4)->second.begin();
+ for (; itTri != (*it4)->second.end(); itTri++)
+ _possibleActions[*itTri].insert(tu);
+ _setQuads.push_back(tu);
}
// Create 'TrianglesUnion' for pattern of 2 triangles
@@ -134,22 +148,25 @@ Recombine2D::Recombine2D(GFace *gf)
*/
std::map<MEdge, std::list<MTriangle*> >::iterator itedge = mapEdge.begin();
for (; itedge != mapEdge.end(); itedge++) {
- if ((*itedge).second.size() == 2) {
+ if (itedge->second.size() == 2) {
std::list<MEdge> listEmbEdge;
- listEmbEdge.push_back((*itedge).first);
+ listEmbEdge.push_back(itedge->first);
std::list<Rec2d_vertex*> listVert;
- listVert.push_back(mapV2N[(*itedge).first.getVertex(0)]);
- listVert.push_back(mapV2N[(*itedge).first.getVertex(1)]);
- TrianglesUnion *tu = new TrianglesUnion(gf, (*itedge).second, listEmbEdge, listVert, mapV2LcUV);
+ listVert.push_back(mapV2N[itedge->first.getVertex(0)]);
+ listVert.push_back(mapV2N[itedge->first.getVertex(1)]);
+ TrianglesUnion *tu = new TrianglesUnion(gf, itedge->second, listEmbEdge, listVert, mapV2LcUV);
_setQuads.push_back(tu);
TrianglesUnion::_NumEdge++;
TrianglesUnion::_ValEdge += tu->getEdgeValue();
+ std::list<MTriangle*>::iterator itTri = itedge->second.begin();
+ for (; itTri != itedge->second.end(); itTri++)
+ _possibleActions[*itTri].insert(tu);
}
else
- ;//Msg::Info("[bord] %d", (*itedge).second.size());
+ ;//Msg::Info("[bord] %d", itedge->second.size());
}
TrianglesUnion::addRec();
- _setQuads.sort(lessTriUnion());
+ //_setQuads.sort(lessTriUnion());
_recombine(true);
_applied = false;
@@ -161,31 +178,324 @@ void Recombine2D::_recombine(bool a)
SetBoundingBox();
int i = 0;
- while (!_setQuads.empty() && i < 1000) {
- std::list<TrianglesUnion*>::iterator it = _setQuads.begin();
- Msg::Info("------------------ %d", _setQuads.size());
- {
- Msg::Info("Expected return %lf", (*it)->getReturn());
+ std::list<TrianglesUnion*> vectTu;
+
+ while (!_setQuads.empty() && i < 2000) {
+ TrianglesUnion *tu;
+ if (_lastQuad.empty()) {
+ if (vectTu.size() < 2)
+ tu = *std::max_element(_setQuads.begin(), _setQuads.end(), lessTriUnion());
+ else {
+ vectTu.pop_front();
+ //tu = *std::max_element(vectTu.begin(),vectTu.end(), lessTriUnion());
+ tu = *vectTu.begin();
+ }
+ vectTu.clear();
+ vectTu.insert(vectTu.begin(), tu);
+
+ Msg::Info("------------------ %d", _setQuads.size());
+ _lookAhead(vectTu, HORIZ2);
+ }
+ else {
+ tu = *_lastQuad.begin();
+ vectTu.clear();
+ vectTu.insert(vectTu.begin(), tu);
+ std::set<MTriangle*> touched;
+ for (int i = 0; i < tu->getNumTriangles(); i++) {
+ touched.insert(tu->getTriangle(i));
+ }
+ std::list<TrianglesUnion*>::iterator it = _lastQuad.begin();
+ while (it != _lastQuad.end()) {
+ bool toBeErased = false;
+ for (int i = 0; i < (*it)->getNumTriangles(); i++) {
+ if (touched.find((*it)->getTriangle(i)) != touched.end())
+ toBeErased = true;
+ }
+ if (toBeErased)
+ _lastQuad.erase(it++);
+ else
+ ++it;
+ }
+ //_lastQuad.erase(_lastQuad.begin());
+ }
+ tu = *vectTu.begin();
+ tu->select();
+ _quads.push_back(tu->createQuad());
+ _removeImpossible(tu);
+ { // draw Mesh
+ Msg::Info("Expected return %lf", tu->getReturn());
_applied = false;
- apply();
+ apply(true);
_applied = false;
CTX::instance()->mesh.changed = (ENT_ALL);
//drawContext::global()->draw();
- //FlGui::instance()->check();
+ FlGui::instance()->check();
drawContext::global()->draw();
- //if (!Msg::GetAnswer("Continue ?", 1, "no", "yes"))
- // Msg::Info("I continue anyway :p");
- }
- (*it)->select();
+ double t = Cpu();
+ while (Cpu() - t < .0001);
+ if (i % 1000 == 0)
+ if (!Msg::GetAnswer("Continue ?", 1, "no", "yes"))
+ Msg::Info("I continue anyway :p");
+ }
+
+ ++i;
+ }
+ Msg::Info("Done %d (%d)", i, _pairs.size());
+}
+
+void printIt(std::vector<list_actions::iterator> &vIt, int size)
+{
+ switch (size) {
+ case 3 :
+ Msg::Info("--Iterators : %d , %d , %d", vIt[0].getPos(), vIt[1].getPos(), vIt[2].getPos());
+ break;
+ case 2 :
+ Msg::Info("--Iterators : %d , %d", vIt[0].getPos(), vIt[1].getPos());
+ break;
+ case 1 :
+ Msg::Info("--Iterators : %d", vIt[0].getPos());
+ break;
+ default :
+ break;
+ }
+}
+
+void Recombine2D::_lookAhead(std::list<TrianglesUnion*> &candidate, int horiz)
+{
+ if (horiz < 1)
+ horiz = 1;
+ _lessActions.clear();
+
+ double maxReturn = .0;
+ //TrianglesUnion *bestAction;
+ int numbers[2] = {0, 0}; // number of embedded verts & edges
+ double values[2] = {.0, .0}; // embedded vert & edge values
+ std::set<MTriangle*> setTri;
+ std::map<Rec2d_vertex*, int> modifiedVert;
+ std::vector<list_actions::iterator> vecIt(horiz);
+ std::vector<std::pair<TrianglesUnion*, int> > bestSequence;
+
+ list_actions la(horiz);
+ {
+ setTriUnion candidates;
+ _getIncompatibles(*candidate.begin(), candidates);
+ la.add(candidates);
+ }
+
+ vecIt[0] = la.begin();
+ (*vecIt[0])->addTri(setTri);
+ (*vecIt[0])->addInfluence(numbers, values, modifiedVert);
+ _addNeighbors(horiz, vecIt, 0, &la);
+
+ int current = 1;
+ int num=0;
+
+ while (current > 0) {
+ bool best = false;
+ // set first acceptable action for each step
+ while (current < horiz && vecIt[current-1] != la.end()) {
+ vecIt[current] = vecIt[current-1];
+ while (++vecIt[current] != la.end() && (*vecIt[current])->isIn(setTri));
+ if (vecIt[current] != la.end()) {
+ (*vecIt[current])->addTri(setTri);
+ (*vecIt[current])->addInfluence(numbers, values, modifiedVert);
+ _addNeighbors(horiz, vecIt, current, &la);
+ }
+ ++current;
+ }
+
+ { // color
+ for (int i = 1; i < current && vecIt[i] != la.end(); ++i) {
+ TrianglesUnion *tu = *vecIt[i];
+ int pos = vecIt[i].getPos();
+ double frac = (double)pos / (double)(horiz-1);
+ unsigned int col;
+ //la.sizes();
+ //Msg::Info("%d / %d -> %d", pos, horiz, (int)(200.*frac));
+ col = CTX::instance()->packColor(50 - (int)(25.*frac), 200 - (int)(25.*frac), (int)(200.*frac), 255);
+ for (int j = 0; j < tu->getNumTriangles(); ++j) {
+ tu->getTriangle(j)->setCol(col);
+ }
+ }
+ }
+
+ double expectedReturn = TrianglesUnion::computeReturn(numbers, values, modifiedVert);
+ if (maxReturn < expectedReturn) {
+ int sizeSequence = current;
+ if (vecIt[sizeSequence - 1] == la.end())
+ --sizeSequence;
+ expectedReturn -= .1 * _checkIsolatedTri(vecIt, sizeSequence, setTri);
+ if (maxReturn < expectedReturn) {
+ maxReturn = expectedReturn;
+ //bestAction = *vecIt[0];
+ bestSequence.resize(sizeSequence);
+ for (int i = 0; i < sizeSequence; ++i)
+ bestSequence[i] = std::make_pair(*vecIt[i], vecIt[i].getPos());
+ { //color
+ best = true;
+ TrianglesUnion *tu = *vecIt[0];
+ unsigned int col;
+ col = CTX::instance()->packColor(190, 0, 190, 255);
+ for (int j = 0; j < tu->getNumTriangles(); ++j) {
+ tu->getTriangle(j)->setCol(col);
+ }
+ }
+ }
+ }
+ else { //color
+ TrianglesUnion *tu = *vecIt[0];
+ unsigned int col;
+ col = CTX::instance()->packColor(190, 0, 0, 255);
+ for (int j = 0; j < tu->getNumTriangles(); ++j) {
+ tu->getTriangle(j)->setCol(col);
+ }
+ }
- _quads.push_back((*it)->createQuad());
+ if (true || best) { // draw mesh
+ apply(false);
+ CTX::instance()->mesh.changed = (ENT_ALL);
+ FlGui::instance()->check();
+ drawContext::global()->draw();
+ double t = Cpu();
+ while (Cpu() - t < .0001);
+ //if (i % 1 == 0)
+ if (/*false && best &&*/ !Msg::GetAnswer("Continue ?", 1, "no", "yes"))
+ Msg::Info("I continue anyway :p");
+ best = false;
+ }
- _removeImpossible(it);
+ { // decolor
+ for (int i = 0; i < current && vecIt[i] != la.end(); ++i) {
+ TrianglesUnion *tu = *vecIt[i];
+ int pos = vecIt[i].getPos();
+ double frac = (double)pos / (double)horiz;
+ unsigned int col;
+ col = CTX::instance()->packColor(255, 255, 0, 255);
+ for (int j = 0; j < tu->getNumTriangles(); ++j) {
+ tu->getTriangle(j)->setCol(col);
+ }
+ }
+ }
- _setQuads.sort(lessTriUnion());
- i++;
+ if (vecIt[current-1] != la.end() || --current > 0) {
+ do {
+ (*vecIt[current-1])->removeTri(setTri);
+ (*vecIt[current-1])->removeInfluence(numbers, values, modifiedVert);
+ _removeNeighbors(horiz, current - 1, &la);
+ while (++vecIt[current-1] != la.end() && (*vecIt[current-1])->isIn(setTri));
+ } while (vecIt[current-1] == la.end() && --current > 0);
+ if (current > 0) {
+ (*vecIt[current-1])->addTri(setTri);
+ (*vecIt[current-1])->addInfluence(numbers, values, modifiedVert);
+ _addNeighbors(horiz, vecIt, current - 1, &la);
+ }
+ }
}
- Msg::Info("Done %d (%d)", i, _pairs.size());
+
+ { // color
+ for (int i = 1; i < bestSequence.size(); ++i) {
+ TrianglesUnion *tu = bestSequence[i].first;
+ unsigned int col = CTX::instance()->packColor(50, 200, 0, 255);
+ for (int j = 0; j < tu->getNumTriangles(); ++j) {
+ tu->getTriangle(j)->setCol(col);
+ }
+ }
+ }
+
+ candidate.clear();
+ for (int i = 0; i < bestSequence.size(); ++i) {
+ candidate.insert(candidate.end(), bestSequence[i].first);
+ }
+}
+
+void Recombine2D::_getIncompatibles(const TrianglesUnion *tu, setTriUnion &set)
+{
+ for (int i = 0; i < tu->getNumTriangles(); i++) {
+ mapTriUnion::iterator it = _lessActions.find(tu->getTriangle(i));
+ if (it == _lessActions.end()) {
+ it = _possibleActions.find(tu->getTriangle(i));
+ if (it == _possibleActions.end()) {
+ Msg::Error("[Rec2D] error no triangle !, stopping with partially filled set");
+ }
+ _lessActions.insert(*it);
+ }
+ setTriUnion::iterator it2 = it->second.begin();
+ for (; it2 != it->second.end(); ++it2) {
+ set.insert(*it2);
+ }
+ }
+}
+
+void Recombine2D::_getNeighbors(const TrianglesUnion *tu, setTriUnion &set)
+{
+ setTriUnion incomp;
+ setTriUnion neighbors;
+ _getIncompatibles(tu, incomp);
+ setTriUnion::iterator it = incomp.begin();
+ for (; it != incomp.end(); ++it)
+ _getIncompatibles(*it, neighbors);
+ it = incomp.begin();
+ for (; it != incomp.end(); ++it)
+ neighbors.erase(*it);
+ it = neighbors.begin();
+ for (; it != neighbors.end(); ++it)
+ set.insert(*it);
+}
+
+void Recombine2D::_addNeighbors(int horiz, std::vector<list_actions::iterator> &vecIt, int current, list_actions *la)
+{
+ if (current < horiz - 1 && current > -1) {
+ //Msg::Info("+ra %d", ++ra);
+ setTriUnion neighbors;
+ _getNeighbors(*vecIt[current], neighbors);
+ la->add(neighbors);
+ }
+}
+
+void Recombine2D::_removeNeighbors(int horiz, int current, list_actions *la)
+{
+ if (current < horiz - 1 && current > -1) {
+ //Msg::Info("-ra %d", --ra);
+ la->pop_back();
+ }
+}
+
+
+int Recombine2D::_checkIsolatedTri(std::vector<list_actions::iterator> &vecIt, int size, std::set<MTriangle*> &seqTri)
+{
+ setTriUnion setTu;
+ for (int i = 0; i < size; ++i) {
+ _getIncompatibles(*vecIt[i], setTu);
+ }
+
+ std::set<MTriangle*> setBoundaryTri;
+ setTriUnion::iterator it = setTu.begin();
+ for (; it != setTu.end(); ++it) {
+ for (int i = 0; i < (*it)->getNumTriangles(); i++)
+ if (seqTri.find((*it)->getTriangle(i)) == seqTri.end())
+ setBoundaryTri.insert((*it)->getTriangle(i));
+ }
+
+ int num = 0;
+ std::set<MTriangle*>::iterator itTri = setBoundaryTri.begin();
+ for (; itTri != setBoundaryTri.end(); ++itTri) {
+ mapTriUnion::iterator it1 = _lessActions.find(*itTri);
+ if (it1 == _lessActions.end()) {
+ it1 = _possibleActions.find(*itTri);
+ if (it1 == _possibleActions.end()) {
+ Msg::Error("[Rec2D] error no triangle !, stopping with partially filled set");
+ }
+ _lessActions.insert(*it1);
+ }
+ int numPossibleRecombinations = 0;
+ setTriUnion::iterator it2 = it1->second.begin();
+ for (; it2 != it1->second.end(); ++it2)
+ if (setTu.find(*it2) == setTu.end())
+ ++numPossibleRecombinations;
+ if (!numPossibleRecombinations)
+ ++num;
+ }
+ return num;
}
MQuadrangle *TrianglesUnion::createQuad() const
@@ -222,24 +532,24 @@ MQuadrangle *TrianglesUnion::createQuad() const
v2 = listBoundEdge.begin()->getMaxVertex();
std::list<MEdge>::iterator it = listBoundEdge.begin();
for (it++; it != listBoundEdge.end(); it++) {
- if (v2 == (*it).getMinVertex()) {
- v3 = (*it).getMaxVertex();
+ if (v2 == it->getMinVertex()) {
+ v3 = it->getMaxVertex();
listBoundEdge.erase(it);
break;
}
- if (v2 == (*it).getMaxVertex()) {
- v3 = (*it).getMinVertex();
+ if (v2 == it->getMaxVertex()) {
+ v3 = it->getMinVertex();
listBoundEdge.erase(it);
break;
}
}
for (it = listBoundEdge.begin(); it != listBoundEdge.end(); it++) {
- if (v3 == (*it).getMinVertex()) {
- v4 = (*it).getMaxVertex();
+ if (v3 == it->getMinVertex()) {
+ v4 = it->getMaxVertex();
break;
}
- if (v3 == (*it).getMaxVertex()) {
- v4 = (*it).getMinVertex();
+ if (v3 == it->getMaxVertex()) {
+ v4 = it->getMinVertex();
break;
}
}
@@ -250,50 +560,96 @@ MQuadrangle *TrianglesUnion::createQuad() const
int Recombine2D::apply(bool a)
{
- if (!_haveParam || _applied) return 0;
-
- std::vector<MTriangle*> triangles2;
- for (int i = 0; i < _gf->triangles.size(); i++) {
- if (_isolated.find(_gf->triangles[i]) != _isolated.end())
- delete _gf->triangles[i];
- else
- triangles2.push_back(_gf->triangles[i]);
+ //Msg::Info("(%d, %d, %d)", _quads.size(), _gf->triangles.size(), _isolated.size());
+ if (a) {
+ std::vector<MTriangle*> triangles2;
+ for (int i = 0; i < _gf->triangles.size(); i++) {
+ if (_isolated.find(_gf->triangles[i]) != _isolated.end())
+ delete _gf->triangles[i];
+ else
+ triangles2.push_back(_gf->triangles[i]);
+ }
+ _gf->triangles = triangles2;
+ _gf->quadrangles = _quads;
+ _isolated.clear();
}
- _gf->triangles = triangles2;
- _gf->quadrangles = _quads;
- _isolated.clear();
_applied = true;
return 1;
}
-void Recombine2D::_removeImpossible(std::list<TrianglesUnion*>::iterator it)
+void Recombine2D::_removeImpossible(TrianglesUnion *tu)
{
- std::set<MTriangle *> touched;
+ for (int i = 0; i < tu->getNumTriangles(); i++) {
+ _isolated.insert(tu->getTriangle(i));
+ }// for test
+
+ std::set<MTriangle*> touched;
+ for (int i = 0; i < tu->getNumTriangles(); i++) {
+ touched.insert(tu->getTriangle(i));
+ }
+
+ setTriUnion incomp;
+ std::set<MTriangle*>::iterator itTri = touched.begin();
+ for (; itTri != touched.end(); ++itTri) {
+ mapTriUnion::iterator it = _possibleActions.find(*itTri);
+ if (it != _possibleActions.end()) {
+ setTriUnion::iterator it2 = it->second.begin();
+ for (; it2 != it->second.end(); ++it2)
+ incomp.insert(*it2);
+ }
+ }
+
+ setTriUnion::iterator itTu = incomp.begin();
+ for (; itTu != incomp.end(); ++itTu) {
+ for (int i = 0; i < (*itTu)->getNumTriangles(); i++) {
+ mapTriUnion::iterator it4 = _possibleActions.find((*itTu)->getTriangle(i));
+ it4->second.erase(*itTu);
+ TrianglesUnion *tu;
+ if (touched.find(it4->first) == touched.end())
+ switch (it4->second.size()) {
+ case 1 :
+ _lastQuad.insert(_lastQuad.begin(), *it4->second.begin());
+ tu = *it4->second.begin();
+ unsigned int col;
+ col = CTX::instance()->packColor(0, 0, 0, 0);
+ for (int j = 0; j < tu->getNumTriangles(); ++j)
+ tu->getTriangle(j)->setCol(col);
+ break;
+ case 0 :
+ default :
+ break;
+ }
+ }
+ //_setQuads.erase(*itTu); //marche pas avec list
+ }
- for (int i = 0; i < (*it)->getNumTriangles(); i++) {
- touched.insert((*it)->getTriangle(i));
+ std::list<TrianglesUnion*>::iterator it = _setQuads.begin();
+ while (it != _setQuads.end()) {
+ if (incomp.find(*it) != incomp.end())
+ it = _setQuads.erase(it); // replacement de _setQuads.erase(*itTu);
+ else
+ ++it;
+ }
+
+ for (int i = 0; i < tu->getNumTriangles(); i++) {
+ _possibleActions.erase(tu->getTriangle(i));
}
- for (it = _setQuads.begin(); it != _setQuads.end();) {
+ /*std::list<TrianglesUnion*>::iterator it = _setQuads.begin();
+ while (it != _setQuads.end()) {
bool b = false;
for (int i = 0; i < (*it)->getNumTriangles(); i++) {
if (touched.find((*it)->getTriangle(i)) != touched.end())
b = true;
}
- if (b) {
- for (int i = 0; i < (*it)->getNumTriangles(); i++) {
- _isolated.insert((*it)->getTriangle(i));
- }
- _setQuads.erase(it++);
- }
+ if (b)
+ it = _setQuads.erase(it);
else
- it++;
- }
+ ++it;
+ }*/
}
-
-
Rec2d_vertex::Rec2d_vertex(MVertex *v,
std::list<MTriangle*> &setTri,
std::map<MVertex*, std::set<GEdge*> > &mapVert)
@@ -304,7 +660,7 @@ Rec2d_vertex::Rec2d_vertex(MVertex *v,
if (_onWhat < 0) {
std::map<MVertex*, std::set<GEdge*> >::iterator it = mapVert.find(v);
if (it != mapVert.end()) {
- _angle = _computeAngle(v, setTri, (*it).second);
+ _angle = _computeAngle(v, setTri, it->second);
}
else {
Msg::Warning("[meshRecombine2D] Can't compute corner angle, setting angle to pi/2");
@@ -324,7 +680,7 @@ void Rec2d_vertex::_initStaticTable()
int nE = 5, nF = 10;
_Vvalues = new double*[2];
- _Vvalues[0] = new double[nE]; //AVERIFIER SI CA MARCHE
+ _Vvalues[0] = new double[nE];
for (int i = 0; i < nE; i++)
_Vvalues[0][i] = 1. - fabs(2. / (i+1) - 1.);
_Vvalues[1] = new double[nF];
@@ -345,7 +701,7 @@ void Rec2d_vertex::_initStaticTable()
double Rec2d_vertex::_computeAngle(MVertex *v,
std::list<MTriangle*> &setTri,
std::set<GEdge*> &setEdge)
-{
+ {
if (setEdge.size() != 2) {
Msg::Warning("[meshRecombine2D] Wrong number of edge : %d, returning pi/2", setEdge.size());
return M_PI / 2.;
@@ -364,9 +720,9 @@ double Rec2d_vertex::_computeAngle(MVertex *v,
vectub -= vectv;
double normlb, normub;
if ((normlb = norm(vectlb)) > prec * (normub = norm(vectub)))
- firstDer1 = vectub;
+ firstDer1 = ge->firstDer(ge->getUpperBound());
else if (normub > prec * normlb)
- firstDer1 = vectlb;
+ firstDer1 = ge->firstDer(ge->getLowerBound());
else {
Msg::Warning("[meshRecombine2D] Bad precision, returning pi/2");
return M_PI / 2.;
@@ -375,12 +731,8 @@ double Rec2d_vertex::_computeAngle(MVertex *v,
firstDer0 = firstDer1;
}
- Msg::Info("-- %lf, %lf, %lf", norm(firstDer0), norm(firstDer1), dot(firstDer0, firstDer1));
- double n = 1 / norm(firstDer0);
- firstDer0 = firstDer0 * n;
- n = 1 / norm(firstDer1);
- firstDer1 = firstDer1 * n;
- Msg::Info("--Angles : %lf, %lf, %lf", norm(firstDer0), norm(firstDer1), dot(firstDer0, firstDer1));
+ firstDer0 = firstDer0 * (1 / norm(firstDer0));
+ firstDer1 = firstDer1 * (1 / norm(firstDer1));
double angle1 = acos(dot(firstDer0, firstDer1));
double angle2 = 2. * M_PI - angle1;
@@ -399,10 +751,6 @@ double Rec2d_vertex::_computeAngle(MVertex *v,
angleMesh += angle3Vertices(t->getVertex((k+1)%3), v, t->getVertex((k+2)%3));
}
- Msg::Info("Angles : %lf, %lf, %lf", angle1, angleMesh, angle2);
-
- return M_PI / 2.;
-
if (angleMesh < M_PI)
return angle1;
return angle2;
@@ -420,10 +768,14 @@ double Rec2d_vertex::getReward(int n)
{
if (n == 0) return .0;
if (_onWhat > -1) {
- if (n > 0)
- return _Vbenefs[_onWhat][_numEl-1];
- else
- return - _Vbenefs[_onWhat][_numEl-2];
+ switch (n) {
+ case 1 :
+ return _Vbenefs[_onWhat][_numEl-1];
+ case -1 :
+ return - _Vbenefs[_onWhat][_numEl-2];
+ default :
+ return _Vvalues[_onWhat][_numEl-1+n] - _Vvalues[_onWhat][_numEl-1];
+ }
}
else
return fabs(2./M_PI*_angle/_numEl - 1.) - fabs(2./M_PI*_angle/(_numEl + n) - 1.);
@@ -452,7 +804,7 @@ TrianglesUnion::TrianglesUnion(GFace *gf,
double sumlc = .0, u[2], v[2];
MVertex *vert[2];
for (int i = 0; i < 2; i++) {
- vert[i] = (*ite).getVertex(i);
+ vert[i] = ite->getVertex(i);
vert[i]->getParameter(0, u[i]);
vert[i]->getParameter(1, v[i]); // Warning : should check if vertex on face or on edge
std::map<MVertex*,double>::iterator itlc;
@@ -460,13 +812,13 @@ TrianglesUnion::TrianglesUnion(GFace *gf,
sumlc += v2lcUV[vert[i]] = BGM_MeshSize(gf, u[i], v[i], vert[i]->x(), vert[i]->y(), vert[i]->z());
}
else
- sumlc += (*itlc).second;
+ sumlc += itlc->second;
}
sumlc = .2; // FIXME BGM_MeshSize returns wrong meshsize
double length = _computeEdgeLength(gf, vert, u, v, 0);
- double a = 0;
+ double a = .0;
_embEdgeValue += length / sumlc * (a + (2-a) *_computeAlignment(*ite, t));
//Msg::Info("Edge a : {%lf/.1 = %lf <-> %lf} => %lf", length, 2 * length / sumlc, 1 - _computeAlignment(*ite, t),length / sumlc * (1 + _computeAlignment(*ite, t)));
}
@@ -536,17 +888,77 @@ void TrianglesUnion::_update()
{
if (_computation >= _RECOMPUTE)
return;
- double r = _globValIfExecuted;
double alpha = _ValVert - _embVertValue;
for (int i = 0; i < _numBoundVert; i++) {
alpha += _vertices[i]->getReward(-1);
}
alpha /= _NumVert - _numEmbVert;
double beta = (_ValEdge - _embEdgeValue) / (_NumEdge - _numEmbEdge);
- _globValIfExecuted = alpha * beta;
+ _globValIfExecuted = alpha * alpha * beta;
_computation = _RECOMPUTE;
}
+
+void TrianglesUnion::addTri(std::set<MTriangle*> &setTri) const
+{
+ for (int i = 0; i < _numTri; ++i)
+ setTri.insert(_triangles[i]);
+}
+
+void TrianglesUnion::removeTri(std::set<MTriangle*> &setTri) const
+{
+ for (int i = 0; i < _numTri; ++i)
+ setTri.erase(_triangles[i]);
+}
+
+bool TrianglesUnion::isIn(std::set<MTriangle*> &setTri) const
+{
+ for (int i = 0; i < _numTri; ++i)
+ if (setTri.find(_triangles[i]) != setTri.end())
+ return true;
+ return false;
+}
+
+void TrianglesUnion::addInfluence(int *num, double *val, std::map<Rec2d_vertex*, int> &mapVert) const
+{
+ num[0] += _numEmbVert;
+ num[1] += _numEmbEdge;
+ val[0] += _embVertValue;
+ val[1] += _embEdgeValue;
+ for (int i = 0; i < _numBoundVert; ++i) {
+ if (mapVert.find(_vertices[i]) == mapVert.end())
+ mapVert[_vertices[i]] = -1;
+ else
+ mapVert[_vertices[i]] -= 1;
+ }
+}
+
+void TrianglesUnion::removeInfluence(int *num, double *val, std::map<Rec2d_vertex*, int> &mapVert) const
+{
+ num[0] -= _numEmbVert;
+ num[1] -= _numEmbEdge;
+ val[0] -= _embVertValue;
+ val[1] -= _embEdgeValue;
+ for (int i = 0; i < _numBoundVert; ++i) {
+ mapVert[_vertices[i]] += 1;
+ }
+}
+
+double TrianglesUnion::computeReturn(const int *num,
+ const double *val,
+ const std::map<Rec2d_vertex*,
+ int> &mapVert)
+{
+ double alpha = _ValVert - val[0];
+ std::map<Rec2d_vertex*, int>::const_iterator it = mapVert.begin();
+ for (; it != mapVert.end(); ++it) {
+ alpha += it->first->getReward(it->second);
+ }
+ alpha /= _NumVert - num[0];
+ double beta = (_ValEdge - val[1]) / (_NumEdge - num[1]);
+ return alpha * beta;
+}
+
bool TrianglesUnion::operator<(TrianglesUnion &other)
{
_update();
@@ -555,6 +967,11 @@ bool TrianglesUnion::operator<(TrianglesUnion &other)
}
bool lessTriUnion::operator()(TrianglesUnion *tu1, TrianglesUnion *tu2) const
+{
+ return *tu1 < *tu2;
+}
+
+bool lessTriUnionInv::operator()(TrianglesUnion *tu1, TrianglesUnion *tu2) const
{
return *tu2 < *tu1;
}