diff --git a/Mesh/Voronoi3D.cpp b/Mesh/Voronoi3D.cpp index 500125efe3a2177c065d18aed535ec84e74ccbda..7f8933fc065987f7a20788bd1cde3726bb270997 100755 --- a/Mesh/Voronoi3D.cpp +++ b/Mesh/Voronoi3D.cpp @@ -1,4 +1,4 @@ -// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle +// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle // // See the LICENSE.txt file for license information. Please report all // bugs and problems to <gmsh@geuz.org>. diff --git a/Mesh/Voronoi3D.h b/Mesh/Voronoi3D.h index d8b28456a27018dc91a734857ec92d7c81ef22e2..861935ebc5a3fa0c2ff189c71d7f5ead41bf6d2d 100755 --- a/Mesh/Voronoi3D.h +++ b/Mesh/Voronoi3D.h @@ -1,4 +1,4 @@ -// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle +// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle // // See the LICENSE.txt file for license information. Please report all // bugs and problems to <gmsh@geuz.org>. @@ -18,4 +18,4 @@ class clip{ double max(double,double); int category(int,int,int,int); void print_segment(SPoint3,SPoint3,std::ofstream&); -}; \ No newline at end of file +}; diff --git a/Mesh/directions3D.cpp b/Mesh/directions3D.cpp index 65d1bd2aa9e317a810dfb0e9a240f62d36a19044..724466562d826a96c8da83b1f367d6da4c520f63 100755 --- a/Mesh/directions3D.cpp +++ b/Mesh/directions3D.cpp @@ -1,4 +1,4 @@ -// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle +// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle // // See the LICENSE.txt file for license information. Please report all // bugs and problems to <gmsh@geuz.org>. diff --git a/Mesh/directions3D.h b/Mesh/directions3D.h index 499cff24351c850c640dac08b5d208444feae935..317724ab7490b3dcac0d2ad67fca2377dce5a0fe 100755 --- a/Mesh/directions3D.h +++ b/Mesh/directions3D.h @@ -1,11 +1,11 @@ -// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle +// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle // // See the LICENSE.txt file for license information. Please report all // bugs and problems to <gmsh@geuz.org>. #include "GFace.h" #include "MElementOctree.h" -#if defined(HAVE_ANN) +#if defined(HAVE_ANN) #include <ANN/ANN.h> #endif @@ -23,14 +23,14 @@ class Matrix{ void set_m32(double); void set_m13(double); void set_m23(double); - void set_m33(double); + void set_m33(double); double get_m11(); double get_m21(); double get_m31(); - double get_m12(); + double get_m12(); double get_m22(); double get_m32(); - double get_m13(); + double get_m13(); double get_m23(); double get_m33(); }; @@ -55,4 +55,4 @@ class Frame_field{ static void print_field1(); static void print_field2(); static void clear(); -}; \ No newline at end of file +}; diff --git a/Mesh/meshRecombine2D.cpp b/Mesh/meshRecombine2D.cpp index 36f2e7a8318efb9b5431ce442f873a3bf0469e8b..9b3bdc1116fb1649411239e656beae5655b7ac8c 100644 --- a/Mesh/meshRecombine2D.cpp +++ b/Mesh/meshRecombine2D.cpp @@ -1,4 +1,4 @@ -// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle +// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle // // See the LICENSE.txt file for license information. Please report all // bugs and problems to <gmsh@geuz.org>. @@ -25,27 +25,27 @@ Recombine2D::Recombine2D(GFace *gf, int horizon) : _horizon(horizon), _gf(gf), _benef(.0), _applied(true) { if (_horizon < 1) _horizon = 1; - + //return; - + //laplaceSmoothing(gf,100); //gf->model()->writeMSH("befSquare.msh"); - + Msg::Info("Branching with horizon %d", _horizon = HORIZ); - + _haveParam = gf->geomType() != GEntity::DiscreteSurface || gf->getCompound(); if(!_haveParam) { Msg::Warning("Surfaces must have prametrization in order to recombine triangles"); //return; _haveParam=true; } - + e2t_cont adj; _buildEdgeToElement(gf->triangles, adj); - + std::map<MElement*, int> boundaryTriangle; std::pair<Set_Recomb::iterator, bool> ret; - + for (e2t_cont::iterator it = adj.begin(); it != adj.end(); ++it) { if (it->second.second) { RecombTriangle *rt = new RecombTriangle(it->first, it->second.first, it->second.second); @@ -64,16 +64,16 @@ Recombine2D::Recombine2D(GFace *gf, int horizon) boundaryTriangle[it->second.first] = 2; } } - + /*std::map<MElement*, int>::iterator itt; - + for (itt = boundaryTriangle.begin(); itt != boundaryTriangle.end(); ++itt) { RecombTriangle *rt; if (itt->second == 1) rt = new RecombTriangle(itt->first, -100); else rt = new RecombTriangle(itt->first, -200); - + ret = _pairs.insert(rt); while (!ret.second) { rt->touch(); @@ -81,13 +81,13 @@ Recombine2D::Recombine2D(GFace *gf, int horizon) } _possibleRecomb[itt->first].insert(rt); }*/ - + Msg::Info("poss %d", _possibleRecomb.size()); Msg::Info("pairs %d", _pairs.size()); Msg::Info("tri %d", gf->triangles.size()); - + _recombine(); - + Msg::Info("new %d ", NEW_NODE); Msg::Info("del %d ", DEL_NODE); _applied = false; @@ -121,18 +121,18 @@ void Recombine2D::_buildEdgeToElement(std::vector<E*> &elements, e2t_cont &adj) void Recombine2D::_recombine() { //SetBoundingBox(); - + /*Map_Tri_Recomb::iterator itt; - + for (itt = _possibleRecomb.begin(); itt != _possibleRecomb.end(); ++itt) { if }*/ - + int i = 0; while ((!_pairs.empty() || !_lastRecomb.empty()) && i < 1000) { Msg::Info("%d",i); Set_Recomb::iterator it = _bestNextRecombination(); - + /*if (it != _pairs.begin()) { std::set<MElement*> isolatedCP(_isolated); _isolated.clear(); @@ -158,21 +158,21 @@ void Recombine2D::_recombine() Msg::Info("I continue anyway :p"); _isolated = isolatedCP; }*/ - + _benef += (*it)->getBenef(); if ((*it)->isQuad()) _quads.push_back((*it)->createQuad()); else for (int i = 0; i < (*it)->numTriangles(); i++) _isolated.insert((*it)->triangle(i)); - + /*Set_Recomb::iterator itmp = _pairs.begin(); while (itmp != it) { _rmRT(*itmp); itmp = _pairs.begin(); }*/ _removeImpossible(it); - + i++; } Msg::Info("Done %d (%d)", i, _pairs.size()); @@ -182,23 +182,23 @@ Set_Recomb::iterator Recombine2D::_bestNextRecombination() { if (!_lastRecomb.empty()) return _lastRecomb.begin(); - - + + Map_Tri_Node touched; std::vector<Recomb2D_Node*> nodes(_horizon); Set_Recomb::iterator it = _pairs.begin(), itmax = _pairs.begin(); - + int h = 0, nSkip = 0; double maxBenef = _horizon * -200.; bool haveOptimal = false; int maxH = 0; - + int numb = 0; while (!haveOptimal) { numb++; while (h < _horizon && it != _pairs.end()) { bool possible = true; - + Map_Tri_Node::iterator itm; for (int i = 0; i < (*it)->numTriangles(); i++) { itm = touched.find((*it)->triangle(i)); @@ -207,7 +207,7 @@ Set_Recomb::iterator Recombine2D::_bestNextRecombination() possible = false; } } - + if (possible) { if (h == 0) nodes[h] = new Recomb2D_Node(it, nSkip, _horizon, touched); else nodes[h] = new Recomb2D_Node(it, nSkip, nodes[h-1]); @@ -217,23 +217,23 @@ NEW_NODE++; if (h >= _horizon || nodes[h - 1]->getBound() < maxBenef) break; } - + it++; } - + if (--h < 0) haveOptimal = true; - + if (nodes[h]->getTotBenef() > maxBenef && h >= maxH) { maxBenef = nodes[h]->getTotBenef(); itmax = nodes[0]->getItRecomb(); maxH = h; } - + /*nSkip = 0; while (h >= 0 && nSkip < 2) { nodes[h]->erase(); nSkip += nodes[h--]->getnSkip(); - + if (h == _horizon - 2 && nodes[h]->getnSkip() >= 3) nodes[h--]->erase(); //probleme si h=-1 while (h >= 0 && nodes[h]->beSkipped() && @@ -245,7 +245,7 @@ NEW_NODE++; DEL_NODE++; delete nodes[h--]; } - + if (h < 0) haveOptimal = true; else { it = ++nodes[h]->getItRecomb(); @@ -254,9 +254,9 @@ DEL_NODE++; DEL_NODE++; } } - + Msg::Info("num test %d", numb); - + return itmax; } @@ -264,10 +264,10 @@ void Recombine2D::_removeImpossible(Set_Recomb::iterator it) { std::vector<MElement *> t; (*it)->triangles(t); - + Map_Tri_Recomb::iterator itm; std::set<RecombTriangle*>::iterator itrt; - + for (int i = 0; i < t.size(); i++) // loop over touched triangles if ((itm = _possibleRecomb.find(t[i])) != _possibleRecomb.end()) { itrt = itm->second.begin(); @@ -278,7 +278,7 @@ void Recombine2D::_removeImpossible(Set_Recomb::iterator it) } _possibleRecomb.erase(itm); } - + RecombTriangle *rt = *it; if (!_lastRecomb.erase(*it)) _pairs.erase(it); @@ -289,7 +289,7 @@ void Recombine2D::_rmRT(RecombTriangle *rt, MElement *el) { MElement *tri; Map_Tri_Recomb::iterator itmtri; - + for (int i = 0; i < rt->numTriangles(); i++) if ((tri = rt->triangle(i)) != el && (itmtri = _possibleRecomb.find(tri)) != _possibleRecomb.end()) { @@ -305,7 +305,7 @@ void Recombine2D::_rmRT(RecombTriangle *rt, MElement *el) default :; } } - + _pairs.erase(rt); _lastRecomb.erase(rt); delete rt; @@ -314,7 +314,7 @@ void Recombine2D::_rmRT(RecombTriangle *rt, MElement *el) int Recombine2D::apply() { if (!_haveParam || _applied) return 0; - + std::vector<MTriangle*> triangles2; for (int i = 0; i < _gf->triangles.size(); i++) { if (_isolated.find(_gf->triangles[i]) == _isolated.end()) @@ -324,7 +324,7 @@ int Recombine2D::apply() } _gf->triangles = triangles2; _gf->quadrangles = _quads; - + _applied = true; //_gf->model()->writeMSH("recSquare.msh"); laplaceSmoothing(_gf,10); @@ -356,10 +356,10 @@ RecombTriangle::RecombTriangle(const MEdge &me, MElement *t1, MElement *t2) _angle = std::max(fabs(90. - a2), _angle); _angle = std::max(fabs(90. - a3), _angle); _angle = std::max(fabs(90. - a4), _angle); - + _benefit = std::max(1. - _angle/90., .0); //_benefit = compute_alignment(me,t1,t2); - + //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 diff --git a/Mesh/meshRecombine2D.h b/Mesh/meshRecombine2D.h index f0be64f1481cd2c7c1e4dee16143a40e618bb240..f21dcf936b1e80ed06a334d7949e8e8db6030635 100644 --- a/Mesh/meshRecombine2D.h +++ b/Mesh/meshRecombine2D.h @@ -1,4 +1,4 @@ -// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle +// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle // // See the LICENSE.txt file for license information. Please report all // bugs and problems to <gmsh@geuz.org>. @@ -45,16 +45,16 @@ typedef std::vector<setTriUnion> vectSetTriUnion; class list_actions { private : vectSetTriUnion _cont; - + public : list_actions(int h) {_cont.reserve(h);} - + class iterator { private : int pos; list_actions *la; setTriUnion::iterator it; - + public : void print() { Msg::Info("iterator(pos %d, la %d, tu %d)", pos, la, *it); @@ -119,7 +119,7 @@ class list_actions { iterator t = *this; int p = pos; ++it; - while (it == la->_end(p) && ++p != la->_size()) { + while (it == la->_end(p) && ++p != la->_size()) { ++pos; it = la->_begin(pos); } @@ -127,7 +127,7 @@ class list_actions { } inline int getPos() {return pos;} }; - + int add(setTriUnion &s) { setTriUnion::iterator it = s.begin(); while (it != s.end()) { @@ -159,12 +159,12 @@ class list_actions { return iterator(this, 0, _begin()); } iterator end() { - return iterator(this, _cont.size() - 1, _end()); + return iterator(this, _cont.size() - 1, _end()); } void pop_back() { - _cont.pop_back(); + _cont.pop_back(); } - + void sizes(){ switch (_cont.size()) { case 3 : @@ -183,7 +183,7 @@ class list_actions { //for (int i = 0; i < _cont.size(); ++i) // Msg::Info("[%d] %d", i+1, _cont[i].size()); } - + private : setTriUnion::iterator _end(int pos = -1) { if (_cont.size() == 0) { @@ -204,7 +204,7 @@ class list_actions { return _cont[pos].begin(); } inline int _size() {return _cont.size();} - + }; @@ -220,23 +220,23 @@ class Recombine2D { std::set<MElement*> _isolated; std::vector<MQuadrangle*> _quads; static int ra; - + template <class E> void _buildEdgeToElement(std::vector<E*> &, e2t_cont &); void _recombine(); Set_Recomb::iterator _bestNextRecombination(); void _removeImpossible(Set_Recomb::iterator); void _rmRT(RecombTriangle *, MElement *e = NULL); - + public : Recombine2D(GFace*, int horizon); Recombine2D(GFace*); ~Recombine2D(); - + int apply(); inline double getBenef() const {return _benef;} inline int numTriangle() const {return _isolated.size();} - + private : //std::set<TrianglesUnion*, lessTriUnion> _setQuads; std::list<TrianglesUnion*> _setQuads; @@ -245,7 +245,7 @@ class Recombine2D { mapTriUnion _lessActions; void _removeImpossible(TrianglesUnion*); void _recombine(bool); - + void _lookAhead(std::list<TrianglesUnion*> &, int horiz); void _lookAhead2(std::list<TrianglesUnion*> &, int horiz); void _getIncompatibles(const TrianglesUnion*, setTriUnion &); @@ -274,16 +274,16 @@ class RecombTriangle { double _benefit; bool _formingQuad; int _sign; - + public : RecombTriangle(const MEdge &, MElement *, MElement *); RecombTriangle(MElement *, double); - + double touch(); - + MQuadrangle *createQuad() const; bool operator<(const RecombTriangle &other) const; - + inline void triangles(std::vector<MElement *> &v) const {v = _t;} inline int numTriangles() const {return (int) _t.size();} MElement *triangle(int i) const { @@ -293,27 +293,27 @@ class RecombTriangle { } inline bool isQuad() const {return _formingQuad;} inline double getBenef() const {return _benefit;} - + double compute_alignment(const MEdge&, MElement*, MElement*); }; class Recomb2D_Node { - private : + private : double _benef, _totBenef; int _nskip, _depth; Set_Recomb::iterator _recomb; Map_Tri_Node *_touch; std::vector<Map_Tri_Node::iterator> _vit; std::set<int> _blocking; - + public : //Node(){} //Recomb2D_Node(Set_Recomb::iterator, Map_Tri_Node &, int, double ben = .0); - + Recomb2D_Node(Set_Recomb::iterator, int nskip, int depth, Map_Tri_Node &); Recomb2D_Node(Set_Recomb::iterator, int nskip, Recomb2D_Node *); ~Recomb2D_Node(); - + void erase(); void blocking(const Map_Tri_Node::iterator &); inline bool isBetter() {return _blocking.size() < 2;} @@ -333,17 +333,17 @@ class Rec2d_vertex { // GEntity : virtual Range<double> parBounds(int i) static double **_Vvalues; static double **_Vbenefs; - + public : Rec2d_vertex(MVertex *, std::list<MTriangle*> &, std::map<MVertex*, std::set<GEdge*> > &, bool); - + inline void changeNumEl(int c) {_numEl += c;} double getReward(); double getReward(int); //int onwhat() {return _onWhat;} //int numel() {return _numEl;} - + private : void _initStaticTable(); double _computeAngle(MVertex *, std::list<MTriangle*> &, std::set<GEdge*> &); @@ -356,17 +356,17 @@ class TrianglesUnion { Rec2d_vertex **_vertices; MTriangle **_triangles; static int _RECOMPUTE; //incremented when a recombination is executed - + public: static int _NumEdge, _NumVert; static double _ValEdge, _ValVert; private: - + public : TrianglesUnion(GFace *, std::list<MTriangle*> &, std::list<MEdge> &, std::list<Rec2d_vertex*> &, std::map<MVertex*,double> &, std::map<MEdge, std::list<MTriangle*>, Less_Edge>&); - + bool operator<(TrianglesUnion &other); void addTri(std::set<MTriangle*>&) const; void removeTri(std::set<MTriangle*>&) const; @@ -406,11 +406,11 @@ class TrianglesUnion { inline static double getActualReturn() { return _ValEdge/_NumEdge * _ValVert/_NumVert * _ValVert/_NumVert; } - + static void clear() {_NumEdge = 0; _NumVert = 0; _ValEdge = .0, _ValVert = .0;} - private: + private: double _computeEdgeLength(GFace*, MVertex**, double *u, double *v, int numIntermedPoint= 1); double _computeAlignment(const MEdge&, std::list<MTriangle*>&, diff --git a/Mesh/meshRecombine2D_2.cpp b/Mesh/meshRecombine2D_2.cpp index 3bd734bcc63420a67fbe79e4b7bc08f8f0e884eb..70a94f3bff98f151f10a9b4474b4bc6470ddd930 100644 --- a/Mesh/meshRecombine2D_2.cpp +++ b/Mesh/meshRecombine2D_2.cpp @@ -1,4 +1,4 @@ -// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle +// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle // // See the LICENSE.txt file for license information. Please report all // bugs and problems to <gmsh@geuz.org>. @@ -47,16 +47,16 @@ Recombine2D::Recombine2D(GFace *gf) { //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 ) @@ -65,7 +65,7 @@ Recombine2D::Recombine2D(GFace *gf) 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; @@ -80,7 +80,7 @@ Recombine2D::Recombine2D(GFace *gf) } } } - + // 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; @@ -100,10 +100,10 @@ Recombine2D::Recombine2D(GFace *gf) TrianglesUnion::_NumVert++; TrianglesUnion::_ValVert += n->getReward(); } - + // store mesh size for better performance std::map<MVertex*,double> mapV2LcUV; - + // Create 'TrianglesUnion' for pattern of 4 triangles /* +-----+ @@ -151,7 +151,7 @@ Recombine2D::Recombine2D(GFace *gf) _possibleActions[*itTri].insert(tu); _setQuads.push_back(tu); } - + // Create 'TrianglesUnion' for pattern of 2 triangles /* +---+ @@ -181,7 +181,7 @@ Recombine2D::Recombine2D(GFace *gf) } TrianglesUnion::addRec(); //_setQuads.sort(lessTriUnion()); - + if (SMOOTH) laplaceSmoothing(_gf,10); Msg::Error("%d", popopo); _recombine(true); @@ -192,11 +192,11 @@ Recombine2D::Recombine2D(GFace *gf) 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); @@ -210,7 +210,7 @@ void Recombine2D::_recombine(bool a) } vectTu.clear(); vectTu.insert(vectTu.begin(), tu); - + Msg::Info("------------------ %d", _setQuads.size()); _lookAhead(vectTu, HORIZ2); } @@ -254,7 +254,7 @@ void Recombine2D::_recombine(bool a) if (!Msg::GetAnswer("Continue ?", 1, "no", "yes")) Msg::Info("I continue anyway :p"); } - + ++i; } Msg::Info("Done %d (%d)", i, _pairs.size()); @@ -282,10 +282,10 @@ 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 @@ -294,22 +294,22 @@ void Recombine2D::_lookAhead(std::list<TrianglesUnion*> &candidate, int horiz) 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 @@ -338,7 +338,7 @@ void Recombine2D::_lookAhead(std::list<TrianglesUnion*> &candidate, int horiz) } } } - + double expectedReturn = TrianglesUnion::computeReturn(numbers, values, modifiedVert); if (maxReturn < expectedReturn) { int sizeSequence = current; @@ -370,7 +370,7 @@ void Recombine2D::_lookAhead(std::list<TrianglesUnion*> &candidate, int horiz) tu->getTriangle(j)->setCol(col); } } - + if (/*true ||*/ best /*|| (((int)(Cpu()*100) - t) % 5 == 0 && newt)*/) { // draw mesh newt = false; apply(false); @@ -387,7 +387,7 @@ void Recombine2D::_lookAhead(std::list<TrianglesUnion*> &candidate, int horiz) } /*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]; @@ -400,7 +400,7 @@ void Recombine2D::_lookAhead(std::list<TrianglesUnion*> &candidate, int horiz) } } } - + if (vecIt[current-1] != la.end() || --current > 0) { do { (*vecIt[current-1])->removeTri(setTri); @@ -415,7 +415,7 @@ void Recombine2D::_lookAhead(std::list<TrianglesUnion*> &candidate, int horiz) } } } - + { // color for (int i = 1; i < bestSequence.size(); ++i) { TrianglesUnion *tu = bestSequence[i].first; @@ -437,10 +437,10 @@ 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 @@ -448,19 +448,19 @@ void Recombine2D::_lookAhead2(std::list<TrianglesUnion*> &candidate, int horiz) 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(); @@ -468,7 +468,7 @@ void Recombine2D::_lookAhead2(std::list<TrianglesUnion*> &candidate, int horiz) bool haveSequence = false; int current = 1, currentLayer = 0, maxLayer = 0; int num=0; - + while (current > 0) { //Msg::Info(" "); bool best = false; @@ -526,7 +526,7 @@ void Recombine2D::_lookAhead2(std::list<TrianglesUnion*> &candidate, int horiz) } 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(); @@ -548,7 +548,7 @@ Msg::Info(" ");*/ } } } - + double expectedReturn = TrianglesUnion::computeReturn(numbers, values, modifiedVert); if (maxReturn < expectedReturn) { int sizeSequence = current; @@ -590,7 +590,7 @@ Msg::Info(" ");*/ tu->getTriangle(j)->setCol(col); } } - + if (/*true ||*/ best /*|| (((int)(Cpu()*100) - t) % 5 == 0 && newt)*/) { // draw mesh //newt = false; apply(false); @@ -610,7 +610,7 @@ Msg::Info(" ");*/ } /*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]; @@ -623,7 +623,7 @@ Msg::Info(" ");*/ } } } - + if (vecIt[current-1] != la.end() || --current > 0) { do { while (current-1 < numLayer[currentLayer]) @@ -678,7 +678,7 @@ vecIt[current-1].print(); Msg::Info(" "); }*/ } - + /*{ // color for (int i = 1; i < bestSequence.size(); ++i) { TrianglesUnion *tu = bestSequence[i].first; @@ -760,17 +760,17 @@ double Recombine2D::_realisticReturn(const int *num, 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()) { @@ -791,9 +791,9 @@ double Recombine2D::_realisticReturn(const int *num, ++numPossible; } } - + setTri.erase(itTri); - + if (numPossible < 1){ ++numIsolatedTri; Msg::Info("'m here"); @@ -809,11 +809,11 @@ double Recombine2D::_realisticReturn(const int *num, 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) { @@ -822,9 +822,9 @@ double Recombine2D::_realisticReturn(const int *num, 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); @@ -847,12 +847,12 @@ double Recombine2D::_realisticReturn(const int *num, } 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; } @@ -863,7 +863,7 @@ int Recombine2D::_checkIsolatedTri(const std::vector<list_actions::iterator> &ve 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) { @@ -871,9 +871,9 @@ int Recombine2D::_checkIsolatedTri(const std::vector<list_actions::iterator> &ve 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) { @@ -914,7 +914,7 @@ int Recombine2D::apply(bool a) _gf->quadrangles = _quads; _isolated.clear(); } - + _applied = true; return 1; } @@ -924,12 +924,12 @@ 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) { @@ -940,7 +940,7 @@ void Recombine2D::_removeImpossible(TrianglesUnion *tu) incomp.insert(*it2); } } - + setTriUnion::iterator itTu = incomp.begin(); for (; itTu != incomp.end(); ++itTu) { for (int i = 0; i < (*itTu)->getNumTriangles(); i++) { @@ -963,8 +963,8 @@ void Recombine2D::_removeImpossible(TrianglesUnion *tu) } } //_setQuads.erase(*itTu); //marche pas avec list - } - + } + std::list<TrianglesUnion*>::iterator it = _setQuads.begin(); while (it != _setQuads.end()) { if (incomp.find(*it) != incomp.end()) @@ -972,11 +972,11 @@ void Recombine2D::_removeImpossible(TrianglesUnion *tu) 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; @@ -998,7 +998,7 @@ Rec2d_vertex::Rec2d_vertex(MVertex *v, { _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()) { @@ -1009,7 +1009,7 @@ Rec2d_vertex::Rec2d_vertex(MVertex *v, _angle = M_PI / 2.; } } - + if (_Vvalues == NULL) _initStaticTable(); } @@ -1020,7 +1020,7 @@ void Rec2d_vertex::_initStaticTable() // _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++) @@ -1028,7 +1028,7 @@ void Rec2d_vertex::_initStaticTable() _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++) @@ -1049,10 +1049,10 @@ double Rec2d_vertex::_computeAngle(MVertex *v, 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; @@ -1072,13 +1072,13 @@ double Rec2d_vertex::_computeAngle(MVertex *v, 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++) { @@ -1092,7 +1092,7 @@ double Rec2d_vertex::_computeAngle(MVertex *v, } angleMesh += angle3Vertices(t->getVertex((k+1)%3), v, t->getVertex((k+2)%3)); } - + if (angleMesh < M_PI) return angle1; return angle2; @@ -1136,7 +1136,7 @@ TrianglesUnion::TrianglesUnion(GFace *gf, _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; @@ -1146,7 +1146,7 @@ TrianglesUnion::TrianglesUnion(GFace *gf, setEdge.insert((*itt)->getEdge(i)); } } - + std::list<MEdge>::iterator ite = e.begin(); for (; ite != e.end(); ite++) { double sumlc = .0, u[2], v[2]; @@ -1162,7 +1162,7 @@ TrianglesUnion::TrianglesUnion(GFace *gf, else sumlc += itlc->second; } - + sumlc = .2; // FIXME BGM_MeshSize returns wrong meshsize double length = _computeEdgeLength(gf, vert, u, v, 0); @@ -1170,7 +1170,7 @@ TrianglesUnion::TrianglesUnion(GFace *gf, _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]; @@ -1186,7 +1186,7 @@ TrianglesUnion::TrianglesUnion(GFace *gf, else sumlc += itlc->second; } - + sumlc = .2; // FIXME BGM_MeshSize returns wrong meshsize double length = _computeEdgeLength(gf, vert, u, v, 0); @@ -1196,14 +1196,14 @@ TrianglesUnion::TrianglesUnion(GFace *gf, } _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; } @@ -1242,7 +1242,7 @@ double TrianglesUnion::_computeAlignment(const MEdge &e, std::list<MTriangle*> & 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]; @@ -1286,7 +1286,7 @@ void TrianglesUnion::_update() alpha /= _NumVert - _numEmbVert; double beta = (_ValEdge - _embEdgeValue + _boundEdgeValue) / (_NumEdge - _numEmbEdge + _numboundEdge); _globValIfExecuted = alpha * alpha * beta; - + _computation = _RECOMPUTE; } @@ -1318,7 +1318,7 @@ MQuadrangle *TrianglesUnion::createQuad() const } if (listBoundEdge.size() != 4) Msg::Warning("[meshRecombine2D] Not 4 edges !"); - + MVertex *v1, *v2, *v3, *v4; v1 = listBoundEdge.begin()->getMinVertex(); v2 = listBoundEdge.begin()->getMaxVertex(); @@ -1345,7 +1345,7 @@ MQuadrangle *TrianglesUnion::createQuad() const break; } } - + return new MQuadrangle(v1, v2, v3, v4); } @@ -1381,7 +1381,7 @@ void TrianglesUnion::addInfluence(int *num, double *val, std::map<Rec2d_vertex*, if (mapVert.find(_vertices[i]) == mapVert.end()) mapVert[_vertices[i]] = -1; else - mapVert[_vertices[i]] -= 1; + mapVert[_vertices[i]] -= 1; } } @@ -1394,7 +1394,7 @@ void TrianglesUnion::removeInfluence(int *num, double *val, std::map<Rec2d_verte val[1] -= _embEdgeValue; val[1] += _boundEdgeValue; for (int i = 0; i < _numBoundVert; ++i) { - mapVert[_vertices[i]] += 1; + mapVert[_vertices[i]] += 1; } } @@ -1416,7 +1416,7 @@ bool TrianglesUnion::operator<(TrianglesUnion &other) { _update(); other._update(); - return _globValIfExecuted < other._globValIfExecuted; + return _globValIfExecuted < other._globValIfExecuted; } bool lessTriUnion::operator()(TrianglesUnion *tu1, TrianglesUnion *tu2) const @@ -1438,10 +1438,10 @@ construction : - list of vertex - set of recomb - map tri to recomb - + destruction : - vertices, recomb - + execution : - take best recomb - determine recomb to execute diff --git a/Mesh/periodical.cpp b/Mesh/periodical.cpp index 6335956b9c8151235ee43e1534f063e7204c221a..9e28645033ef8be07d1b1ab209895b83b4a8fbca 100755 --- a/Mesh/periodical.cpp +++ b/Mesh/periodical.cpp @@ -1,4 +1,4 @@ -// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle +// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle // // See the LICENSE.txt file for license information. Please report all // bugs and problems to <gmsh@geuz.org>. diff --git a/Mesh/periodical.h b/Mesh/periodical.h index f40fdb84239d46de39259d710c8542f373838b76..2a167342d5703ad1bbad110f898c7e45e0ec86a5 100755 --- a/Mesh/periodical.h +++ b/Mesh/periodical.h @@ -1,4 +1,4 @@ -// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle +// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle // // See the LICENSE.txt file for license information. Please report all // bugs and problems to <gmsh@geuz.org>. @@ -24,4 +24,4 @@ class voroMetal3D{ void print_geo_volume(int,int,std::ofstream&); void print_geo_line_loop(int,std::vector<int>&,std::vector<int>&,std::ofstream&); void print_geo_face_loop(int,std::vector<int>&,std::ofstream&); -}; \ No newline at end of file +}; diff --git a/Mesh/simple3D.cpp b/Mesh/simple3D.cpp index 4fd499293cdbc1d5385813f295ed979b5b564ab6..80d2a4e3d82a9b5e4442f3ed7621aa0bb0a945ed 100755 --- a/Mesh/simple3D.cpp +++ b/Mesh/simple3D.cpp @@ -1,4 +1,4 @@ -// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle +// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle // // See the LICENSE.txt file for license information. Please report all // bugs and problems to <gmsh@geuz.org>. diff --git a/Mesh/simple3D.h b/Mesh/simple3D.h index 8d0538aae438450812584533db415666f7891f92..1086e518694d204af67c83516c716ac63bf0c6c9 100755 --- a/Mesh/simple3D.h +++ b/Mesh/simple3D.h @@ -1,4 +1,4 @@ -// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle +// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle // // See the LICENSE.txt file for license information. Please report all // bugs and problems to <gmsh@geuz.org>. @@ -32,4 +32,4 @@ class Filler{ void treat_region(GRegion*); static int get_nbr_new_vertices(); static MVertex* get_new_vertex(int); -}; \ No newline at end of file +};