diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp index 9a4f7d7ceb7f0010a5de48c7274efbebf851e6af..48ab19d1bfef3a6cd1942fdee952dead1fdbcdd6 100644 --- a/Mesh/meshGFaceOptimize.cpp +++ b/Mesh/meshGFaceOptimize.cpp @@ -964,7 +964,7 @@ static void _relocateVertex (GFace *gf, factor *= 0.5; } - if (success){ + if (success){ // ne devrait-on pas utiliser newp ici au lieu de after ??? (vu la boucle precedente) ver->setParameter(0, after.x()); ver->setParameter(1, after.y()); GPoint pt = gf->point(after); @@ -1522,7 +1522,7 @@ void _triangleSplit (GFace *gf, MElement *t, bool swop = false) { } -/*struct RecombineTriangle +struct RecombineTriangle { MElement *t1, *t2; double angle; @@ -1561,7 +1561,7 @@ void _triangleSplit (GFace *gf, MElement *t, bool swop = false) { //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) { @@ -1863,12 +1863,6 @@ 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; @@ -1928,204 +1922,6 @@ 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; @@ -2339,136 +2135,3 @@ 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 93bcf5c92429a01a6e14a26212e67dff787f39bd..55b2d34e371aef4b2dcfe9d4c7a66723f4a8af78 100644 --- a/Mesh/meshGFaceOptimize.h +++ b/Mesh/meshGFaceOptimize.h @@ -147,111 +147,4 @@ 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