Skip to content
Snippets Groups Projects
Commit 7742ab29 authored by Amaury Johnen's avatar Amaury Johnen
Browse files

remove old code for 2D recombinations

parent 382d7d40
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Mesh/meshRecombine2D.cpp deleted 100644 → 0
+ 0
− 490
View file @ 382d7d40
// 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);
}
}
}
//
Mesh/meshRecombine2D.h deleted 100644 → 0
+ 0
− 420
View file @ 382d7d40
// 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
Mesh/meshRecombine2D_2.cpp deleted 100644 → 0
+ 0
− 1450
View file @ 382d7d40
// 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*> &lt,
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*/
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment