diff --git a/Mesh/meshRecombine2D.cpp b/Mesh/meshRecombine2D.cpp
index ac9a18f2e332f38115e1cf6fe202d357dcfd9e3b..c28768c0b3d22998c42f669faeee1801fcb8dfda 100644
--- a/Mesh/meshRecombine2D.cpp
+++ b/Mesh/meshRecombine2D.cpp
@@ -28,8 +28,8 @@ Recombine2D::Recombine2D(GFace *gf, int horizon)
//return;
- laplaceSmoothing(gf,100);
- gf->model()->writeMSH("befSquare.msh");
+ //laplaceSmoothing(gf,100);
+ //gf->model()->writeMSH("befSquare.msh");
Msg::Info("Branching with horizon %d", _horizon = HORIZ);
@@ -120,7 +120,7 @@ void Recombine2D::_buildEdgeToElement(std::vector<E*> &elements, e2t_cont &adj)
void Recombine2D::_recombine()
{
- SetBoundingBox();
+ //SetBoundingBox();
/*Map_Tri_Recomb::iterator itt;
diff --git a/Mesh/meshRecombine2D.h b/Mesh/meshRecombine2D.h
index c441b108db46543422de546ab38b5d1b4217902f..d80ca7e9950b2aae97331bb09ce2627f6f6d7d23 100644
--- a/Mesh/meshRecombine2D.h
+++ b/Mesh/meshRecombine2D.h
@@ -20,12 +20,14 @@
class RecombTriangle;
class Recomb2D_Node;
-class Rec2d_node;
-class Rec2d_edge;
-class Rec2d_cycle;
+class TrianglesUnion;
+class Rec2d_vertex;
struct lessRecombTri {
bool operator()(RecombTriangle *rt1, RecombTriangle *rt2) const;
};
+struct lessTriUnion {
+ bool operator()(TrianglesUnion *, TrianglesUnion *) const;
+};
typedef std::set<RecombTriangle*,lessRecombTri> Set_Recomb;
typedef std::map<MElement*,std::set<RecombTriangle*> > Map_Tri_Recomb;
@@ -55,13 +57,15 @@ class Recombine2D {
~Recombine2D();
int apply();
- double getBenef() const {return _benef;}
- int numTriangle() const {return _isolated.size();}
+ int apply(bool);
+ inline double getBenef() const {return _benef;}
+ inline int numTriangle() const {return _isolated.size();}
private :
- std::set<Rec2d_node*> _nodes;
- std::set<Rec2d_edge*> _edges;
- std::set<Rec2d_cycle*> _cycles;
+ //std::set<TrianglesUnion*, lessTriUnion> _setQuads;
+ std::list<TrianglesUnion*> _setQuads;
+ void _removeImpossible(std::list<TrianglesUnion*>::iterator);
+ void _recombine(bool);
};
class RecombTriangle {
@@ -82,15 +86,15 @@ class RecombTriangle {
MQuadrangle *createQuad() const;
bool operator<(const RecombTriangle &other) const;
- void triangles(std::vector<MElement *> &v) const {v = _t;}
- int numTriangles() const {return (int) _t.size();}
+ 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 _t[i];
return NULL;
}
- bool isQuad() const {return _formingQuad;}
- double getBenef() const {return _benefit;}
+ inline bool isQuad() const {return _formingQuad;}
+ inline double getBenef() const {return _benefit;}
double compute_alignment(const MEdge&, MElement*, MElement*);
};
@@ -114,51 +118,71 @@ class Recomb2D_Node {
void erase();
void blocking(const Map_Tri_Node::iterator &);
- bool isBetter() {return _blocking.size() < 2;}
- Set_Recomb::iterator getItRecomb() {return _recomb;}
- double getTotBenef() {return _totBenef;}
- int getnSkip() {return _nskip;}
- double getBound() {return _totBenef + _benef * _depth;}
+ 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_node {
+class Rec2d_vertex {
private :
- std::set<Rec2d_edge*> _freeEdges;
- std::set<Rec2d_cycle*> _cycle3;
- std::set<Rec2d_cycle*> _cycle4;
-
- public :
- Rec2d_node(){}
-
- void addCycle(Rec2d_cycle*, int n = 0);
- void addEdge(Rec2d_edge*);
+ // _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;
- void print();
-};
-class Rec2d_edge {
- private :
- Rec2d_node *_nodes[2];
- std::set<Rec2d_cycle*> _cycles;
+ void _initStaticTable();
+ double _computeAngle(MVertex *, std::list<MTriangle*> &, std::set<GEdge*> &);
public :
- Rec2d_edge(Rec2d_node*, Rec2d_node*);
-
- void addCycle(Rec2d_cycle*);
+ Rec2d_vertex(MVertex *, std::list<MTriangle*> &,
+ std::map<MVertex*, std::set<GEdge*> > &);
- void print();
+ inline void changeNumEl(int c) {_numEl += c;}
+ double getReward();
+ double getReward(int);
};
-class Rec2d_cycle {
+
+class TrianglesUnion {
private :
- std::set<Rec2d_edge*> _edges;
-
- public :
- Rec2d_cycle(){}
- void addEdge(Rec2d_edge*);
- int size() {return _edges.size();}
+ int _numEmbEdge, _numEmbVert, _numBoundVert, _numTri, _computation;
+ double _embEdgeValue, _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:
+
+ double _computeEdgeLength(GFace *, MVertex **, double *u, double *v,
+ int numIntermedPoint= 1);
+ double _computeAlignment(MEdge &, std::list<MTriangle*> &);
+ void _update();
- void print();
+ public :
+ TrianglesUnion(GFace *, std::list<MTriangle*> &, std::list<MEdge> &,
+ std::list<Rec2d_vertex*> &, std::map<MVertex*,double*> &);
+
+ bool operator<(TrianglesUnion &other);
+ inline double getEdgeValue() {return _embEdgeValue;}
+ inline double getNumTriangles() {return _numTri;}
+ inline MTriangle *getTriangle(int num) {return _triangles[num];}
+ void select() {
+ _ValEdge -= _embEdgeValue;
+ _NumEdge -= _numEmbEdge;
+ _ValVert -= _embVertValue;
+ _NumVert -= _numEmbVert;
+ for (int i = 0; i < _numBoundVert; i++) {
+ _ValVert += _vertices[i]->getReward(-1);
+ }
+ _RECOMPUTE++;
+ }
+ MQuadrangle *createQuad() const;
};
-
-
#endif
diff --git a/Mesh/meshRecombine2D_2.cpp b/Mesh/meshRecombine2D_2.cpp
index 8ed74c74205b0bead78408e1962c8bc3582412d5..15ffe9716dd3d72a86325abe4cd247c6e473c781 100644
--- a/Mesh/meshRecombine2D_2.cpp
+++ b/Mesh/meshRecombine2D_2.cpp
@@ -8,6 +8,7 @@
//
#include "meshRecombine2D.h"
+#include "BackgroundMesh.h"
#include "GFace.h"
#include <cmath>
#include <FL/Fl.H>
@@ -18,102 +19,541 @@
static int HORIZ = 20;
+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;
+
Recombine2D::Recombine2D(GFace *gf)
: _horizon(0), _gf(gf), _benef(.0), _applied(true)
{
- std::map<MEdge, Rec2d_edge*, Less_Edge> mapEdge;
- std::map<MVertex*, Rec2d_node*> mapNode;
+ 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;
for (unsigned int i = 0; i < gf->triangles.size(); i++) {
- Rec2d_cycle *rc = new Rec2d_cycle();
- _cycles.insert(rc);
MTriangle *t = gf->triangles[i];
-
for (int j = 0; j < 3; j++) {
- MVertex *v = t->getVertex(j);
- std::map<MVertex*, Rec2d_node*>::iterator itn = mapNode.find(v);
- if (itn == mapNode.end()) {
- Rec2d_node *rn = new Rec2d_node();
- mapNode[v] = rn;
- _nodes.insert(rn);
- rn->addCycle(rc, 3);
+ mapVert[t->getVertex(j)].push_back(t);
+ mapEdge[t->getEdge(j)].push_back(t);
+ }
+ }
+
+ // 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)
+ fourTri.push_back(itvert);
+ Rec2d_vertex *n = new Rec2d_vertex((*itvert).first, (*itvert).second, mapCornerVert);
+ itV2N = mapV2N.insert(itV2N, std::make_pair((*itvert).first,n));
+ TrianglesUnion::_NumVert++;
+ TrianglesUnion::_ValVert += n->getReward();
+ }
+
+ // store mesh size and parametric coordinates 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);
+ }
}
- else
- (*itn).second->addCycle(rc, 3);
+ 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);
+ _setQuads.push_back(new TrianglesUnion(gf, (*(*it4)).second, listEmbEdge, listVert, mapV2LcUV));
+ }
+
+ // 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);
+ _setQuads.push_back(tu);
+ TrianglesUnion::_NumEdge++;
+ TrianglesUnion::_ValEdge += tu->getEdgeValue();
+ }
+ else
+ Msg::Info("[bord] %d", (*itedge).second.size());
+ }
+ _setQuads.sort(lessTriUnion());
+
+ _recombine(true);
+ _applied = true;*/
+}
+
+
+void Recombine2D::_recombine(bool a)
+{
+ int i = 0;
+ while (!_setQuads.empty()) {
+ Msg::Info("%d",i);
+ std::list<TrianglesUnion*>::iterator it = _setQuads.begin();
- for (int j = 0; j < 3; j++) {
- MEdge e = t->getEdge(j);
- std::map<MEdge, Rec2d_edge*>::iterator ite = mapEdge.find(e);
- if (ite == mapEdge.end()) {
- Rec2d_node *n0 = mapNode[e.getVertex(0)];
- Rec2d_node *n1 = mapNode[e.getVertex(1)];
- Rec2d_edge *re = new Rec2d_edge(n0, n1);
- n0->addEdge(re);
- n1->addEdge(re);
- mapEdge[e] = re;
- _edges.insert(re);
- rc->addEdge(re);
- re->addCycle(rc);
- }
- else {
- rc->addEdge((*ite).second);
- (*ite).second->addCycle(rc);
+ (*it)->select();
+
+ _quads.push_back((*it)->createQuad());
+
+ _removeImpossible(it);
+
+ _setQuads.sort(lessTriUnion());
+ i++;
+ }
+ Msg::Info("Done %d (%d)", i, _pairs.size());
+}
+
+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);
}
- }ZzZ
+ if (listBoundEdge.size() != 4)
+ Msg::Warning("[meshRecombine2D] Not 4 edges !");
+
+ MVertex *v1, *v2, *v3, *v4;
+ v1 = listBoundEdge.begin()->getMinVertex();
+ v2 = listBoundEdge.begin()->getMinVertex();
+ 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 Rec2d_node::addCycle(Rec2d_cycle *c, int n)
+
+int Recombine2D::apply(bool a)
{
- if (n == 3) {
- _cycle3.insert(c);
- return;
+ if (!_haveParam || _applied) return 0;
+
+ std::vector<MTriangle*> triangles2;
+ for (int i = 0; i < _gf->triangles.size(); i++) {
+ delete _gf->triangles[i];
}
- if (n == 4) {
- _cycle4.insert(c);
- return;
+ _gf->triangles = triangles2;
+ _gf->quadrangles = _quads;
+
+ _applied = true;
+ return 1;
+}
+
+void Recombine2D::_removeImpossible(std::list<TrianglesUnion*>::iterator it)
+{
+ std::set<MTriangle *> touched;
+
+ for (int i = 0; i < (*it)->getNumTriangles(); i++) {
+ touched.insert((*it)->getTriangle(i));
}
- if (c->size() == 3) {
- _cycle3.insert(c);
- return;
+
+ for (it = _setQuads.begin(); 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)
+ _setQuads.erase(it++);
+ else
+ it++;
}
- if (c->size() == 4) {
- _cycle4.insert(c);
- return;
+}
+
+
+
+Rec2d_vertex::Rec2d_vertex(MVertex *v,
+ std::list<MTriangle*> &setTri,
+ std::map<MVertex*, std::set<GEdge*> > &mapVert)
+: _numEl(setTri.size()), _angle(.0)
+{
+ _onWhat = v->onWhat()->dim() - 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.;
+ }
}
- Msg::Error("[Recombine2D] Can't tell if it's a cycle3 or a cycle4");
+
+ if (_Vvalues == NULL)
+ _initStaticTable();
}
-void Rec2d_node::addEdge(Rec2d_edge *e)
+
+void Rec2d_vertex::_initStaticTable()
{
- _freeEdges.insert(e); //que faire ?
+ 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]; //AVERIFIER SI CA MARCHE
+ 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];
+ }
}
-void Rec2d_node::print()
+
+
+double Rec2d_vertex::_computeAngle(MVertex *v,
+ std::list<MTriangle*> &setTri,
+ std::set<GEdge*> &setEdge)
{
- Msg::Info("Node : %d cycle3, %d cycle4, %d edges", _cycle3.size(), _cycle4.size(), _freeEdges.size());
+ 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 = vectub;
+ else if (normub > prec * normlb)
+ firstDer1 = vectlb;
+ else {
+ Msg::Warning("[meshRecombine2D] Bad precision, returning pi/2");
+ return M_PI / 2.;
+ }
+ if (k == 0)
+ firstDer0 = firstDer1;
+ }
+
+ Msg::Info("-- %lf, %lf, %lf", norm(firstDer0), norm(firstDer1), dot(firstDer0, firstDer1));
+ double n = 1 / norm(firstDer0);
+ firstDer0 = firstDer0 * n;
+ n = 1 / norm(firstDer1);
+ firstDer1 = firstDer1 * n;
+ Msg::Info("--Angles : %lf, %lf, %lf", norm(firstDer0), norm(firstDer1), dot(firstDer0, firstDer1));
+
+ 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));
+ }
+
+ Msg::Info("Angles : %lf, %lf, %lf", angle1, angleMesh, angle2);
+
+ return M_PI / 2.;
+
+ if (angleMesh < M_PI)
+ return angle1;
+ return angle2;
}
-Rec2d_edge::Rec2d_edge(Rec2d_node *n0, Rec2d_node *n1)
+double Rec2d_vertex::getReward()
{
- _nodes[0] = n0;
- _nodes[1] = n1;
+ if (_onWhat > -1)
+ return _Vvalues[_onWhat][_numEl-1];
+ else
+ return 1. - fabs(2. / M_PI * _angle / _numEl - 1.);
}
-void Rec2d_edge::addCycle(Rec2d_cycle *c)
+
+double Rec2d_vertex::getReward(int n)
{
- _cycles.insert(c);
+ if (n == 0) return .0;
+ if (_onWhat > -1) {
+ if (n > 0)
+ return _Vbenefs[_onWhat][_numEl-1];
+ else
+ return - _Vbenefs[_onWhat][_numEl-2];
+ }
+ else
+ return fabs(2./M_PI*_angle/_numEl - 1.) - fabs(2./M_PI*_angle/(_numEl + n) - 1.);
}
-void Rec2d_edge::print()
+
+
+TrianglesUnion::TrianglesUnion(GFace *gf,
+ std::list<MTriangle*> &t,
+ std::list<MEdge> &e,
+ std::list<Rec2d_vertex*> &v,
+ std::map<MVertex*,double*> &v2lcUV)
{
- Msg::Info("Edge : %d cycles", _cycles.size());
+ _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();
+ for (int k = 0; itt != t.end(); itt++, k++) {
+ _triangles[k] = *itt;
+ }
+
+ 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);
+ std::map<MVertex*,double*>::iterator itlc;
+ if ((itlc = v2lcUV.find(vert[i])) == v2lcUV.end()) {
+ double *a = new double[3];
+ gf->XYZtoUV(vert[i]->x(), vert[i]->y(), vert[i]->z(), a[1], a[2], 1.);
+ sumlc += a[0] = BGM_MeshSize(gf, a[1], a[2], vert[i]->x(), vert[i]->y(), vert[i]->z());
+ u[i] = a[1];
+ v[i] = a[2];
+ v2lcUV[vert[i]] = a;
+ }
+ else {
+ sumlc += (*itlc).second[0];
+ u[i] = (*itlc).second[1];
+ v[i] = (*itlc).second[2];
+ }
+ }
+
+ double length = _computeEdgeLength(gf, vert, u, v, 0);
+ _embEdgeValue += length / sumlc * (1 + _computeAlignment(*ite, t));
+ Msg::Info("Edge a : %lf/%lf = %lf <-> %lf", length, sumlc / 2, 2 * length / sumlc, _computeAlignment(*ite, t));
+ }
+
+ _vertices = new Rec2d_vertex*[_numBoundVert];
+ std::list<Rec2d_vertex*>::iterator itv = v.begin();
+ for (int k = 0; itv != v.end() && k < _numBoundVert; itv++, k++) {
+ _globValIfExecuted += (*itv)->getReward(-1);
+ _vertices[k] = *itv;
+ }
+ for (_numEmbVert = 0; itv != v.end(); itv++, _numEmbVert++)
+ _embVertValue += (*itv)->getReward();
+
+ _computation = _RECOMPUTE - 1;
+ _update();
}
-void Rec2d_cycle::addEdge(Rec2d_edge *e)
+
+double TrianglesUnion::_computeEdgeLength(GFace *gf, MVertex **vert,
+ double *u, double *v, int n)
{
- _edges.insert(e);
+ 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;
}
-void Rec2d_cycle::print()
+double TrianglesUnion::_computeAlignment(MEdge &e, std::list<MTriangle*> <)
{
- Msg::Info("Cycle : %d edges", _edges.size());
+ 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 > 2) {
+ Msg::Warning("[meshRecombine2D] error in alignment computation, returning 0");
+ return .0;
+ }
+ 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) / (_NumEdge - _numEmbEdge);
+ _globValIfExecuted = alpha * beta;
+
+ _computation = _RECOMPUTE;
+}
+bool TrianglesUnion::operator<(TrianglesUnion &other)
+{
+ _update();
+ other._update();
+ return _globValIfExecuted < other._globValIfExecuted;
}
+bool lessTriUnion::operator()(TrianglesUnion *tu1, TrianglesUnion *tu2) const
+{
+ return *tu1 < *tu2;
+}
+
+/*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*/