From ea2ec757c7db07f387feeb461bf7a43ccaac8e01 Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Mon, 16 Jan 2012 07:44:08 +0000
Subject: [PATCH]
---
Mesh/meshGFaceRecombine.cpp | 734 ++++++++++++++++++++++++++++++++++--
Mesh/meshGFaceRecombine.h | 132 ++++++-
2 files changed, 814 insertions(+), 52 deletions(-)
diff --git a/Mesh/meshGFaceRecombine.cpp b/Mesh/meshGFaceRecombine.cpp
index 53abddba3d..af69ab0b1a 100644
--- a/Mesh/meshGFaceRecombine.cpp
+++ b/Mesh/meshGFaceRecombine.cpp
@@ -7,21 +7,49 @@
// Amaury Johnen (a.johnen@ulg.ac.be)
//
-// TODO :
-// - recombine w. take care of parity, update of quality...
-
#include "meshGFaceRecombine.h"
#include "MElement.h"
#include "MTriangle.h"
+#include "MQuadrangle.h"
+#include "meshGFaceOptimize.h"
+
+#include "drawContext.h"
+#include "FlGui.h"
+#include "Context.h"
+#include "OS.h"
-#define EDGE_BASE 2
-#define EDGE_QUAD 1
+#define REC2D_EDGE_BASE 1
+#define REC2D_EDGE_QUAD 1
+#define REC2D_WAIT_TIME 1
+
+// #define REC2D_SMOOTH
+ #define REC2D_DRAW
Recombine2D *Recombine2D::_current = NULL;
bool Recombine2D::bo = 0;
double **Rec2DVertex::_qualVSnum = NULL;
double **Rec2DVertex::_gains = NULL;
+void Recombine2D::associateParity(int pOld, int pNew)
+{
+ std::set<Rec2DVertex*>::iterator it = _current->_parities[pOld].begin();
+ for (; it != _current->_parities[pOld].end(); ++it)
+ (*it)->setParity(pNew);
+ _current->_parities.erase(pOld);
+ if (pOld % 2)
+ pOld -= 1;
+ else
+ pOld +=1;
+ if (pNew % 2)
+ pNew -= 1;
+ else
+ pNew +=1;
+ it = _current->_parities[pOld].begin();
+ for (; it != _current->_parities[pOld].end(); ++it)
+ (*it)->setParity(pNew);
+ _current->_parities.erase(pOld);
+}
+
/** Recombine2D **/
/*******************/
Recombine2D::Recombine2D(GFace *gf) : _gf(gf)
@@ -38,6 +66,8 @@ Recombine2D::Recombine2D(GFace *gf) : _gf(gf)
_numEdge = _numVert = 0;
_valEdge = _valVert = .0;
+ _tri = _gf->triangles;
+
// Be able to compute geometrical angle at corners
std::map<MVertex*, std::set<GEdge*> > mapCornerVert;
{
@@ -91,11 +121,11 @@ Recombine2D::Recombine2D(GFace *gf) : _gf(gf)
itE2E = _edges.insert(itE2E, std::make_pair(e,re));
_vertices[e.getVertex(0)]->add(re);
_vertices[e.getVertex(1)]->add(re);
- _numEdge += EDGE_BASE;
- _valEdge += EDGE_BASE * re->getQual();
+ _numEdge += REC2D_EDGE_BASE;
+ _valEdge += REC2D_EDGE_BASE * re->getQual();
if (itedge->second.size() == 1) {
- _numEdge += EDGE_QUAD;
- _valEdge += EDGE_QUAD * re->getQual();
+ _numEdge += REC2D_EDGE_QUAD;
+ _valEdge += REC2D_EDGE_QUAD * re->getQual();
_vertices[e.getVertex(0)]->setOnBoundary();
_vertices[e.getVertex(1)]->setOnBoundary();
boundV2E[e.getVertex(0)].push_back(e);
@@ -145,7 +175,7 @@ Recombine2D::Recombine2D(GFace *gf) : _gf(gf)
boundV2E.erase(itBound);
}
if (_vertices[actualV]->getParity() != actualParity+p)
- Msg::Error("Not even number of mesh edges on boundary !");
+ Msg::Error("[Recombine2D] Not even number of mesh edges on boundary !");
boundV2E.erase(boundV2E.begin());
itBound = boundV2E.begin();
@@ -159,19 +189,68 @@ Recombine2D::Recombine2D(GFace *gf) : _gf(gf)
std::map<MEdge, std::list<MTriangle*> >::iterator ite = mapEdge.begin();
for (; ite != mapEdge.end(); ++ite)
if (ite->second.size() == 2) {
- _actions.insert(new Rec2DTwoTri2Quad(_edges[ite->first], ite->second));
- _actions.insert(new Rec2DCollapseTri(_edges[ite->first], ite->second));
+ Rec2DTwoTri2Quad *q = new Rec2DTwoTri2Quad(_edges[ite->first], ite->second);
+ //Rec2DCollapseTri *c = new Rec2DCollapseTri(_edges[ite->first], ite->second);
+ _actions.push_back(q);
+ //_actions.insert(c);
+ std::list<MTriangle*>::iterator it = ite->second.begin();
+ _mea[*it].insert(q);
+ //_mea[*it].insert(c);
+ ++it;
+ _mea[*it].insert(q);
+ //_mea[*it].insert(c);
}
}
+ // Rec2DTwoTri2Quad
+ /*{
+ std::list<std::map<MVertex*, std::list<MTriangle*> >::iterator>::iterator it4;
+ for (it4 = fourTri.begin(); it4 != fourTri.end(); it4++) {
+ Rec2DVertex *rv = _vertices[(*it4)->first];
+ std::set<MEdge, Less_Edge> setEdge;
+ {
+ std::list<MTriangle*>::iterator it = (*it4)->second.begin();
+ for (; it != (*it4)->second.end(); ++it) {
+ setEdge.insert((*it)->getEdge(0));
+ setEdge.insert((*it)->getEdge(1));
+ setEdge.insert((*it)->getEdge(2));
+ }
+ }
+
+ std::set<MEdge>::iterator it = setEdge.begin();
+ for (; it != setEdge.end(); ++it) {
+ std::map<MEdge, std::list<MTriangle*> >::iterator it2;
+ it2 = mapEdge.find(*it);
+ if (it2 != mapEdge.end())
+ if (it2->second.size() == 2) {
+ Rec2DTwoTri2Quad *q = new Rec2DTwoTri2Quad(_edges[it2->first], it2->second);
+ _actions.insert(q);
+ std::list<MTriangle*>::iterator it = it2->second.begin();
+ _mea[*it].insert(q);
+ ++it;
+ _mea[*it].insert(q);
+ }
+ else
+ Msg::Error("fjdaks;lfjkdsa;");
+ }
+ }
+ }*/
+
// Create the 'Rec2DFourTri2Quad' and 'Rec2DCollapseTri'
{
std::list<std::map<MVertex*, std::list<MTriangle*> >::iterator>::iterator it4;
for (it4 = fourTri.begin(); it4 != fourTri.end(); it4++) {
Rec2DVertex *rv = _vertices[(*it4)->first];
if (!rv->getIsOnBoundary()) {
- _actions.insert(new Rec2DFourTri2Quad(rv, (*it4)->second));
- _actions.insert(new Rec2DCollapseTri(rv, (*it4)->second));
+ Rec2DFourTri2Quad *q = new Rec2DFourTri2Quad(rv, (*it4)->second);
+ //Rec2DCollapseTri *c = new Rec2DCollapseTri(rv, (*it4)->second);
+ _actions.push_back(q);
+ //_actions.insert(c);
+ std::list<MTriangle*>::iterator it = (*it4)->second.begin();
+ for (; it != (*it4)->second.end(); ++it) {
+ _mea[*it].insert(q);
+ //_mea[*it].insert(c);
+ }
}
}
}
@@ -202,11 +281,340 @@ Recombine2D::~Recombine2D()
//delete vertices, edges, actions;
}
+bool Recombine2D::recombine()
+{
+ Rec2DAction *nextAction;
+ while (_actions.size() > 0) {
+ double t = Cpu();
+ nextAction = *std::max_element(_actions.begin(), _actions.end(),
+ lessRec2DAction() );
+ if (nextAction->getReward() > 0) {
+ if (_isAllQuad(nextAction))
+ _choose(nextAction);
+ else {
+ nextAction->color(190, 0, 0);
+ _removeReference(nextAction);
+ }
+ }
+ else {
+ nextAction->color(0, 100, 0);
+ _removeReference(nextAction);
+ }
+
+ _gf->triangles = _tri;
+ _gf->quadrangles = _quad;
+#ifdef REC2D_DRAW
+ CTX::instance()->mesh.changed = (ENT_ALL);
+ FlGui::instance()->check();
+ drawContext::global()->draw();
+ while (Cpu()-t < REC2D_WAIT_TIME)
+ FlGui::instance()->check();
+#endif
+#ifdef REC2D_SMOOTH
+ laplaceSmoothing(_gf,1);
+#endif
+
+ std::map<int, std::set<Rec2DVertex*> >::iterator it = _parities.begin();
+ for (; it != _parities.end(); ++it) {
+ Msg::Info("par %d, #%d", it->first, it->second.size());
+ }
+ }
+
+ Msg::Info("State");
+ Msg::Info("-----");
+ Msg::Info("numEdge %d (%d), valEdge %g", _numEdge, _edges.size(), _valEdge);
+ Msg::Info("numVert %d (%d), valVert %g", _numVert, _vertices.size(), _valVert);
+ Msg::Info("global Value %g", Recombine2D::getGlobalValue(0,0,0,0));
+ Msg::Info("num action %d", _actions.size());
+ std::map<int, std::set<Rec2DVertex*> >::iterator it = _parities.begin();
+ for (; it != _parities.end(); ++it) {
+ Msg::Info("par %d, #%d", it->first, it->second.size());
+ }
+ mapofElementActions::iterator it2 = _mea.begin();
+ int a = 0;
+ for (; it2 != _mea.end(); ++it2) {
+ a += it2->second.size();
+ }
+ if (a > 0)
+ Msg::Info("%d elements in mapof elem to action, with average action %g", _mea.size(), (double)a/_mea.size());
+ Msg::Info(" ");
+
+ return 1;
+}
+
+bool Recombine2D::_isAllQuad(Rec2DAction *action)
+{
+ int p[4];
+ action->getParities(p);
+
+ if (p[0] < 0 && p[1] < 0 && p[2] < 0 && p[3] < 0)
+ return true;
+
+ if (action->isObsolete()) {
+ Msg::Error("[Recombine2D] obsolete action");
+ return false;
+ }
+
+ int min = 100, index;
+ for (int i = 0; i < 4; ++i) {
+ if (p[i] > -1 && min > p[i]) {
+ min = p[i];
+ index = i;
+ }
+ }
+
+ //Quid si p[0..4] = -1 ?
+
+ std::map<Rec2DVertex*, int> newParVert;
+ std::map<int, int> parAssoc; // first parity value become second parity value
+
+ for (int i = 0; i < 4; i += 2) {
+ int par;
+ if ((index/2) * 2 == i)
+ par = min;
+ else if (min % 2)
+ par = min - 1;
+ else
+ par = min + 1;
+ for (int j = 0; j < 2; ++j) {
+ if (p[i+j] == -1)
+ newParVert.insert(std::make_pair(action->getRVertex(i+j), par));
+ else if (p[i+j] != par)
+ parAssoc.insert(std::make_pair(p[i+j], par));
+ }
+ }
+
+ std::set<Rec2DAction*> neighbours;
+ _getNeighbours(action, neighbours);
+
+ while (neighbours.size() > 0) {
+ std::set<Rec2DAction*>::iterator it = neighbours.begin();
+ if ((*it)->getNum() != 1) {
+ neighbours.erase(it);
+ continue;
+ }
+ Rec2DVertex *v0 = (*it)->getRVertex(0);
+ Rec2DVertex *v1 = (*it)->getRVertex(1);
+ std::map<Rec2DVertex*, int>::iterator it2;
+ int p0 = v0->getParity(), p1 = v1->getParity();
+ if (p0 < 0 && (it2 = newParVert.find(v0)) == newParVert.end()) {
+ neighbours.erase(it);
+ continue;
+ }
+ if (p0 < 0)
+ p0 = it2->second;
+
+ if (p1 < 0 && (it2 = newParVert.find(v1)) == newParVert.end()) {
+ neighbours.erase(it);
+ continue;
+ }
+ if (p1 < 0)
+ p1 = it2->second;
+
+ std::map<int, int>::iterator it3;
+ if ((it3 = parAssoc.find(p0)) != parAssoc.end())
+ p0 = it3->second;
+ if ((it3 = parAssoc.find(p1)) != parAssoc.end())
+ p1 = it3->second;
+
+ if (p0 != p1) {
+ neighbours.erase(it);
+ continue;
+ }
+
+ int par;
+ bool hasChange = false;
+ if (p0 % 2)
+ par = p0 - 1;
+ else
+ par = p0 + 1;
+ (*it)->getParities(p);
+ for (int j = 2; j < 4; ++j) {
+ if (p[j] == -1) {
+ if ((it2 = newParVert.find((*it)->getRVertex(j))) == newParVert.end()) {
+ newParVert.insert(std::make_pair(action->getRVertex(j), par));
+ hasChange = true;
+ }
+ else if (it2->second != par)
+ return false;
+ }
+ else if (p[j] != par) {
+ if ((it3 = parAssoc.find(p[j])) == parAssoc.end()) {
+ parAssoc.insert(std::make_pair(p[j], par));
+ hasChange = true;
+ }
+ else {
+ int a = it3->second;
+ while (a != par && (it3 = parAssoc.find(a)) != parAssoc.end()) {
+ a = it3->second;
+ }
+ if (a != par)
+ return false;
+ }
+ }
+ }
+ Rec2DAction *ra = *it;
+ neighbours.erase(it);
+ if (hasChange)
+ _getNeighbours(ra, neighbours);
+ }
+ return true;
+}
+
+void Recombine2D::_choose(Rec2DAction *action)
+{
+ action->apply();
+ _removeIncompatible(action);
+ ++_numChange;
+}
+
+void Recombine2D::_removeIncompatible(Rec2DAction *action)
+{
+ std::set<MTriangle*> touchedTri;
+ action->getTriangles(touchedTri);
+ std::set<MTriangle*>::iterator itTouched;
+
+ std::set<Rec2DAction*> incomp;
+ {
+ itTouched = touchedTri.begin();
+ for (; itTouched != touchedTri.end(); ++itTouched) {
+ mapofElementActions::iterator ite2a = _mea.find(*itTouched);
+ if (ite2a != _mea.end()) {
+ std::set<Rec2DAction*>::iterator ita = ite2a->second.begin();
+ for (; ita != ite2a->second.end(); ++ita)
+ incomp.insert(*ita);
+ }
+ else
+ Msg::Error("[Recombine2D] didn't found triangle %d (rmIncomp1)", (*itTouched)->getNum());
+ }
+ }
+
+ std::set<Rec2DAction*>::iterator itIncomp = incomp.begin();
+ for (; itIncomp != incomp.end(); ++itIncomp) {
+ std::set<MTriangle*> touched;
+ (*itIncomp)->getTriangles(touched);
+ for (itTouched = touched.begin(); itTouched != touched.end(); ++itTouched) {
+ mapofElementActions::iterator ite2a = _mea.find(*itTouched);
+ if (ite2a != _mea.end()) {
+ bool b = ite2a->second.erase(*itIncomp);
+ if (!b)
+ Msg::Error("[Recombine2D] Tri->Action not found (rmIncomp)");
+ }
+ else
+ Msg::Error("[Recombine2D] didn't found triangle %d (rmIncomp)", (*itTouched)->getNum());
+ }
+ //bool b = _actions.erase(*itIncomp);
+ setofRec2DAction::iterator it2 = _actions.begin();
+ bool b = false;
+ while (it2 != _actions.end()) {
+ if (*it2 == *itIncomp) {
+ _actions.erase(it2);
+ b = true;
+ break;
+ }
+ ++it2;
+ }
+ if (!b)
+ Msg::Error("[Recombine2D] action not found (rmIncomp2)");
+ delete *itIncomp;
+ }
+
+ itTouched = touchedTri.begin();
+ for (; itTouched != touchedTri.end(); ++itTouched) {
+ bool b = _mea.erase(*itTouched);
+ if (!b)
+ Msg::Error("[Recombine2D] triangle not found %d", (*itTouched)->getNum());
+ }
+}
+
+void Recombine2D::_removeReference(Rec2DAction *action)
+{
+ std::set<MTriangle*> touchedTri;
+ action->getTriangles(touchedTri);
+ std::set<MTriangle*>::iterator it = touchedTri.begin();
+ for (; it != touchedTri.end(); ++it) {
+ mapofElementActions::iterator ite2a = _mea.find(*it);
+ if (ite2a != _mea.end()) {
+ bool b = ite2a->second.erase(action);
+ if (!b)
+ Msg::Error("[Recombine2D] Tri->Action not found (rmRef)");
+ }
+ else
+ Msg::Error("[Recombine2D] didn't found triangle %d (rmRef)", (*it)->getNum());
+ }
+ //bool b = _actions.erase(action);
+ setofRec2DAction::iterator it2 = _actions.begin();
+ bool b = false;
+ while (it2 != _actions.end()) {
+ if (*it2 == action) {
+ _actions.erase(it2);
+ b = true;
+ break;
+ }
+ ++it2;
+ }
+ if (!b)
+ Msg::Error("[Recombine2D] action not found (rmRef)");
+ delete action;
+}
+
+void Recombine2D::_getNeighbours(Rec2DAction *ra, std::set<Rec2DAction*> &actions)
+{
+ std::set<MTriangle*> touchedTri;
+ ra->getTriangles(touchedTri);
+ std::set<MTriangle*>::iterator itTouched;
+
+ std::set<Rec2DAction*> incomp;
+ {
+ itTouched = touchedTri.begin();
+ for (; itTouched != touchedTri.end(); ++itTouched) {
+ mapofElementActions::iterator ite2a = _mea.find(*itTouched);
+ if (ite2a != _mea.end()) {
+ std::set<Rec2DAction*>::iterator ita = ite2a->second.begin();
+ for (; ita != ite2a->second.end(); ++ita)
+ incomp.insert(*ita);
+ }
+ else
+ Msg::Error("[Recombine2D.] didn't found triangle %d (rmIncomp1)", (*itTouched)->getNum());
+ }
+ }
+
+ std::set<Rec2DAction*>::iterator it = incomp.begin();
+ for (; it != incomp.end(); ++it) {
+ (*it)->getTriangles(touchedTri);
+ }
+
+ std::set<Rec2DAction*> neighbours;
+ {
+ itTouched = touchedTri.begin();
+ for (; itTouched != touchedTri.end(); ++itTouched) {
+ mapofElementActions::iterator ite2a = _mea.find(*itTouched);
+ if (ite2a != _mea.end()) {
+ std::set<Rec2DAction*>::iterator ita = ite2a->second.begin();
+ for (; ita != ite2a->second.end(); ++ita)
+ neighbours.insert(*ita);
+ }
+ else
+ Msg::Error("[Recombine2D.] didn't found triangle %d (rmIncomp1)", (*itTouched)->getNum());
+ }
+ std::set<Rec2DAction*>::iterator itIncomp = incomp.begin();
+ for (; itIncomp != incomp.end(); ++itIncomp) {
+ bool b = neighbours.erase(*itIncomp);
+ if (!b)
+ Msg::Error("[Recombine2D] didn't found incomp action");
+ }
+ }
+
+ for (it = neighbours.begin(); it != neighbours.end(); ++it) {
+ actions.insert(*it);
+ }
+}
+
double Recombine2D::getGlobalValue(int numEdge, double valEdge,
int numVert, double valVert )
{
- return (_current->_valEdge + valEdge) / (_current->_numEdge + numEdge) *
- (_current->_valVert + valVert) / (_current->_numVert + numVert);
+ double a = (_current->_valVert + valVert) / (_current->_numVert + numVert);
+ return (_current->_valEdge + valEdge) / (_current->_numEdge + numEdge) * a * a;
}
Rec2DEdge* Recombine2D::getREdge(MEdge e)
@@ -227,6 +635,85 @@ Rec2DVertex* Recombine2D::getRVertex(MVertex *v)
return NULL;
}
+void Recombine2D::remove(MEdge e)
+{
+ mapofEdges::iterator it = _current->_edges.find(e);
+ if (it != _current->_edges.end())
+ _current->_edges.erase(it);
+ else
+ Msg::Error("[Recombine2D::getREdge] edge not found");
+}
+
+void Recombine2D::remove(MVertex *v)
+{
+ mapofVertices::iterator it = _current->_vertices.find(v);
+ if (it != _current->_vertices.end())
+ _current->_vertices.erase(it);
+ else
+ Msg::Error("[Recombine2D::getRVertex] vertex not found");
+}
+
+void Recombine2D::remove(MTriangle *tri)
+{
+ std::vector<MTriangle*>::iterator it = _current->_tri.begin();
+ for (; it != _current->_tri.end(); ++it) {
+ if (*it == tri) {
+ _current->_tri.erase(it);
+ return;
+ }
+ }
+ Msg::Error("[Recombine2D] triangle %d was not there", tri->getNum());
+}
+
+void Recombine2D::add(MQuadrangle *quad)
+{
+ _current->_quad.push_back(quad);
+}
+
+
+/** Rec2DAction **/
+/*******************/
+bool Rec2DAction::operator<(Rec2DAction &other)
+{
+ return getReward() < other.getReward();
+}
+
+MQuadrangle* Rec2DAction::_createQuad(std::vector<Rec2DEdge*> &boundEdge) const
+{
+ if (boundEdge.size() != 4) {
+ Msg::Error("[Rec2DAction] Not 4 edges for quad creation");
+ return new MQuadrangle(NULL, NULL, NULL, NULL);
+ }
+ MVertex *v1, *v2, *v3, *v4;
+ v1 = boundEdge[0]->getRVertex(0)->getMVertex();
+ v2 = boundEdge[0]->getRVertex(1)->getMVertex();
+ int I;
+ for (unsigned int i = 1; i < 4; ++i) {
+ if (v2 == boundEdge[i]->getRVertex(0)->getMVertex()) {
+ v3 = boundEdge[i]->getRVertex(1)->getMVertex();
+ I = i;
+ break;
+ }
+ if (v2 == boundEdge[i]->getRVertex(1)->getMVertex()) {
+ v3 = boundEdge[i]->getRVertex(0)->getMVertex();
+ I = i;
+ break;
+ }
+ }
+ for (unsigned int i = 1; i < 4; ++i) {
+ if (i == I) continue;
+ if (v3 == boundEdge[i]->getRVertex(0)->getMVertex()) {
+ v4 = boundEdge[i]->getRVertex(1)->getMVertex();
+ break;
+ }
+ if (v3 == boundEdge[i]->getRVertex(1)->getMVertex()) {
+ v4 = boundEdge[i]->getRVertex(0)->getMVertex();
+ break;
+ }
+ }
+ return new MQuadrangle(v1, v2, v3, v4);
+}
+
/** Rec2DTwoTri2Quad **/
/************************/
@@ -282,19 +769,97 @@ bool Rec2DTwoTri2Quad::isObsolete()
return false;
}
+void Rec2DTwoTri2Quad::apply()
+{
+ if (isObsolete()) {
+ Msg::Error("[Rec2DTwoTri2Quad] Applying obsolet action");
+ return;
+ }
+
+ int min = 100, index = -1;
+ for (int i = 0; i < 4; ++i) {
+ if (_vertices[0]->getParity() > -1 && min > _vertices[0]->getParity()) {
+ min = _vertices[0]->getParity();
+ index = i;
+ }
+ }
+ if (index == -1) {
+ int par = Recombine2D::highPar();
+ _vertices[0]->setParity(par);
+ _vertices[1]->setParity(par);
+ _vertices[2]->setParity(par+1);
+ _vertices[3]->setParity(par+1);
+ }
+ else {
+ for (int i = 0; i < 4; i += 2) {
+ int par;
+ if ((index/2) * 2 == i)
+ par = min;
+ else if (min % 2)
+ par = min - 1;
+ else
+ par = min + 1;
+ Msg::Info("%d", par);
+ for (int j = 0; j < 2; ++j) {
+ Msg::Info(" %d %d", _vertices[i+j]->getParity(), par);
+ if (_vertices[i+j]->getParity() == -1)
+ _vertices[i+j]->setParity(par);
+ else if (_vertices[i+j]->getParity() != par)
+ Recombine2D::associateParity(_vertices[0]->getParity(), par);
+ }
+ }
+ }
+
+ double valEdge = - REC2D_EDGE_BASE * _edges[0]->getQual();
+ for (int i = 1; i < 5; ++i)
+ valEdge += REC2D_EDGE_QUAD * _edges[i]->getQual();
+
+ Recombine2D::addNumEdge(4*REC2D_EDGE_QUAD - REC2D_EDGE_BASE);
+ Recombine2D::addValEdge(valEdge);
+ //Recombine2D::addValVert(_vertices[0]->getGain(-1) + _vertices[1]->getGain(-1));
+
+ Recombine2D::remove(_triangles[0]); // FIXME not performant at all !!!
+ Recombine2D::remove(_triangles[1]);
+ std::vector<Rec2DEdge*> edges;
+ for (int i = 1; i < 5; ++i)
+ edges.push_back(_edges[i]);
+ Recombine2D::add(_createQuad(edges));
+
+ _vertices[0]->remove(_edges[0]);
+ _vertices[1]->remove(_edges[0]);
+ Recombine2D::remove(_edges[0]->getMEdge()); // could also delete REdge
+ delete _edges[0];
+}
+
+void Rec2DTwoTri2Quad::getParities(int *par)
+{
+ par[0] = _vertices[0]->getParity();
+ par[1] = _vertices[1]->getParity();
+ par[2] = _vertices[2]->getParity();
+ par[3] = _vertices[3]->getParity();
+}
+
void Rec2DTwoTri2Quad::_computeGlobVal()
{
- double valEdge = - EDGE_BASE * _edges[0]->getQual();
+ double valEdge = - REC2D_EDGE_BASE * _edges[0]->getQual();
for (int i = 1; i < 5; ++i)
- valEdge += EDGE_QUAD * _edges[i]->getQual();
+ valEdge += REC2D_EDGE_QUAD * _edges[i]->getQual();
double valVert = _vertices[0]->getGain(-1) + _vertices[1]->getGain(-1);
- _globValIfExecuted = Recombine2D::getGlobalValue(4*EDGE_QUAD - EDGE_BASE,
+ _globValIfExecuted = Recombine2D::getGlobalValue(4*REC2D_EDGE_QUAD - REC2D_EDGE_BASE,
valEdge, 0, valVert );
_lastUpdate = Recombine2D::getNumChange();
}
+void Rec2DTwoTri2Quad::color(int a, int b, int c)
+{
+ unsigned int col = CTX::instance()->packColor(a, b, c, 255);
+ _triangles[0]->setCol(col);
+ _triangles[1]->setCol(col);
+}
+
+
/** Rec2DFourTri2Quad **/
/*************************/
Rec2DFourTri2Quad::Rec2DFourTri2Quad(Rec2DVertex *rv, std::list<MTriangle*> &tri)
@@ -359,23 +924,77 @@ bool Rec2DFourTri2Quad::isObsolete()
return false;
}
+void Rec2DFourTri2Quad::apply()
+{
+ double valEdge = 0;
+ for (int i = 0; i < 4; ++i)
+ valEdge -= REC2D_EDGE_BASE * _edges[i]->getQual();
+ for (int i = 4; i < 8; ++i)
+ valEdge += REC2D_EDGE_QUAD * _edges[i]->getQual();
+
+ double valVert = - _vertices[4]->getQual();
+ for (int i = 0; i < 4; ++i)
+ valEdge += _vertices[i]->getGain(-1);
+
+ Recombine2D::addNumEdge(4*REC2D_EDGE_QUAD - 4*REC2D_EDGE_BASE);
+ Recombine2D::addNumVert(-1);
+ Recombine2D::addValEdge(valEdge);
+ Recombine2D::addValVert(valVert);
+
+ Recombine2D::remove(_triangles[0]); // FIXME not performant at all !!!
+ Recombine2D::remove(_triangles[1]);
+ Recombine2D::remove(_triangles[2]);
+ Recombine2D::remove(_triangles[3]);
+ std::vector<Rec2DEdge*> edges;
+ for (int i = 4; i < 8; ++i)
+ edges.push_back(_edges[i]);
+ Recombine2D::add(_createQuad(edges));
+
+ for (int i = 0; i < 4; ++i) {
+ _edges[i]->getRVertex(0)->remove(_edges[i]);
+ _edges[i]->getRVertex(1)->remove(_edges[i]);
+ Recombine2D::remove(_edges[i]->getMEdge());
+ delete _edges[i];
+ }
+ Recombine2D::remove(_vertices[4]->getMVertex());
+ delete _vertices[4];
+}
+
+void Rec2DFourTri2Quad::getParities(int *par)
+{
+ par[0] = _vertices[0]->getParity();
+ par[1] = _vertices[1]->getParity();
+ par[2] = _vertices[2]->getParity();
+ par[3] = _vertices[3]->getParity();
+}
+
void Rec2DFourTri2Quad::_computeGlobVal()
{
double valEdge = 0;
for (int i = 0; i < 4; ++i)
- valEdge -= EDGE_BASE * _edges[i]->getQual();
+ valEdge -= REC2D_EDGE_BASE * _edges[i]->getQual();
for (int i = 4; i < 8; ++i)
- valEdge += EDGE_QUAD * _edges[i]->getQual();
+ valEdge += REC2D_EDGE_QUAD * _edges[i]->getQual();
- double valVert = - _vertices[0]->getQual();
- for (int i = 1; i < 5; ++i)
+ double valVert = - _vertices[4]->getQual();
+ for (int i = 0; i < 4; ++i)
valEdge += _vertices[i]->getGain(-1);
- _globValIfExecuted = Recombine2D::getGlobalValue(4*EDGE_QUAD - 4*EDGE_BASE,
+ _globValIfExecuted = Recombine2D::getGlobalValue(4*REC2D_EDGE_QUAD - 4*REC2D_EDGE_BASE,
valEdge, -1, valVert );
_lastUpdate = Recombine2D::getNumChange();
}
+void Rec2DFourTri2Quad::color(int a, int b, int c)
+{
+ unsigned int col = CTX::instance()->packColor(a, b, c, 255);
+ _triangles[0]->setCol(col);
+ _triangles[1]->setCol(col);
+ _triangles[2]->setCol(col);
+ _triangles[3]->setCol(col);
+}
+
+
/** Rec2DCollapseTri **/
/*********************/
Rec2DCollapseTri::Rec2DCollapseTri(Rec2DVertex *rv, std::list<MTriangle*> &tri)
@@ -464,6 +1083,19 @@ bool Rec2DCollapseTri::isObsolete()
return false;
}
+void Rec2DCollapseTri::apply()
+{
+ return;
+}
+
+void Rec2DCollapseTri::getParities(int *par)
+{
+ par[0] = _vertices[0]->getParity();
+ par[1] = _vertices[1]->getParity();
+ par[2] = _vertices[2]->getParity();
+ par[3] = _vertices[3]->getParity();
+}
+
void Rec2DCollapseTri::_computeGlobVal()
{
int b0 = _vertices[0]->getIsOnBoundary();
@@ -602,6 +1234,18 @@ void Rec2DCollapseTri::_qualCavity(double &valVert, double &valEdge,
}
}
+void Rec2DCollapseTri::color(int a, int b, int c)
+{
+ unsigned int col = CTX::instance()->packColor(a, b, c, 255);
+ _triangles[0]->setCol(col);
+ _triangles[1]->setCol(col);
+ if (_triangles[2]) {
+ _triangles[2]->setCol(col);
+ _triangles[3]->setCol(col);
+ }
+}
+
+
/** Rec2DVertex **/
/*******************/
Rec2DVertex::Rec2DVertex(MVertex *v, std::list<MTriangle*> tri, int onWhat,
@@ -632,7 +1276,7 @@ Rec2DVertex::Rec2DVertex(double u, double v)
_v = new MVertex(p.x(), p.y(), p.z(), Recombine2D::getGFace());
static int i = 0;
if (++i < 20)
- Msg::Info("resulting point = [%g,%g,%g]", p.x(), p.y(), p.z());
+ Msg::Info("resulting point = [%g, %g, %g]", p.x(), p.y(), p.z());
}
double Rec2DVertex::_computeAngle(MVertex *v,
@@ -751,7 +1395,7 @@ void Rec2DVertex::getNeighbours(std::map<Rec2DVertex*, int> &vert,
else
rv = _edges[i]->getRVertex(1);
if ((it = vert.find(rv)) != vert.end())
- it->second += _edges[i]->getWeight() - EDGE_BASE;
+ it->second += _edges[i]->getWeight() - REC2D_EDGE_BASE;
else
vert[rv] = _edges[i]->getWeight();
}
@@ -788,10 +1432,29 @@ void Rec2DVertex::initStaticTable()
}
}
+void Rec2DVertex::remove(Rec2DEdge *re)
+{
+ std::vector<Rec2DEdge*>::iterator it = _edges.begin();
+ bool b = false;
+ for (; it != _edges.end() ; ++it)
+ if (*it == re) {
+ _edges.erase(it--);
+ b = true;
+ break;
+ }
+ if (b) {
+ Recombine2D::addValVert(getGain(-1));
+ --_numEl;
+ }
+ else
+ Msg::Error("[Rec2DVertex] Edge not found");
+}
+
+
/** Rec2DEdge **/
/*****************/
Rec2DEdge::Rec2DEdge(MEdge e, mapofVertices &vert, std::list<MTriangle*> &tri)
-: _lastChange(-1), _qual(.0), _weight(EDGE_BASE),
+: _lastChange(-1), _qual(.0), _weight(REC2D_EDGE_BASE),
_rv1(vert[e.getVertex(0)]), _rv2(vert[e.getVertex(1)])
{
_rv1->add(this);
@@ -799,7 +1462,7 @@ Rec2DEdge::Rec2DEdge(MEdge e, mapofVertices &vert, std::list<MTriangle*> &tri)
}
Rec2DEdge::Rec2DEdge(Rec2DVertex *rv1, Rec2DVertex *rv2)
-: _lastChange(-1), _qual(.0), _weight(EDGE_BASE), _rv1(rv1), _rv2(rv2)
+: _lastChange(-1), _qual(.0), _weight(REC2D_EDGE_BASE), _rv1(rv1), _rv2(rv2)
{
}
@@ -810,13 +1473,15 @@ double Rec2DEdge::getQual()
return _qual;
}
-void Rec2DEdge::_computeQual()
+void Rec2DEdge::_computeQual() //*
{
static int i = 0; if (++i < 2) Msg::Warning("FIXME compute qual edge and update of vertex");
//_rv1->update();
//_rv2->update();
double adimLength = _straightAdimLength();
double alignment = _straightAlignment();
+ if (adimLength > 1)
+ adimLength = 1/adimLength;
_qual = Recombine2D::edgeReward(adimLength, alignment);
_lastChange = Recombine2D::getNumChange();
}
@@ -853,11 +1518,16 @@ double Rec2DEdge::_straightAlignment()
double alpha2 = angleEdge - angle2;
crossField2d::normalizeAngle(alpha1);
crossField2d::normalizeAngle(alpha2);
- alpha1 = std::min(alpha1, .5 * M_PI - alpha1);
- alpha2 = std::min(alpha2, .5 * M_PI - alpha2);
+ alpha1 = 1 - 4 * std::min(alpha1, .5 * M_PI - alpha1) / M_PI;
+ alpha2 = 1 - 4 * std::min(alpha2, .5 * M_PI - alpha2) / M_PI;
double lc1 = (*Recombine2D::bgm())(_rv1->u(), _rv1->v(), .0);
double lc2 = (*Recombine2D::bgm())(_rv1->u(), _rv1->v(), .0);
return (alpha1/lc1 + alpha2/lc2) / (1/lc1 + 1/lc2);
-}
\ No newline at end of file
+}
+
+bool lessRec2DAction::operator()(Rec2DAction *ra1, Rec2DAction *ra2) const
+{
+ return *ra1 < *ra2;
+}
diff --git a/Mesh/meshGFaceRecombine.h b/Mesh/meshGFaceRecombine.h
index 619221cec5..868fff33f4 100644
--- a/Mesh/meshGFaceRecombine.h
+++ b/Mesh/meshGFaceRecombine.h
@@ -19,12 +19,16 @@ class Rec2DEdge;
class Rec2DVertex;
class Rec2DAction;
-typedef std::set<Rec2DAction*/*, lessRec2DAction*/> setofRec2DAction;
+typedef std::list<Rec2DAction*> setofRec2DAction;
typedef std::map<MVertex*, Rec2DVertex*> mapofVertices;
typedef std::map<MEdge, Rec2DEdge*, Less_Edge> mapofEdges;
typedef std::map<MElement*, std::set<Rec2DAction*> > mapofElementActions;
typedef std::map<MQuadrangle*, std::set<MElement*> > mapofAdjacencies;
+struct lessRec2DAction {
+ bool operator()(Rec2DAction*, Rec2DAction*) const;
+};
+
class Recombine2D {
private :
int _numEdge, _numVert, _numChange;
@@ -40,6 +44,9 @@ class Recombine2D {
std::vector<Rec2DAction*> _obsoleteAction;
static Recombine2D *_current;
+ std::vector<MTriangle*> _tri;
+ std::vector<MQuadrangle*> _quad;
+
public :
static bool bo;
@@ -47,6 +54,7 @@ class Recombine2D {
Recombine2D(GModel*);
~Recombine2D();
+ bool recombine();
void apply(bool){return;}
static inline double edgeReward(double adimLength, double alignment)
@@ -57,11 +65,36 @@ class Recombine2D {
int numVert, double valVert );
static Rec2DEdge* getREdge(MEdge);
static Rec2DVertex* getRVertex(MVertex*);
+ static void remove(MEdge);
+ static void remove(MVertex*);
static inline GFace* getGFace() {return _current->_gf;}
static inline void addObsleteAction(Rec2DAction *a) {
_current->_obsoleteAction.push_back(a);
+ } // correct ?
+ static inline void addNumEdge(int a) {_current->_numEdge += a;}
+ static inline void addNumVert(int a) {_current->_numVert += a;}
+ static inline void addValEdge(double a) {_current->_valEdge += a;}
+ static inline void addValVert(double a) {_current->_valVert += a;}
+
+ static void remove(MTriangle*);
+ static void add(MQuadrangle*);
+
+ static inline void addVert2Par(Rec2DVertex *rv, int a)
+ {_current->_parities[a].insert(rv);}
+ static void associateParity(int pOld, int pNew);
+ static inline int highPar() {
+ std::map<int, std::set<Rec2DVertex*> >::iterator it;
+ it = _current->_parities.end();
+ --it;
+ return it->first;
}
- //static inline int* getPtrNumChange() {return &_numChange;}
+
+ private :
+ bool _isAllQuad(Rec2DAction *action);
+ void _choose(Rec2DAction*);
+ void _removeIncompatible(Rec2DAction*);
+ void _removeReference(Rec2DAction*);
+ void _getNeighbours(Rec2DAction*, std::set<Rec2DAction*>&);
};
class Rec2DAction {
@@ -73,10 +106,26 @@ class Rec2DAction {
double _globValIfExecuted;
public :
- inline int getVectPosition() {return _position;}
+ bool operator<(Rec2DAction&);
+ inline int getVectPosition() const {return _position;}
inline void setVectPosition(int p) {_position = p;}
- //virtual double getReward() = 0;
+
+ virtual int getNum() = 0;
+ virtual double getReward() {
+ if (_lastUpdate < Recombine2D::getNumChange())
+ _computeGlobVal();
+ return _globValIfExecuted - Recombine2D::getGlobalValue(0, .0, 0, .0);
+ }
virtual bool isObsolete() = 0;
+ virtual void getTriangles(std::set<MTriangle*>&) = 0;
+ virtual MTriangle* getTriangle(int) = 0;
+ virtual void apply() = 0;
+ virtual void getParities(int*) = 0;
+ virtual Rec2DVertex* getRVertex(int) = 0;
+ virtual void color(int, int, int) = 0;
+
+ protected :
+ MQuadrangle* _createQuad(std::vector<Rec2DEdge*>&) const;
private :
virtual void _computeGlobVal() = 0;
@@ -90,9 +139,20 @@ class Rec2DTwoTri2Quad : public Rec2DAction {
public :
Rec2DTwoTri2Quad(Rec2DEdge*, std::list<MTriangle*>&);
- //double getReward();
bool isObsolete();
-
+ void getTriangles(std::set<MTriangle*> &tri) {
+ tri.insert(_triangles[0]);
+ tri.insert(_triangles[1]);
+ }
+ MTriangle* getTriangle(int a) {
+ return _triangles[a];
+ }
+ void apply();
+ void getParities(int*);
+ Rec2DVertex* getRVertex(int i) {return _vertices[i];}
+ void color(int, int, int);
+
+ int getNum() {return 1;}
private :
void _computeGlobVal();
};
@@ -105,9 +165,22 @@ class Rec2DFourTri2Quad : public Rec2DAction {
public :
Rec2DFourTri2Quad(Rec2DVertex*, std::list<MTriangle*>&);
- //double getReward();
bool isObsolete();
-
+ void getTriangles(std::set<MTriangle*> &tri) {
+ tri.insert(_triangles[0]);
+ tri.insert(_triangles[1]);
+ tri.insert(_triangles[2]);
+ tri.insert(_triangles[3]);
+ }
+ MTriangle* getTriangle(int a) {
+ return _triangles[a];
+ }
+ void apply();
+ void getParities(int*);
+ Rec2DVertex* getRVertex(int i) {return _vertices[i];}
+ void color(int, int, int);
+
+ int getNum() {return 2;}
private :
void _computeGlobVal();
};
@@ -122,15 +195,33 @@ class Rec2DCollapseTri : public Rec2DAction {
public :
Rec2DCollapseTri(Rec2DVertex*, std::list<MTriangle*>&);
Rec2DCollapseTri(Rec2DEdge*, std::list<MTriangle*>&);
- //double getReward();
+ double getReward() {
+ return std::max(_globValIfExecuted, _globValIfExecuted2)
+ - Recombine2D::getGlobalValue(0, .0, 0, .0);
+ }
bool isObsolete();
-
+ void getTriangles(std::set<MTriangle*> &tri) {
+ tri.insert(_triangles[0]);
+ tri.insert(_triangles[1]);
+ if (_triangles[2]) {
+ tri.insert(_triangles[2]);
+ tri.insert(_triangles[3]);
+ }
+ }
+ MTriangle* getTriangle(int a) {
+ return _triangles[a];
+ }
+ void apply();
+ void getParities(int*);
+ Rec2DVertex* getRVertex(int i) {return _vertices[i];}
+ void color(int, int, int);
+
+ int getNum() {return 3;}
private :
void _computeGlobVal();
void _qualCavity(double &valVert, double &valEdge, int &numEdge,
std::map<Rec2DVertex*, int> &vertices,
Rec2DVertex *imposedVert = NULL );
-
};
class Rec2DListAction {
@@ -158,13 +249,17 @@ class Rec2DVertex {
std::map<MVertex*, std::set<GEdge*> >);
Rec2DVertex(double u, double v);
- void add(Rec2DEdge *re) {_edges.push_back(re);}
- void setOnBoundary() {if (_onWhat > 0) _onWhat = -1;
- static bool a=1; if(a){Msg::Error("FIXME boundary");a=0;}}
- bool getIsOnBoundary() {return _onWhat < 1;}
+ inline void add(Rec2DEdge *re) {_edges.push_back(re);}
+ void remove(Rec2DEdge*);
+ inline void setOnBoundary() {if (_onWhat > 0) _onWhat = 0;
+ static bool a=1; if(a){Msg::Warning("FIXME boundary");a=0;}}
+ inline bool getIsOnBoundary() {return _onWhat < 1;}
double getQual(int numEdge = -1);
inline int getParity() {return _parity;}
- inline void setParity(int p) {_parity = p;}
+ inline void setParity(int p) {
+ _parity = p;
+ Recombine2D::addVert2Par(this, p);
+ }
inline MVertex* getMVertex() {return _v;}
void getxyz(double *xyz) {
xyz[0] = _v->x(); xyz[1] = _v->y(); xyz[2] = _v->z();}
@@ -205,7 +300,4 @@ class Rec2DEdge {
double _straightAlignment();
};
-#endif
-
-//idee : lors du parcours de larbre avec un horizon : copier les Rec2dVertex et les Rec2DEdge.
-//idee : faire carrŽment : arrete d'un triangle = 0, arrete d'un quad = 1/2
+#endif
\ No newline at end of file
--
GitLab