Select Git revision
meshRecombine2D_2.cpp
Forked from
gmsh / gmsh
8941 commits behind the upstream repository.
Christophe Geuzaine authored
meshRecombine2D_2.cpp 45.74 KiB
// 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*/