diff --git a/Mesh/meshRecombine2D.cpp b/Mesh/meshRecombine2D.cpp
deleted file mode 100644
index 6df5959ff2b4560b7281d69395c9d203fc4cf1d6..0000000000000000000000000000000000000000
--- a/Mesh/meshRecombine2D.cpp
+++ /dev/null
@@ -1,490 +0,0 @@
-// Gmsh - Copyright (C) 1997-2016 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to the public mailing list <gmsh@onelab.info>.
-//
-// Contributor(s):
-// Amaury Johnen, adapted from meshGFaceOptimize
-//
-
-#include "meshRecombine2D.h"
-#include "GFace.h"
-#include <cmath>
-#include <FL/Fl.H>
-#include "drawContext.h"
-#include "FlGui.h"
-#include "Context.h"
-#include "OpenFile.h"//pas propre
-
-static int HORIZ = 20;
-
-static int NEW_NODE = 0;
-static int DEL_NODE = 0;
-
-Recombine2D::Recombine2D(GFace *gf, int horizon)
-: _horizon(horizon), _gf(gf), _benef(.0), _applied(true)
-{
- if (_horizon < 1) _horizon = 1;
-
- //return;
-
- //laplaceSmoothing(gf,100);
- //gf->model()->writeMSH("befSquare.msh");
-
- Msg::Info("Branching with horizon %d", _horizon = HORIZ);
-
- _haveParam = gf->geomType() != GEntity::DiscreteSurface || gf->getCompound();
- if(!_haveParam) {
- Msg::Warning("Surfaces must have prametrization in order to recombine triangles");
- //return;
- _haveParam=true;
- }
-
- e2t_cont adj;
- _buildEdgeToElement(gf->triangles, adj);
-
- std::map<MElement*, int> boundaryTriangle;
- std::pair<Set_Recomb::iterator, bool> ret;
-
- for (e2t_cont::iterator it = adj.begin(); it != adj.end(); ++it) {
- if (it->second.second) {
- RecombTriangle *rt = new RecombTriangle(it->first, it->second.first, it->second.second);
- ret = _pairs.insert(rt);
- while (ret.second == false) {
- rt->touch();
- ret = _pairs.insert(rt);
- }
- _possibleRecomb[it->second.first].insert(rt);
- _possibleRecomb[it->second.second].insert(rt);
- }
- else {
- if (boundaryTriangle.find(it->second.first) == boundaryTriangle.end())
- boundaryTriangle[it->second.first] = 1;
- else
- boundaryTriangle[it->second.first] = 2;
- }
- }
-
- /*std::map<MElement*, int>::iterator itt;
-
- for (itt = boundaryTriangle.begin(); itt != boundaryTriangle.end(); ++itt) {
- RecombTriangle *rt;
- if (itt->second == 1)
- rt = new RecombTriangle(itt->first, -100);
- else
- rt = new RecombTriangle(itt->first, -200);
-
- ret = _pairs.insert(rt);
- while (!ret.second) {
- rt->touch();
- ret = _pairs.insert(rt);
- }
- _possibleRecomb[itt->first].insert(rt);
- }*/
-
- Msg::Info("poss %d", _possibleRecomb.size());
- Msg::Info("pairs %d", _pairs.size());
- Msg::Info("tri %d", gf->triangles.size());
-
- _recombine();
-
- Msg::Info("new %d ", NEW_NODE);
- Msg::Info("del %d ", DEL_NODE);
- _applied = false;
-}
-
-Recombine2D::~Recombine2D()
-{
- if (!_applied)
- for (int i = 0; i < _quads.size(); i++)
- delete _quads[i];
-}
-
-template <class E>
-void Recombine2D::_buildEdgeToElement(std::vector<E*> &elements, e2t_cont &adj)
-{
- for (unsigned int i = 0; i < elements.size(); i++) {
- E *el = elements[i];
- for (int j = 0; j < el->getNumEdges(); j++) {
- MEdge e = el->getEdge(j);
- e2t_cont::iterator it = adj.find(e);
- if (it == adj.end()) {
- std::pair<MElement*, MElement*> one = std::make_pair(el, (MElement*)0);
- adj[e] = one;
- }
- else
- it->second.second = el;
- }
- }
-}
-
-void Recombine2D::_recombine()
-{
- //SetBoundingBox();
-
- /*Map_Tri_Recomb::iterator itt;
-
- for (itt = _possibleRecomb.begin(); itt != _possibleRecomb.end(); ++itt) {
- if
- }*/
-
- int i = 0;
- while ((!_pairs.empty() || !_lastRecomb.empty()) && i < 1000) {
- Msg::Info("%d",i);
- Set_Recomb::iterator it = _bestNextRecombination();
-
- /*if (it != _pairs.begin()) {
- std::set<MElement*> isolatedCP(_isolated);
- _isolated.clear();
- Set_Recomb::iterator itnp = _pairs.begin();
- int n = 0;
- while (itnp != it) {
- n++;
- for (int k = 0; k < (*itnp)->numTriangles(); k++){
- Msg::Info("recomb[%d] %d : %d", n, k+1, (*itnp)->triangle(k)->getNum());
- _isolated.insert((*itnp)->triangle(k));
- }
- itnp++;
- }
- Msg::Info("%d triangles et %d quads", _gf->triangles.size(), _gf->quadrangles.size());
- _applied = false;
- apply();
- _applied = false;
- CTX::instance()->mesh.changed = (ENT_ALL);
- drawContext::global()->draw();
- FlGui::instance()->check();
- drawContext::global()->draw();
- if (!Msg::GetAnswer("Continue ?", 1, "no", "yes"))
- Msg::Info("I continue anyway :p");
- _isolated = isolatedCP;
- }*/
-
- _benef += (*it)->getBenef();
- if ((*it)->isQuad())
- _quads.push_back((*it)->createQuad());
- else
- for (int i = 0; i < (*it)->numTriangles(); i++)
- _isolated.insert((*it)->triangle(i));
-
- /*Set_Recomb::iterator itmp = _pairs.begin();
- while (itmp != it) {
- _rmRT(*itmp);
- itmp = _pairs.begin();
- }*/
- _removeImpossible(it);
-
- i++;
- }
- Msg::Info("Done %d (%d)", i, _pairs.size());
-}
-
-Set_Recomb::iterator Recombine2D::_bestNextRecombination()
-{
- if (!_lastRecomb.empty())
- return _lastRecomb.begin();
-
-
- Map_Tri_Node touched;
- std::vector<Recomb2D_Node*> nodes(_horizon);
- Set_Recomb::iterator it = _pairs.begin(), itmax = _pairs.begin();
-
- int h = 0, nSkip = 0;
- double maxBenef = _horizon * -200.;
- bool haveOptimal = false;
- int maxH = 0;
-
- int numb = 0;
- while (!haveOptimal) {
- numb++;
- while (h < _horizon && it != _pairs.end()) {
- bool possible = true;
-
- Map_Tri_Node::iterator itm;
- for (int i = 0; i < (*it)->numTriangles(); i++) {
- itm = touched.find((*it)->triangle(i));
- if (itm != touched.end()) {
- itm->second->blocking(itm);
- possible = false;
- }
- }
-
- if (possible) {
- if (h == 0) nodes[h] = new Recomb2D_Node(it, nSkip, _horizon, touched);
- else nodes[h] = new Recomb2D_Node(it, nSkip, nodes[h-1]);
-NEW_NODE++;
- h++;
- nSkip = 0;
- if (h >= _horizon || nodes[h - 1]->getBound() < maxBenef)
- break;
- }
-
- it++;
- }
-
- if (--h < 0) haveOptimal = true;
-
- if (nodes[h]->getTotBenef() > maxBenef && h >= maxH) {
- maxBenef = nodes[h]->getTotBenef();
- itmax = nodes[0]->getItRecomb();
- maxH = h;
- }
-
- /*nSkip = 0;
- while (h >= 0 && nSkip < 2) {
- nodes[h]->erase();
- nSkip += nodes[h--]->getnSkip();
-
- if (h == _horizon - 2 && nodes[h]->getnSkip() >= 3) nodes[h--]->erase(); //probleme si h=-1
- while (h >= 0 &&
- nodes[h]->beSkipped() &&
- nodes[h]->getnSkip() >= 2 * (_horizon - h))
- nodes[h--]->erase();*/
- while (h >= 0 &&
- (nodes[h]->isBetter() ||
- nodes[h]->getnSkip() >= 2 * (_horizon - h) - 1)){
-DEL_NODE++;
- delete nodes[h--];
- }
-
- if (h < 0) haveOptimal = true;
- else {
- it = ++nodes[h]->getItRecomb();
- nSkip = nodes[h]->getnSkip() + 1;
- delete nodes[h];
-DEL_NODE++;
- }
- }
-
- Msg::Info("num test %d", numb);
-
- return itmax;
-}
-
-void Recombine2D::_removeImpossible(Set_Recomb::iterator it)
-{
- std::vector<MElement *> t;
- (*it)->triangles(t);
-
- Map_Tri_Recomb::iterator itm;
- std::set<RecombTriangle*>::iterator itrt;
-
- for (int i = 0; i < t.size(); i++) // loop over touched triangles
- if ((itm = _possibleRecomb.find(t[i])) != _possibleRecomb.end()) {
- itrt = itm->second.begin();
- while (itrt != itm->second.end()) { // loop over incompatible RecombTriangle
- if (*itrt != *it)
- _rmRT(*itrt, t[i]);
- itrt++;
- }
- _possibleRecomb.erase(itm);
- }
-
- RecombTriangle *rt = *it;
- if (!_lastRecomb.erase(*it))
- _pairs.erase(it);
- delete rt;
-}
-
-void Recombine2D::_rmRT(RecombTriangle *rt, MElement *el)
-{
- MElement *tri;
- Map_Tri_Recomb::iterator itmtri;
-
- for (int i = 0; i < rt->numTriangles(); i++)
- if ((tri = rt->triangle(i)) != el &&
- (itmtri = _possibleRecomb.find(tri)) != _possibleRecomb.end()) {
- itmtri->second.erase(rt);
- switch (itmtri->second.size()) {
- case 1 :
- _pairs.erase(*itmtri->second.begin());
- _lastRecomb.insert(*itmtri->second.begin());
- break;
- case 0 :
- _isolated.insert(tri);
- _possibleRecomb.erase(itmtri);
- default :;
- }
- }
-
- _pairs.erase(rt);
- _lastRecomb.erase(rt);
- delete rt;
-}
-
-int Recombine2D::apply()
-{
- 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]);
- }
- _gf->triangles = triangles2;
- _gf->quadrangles = _quads;
-
- _applied = true;
- //_gf->model()->writeMSH("recSquare.msh");
- laplaceSmoothing(_gf,10);
- //_gf->model()->writeMSH("aftSquare.msh");
- return 1;
-}
-
-RecombTriangle::RecombTriangle(const MEdge &me, MElement *t1, MElement *t2)
-: _formingQuad(true)
-{
- _t.push_back(t1);
- _t.push_back(t2);
- _n1 = me.getVertex(0);
- _n2 = me.getVertex(1);
- _sign = (std::rand() % 2) * 2 - 1;
-
- if(t1->getVertex(0) != _n1 && t1->getVertex(0) != _n2) _n3 = t1->getVertex(0);
- else if(t1->getVertex(1) != _n1 && t1->getVertex(1) != _n2) _n3 = t1->getVertex(1);
- else if(t1->getVertex(2) != _n1 && t1->getVertex(2) != _n2) _n3 = t1->getVertex(2);
- if(t2->getVertex(0) != _n1 && t2->getVertex(0) != _n2) _n4 = t2->getVertex(0);
- else if(t2->getVertex(1) != _n1 && t2->getVertex(1) != _n2) _n4 = t2->getVertex(1);
- else if(t2->getVertex(2) != _n1 && t2->getVertex(2) != _n2) _n4 = t2->getVertex(2);
-
- double a1 = 180 * angle3Vertices(_n1, _n4, _n2) / M_PI;
- double a2 = 180 * angle3Vertices(_n4, _n2, _n3) / M_PI;
- double a3 = 180 * angle3Vertices(_n2, _n3, _n1) / M_PI;
- double a4 = 180 * angle3Vertices(_n3, _n1, _n4) / M_PI;
- _angle = fabs(90. - a1);
- _angle = std::max(fabs(90. - a2), _angle);
- _angle = std::max(fabs(90. - a3), _angle);
- _angle = std::max(fabs(90. - a4), _angle);
-
- _benefit = std::max(1. - _angle/90., .0);
- //_benefit = compute_alignment(me,t1,t2);
-
- //cost_alignment = Temporary::compute_alignment(me,_t1,_t2); //addition for class Temporary
- //total_cost = Temporary::compute_total_cost(cost_quality,cost_alignment); //addition for class Temporary
- //total_cost = 100.0*cost_quality; //addition for class Temporary
-}
-
-RecombTriangle::RecombTriangle(MElement *t, double a)
-: _benefit(a), _formingQuad(false)
-{
- _t.push_back(t);
- _sign = (std::rand() % 2) * 2 - 1;
-}
-
-double RecombTriangle::touch()
-{
- if (_benefit == .0)
- return _benefit = _sign * std::pow(2., -1022);
- else // add smaller possible value
- return _benefit += _sign * std::pow(2., std::floor(std::log(std::fabs(_benefit))/std::log(2.))-52);
-}
-
-MQuadrangle *RecombTriangle::createQuad() const
-{
- switch (_t.size()) {
- case 1 :
- {
- Msg::Error("Can't make a quad from only one triangle");
- return new MQuadrangle(_t[0]->getVertex(0),
- _t[0]->getVertex(1),
- _t[0]->getVertex(2),
- _t[0]->getVertex(0));
- }
- case 2 :
- {
- int orientation = 0;
- for (int i = 0; i < 3; i++) {
- if (_t[0]->getVertex(i) == _n1) {
- if (_t[0]->getVertex((i + 1) % 3) == _n2)
- orientation = 1;
- else
- orientation = -1;
- break;
- }
- }
- if (orientation < 0)
- return new MQuadrangle(_n1, _n3, _n2, _n4);
- else
- return new MQuadrangle(_n1, _n4, _n2, _n3);
- }
- default :
- {
- Msg::Warning("Quad creation not yet implemented for more than two triangles");
- return new MQuadrangle(_n1, _n2, _n3, _n4);
- }
- }
-}
-
-double RecombTriangle::compute_alignment(const MEdge &edge, MElement *e1, MElement* e2)
-{
- return 1 - Temporary::compute_alignment(edge, e1, e2);
-}
-
-bool RecombTriangle::operator<(const RecombTriangle &other) const
-{
- return _benefit < other._benefit;
-}
-
-bool lessRecombTri::operator()(RecombTriangle *rt1, RecombTriangle *rt2) const
-{
- return *rt2 < *rt1;
-}
-
-//{Recomb2D_Node::Recomb2D_Node(Set_Recomb::iterator it,
-// Map_Tri_Node &touched,
-// int nskp, double ben) : _nskip(nskp), _recomb(it)
-//{
-// _benef = (*it)->getBenef();
-// _totBenef = ben + _benef;
-// for (int i = 0; i < (*it)->numTriangles(); i++)
-// _vit.push_back(touched.insert( touched.begin(), std::make_pair((*it)->triangle(i), this) ));
-// _touch = &touched;
-//}
-//}
-
-Recomb2D_Node::Recomb2D_Node(Set_Recomb::iterator it, int nskp, int horiz,
- Map_Tri_Node &touched)
-: _nskip(nskp), _recomb(it)
-{
- _benef = (*it)->getBenef();
- _totBenef = _benef;
- _depth = horiz - 1;
- for (int i = 0; i < (*it)->numTriangles(); i++)
- _vit.push_back(touched.insert( touched.begin(), std::make_pair((*it)->triangle(i), this) ));
- _touch = &touched;
-}
-
-Recomb2D_Node::Recomb2D_Node(Set_Recomb::iterator it, int nskp, Recomb2D_Node *n)
-: _nskip(nskp), _recomb(it)
-{
- _benef = (*it)->getBenef();
- _totBenef = n->_benef + _benef;
- _depth = n->_depth - 1;
- _touch = n->_touch;
- for (int i = 0; i < (*it)->numTriangles(); i++)
- _vit.push_back(_touch->insert( _touch->begin(), std::make_pair((*it)->triangle(i), this) ));
-}
-
-Recomb2D_Node::~Recomb2D_Node()
-{
- for (int i = 0; i < _vit.size(); i++)
- _touch->erase(_vit[i]);
-}
-
-void Recomb2D_Node::erase()
-{
- for (int i = 0; i < _vit.size(); i++)
- _touch->erase(_vit[i]);
-}
-
-void Recomb2D_Node::blocking(const Map_Tri_Node::iterator &it)
-{
- for (int i = 0; i < _vit.size(); i++) {
- if (it == _vit[i]) {
- _blocking.insert(i);
- }
- }
-}
-
-//
diff --git a/Mesh/meshRecombine2D.h b/Mesh/meshRecombine2D.h
deleted file mode 100644
index 1b3faa39640c5a231e3243904fd5b895c354a456..0000000000000000000000000000000000000000
--- a/Mesh/meshRecombine2D.h
+++ /dev/null
@@ -1,420 +0,0 @@
-// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
-//
-// Contributor(s):
-// Amaury Johnen
-//
-
-#ifndef _MESH_RECOMBINE_2D_H_
-#define _MESH_RECOMBINE_2D_H_
-
-#include <map>
-#include <set>
-#include <vector>
-#include "MVertex.h"
-#include "MTriangle.h"
-#include "MQuadrangle.h"
-#include "meshGFaceOptimize.h"
-
-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);
- }
- return *this;
- }
- iterator& operator--() {
- int p = pos;
- while (it == la->_begin(p) && --p > -1) {
- --pos;
- it = la->_end(pos);
- }
- if (p > -1)
- --it;
- 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 :
- GFace *_gf;
- int _horizon;
- double _benef;
- bool _haveParam, _applied;
- Set_Recomb _pairs, _lastRecomb;
- 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 &);
- void _recombine();
- Set_Recomb::iterator _bestNextRecombination();
- void _removeImpossible(Set_Recomb::iterator);
- void _rmRT(RecombTriangle *, MElement *e = NULL);
-
- public :
- Recombine2D(GFace*, int horizon);
- Recombine2D(GFace*);
- ~Recombine2D();
-
- int apply();
- inline double getBenef() const {return _benef;}
- inline int numTriangle() const {return _isolated.size();}
-
- private :
- //std::set<TrianglesUnion*, lessTriUnion> _setQuads;
- std::list<TrianglesUnion*> _setQuads;
- std::list<TrianglesUnion*> _lastQuad;
- mapTriUnion _possibleActions;
- mapTriUnion _lessActions;
- void _removeImpossible(TrianglesUnion*);
- void _recombine(bool);
-
- void _lookAhead(std::list<TrianglesUnion*> &, int horiz);
- void _lookAhead2(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(const std::vector<list_actions::iterator> &,
- int size, std::set<MTriangle*> &);
- double _realisticReturn(const int *num,
- const double *val,
- const std::map<Rec2d_vertex*, int> &,
- const std::vector<list_actions::iterator> &,
- int size,
- const std::set<MTriangle*> &);
-
- public :
- int apply(bool);
-};
-
-class RecombTriangle {
- private :
- std::vector<MElement *> _t;
- MVertex *_n1, *_n2, *_n3, *_n4;
- double _angle;
- double _benefit;
- bool _formingQuad;
- int _sign;
-
- public :
- RecombTriangle(const MEdge &, MElement *, MElement *);
- RecombTriangle(MElement *, double);
-
- double touch();
-
- MQuadrangle *createQuad() const;
- bool operator<(const RecombTriangle &other) const;
-
- inline void triangles(std::vector<MElement *> &v) const {v = _t;}
- inline int numTriangles() const {return (int) _t.size();}
- MElement *triangle(int i) const {
- if (i >= 0 && i < _t.size())
- return _t[i];
- return NULL;
- }
- inline bool isQuad() const {return _formingQuad;}
- inline double getBenef() const {return _benefit;}
-
- double compute_alignment(const MEdge&, MElement*, MElement*);
-};
-
-class Recomb2D_Node {
- private :
- double _benef, _totBenef;
- int _nskip, _depth;
- Set_Recomb::iterator _recomb;
- Map_Tri_Node *_touch;
- std::vector<Map_Tri_Node::iterator> _vit;
- std::set<int> _blocking;
-
- public :
- //Node(){}
- //Recomb2D_Node(Set_Recomb::iterator, Map_Tri_Node &, int, double ben = .0);
-
- Recomb2D_Node(Set_Recomb::iterator, int nskip, int depth, Map_Tri_Node &);
- Recomb2D_Node(Set_Recomb::iterator, int nskip, Recomb2D_Node *);
- ~Recomb2D_Node();
-
- void erase();
- void blocking(const Map_Tri_Node::iterator &);
- inline bool isBetter() {return _blocking.size() < 2;}
- inline Set_Recomb::iterator getItRecomb() {return _recomb;}
- inline double getTotBenef() {return _totBenef;}
- inline int getnSkip() {return _nskip;}
- inline double getBound() {return _totBenef + _benef * _depth;}
-};
-
-
-class Rec2d_vertex {
- private :
- // _onWhat={-1:corner,0:edge,1:face,(2:volume)}
- int _onWhat, _numEl;
- double _angle;
- // GEdge : virtual SVector3 firstDer(double par) const = 0;
- // GEntity : virtual Range<double> parBounds(int i)
- static double **_Vvalues;
- static double **_Vbenefs;
-
- public :
- Rec2d_vertex(MVertex *, std::list<MTriangle*> &,
- std::map<MVertex*, std::set<GEdge*> > &, bool);
-
- 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 {
- private :
- int _numEmbEdge, _numboundEdge, _numEmbVert, _numBoundVert, _numTri, _computation;
- double _embEdgeValue, _boundEdgeValue, _embVertValue, _globValIfExecuted;
- Rec2d_vertex **_vertices;
- MTriangle **_triangles;
- static int _RECOMPUTE; //incremented when a recombination is executed
-
- public:
- static int _NumEdge, _NumVert;
- static double _ValEdge, _ValVert;
- private:
-
- public :
- TrianglesUnion(GFace *, std::list<MTriangle*> &, std::list<MEdge> &,
- std::list<Rec2d_vertex*> &, std::map<MVertex*,double> &,
- std::map<MEdge, std::list<MTriangle*>, Less_Edge>&);
-
- bool operator<(TrianglesUnion &other);
- 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;
- _ValEdge += _boundEdgeValue;
- _NumEdge += _numboundEdge;
- _ValVert -= _embVertValue;
- _NumVert -= _numEmbVert;
- for (int i = 0; i < _numBoundVert; i++) {
- _ValVert += _vertices[i]->getReward(-1);
- _vertices[i]->changeNumEl(-1);
- }
- _RECOMPUTE++;
- }
- MQuadrangle* createQuad() const;
- void print() {
- Msg::Info("Printing TU (%d,%d,%d,%d)", _numEmbVert, _numBoundVert, _numEmbEdge, _numTri);
- for (int i = 0; i < _numTri; i++)
- Msg::Info("Triangle %d", _triangles[i]->getNum());
- }
- static inline void addRec() {_RECOMPUTE++;}
- inline double getReturn() {return _globValIfExecuted;}
- inline static double getActualReturn() {
- return _ValEdge/_NumEdge * _ValVert/_NumVert * _ValVert/_NumVert;
- }
-
- static void clear() {_NumEdge = 0; _NumVert = 0;
- _ValEdge = .0, _ValVert = .0;}
-
- private:
- double _computeEdgeLength(GFace*, MVertex**, double *u, double *v,
- int numIntermedPoint= 1);
- double _computeAlignment(const MEdge&, std::list<MTriangle*>&,
- std::map<MEdge, std::list<MTriangle*>, Less_Edge>&);
- void _update();
-};
-#endif
diff --git a/Mesh/meshRecombine2D_2.cpp b/Mesh/meshRecombine2D_2.cpp
deleted file mode 100644
index 20abf4208379c64c8af2a194eee54a15f9557712..0000000000000000000000000000000000000000
--- a/Mesh/meshRecombine2D_2.cpp
+++ /dev/null
@@ -1,1450 +0,0 @@
-// Gmsh - Copyright (C) 1997-2016 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to the public mailing list <gmsh@onelab.info>.
-//
-// Contributor(s):
-// Amaury Johnen, adapted from meshGFaceOptimize
-//
-
-#include "meshRecombine2D.h"
-#include "BackgroundMeshTools.h"
-#include "GFace.h"
-#include <cmath>
-#include <FL/Fl.H>
-#include "drawContext.h"
-#include "FlGui.h"
-#include "Context.h"
-#include "OpenFile.h"//pas propre
-#include "Field.h"
-#include "OS.h"
-#define HORIZ2 1
-#define SMOOTH 0
-
-double **Rec2d_vertex::_Vvalues = NULL;
-double **Rec2d_vertex::_Vbenefs = NULL;
-int TrianglesUnion::_RECOMPUTE = 0;
-int TrianglesUnion::_NumEdge = 0;
-int TrianglesUnion::_NumVert = 0;
-double TrianglesUnion::_ValEdge = .0;
-double TrianglesUnion::_ValVert = .0;
-int Recombine2D::ra = 0;
-
-/*
- a faire :
- - fonction retournant qualite
- - verifier la qualite
- - action supplementaire (faire classe de base)
- - Faire un super conteneur pour les actions
- - pourquoi infini
- - pourquoi premier pas meilleur
- - implementer field
- ~ limiter sequence
-*/
-
-Recombine2D::Recombine2D(GFace *gf)
-: _horizon(0), _gf(gf), _benef(.0), _applied(true)
-{
- //GModel::current()->createTopologyFromMesh(1);
- int popopo = 0;
-
- TrianglesUnion::clear();
- Msg::Warning("[meshRecombine2D] Alignement computation ok uniquely for xy or yz plane mesh !");
-
- // be able to compute geometrical angle at corners
- std::map<MVertex*, std::set<GEdge*> > mapCornerVert;
- {
- std::list<GEdge*> listge = gf->edges();
- std::list<GEdge*>::iterator itge = listge.begin();
-
- for (; itge != listge.end(); itge++) {
- if ((*itge)->getBeginVertex()->getNumMeshElements() > 1 ||
- (*itge)->getEndVertex()->getNumMeshElements() > 1 )
- Msg::Warning("[meshRecombine2D] Why more than one MPoint, what should I do ?");
- mapCornerVert[(*itge)->getBeginVertex()->getMeshElement(0)->getVertex(0)].insert(*itge);
- mapCornerVert[(*itge)->getEndVertex()->getMeshElement(0)->getVertex(0)].insert(*itge);
- }
- }
-
- // Find all Vertices and edges
- std::map<MVertex*, std::list<MTriangle*> > mapVert;
- std::map<MEdge, std::list<MTriangle*>, Less_Edge> mapEdge;
- std::map<MVertex*, std::list<MEdge> > mapVtoE;
- for (unsigned int i = 0; i < gf->triangles.size(); i++) {
- MTriangle *t = gf->triangles[i];
- for (int j = 0; j < 3; j++) {
- mapVert[t->getVertex(j)].push_back(t);
- mapEdge[t->getEdge(j)].push_back(t);
- for (int k = 0; k < 3; k++) {
- mapVtoE[t->getVertex(j)].push_back(t->getEdge(k));
- }
- }
- }
-
- // Create the 'Rec2d_vertex' and store iterator to vertices which have degree 4
- std::map<MVertex*, std::list<MTriangle*> >::iterator itvert = mapVert.begin();
- std::map<MVertex*, Rec2d_vertex*> mapV2N;
- 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)
- fourTri.push_back(itvert);
- std::list<MEdge>::iterator itti = mapVtoE[itvert->first].begin();
- bool ab = false;
- for (; itti != mapVtoE[itvert->first].end(); ++itti) {
- if (mapEdge[*itti].size() == 1) ab = true;
- }
- if (ab) ++popopo;
- Rec2d_vertex *n = new Rec2d_vertex(itvert->first, itvert->second, mapCornerVert, ab);
- itV2N = mapV2N.insert(itV2N, std::make_pair(itvert->first,n));
- TrianglesUnion::_NumVert++;
- TrianglesUnion::_ValVert += n->getReward();
- }
-
- // store mesh size for better performance
- std::map<MVertex*,double> mapV2LcUV;
-
- // Create 'TrianglesUnion' for pattern of 4 triangles
-/*
- +-----+
- |\ /|
- | \ / |
- | + |
- | / \ |
- |/ \|
- +-----+
-*/
- std::list<std::map<MVertex*, std::list<MTriangle*> >::iterator>::iterator it4;
- int j = 0;
- for (it4 = fourTri.begin(); it4 != fourTri.end(); it4++) {
- std::list<MEdge> listEmbEdge;
- std::list<Rec2d_vertex*> listVert;
- {
- std::set<MVertex*> setVert;
- std::multiset<MEdge, Less_Edge> msetEdge;
- 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)
- setVert.insert(vt);
- }
- }
- std::multiset<MEdge>::iterator itEdge = msetEdge.begin();
- MEdge me = *(itEdge++);
- for (; itEdge != msetEdge.end(); itEdge++) {
- if (*itEdge == me)
- listEmbEdge.push_back(*itEdge);
- me = *itEdge;
- }
- std::set<MVertex*>::iterator itsVert = setVert.begin();
- for (; itsVert != setVert.end(); itsVert++)
- listVert.push_back(mapV2N[*itsVert]);
- listVert.push_back(mapV2N[(*it4)->first]);
- }
- Msg::Info("%d",++j);
- TrianglesUnion *tu = new TrianglesUnion(gf, (*it4)->second, listEmbEdge, listVert, mapV2LcUV, mapEdge);
- 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
-/*
- +---+
- |\ |
- | \ |
- | \|
- +---+
-*/
- std::map<MEdge, std::list<MTriangle*> >::iterator itedge = mapEdge.begin();
- for (; itedge != mapEdge.end(); itedge++) {
- if (itedge->second.size() == 2) {
- std::list<MEdge> listEmbEdge;
- 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, mapEdge);
- _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());
- }
- TrianglesUnion::addRec();
- //_setQuads.sort(lessTriUnion());
-
- if (SMOOTH) laplaceSmoothing(_gf,10);
- Msg::Error("%d", popopo);
- _recombine(true);
- _applied = false;
-}
-
-
-void Recombine2D::_recombine(bool a)
-{
- SetBoundingBox();
-
- int i = 0;
- std::list<TrianglesUnion*> vectTu;
- TrianglesUnion::addRec();
-
- while (!_setQuads.empty() && i < 2000) {
- TrianglesUnion *tu;
- if (SMOOTH) laplaceSmoothing(_gf,1);
- if (_lastQuad.empty()) {
- if (/*true ||*/ vectTu.size() < 2 || i % 1 == 0)
- 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;
- }
- }
- tu = *vectTu.begin();
- tu->select();
- _quads.push_back(tu->createQuad());
- _removeImpossible(tu);
- if (i % 5 == 0) { // draw Mesh
- Msg::Info("Expected return %lf", tu->getReturn());
- _applied = false;
- apply(true);
- _applied = false;
- CTX::instance()->mesh.changed = (ENT_ALL);
- //drawContext::global()->draw();
- FlGui::instance()->check();
- drawContext::global()->draw();
- 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();
-
- int t = (int) Cpu() * 100;
- bool newt = true;
-
- double maxReturn = -100.;
- //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 = _realisticReturn(numbers, values, modifiedVert, 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);
- }
- }
-
- if (/*true ||*/ best /*|| (((int)(Cpu()*100) - t) % 5 == 0 && newt)*/) { // draw mesh
- newt = false;
- apply(false);
- if (SMOOTH) laplaceSmoothing(_gf,1);
- 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;
- }
- /*if (((int)(Cpu()*100) - t) % 5 != 0)
- newt = true;*/
-
- { // 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);
- }
- }
- }
-
- 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);
- }
- }
- }
-
- { // 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::_lookAhead2(std::list<TrianglesUnion*> &candidate, int horiz)
-{
- if (horiz < 1)
- horiz = 1;
- _lessActions.clear();
-
- //int t = (int) Cpu() * 100;
- //bool newt = true;
-
- double maxReturn = -100.;
- 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);
-
- std::vector<list_actions::iterator> lastLayer(horiz);
- std::vector<int> numLayer(horiz);
- lastLayer[0] = la.end();
- numLayer[0] = 0;
- bool haveSequence = false;
- int current = 1, currentLayer = 0, maxLayer = 0;
- int num=0;
-
- while (current > 0) {
- //Msg::Info(" ");
- bool best = false;
- // add max of actions in sequence
- while (current < horiz && vecIt[current-1] != la.end()) {
- vecIt[current] = vecIt[current-1];
- //vecIt[current].print();
- while (++vecIt[current] != lastLayer[currentLayer]
- && vecIt[current] != la.end()
- && (*vecIt[current])->isIn(setTri) );
- if (vecIt[current] != lastLayer[currentLayer] && vecIt[current] == la.end())
- {Msg::Error("Ca a foir� mec");Msg::Error(" ");}
- if (vecIt[current] != lastLayer[currentLayer] && vecIt[current] != la.end()) {
- (*vecIt[current])->addTri(setTri);
- (*vecIt[current])->addInfluence(numbers, values, modifiedVert);
- if (currentLayer < maxLayer) {
- if (current == numLayer[currentLayer])
- lastLayer[currentLayer] = --la.end();
- _addNeighbors(horiz, vecIt, current, &la);
- if (current == numLayer[currentLayer])
- ++lastLayer[currentLayer];
- lastLayer[currentLayer+1] = la.end();
- }
- ++current;
- }
- else {
- //Msg::Info("in else");
- if (!haveSequence) {
- ++maxLayer;
- --lastLayer[currentLayer];
- for (int i = numLayer[currentLayer]; i < current; ++i) {
- _addNeighbors(horiz, vecIt, i, &la);
- }
- ++lastLayer[currentLayer];
- ++currentLayer;
- lastLayer[currentLayer] = la.end();
- numLayer[currentLayer] = current;
- }
- else if (currentLayer < maxLayer) {
- /*for (int i = numLayer[currentLayer]; i < current; ++i) {
- _addNeighbors(horiz, vecIt, i, &la);
- }*/
- //Msg::Error("++layer'");
- ++currentLayer;
- lastLayer[currentLayer] = la.end();
- numLayer[currentLayer] = current;
- }
- else{
- //Msg::Error("++current : %d/%d", current, horiz);
- //if(lastLayer[currentLayer] != la.end())
- //Msg::Error(" ");
- ++current;
- }
- }
- }
- haveSequence = true;
- //Msg::Info("maxLayer %d, current %d", maxLayer, current);
-
-/*Msg::Info("LAYER %d - num{%d,%d} - last :", currentLayer, numLayer[0], numLayer[1]);
-lastLayer[0].print();
-lastLayer[1].print();
-Msg::Info("current %d - iterator : ", current-1);
-vecIt[current-1].print();
-Msg::Info(" ");*/
-
- { // 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 = _realisticReturn(numbers, values, modifiedVert, vecIt, sizeSequence, setTri);
- if (maxReturn < expectedReturn) {
- //Msg::Info("best");
- //Msg::Info(" ");
- 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, 130, 0, 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);
- }
- }
-
- if (/*true ||*/ best /*|| (((int)(Cpu()*100) - t) % 5 == 0 && newt)*/) { // draw mesh
- //newt = false;
- apply(false);
- CTX::instance()->mesh.changed = (ENT_ALL);
- FlGui::instance()->check();
- //Msg::Info("before");
- //Msg::Info(" ");
- drawContext::global()->draw();
- //Msg::Info("after");
- //Msg::Info(" ");
- 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;
- }
- /*if (((int)(Cpu()*100) - t) % 5 != 0)
- newt = true;*/
-
- { // 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);
- }
- }
- }
-
- if (vecIt[current-1] != la.end() || --current > 0) {
- do {
- while (current-1 < numLayer[currentLayer])
- --currentLayer;
- (*vecIt[current-1])->removeTri(setTri);
- (*vecIt[current-1])->removeInfluence(numbers, values, modifiedVert);
- //Msg::Info(" - current %d, currentLayer %d", current, currentLayer);
- if (currentLayer < maxLayer) {
- _removeNeighbors(horiz, current - 1, &la);
- if (current-1 == numLayer[currentLayer])
- lastLayer[currentLayer] = la.end();
- else
- lastLayer[currentLayer+1] = la.end();
- }
- while (++vecIt[current-1] != la.end() && (*vecIt[current-1])->isIn(setTri)) {
- //Msg::Info(" - - curent %d & layer %d & numLayer %d", current, currentLayer, numLayer[currentLayer]);
- //vecIt[current-1].print(); lastLayer[currentLayer].print();
- if (vecIt[current-1] == lastLayer[currentLayer]) {
- //Msg::Warning("++layer");
- ++currentLayer;
- numLayer[currentLayer] = current - 1;
- }
- }
- if (currentLayer < maxLayer && vecIt[current-1] == lastLayer[currentLayer]) {
- //Msg::Warning("++ layer");
- //vecIt[current-1].print(); lastLayer[currentLayer].print();
- ++currentLayer;
- numLayer[currentLayer] = current - 1;
- }
- //Msg::Info("b %d, %d - %d, %d", vecIt[current-1] == lastLayer[currentLayer], current-2, numLayer[currentLayer], currentLayer);
- } while (vecIt[current-1] == la.end() && --current > 0);
- //Msg::Info("current %d, currentLayer %d", current, currentLayer);
- if (current > 0) {
- (*vecIt[current-1])->addTri(setTri);
- (*vecIt[current-1])->addInfluence(numbers, values, modifiedVert);
- if (currentLayer < maxLayer) {
- if (current-1 == numLayer[currentLayer])
- --lastLayer[currentLayer];
- _addNeighbors(horiz, vecIt, current - 1, &la);
- if (current-1 == numLayer[currentLayer])
- ++lastLayer[currentLayer];
- lastLayer[currentLayer+1] = la.end();
- }
- }
- }
-/* if (current > 0) {
-Msg::Info("layer %d - num{%d,%d} - last :", currentLayer, numLayer[0], numLayer[1]);
-lastLayer[0].print();
-lastLayer[1].print();
-Msg::Info("current %d - iterator : ", current-1);
-vecIt[current-1].print();
-Msg::Info(" ");
-}*/
- }
-
- /*{ // 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");
- Msg::Error(" ");
- }
- _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();
- }
-}
-
-
-double Recombine2D::_realisticReturn(const int *num,
- const double *val,
- const std::map<Rec2d_vertex*, int> &mapVert,
- const std::vector<list_actions::iterator> &vecIt,
- int size,
- const std::set<MTriangle*> &seqTri)
-{
- //Msg::Info("before (%d) : %lf -> %lf", size, TrianglesUnion::getActualReturn(), TrianglesUnion::computeReturn(num, val, mapVert));
- int newNum[2] = {num[0], num[1]}, numIsolatedTri, numRecomb = size;
- double newVal[2] = {val[0], val[1]};
- std::map<Rec2d_vertex*, int> newMapVert = mapVert;
-
- setTriUnion setIncompTu;
- for (int i = 0; i < size; ++i) {
- _getIncompatibles(*vecIt[i], setIncompTu);
- }
-
- std::set<MTriangle*> setTri = seqTri, touched;
-
- numIsolatedTri = _checkIsolatedTri(vecIt, size, setTri);
- //Msg::Info("numisolated %d, last %d", numIsolatedTri, setTri.size());
-
- std::set<MTriangle*>::iterator itTri = setTri.begin();
- int kl = 0;
- while (itTri != setTri.end()) {
- //Msg::Info("%d.", ++kl);
- mapTriUnion::iterator itTmp = _possibleActions.find(*itTri);
- setTriUnion::iterator it_tu = itTmp->second.begin();
- TrianglesUnion *tu;
- int numPossible = 0;
- for (; it_tu != itTmp->second.end(); ++it_tu) {
- bool possible = true;
- TrianglesUnion *tu1 = *it_tu;
- for (int i = 0; i < tu1->getNumTriangles(); i++)
- if (seqTri.find(tu1->getTriangle(i)) != seqTri.end() ||
- touched.find(tu1->getTriangle(i)) != touched.end() )
- possible = false;
- if (possible) {
- tu = tu1;
- ++numPossible;
- }
- }
-
- setTri.erase(itTri);
-
- if (numPossible < 1){
- ++numIsolatedTri;
- Msg::Info("'m here");
- }
- else if (numPossible > 1)
- Msg::Info("Oh oh, trouble");
- else {
- //Msg::Info("possible");
- ++numRecomb;
- tu->addInfluence(newNum, newVal, newMapVert);
- for (int i = 0; i < tu->getNumTriangles(); i++) {
- touched.insert(tu->getTriangle(i));
- setTri.erase(tu->getTriangle(i));
- }
- _getIncompatibles(tu, setIncompTu);
-
- setTriUnion setTu;
- setTu.clear();
- _getIncompatibles(tu, 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() &&
- touched.find((*it)->getTriangle(i)) == touched.end() )
- setBoundaryTri.insert((*it)->getTriangle(i));
- }
-
- //Msg::Info("size boundary : %d", setBoundaryTri.size());
-
- std::set<MTriangle*>::iterator itTri2 = setBoundaryTri.begin();
- for (; itTri2 != setBoundaryTri.end(); ++itTri2) {
- mapTriUnion::iterator it1 = _possibleActions.find(*itTri2);
- if (it1 == _possibleActions.end()) {
- Msg::Error("[Rec2D] error no triangle !, stopping with partially filled set");
- Msg::Error(" ");
- }
- int numPossibleRecombinations = 0;
- setTriUnion::iterator it2 = it1->second.begin();
- for (; it2 != it1->second.end(); ++it2)
- if (setIncompTu.find(*it2) == setIncompTu.end())
- ++numPossibleRecombinations;
- if (!numPossibleRecombinations) {
- ++numIsolatedTri;
- setTri.erase(*itTri2);
- }
- if (numPossibleRecombinations == 1)
- setTri.insert(*itTri2);
- }
- }
- itTri = setTri.begin();
- }
-
- double oldReturn = TrianglesUnion::getActualReturn();
- double newReturn = TrianglesUnion::computeReturn(newNum, newVal, newMapVert);
-
- //Msg::Info("after (%d) : %lf -> %lf - %d", numRecomb, TrianglesUnion::getActualReturn(), (newReturn-oldReturn) * size/numRecomb + oldReturn, numIsolatedTri);
-
- return (newReturn-oldReturn) * size/numRecomb + oldReturn - numIsolatedTri;
-}
-
-int Recombine2D::_checkIsolatedTri(const 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));
- }
-
- seqTri.clear();
-
- 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;
- if (numPossibleRecombinations == 1)
- seqTri.insert(*itTri);
- }
- return num;
-}
-
-int Recombine2D::apply(bool a)
-{
- //Msg::Info("(%d, %d, %d)", _quads.size(), _gf->triangles.size(), _isolated.size());
- if (a) {
- _gf->triangles.clear();
- //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();
- }
-
- _applied = true;
- return 1;
-}
-
-void Recombine2D::_removeImpossible(TrianglesUnion *tu)
-{
- 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
- }
-
- 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));
- }
-
- /*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)
- it = _setQuads.erase(it);
- else
- ++it;
- }*/
-}
-
-Rec2d_vertex::Rec2d_vertex(MVertex *v,
- std::list<MTriangle*> &setTri,
- std::map<MVertex*, std::set<GEdge*> > &mapVert, bool a)
-: _numEl(setTri.size()), _angle(.0)
-{
- _onWhat = v->onWhat()->dim() - 1;
- if (a) _onWhat = 1;
-
- if (_onWhat < 0) {
- std::map<MVertex*, std::set<GEdge*> >::iterator it = mapVert.find(v);
- if (it != mapVert.end()) {
- _angle = _computeAngle(v, setTri, it->second);
- }
- else {
- Msg::Warning("[meshRecombine2D] Can't compute corner angle, setting angle to pi/2");
- _angle = M_PI / 2.;
- }
- }
-
- if (_Vvalues == NULL)
- _initStaticTable();
-}
-
-void Rec2d_vertex::_initStaticTable()
-{
- if (_Vvalues == NULL) {
- // _Vvalues[onWhat][nEl]; onWhat={0:edge,1:face} / nEl=_numEl-1
- // _Vbenefs[onWhat][nEl]; onWhat={0:edge,1:face} / nEl=_numEl-1 (benef of adding an element)
- int nE = 5, nF = 10;
-
- _Vvalues = new double*[2];
- _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];
- for (int i = 0; i < nF; i++)
- _Vvalues[1][i] = 1. - fabs(4. / (i+1) - 1.);
-
- _Vbenefs = new double*[2];
- _Vbenefs[0] = new double[nE-1];
- for (int i = 0; i < nE-1; i++)
- _Vbenefs[0][i] = _Vvalues[0][i+1] - _Vvalues[0][i];
- _Vbenefs[1] = new double[nF-1];
- for (int i = 0; i < nF-1; i++)
- _Vbenefs[1][i] = _Vvalues[1][i+1] - _Vvalues[1][i];
- }
-}
-
-
-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.;
- }
- const double prec = 100.;
-
- SVector3 vectv = SVector3(v->x(), v->y(), v->z());
- SVector3 firstDer0, firstDer1;
-
- std::set<GEdge*>::iterator it = setEdge.begin();
- for (int k = 0; k < 2; it++, k++) {
- GEdge *ge = *it;
- SVector3 vectlb = ge->position(ge->getLowerBound());
- SVector3 vectub = ge->position(ge->getUpperBound());
- vectlb -= vectv;
- vectub -= vectv;
- double normlb, normub;
- if ((normlb = norm(vectlb)) > prec * (normub = norm(vectub)))
- firstDer1 = ge->firstDer(ge->getUpperBound());
- else if (normub > prec * normlb)
- firstDer1 = ge->firstDer(ge->getLowerBound());
- else {
- Msg::Warning("[meshRecombine2D] Bad precision, returning pi/2");
- return M_PI / 2.;
- }
- if (k == 0)
- firstDer0 = firstDer1;
- }
-
- firstDer0 = firstDer0 * (1 / norm(firstDer0));
- firstDer1 = firstDer1 * (1 / norm(firstDer1));
-
- double angle1 = acos(dot(firstDer0, firstDer1));
- double angle2 = 2. * M_PI - angle1;
-
- double angleMesh = .0;
- std::list<MTriangle*>::iterator it2 = setTri.begin();
- for (; it2 != setTri.end(); it2++) {
- MTriangle *t = *it2;
- int k = 0;
- while (t->getVertex(k) != v && k < 3)
- k++;
- if (k == 3) {
- Msg::Warning("[meshRecombine2D] Wrong triangle, returning pi/2");
- return M_PI / 2.;
- }
- angleMesh += angle3Vertices(t->getVertex((k+1)%3), v, t->getVertex((k+2)%3));
- }
-
- if (angleMesh < M_PI)
- return angle1;
- return angle2;
-}
-
-double Rec2d_vertex::getReward()
-{
- if (_onWhat > -1)
- return _Vvalues[_onWhat][_numEl-1];
- else
- return 1. - fabs(2. / M_PI * _angle / _numEl - 1.);
-}
-
-double Rec2d_vertex::getReward(int n)
-{
- static double t= Cpu();
- if (_onWhat > -1) {
- 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.);
-}
-
-
-TrianglesUnion::TrianglesUnion(GFace *gf,
- std::list<MTriangle*> &t,
- std::list<MEdge> &e,
- std::list<Rec2d_vertex*> &v,
- std::map<MVertex*,double> &v2lcUV,
- std::map<MEdge, std::list<MTriangle*>, Less_Edge> &mapEdge)
-{
- _numTri = t.size();
- _numEmbEdge = e.size();
- _numBoundVert = v.size() == 2 ? 2 : 4;
- _numEmbVert = v.size() - _numBoundVert;
- _globValIfExecuted = _embEdgeValue = _embVertValue = .0;
-
- _triangles = new MTriangle*[_numTri];
- std::list<MTriangle*>::iterator itt = t.begin();
- std::set<MEdge, Less_Edge> setEdge;
- for (int k = 0; itt != t.end(); itt++, k++) {
- _triangles[k] = *itt;
- for (int i = 0; i < 3; ++i) {
- setEdge.insert((*itt)->getEdge(i));
- }
- }
-
- std::list<MEdge>::iterator ite = e.begin();
- for (; ite != e.end(); ite++) {
- double sumlc = .0, u[2], v[2];
- MVertex *vert[2];
- for (int i = 0; i < 2; 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 = v2lcUV.find(vert[i]);
- if (itlc == v2lcUV.end()) {
- 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 = .2; // FIXME BGM_MeshSize returns wrong meshsize
-
- double length = _computeEdgeLength(gf, vert, u, v, 0);
- double a = .0;
- _embEdgeValue += length / sumlc * (a + (2-a)*_computeAlignment(*ite, t, mapEdge));
- setEdge.erase(*ite);
- }
-
- std::set<MEdge, Less_Edge>::iterator ite2 = setEdge.begin();
- for (; ite2 != setEdge.end(); ite2++) {
- double sumlc = .0, u[2], v[2];
- MVertex *vert[2];
- for (int i = 0; i < 2; i++) {
- vert[i] = ite2->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 = v2lcUV.find(vert[i]);
- if (itlc == v2lcUV.end()) {
- 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 = .2; // FIXME BGM_MeshSize returns wrong meshsize
-
- double length = _computeEdgeLength(gf, vert, u, v, 0);
- double a = .0;
- _boundEdgeValue += length / sumlc * (a + (2-a)*_computeAlignment(*ite2, t, mapEdge));
- _boundEdgeValue = 1;
- }
- _boundEdgeValue /= 2.;
- _numboundEdge = 2;
-
- _vertices = new Rec2d_vertex*[_numBoundVert];
- std::list<Rec2d_vertex*>::iterator itv = v.begin();
- for (int k = 0; itv != v.end() && k < _numBoundVert; itv++, k++)
- _vertices[k] = *itv;
- for (_numEmbVert = 0; itv != v.end(); itv++, _numEmbVert++)
- _embVertValue += (*itv)->getReward();
-
- _computation = _RECOMPUTE - 1;
-}
-
-
-double TrianglesUnion::_computeEdgeLength(GFace *gf, MVertex **vert,
- double *u, double *v, int n)
-{
- double dx, dy, dz, l = .0;
- if (n == 0) {
- dx = vert[1]->x() - vert[0]->x();
- dy = vert[1]->y() - vert[0]->y();
- dz = vert[1]->z() - vert[0]->z();
- l = sqrt(dx * dx + dy * dy + dz * dz);
- }
- else if (n == 1) {
- double edgeCenter[2] ={(u[0] + u[1]) * .5, (v[0] + v[1]) * .5};
- GPoint GP = gf->point (edgeCenter[0], edgeCenter[1]);
- dx = vert[1]->x() - GP.x();
- dy = vert[1]->y() - GP.y();
- dz = vert[1]->z() - GP.z();
- l += sqrt(dx * dx + dy * dy + dz * dz);
- dx = vert[0]->x() - GP.x();
- dy = vert[0]->y() - GP.y();
- dz = vert[0]->z() - GP.z();
- l += sqrt(dx * dx + dy * dy + dz * dz);
- }
- else {
- Msg::Warning("[meshRecombine2D] edge length computation not implemented for n=%d, returning 0", n);
- return .0;
- }
- return l;
-}
-double TrianglesUnion::_computeAlignment(const MEdge &e, std::list<MTriangle*> <,
- std::map<MEdge, std::list<MTriangle*>, Less_Edge> &mapEdge)
-{
- std::list<MTriangle*>::iterator itlt = lt.begin();
- //if (lt.size() == 2)
- // return Temporary::compute_alignment(e, *itlt, *(++itlt));
-
- MVertex *v0 = e.getVertex(0);
- MVertex *v1 = e.getVertex(1);
- MTriangle *tris[2];
- int k = 0;
- for (; itlt != lt.end(); itlt++) {
- MTriangle *t = *itlt;
- if ((t->getVertex(0) == v0 || t->getVertex(1) == v0 || t->getVertex(2) == v0) &&
- (t->getVertex(0) == v1 || t->getVertex(1) == v1 || t->getVertex(2) == v1) )
- tris[k++] = t;
- }
- if (k < 2) {
- k = 0;
- std::map<MEdge, std::list<MTriangle*> >::iterator itmEdge = mapEdge.find(e);
- if (itmEdge == mapEdge.end()) {
- Msg::Info("g po :( (szie %d)", mapEdge.size());
- return .0;
- }
- itlt = itmEdge->second.begin();
- int num = 0;
- for (; itlt != itmEdge->second.end(); itlt++) {
- MTriangle *t = *itlt;
- if ((t->getVertex(0) == v0 || t->getVertex(1) == v0 || t->getVertex(2) == v0) &&
- (t->getVertex(0) == v1 || t->getVertex(1) == v1 || t->getVertex(2) == v1) )
- tris[k++] = t;
- }
- }
- if (k < 2 || k > 2) {
- //Msg::Warning("[meshRecombine2D] error in alignment computation, returning 1");
- return 1.;
- }
- return Temporary::compute_alignment(e, tris[0], tris[1]);
-}
-void TrianglesUnion::_update()
-{
- if (_computation >= _RECOMPUTE)
- return;
- double alpha = _ValVert - _embVertValue;
- for (int i = 0; i < _numBoundVert; i++) {
- alpha += _vertices[i]->getReward(-1);
- }
- alpha /= _NumVert - _numEmbVert;
- double beta = (_ValEdge - _embEdgeValue + _boundEdgeValue) / (_NumEdge - _numEmbEdge + _numboundEdge);
- _globValIfExecuted = alpha * alpha * beta;
-
- _computation = _RECOMPUTE;
-}
-
-MQuadrangle *TrianglesUnion::createQuad() const
-{
- std::list<MEdge> listBoundEdge;
- {
- std::multiset<MEdge, Less_Edge> msetEdge;
- for (int i = 0; i < _numTri; i++) {
- MTriangle *t = _triangles[i];
- for (int i = 0; i < 3; i++) {
- msetEdge.insert(t->getEdge(i));
- }
- }
- std::multiset<MEdge>::iterator itEdge = msetEdge.begin();
- MEdge me = *(itEdge++);
- bool precWasSame = false;
- for (; itEdge != msetEdge.end(); itEdge++) {
- if (*itEdge == me)
- precWasSame = true;
- else if (precWasSame)
- precWasSame = false;
- else
- listBoundEdge.push_back(me);
- me = *itEdge;
- }
- if (!precWasSame)
- listBoundEdge.push_back(me);
- }
- if (listBoundEdge.size() != 4)
- Msg::Warning("[meshRecombine2D] Not 4 edges !");
-
- MVertex *v1, *v2, *v3, *v4;
- v1 = listBoundEdge.begin()->getMinVertex();
- v2 = listBoundEdge.begin()->getMaxVertex();
- std::list<MEdge>::iterator it = listBoundEdge.begin();
- for (it++; it != listBoundEdge.end(); it++) {
- if (v2 == it->getMinVertex()) {
- v3 = it->getMaxVertex();
- listBoundEdge.erase(it);
- break;
- }
- 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();
- break;
- }
- if (v3 == it->getMaxVertex()) {
- v4 = it->getMinVertex();
- break;
- }
- }
-
- return new MQuadrangle(v1, v2, v3, v4);
-}
-
-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;
- num[1] -= _numboundEdge;
- val[0] += _embVertValue;
- val[1] += _embEdgeValue;
- val[1] -= _boundEdgeValue;
- 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;
- num[1] += _numboundEdge;
- val[0] -= _embVertValue;
- val[1] -= _embEdgeValue;
- val[1] += _boundEdgeValue;
- 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 * alpha * beta;
-}
-
-bool TrianglesUnion::operator<(TrianglesUnion &other)
-{
- _update();
- other._update();
- return _globValIfExecuted < other._globValIfExecuted;
-}
-
-bool lessTriUnion::operator()(TrianglesUnion *tu1, TrianglesUnion *tu2) const
-{
- return *tu1 < *tu2;
-}
-
-bool lessTriUnionInv::operator()(TrianglesUnion *tu1, TrianglesUnion *tu2) const
-{
- return *tu2 < *tu1;
-}
-
-/*map tri to recomb
-//map edge to double value (constructor)
-set Recomb sorted
-function best(copy vertex*[], copy set Recomb sorted);
-
-construction :
- - list of vertex
- - set of recomb
- - map tri to recomb
-
-destruction :
- - vertices, recomb
-
-execution :
- - take best recomb
- - determine recomb to execute
- - maj E N e n
- - reduce maptrirecomb, setrecomb
- - sort setrecomb*/