diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 47cf062de496ef1963a8708d89d1f0f03a6f3df1..2d5e97f62e128260ba1dab86db23c2e6481570c6 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -1519,7 +1519,7 @@ void _triangleSplit (GFace *gf, MElement *t, bool swop = false) {
}
-struct RecombineTriangle
+/*struct RecombineTriangle
{
MElement *t1, *t2;
double angle;
@@ -1558,7 +1558,7 @@ struct RecombineTriangle
//return angle < other.angle;
return total_cost < other.total_cost; //addition for class Temporary
}
-};
+};*///DEV amaury
static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
{
@@ -1602,17 +1602,15 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
e2t_cont adj;
buildEdgeToElement(gf->triangles, adj);
- std::map<MVertex*,std::pair<MElement*,MElement*> > makeGraphPeriodic;
-
std::vector<RecombineTriangle> pairs;
- int nbValidPairs = 0;
+ std::map<MVertex*,std::pair<MElement*,MElement*> > makeGraphPeriodic;
+
for(e2t_cont::iterator it = adj.begin(); it!= adj.end(); ++it){
if(it->second.second &&
it->second.first->getNumVertices() == 3 &&
it->second.second->getNumVertices() == 3 &&
(emb_edgeverts.find(it->first.getVertex(0)) == emb_edgeverts.end() ||
emb_edgeverts.find(it->first.getVertex(1)) == emb_edgeverts.end())){
-
pairs.push_back(RecombineTriangle(it->first,
it->second.first,
it->second.second));
@@ -1638,7 +1636,6 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
std::sort(pairs.begin(),pairs.end());
std::set<MElement*> touched;
-
if(CTX::instance()->mesh.algoRecombine == 1){
#ifdef HAVE_MATCH
int ncount = gf->triangles.size();
@@ -1812,7 +1809,7 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
Msg::Warning("Gmsh should be compiled with the Blossom IV code and CONCORDE in order to allow the Blossom optimization");
#endif
}
-
+
std::vector<RecombineTriangle>::iterator itp = pairs.begin();
while(itp != pairs.end()){
// recombine if difference between max quad angle and right
@@ -1863,6 +1860,12 @@ void recombineIntoQuads(GFace *gf,
bool topologicalOpti,
bool nodeRepositioning)
{
+ if (false) { //|| CTX::instance()->mesh.algoRecombine == 2) {
+ Msg::Warning("New technic of recombination in development (Amaury)");
+ devRecombine_Amaury(gf,topologicalOpti,nodeRepositioning);
+ return;
+ }
+
double t1 = Cpu();
bool haveParam = true;
@@ -1920,6 +1923,204 @@ void recombineIntoQuads(GFace *gf,
Msg::Info("Simple recombination algorithm completed (%g s)", t2 - t1);
}
+static void _makeQuad_devAmaury(GFace *gf, SetRecombTri::iterator it)
+{
+ MElement *t1 = (*it)->t1;
+ MElement *t2 = (*it)->t2;
+ int orientation = 0;
+ for(int i = 0; i < 3; i++) {
+ if(t1->getVertex(i) == (*it)->n1) {
+ if(t1->getVertex((i + 1) % 3) == (*it)->n2)
+ orientation = 1;
+ else
+ orientation = -1;
+ break;
+ }
+ }
+ MQuadrangle *q;
+ if(orientation < 0)
+ q = new MQuadrangle((*it)->n1, (*it)->n3, (*it)->n2, (*it)->n4);
+ else
+ q = new MQuadrangle((*it)->n1, (*it)->n4, (*it)->n2, (*it)->n3);
+ gf->quadrangles.push_back(q);
+}
+
+static double *_devRecombine_Amaury(GFace *gf, int depth, bool cubicGraph = 1)
+{
+ int success = 1;
+ double glob_benef = 0;
+
+ std::set<MVertex*> emb_edgeverts;
+ {
+ std::list<GEdge*> emb_edges = gf->embeddedEdges();
+ std::list<GEdge*>::iterator ite = emb_edges.begin();
+ while(ite != emb_edges.end()){
+ if(!(*ite)->isMeshDegenerated()){
+ emb_edgeverts.insert((*ite)->mesh_vertices.begin(),
+ (*ite)->mesh_vertices.end() );
+ emb_edgeverts.insert((*ite)->getBeginVertex()->mesh_vertices.begin(),
+ (*ite)->getBeginVertex()->mesh_vertices.end());
+ emb_edgeverts.insert((*ite)->getEndVertex()->mesh_vertices.begin(),
+ (*ite)->getEndVertex()->mesh_vertices.end());
+ }
+ ++ite;
+ }
+ }
+
+ {
+ std::list<GEdge*> _edges = gf->edges();
+ std::list<GEdge*>::iterator ite = _edges.begin();
+ while(ite != _edges.end()){
+ if(!(*ite)->isMeshDegenerated()){
+ if ((*ite)->isSeam(gf)){
+ emb_edgeverts.insert((*ite)->mesh_vertices.begin(),
+ (*ite)->mesh_vertices.end() );
+ emb_edgeverts.insert((*ite)->getBeginVertex()->mesh_vertices.begin(),
+ (*ite)->getBeginVertex()->mesh_vertices.end());
+ emb_edgeverts.insert((*ite)->getEndVertex()->mesh_vertices.begin(),
+ (*ite)->getEndVertex()->mesh_vertices.end());
+ }
+ }
+ ++ite;
+ }
+ }
+
+ e2t_cont adj;
+ buildEdgeToElement(gf->triangles, adj);
+
+ SetRecombTri pairs;
+ MapElemToIncomp incomp;
+
+ for (e2t_cont::iterator it = adj.begin(); it!= adj.end(); ++it) {
+ if (it->second.second &&
+ it->second.first->getNumVertices() == 3 &&
+ it->second.second->getNumVertices() == 3 &&
+ (emb_edgeverts.find(it->first.getVertex(0)) == emb_edgeverts.end() ||
+ emb_edgeverts.find(it->first.getVertex(1)) == emb_edgeverts.end())) {
+ RecombineTriangle *rt = new RecombineTriangle(it->first, it->second.first, it->second.second);
+ pairs.insert(rt);
+ incomp[it->second.first].push_back(rt);
+ incomp[it->second.second].push_back(rt);
+ }
+ }
+
+ std::set<MElement*> touched;
+ SetRecombTri::iterator itp = pairs.begin();
+ while (itp != pairs.end()) {
+ RecombSeq rs(itp, pairs.end(), depth);
+
+ if (rs.optimal() || !rs.shouldBeSkipped()) {
+ _makeQuad_devAmaury(gf, itp);
+ glob_benef += 1 - (*itp)->cost_quality;
+
+ MElement *t1 = (*itp)->t1;
+ MElement *t2 = (*itp)->t2;
+ touched.insert(t1);
+ touched.insert(t2);
+
+ MapElemToIncomp::iterator iti;
+ SetRecombTri::iterator itmp;
+
+ iti = incomp.find(t1);
+ for (int i = 0; i < iti->second.size(); i++)
+ if ((itmp = pairs.find(iti->second[i])) != pairs.end())
+ pairs.erase(itmp);
+ incomp.erase(iti);
+
+ iti = incomp.find(t2);
+ for (int i = 0; i < iti->second.size(); i++)
+ if ((itmp = pairs.find(iti->second[i])) != pairs.end())
+ pairs.erase(itmp);
+ incomp.erase(iti);
+ }
+ else {
+ pairs.erase(itp);
+ }
+ itp = pairs.begin();
+ }
+
+ std::vector<MTriangle*> triangles2;
+ for(unsigned int i = 0; i < gf->triangles.size(); i++){
+ if(touched.find(gf->triangles[i]) == touched.end()){
+ triangles2.push_back(gf->triangles[i]);
+ }
+ else {
+ delete gf->triangles[i];
+ }
+ }
+ gf->triangles = triangles2;
+
+ double *b = new double[2];
+ *b = glob_benef;
+ b[1] = touched.size()/2;
+ return b;
+}
+
+void devRecombine_Amaury(GFace *gf, bool topologicalOpti, bool nodeRepositioning)
+{
+ double t1 = Cpu();
+
+ bool haveParam = true;
+ if(gf->geomType() == GEntity::DiscreteSurface && !gf->getCompound())
+ haveParam = false;
+
+ if(haveParam && topologicalOpti)
+ removeFourTrianglesNodes(gf, false);
+
+ gf->model()->writeMSH("before.msh");
+ Msg::Info("Nombre de triangles : %d", gf->triangles.size());
+ double tt1 = Cpu();
+ int d = 700;
+ double max = 0, maxQuad = 0;
+ while (Cpu() - tt1 < 30) {
+ double tt2 = Cpu();
+ double *b = _devRecombine_Amaury(gf, d);
+ if (*b > max) max = *b;
+ if (b[1] > maxQuad) maxQuad = b[1];
+ Msg::Info("depth %d in (%g s) : benef = %lf / numQuad = %f", d++, Cpu() - tt2, *b, b[1]);
+ if (d > 200 && d < 100) d = 620;
+ break;
+ }
+ Msg::Info("Benef max = %lf", max);
+ Msg::Info("Max quad = %f sur %d", maxQuad, gf->triangles.size()/2);
+
+ gf->model()->writeMSH("raw.msh");
+ /*if(haveParam && nodeRepositioning)
+ laplaceSmoothing(gf, CTX::instance()->mesh.nbSmoothing);
+
+ // blossom-quad algo
+ if(success){
+ if(topologicalOpti){
+ if(haveParam){
+ gf->model()->writeMSH("smoothed.msh");
+ int COUNT = 0;
+ char NAME[256];
+ while(1){
+ int x = removeTwoQuadsNodes(gf);
+ if(x){ sprintf(NAME,"iter%dTQ.msh",COUNT++); gf->model()->writeMSH(NAME);}
+ int y = removeDiamonds(gf);
+ if(y){ sprintf(NAME,"iter%dD.msh",COUNT++); gf->model()->writeMSH(NAME); }
+ laplaceSmoothing(gf);
+ int z = 0; //edgeSwapQuadsForBetterQuality(gf);
+ if(z){ sprintf(NAME,"iter%dS.msh",COUNT++); gf->model()->writeMSH(NAME); }
+ if (!(x+y+z)) break;
+ }
+ }
+ edgeSwapQuadsForBetterQuality(gf);
+ //if(z){ sprintf(NAME,"iter%dS.msh",COUNT++); gf->model()->writeMSH(NAME); }
+ if(haveParam) laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing);
+ }
+ double t2 = Cpu();
+ Msg::Info("Tree recombination algorithm completed (%g s)", t2 - t1);
+ return;
+ }
+
+ gf->model()->writeMSH("after.msh");*/
+
+ double t2 = Cpu();
+ Msg::Info("Simple recombination algorithm completed (%g s)", t2 - t1);
+}
+
void quadsToTriangles(GFace *gf, double minqual = -10000.)
{
std::vector<MQuadrangle*> qds;
@@ -2130,3 +2331,137 @@ void Temporary::quadrilaterize(std::string file_name,double _w1,double _w2,doubl
}
}
+
+bool RecombineTriangle::operator < (const RecombineTriangle &other) const
+{
+ //return angle < other.angle;
+ //return total_cost < other.total_cost; //addition for class Temporary
+ return cost_quality < other.cost_quality; // DEV AMAURY !!!!!!!!!!
+}
+bool RecombSeq::parcours(int d, int countdown)
+{
+ if (d < 1 || countdown < 1) return true;
+
+ if (_toBeSkipped.find(*_it2) != _toBeSkipped.end() ||
+ _toBeKept.find(*_it2) != _toBeKept.end() ) {
+ if (++_it2 == _ite) return true;
+ return parcours(d, countdown);
+ }
+
+ MElement *t1 = (*_it2)->t1;
+ MElement *t2 = (*_it2)->t2;
+
+ bool a, b;
+ std::map< MElement*, std::pair<RecombineTriangle*, bool> >::iterator m1, m2;
+ a = (m1 = _touched.find(t1)) != _touched.end();
+ b = (m2 = _touched.find(t2)) != _touched.end();
+
+ if (a || b) {
+ MElement *t3;
+ if (a) {
+ (*m1).second.second = true;
+ if((*m1).second.first->t1 == t1)
+ t3 = (*m1).second.first->t2;
+ else
+ t3 = (*m1).second.first->t1;
+ if (_touched[t3].second == true) _canBeSkipped.insert((*m1).second.first);
+ }
+ if (b) {
+ (*m2).second.second = true;
+ if((*m2).second.first->t1 == t2)
+ t3 = (*m2).second.first->t2;
+ else
+ t3 = (*m2).second.first->t1;
+ if (_touched[t3].second == true) _canBeSkipped.insert((*m2).second.first);
+ }
+ if (++_it2 == _ite) return true;
+ return parcours(d, --countdown);
+ }
+
+ _touched[t1] = std::make_pair(*_it2,false);
+ _touched[t2] = std::make_pair(*_it2,false);
+ _benef += 1 - (*_it2)->cost_quality;
+ if (++_it2 == _ite) return true;
+ d--;
+ //_perspect = (1 - (*(_it2))->cost_quality) * d;
+ //if (_benef+_perspect < best) return false;
+
+ return parcours(d);
+}
+
+bool RecombSeq::shouldBeSkipped()
+{
+ //if (_canBeSkipped.find(*_it1) == _canBeSkipped.end())
+ // return false;
+ SetRecombTri::iterator itmp = _it1;
+ itmp++;
+ RecombSeq rs(itmp, _ite, _depth);
+ if(rs.getBenef() > getBenef())
+ return true;
+ return false;
+ if (rs.getOptimalBenef(true) > getOptimalBenef(false))
+ return true;
+ return false;
+}
+
+RecombSeq::RecombSeq(SetRecombTri::iterator it,
+ SetRecombTri::iterator itend,
+ int depth,
+ SetRecombTri toBeSkipped,
+ SetRecombTri toBeKept)
+ : _it1(it), _ite(itend), _depth(depth),
+ _toBeSkipped(toBeSkipped), _toBeKept(toBeKept)
+{
+ _benef = 0;
+ _it2 = it;
+
+ if (depth < _toBeKept.size()) return;
+
+ SetRecombTri::iterator itmp = _toBeKept.begin();
+ while (itmp != _toBeKept.end()) {
+ MElement *t1 = (*itmp)->t1;
+ MElement *t2 = (*itmp)->t2;
+ _touched[t1] = std::make_pair(*itmp,false);
+ _touched[t2] = std::make_pair(*itmp,false);
+ _benef += 1 - (*itmp++)->cost_quality;
+ }
+
+ parcours(depth - _toBeKept.size());
+}
+
+double RecombSeq::getOptimalBenef(bool firstSkipable)
+{
+ if (!firstSkipable)
+ if (_canBeSkipped.find(*_it1) != _canBeSkipped.end())
+ _canBeSkipped.erase(*_it1);
+
+ if (_canBeSkipped.empty())
+ return getBenef();
+
+ double maxBenef = getBenef();
+ SetRecombTri toBeKept(_toBeKept);
+ SetRecombTri toBeSkipped(_toBeSkipped);
+ SetRecombTri::iterator it = _canBeSkipped.begin();
+
+ if (_toBeKept.find(*it) == _toBeKept.end()) {
+ toBeSkipped.insert(*it);
+ RecombSeq rs(_it1, _it2, _depth, toBeSkipped, toBeKept);
+ double b;
+ if ((b = rs.getOptimalBenef(firstSkipable)) > maxBenef)
+ maxBenef = b;
+ }
+
+ while(it != --_canBeSkipped.end()) {
+ toBeSkipped.erase(*it);
+ toBeKept.insert(*it++);
+ if (_toBeKept.find(*it) == _toBeKept.end()) {
+ toBeSkipped.insert(*it);
+ RecombSeq rs(_it1, _it2, _depth, toBeSkipped, toBeKept);
+ double b;
+ if ((b = rs.getOptimalBenef(firstSkipable)) > maxBenef)
+ maxBenef = b;
+ }
+ }
+
+ return maxBenef;
+}
diff --git a/Mesh/meshGFaceOptimize.h b/Mesh/meshGFaceOptimize.h
index 55b2d34e371aef4b2dcfe9d4c7a66723f4a8af78..93bcf5c92429a01a6e14a26212e67dff787f39bd 100644
--- a/Mesh/meshGFaceOptimize.h
+++ b/Mesh/meshGFaceOptimize.h
@@ -147,4 +147,111 @@ class Temporary{
static double compute_alignment(const MEdge&,MElement*,MElement*);
};
+
+//DEV Amaury
+void devRecombine_Amaury(GFace *gf,
+ bool topologicalOpti = true,
+ bool nodeRepositioning = true);
+struct RecombineTriangle
+{
+ MElement *t1, *t2;
+ double angle;
+ double cost_quality; //addition for class Temporary
+ double cost_alignment; //addition for class Temporary
+ double total_cost; //addition for class Temporary
+ MVertex *n1, *n2, *n3, *n4;
+ RecombineTriangle(const MEdge &me, MElement *_t1, MElement *_t2)
+ : t1(_t1), t2(_t2)
+ {
+ n1 = me.getVertex(0);
+ n2 = me.getVertex(1);
+
+ if(t1->getVertex(0) != n1 && t1->getVertex(0) != n2) n3 = t1->getVertex(0);
+ else if(t1->getVertex(1) != n1 && t1->getVertex(1) != n2) n3 = t1->getVertex(1);
+ else if(t1->getVertex(2) != n1 && t1->getVertex(2) != n2) n3 = t1->getVertex(2);
+ if(t2->getVertex(0) != n1 && t2->getVertex(0) != n2) n4 = t2->getVertex(0);
+ else if(t2->getVertex(1) != n1 && t2->getVertex(1) != n2) n4 = t2->getVertex(1);
+ else if(t2->getVertex(2) != n1 && t2->getVertex(2) != n2) n4 = t2->getVertex(2);
+
+ double a1 = 180 * angle3Vertices(n1, n4, n2) / M_PI;
+ double a2 = 180 * angle3Vertices(n4, n2, n3) / M_PI;
+ double a3 = 180 * angle3Vertices(n2, n3, n1) / M_PI;
+ double a4 = 180 * angle3Vertices(n3, n1, n4) / M_PI;
+ angle = fabs(90. - a1);
+ angle = std::max(fabs(90. - a2),angle);
+ angle = std::max(fabs(90. - a3),angle);
+ angle = std::max(fabs(90. - a4),angle);
+ cost_quality = 1.0 - std::max(1.0 - angle/90.0,0.0); //addition for class Temporary
+ cost_alignment = Temporary::compute_alignment(me,_t1,_t2); //addition for class Temporary
+ total_cost = Temporary::compute_total_cost(cost_quality,cost_alignment); //addition for class Temporary
+ total_cost = 100.0*cost_quality; //addition for class Temporary
+ }
+ bool operator < (const RecombineTriangle &other) const;
+ /*{
+ //return angle < other.angle;
+ //return total_cost < other.total_cost; //addition for class Temporary
+ }*/
+};
+struct lessRecombTri
+{
+ bool operator()(RecombineTriangle *rt1, RecombineTriangle *rt2)
+ {
+ return *rt1 < *rt2;
+ }
+};
+typedef std::set<RecombineTriangle*,lessRecombTri> SetRecombTri;
+typedef std::map< MElement*, std::vector<RecombineTriangle*> > MapElemToIncomp;
+class RecombSeq {
+ private:
+ //std::vector<int> _skipped;
+ //std::vector<int> _toBeSkipped;
+ //std::map<RecombineTriangle*, int> _numSkip;
+ SetRecombTri _canBeSkipped, _toBeSkipped, _toBeKept;
+ SetRecombTri::iterator _it1, _it2, _ite;
+ //std::set<MElement*> _touched;
+ //std::map< MElement*, std::pair<MElement*, bool> > _touched;
+ std::map< MElement*, std::pair<RecombineTriangle*, bool> > _touched;
+ double _benef, _perspect;
+ int _depth;
+
+ public:
+ RecombSeq(SetRecombTri::iterator it,
+ SetRecombTri::iterator itend,
+ int depth) : _it1(it), _ite(itend), _depth(depth)
+ {
+ _benef = 0;
+ _it2 = it;
+ parcours(depth);
+ }
+
+ RecombSeq(SetRecombTri::iterator it,
+ SetRecombTri::iterator itend,
+ int depth,
+ SetRecombTri toBeSkipped,
+ SetRecombTri toBeKept);
+ bool parcours(int d, int a = 10);
+ double getOptimalBenef(bool firstSkipable);
+
+ bool optimal() const {return _canBeSkipped.empty();}
+ bool shouldBeSkipped();
+ /*void getOptimal(SetRecombTri::iterator &it)
+ {
+ SetRecombTri::iterator itmp = _it1;
+ itmp++;
+ MElement *t1 = (*_it2)->t1;
+ MElement *t2 = (*_it2)->t2;
+ while (_touched.find(t1) != _touched.end() ||
+ _touched.find(t2) != _touched.end()) {
+ if (++itmp == _ite) return;
+ t1 = (*_it2)->t1;
+ t2 = (*_it2)->t2;
+ }
+
+ RecombSeq rs(itmp, _ite, _depth);
+
+ }*/
+
+ double getBenef() const {return _benef;}
+ //std::vector<int> getSkipped() const {return _skipped;}
+};
#endif