diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 2889e567e615565eab25f619d065621d4ea7fcc2..d3063533e87ab485457069699c244d8374260d39 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -1829,10 +1829,6 @@ bool GFaceCompound::checkTopology() const
Msg::Info("-----------------------------------------------------------");
Msg::Info("--- Split surface %d in 2 parts with Laplacian Mesh partitioner", tag());
}
-
-// correctTopo = true;
-// _mapping = MULTISCALE;
-// Msg::Warning("Aspect Ratio (AR=%d) is too high: using multiscale Laplace", AR);
}
else{
Msg::Debug("Correct topology: Genus=%d and Nb boundaries=%d, AR=%g", G, Nb, H/D);
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 738a268894adb8db36184c8cfb6f5ef8055f066b..04cc9e13f0a9fd3af94949671135d9262ff336d4 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1141,8 +1141,6 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces)
it != discFaces.end(); it++)
(*it)->findEdges(map_edges);
-
-
//return if no boundary edges (torus, sphere, ...)
if (map_edges.empty()) return;
@@ -1266,7 +1264,6 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces)
// set boundary edges for each face
for (std::vector<discreteFace*>::iterator it = discFaces.begin(); it != discFaces.end(); it++){
- //printf("set boundary edge for face =%d \n", (*it)->tag());
std::map<int, std::vector<int> >::iterator ite = face2Edges.find((*it)->tag());
if (ite != face2Edges.end()){
std::vector<int> bcEdges = ite->second;
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 63ce4c731a90e442bbf7035f78fe6da0f4f2ed89..4ddf81c4d6b3bc274b85a45305619d0fdd41d694 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -34,6 +34,7 @@ void discreteFace::findEdges(std::map<MEdge, std::vector<int>, Less_Edge> &map_e
else bound_edges.erase(itset);
}
}
+
// for the boundary edges, associate the tag of the current discrete face
for (std::set<MEdge, Less_Edge>::iterator itv = bound_edges.begin();
itv != bound_edges.end(); ++itv){
@@ -49,6 +50,7 @@ void discreteFace::findEdges(std::map<MEdge, std::vector<int>, Less_Edge> &map_e
itmap->second = tagFaces;
}
}
+
}
void discreteFace::setBoundEdges(std::vector<int> tagEdges)
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 857ca25cd78e1331c4e71f4d1ab35c1740008253..ec602f3363c656e30cffff1a73fe839d0dd2185d 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -253,8 +253,7 @@ BDS_Edge *BDS_Mesh::recover_edge_fast(BDS_Point *p1, BDS_Point *p2){
BDS_Point *p2b = f->oppositeVertex(e);
if (p2 == p2b){
if (swap_edge(e, BDS_SwapEdgeTestQuality(false,false))){
- printf("coucou\n");
- return find_edge (p1,p2->iD);
+ return find_edge (p1,p2->iD);
}
}
}
diff --git a/Mesh/highOrderSmoother.cpp b/Mesh/highOrderSmoother.cpp
index ee7a2f9dba10f7ffa9e1bfc53be2cce8292bd550..a18006c36cfa923297319a94b41b841279e4af3c 100644
--- a/Mesh/highOrderSmoother.cpp
+++ b/Mesh/highOrderSmoother.cpp
@@ -1270,18 +1270,16 @@ static int findOptimalLocationsP2(GFace *gf, highOrderSmoother *s)
static int findOptimalLocationsPN(GFace *gf,highOrderSmoother *s)
{
- printf("coucou1\n");
+
e2t_cont adj;
buildEdgeToTriangle(gf->triangles, adj);
int N=0;
- printf("coucou2\n");
for (e2t_cont::iterator it = adj.begin(); it!= adj.end(); ++it){
if (it->second.second)
N += optimalLocationPN_(gf,it->first, dynamic_cast<MTriangle*>(it->second.first),
dynamic_cast<MTriangle*>(it->second.second),s);
}
- printf("coucou3\n");
return N;
}
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index c8989a61abd510bf90c3378da91c801260807500..541f4f54122fa72a071330e1fb91c491c4a27d4c 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -252,25 +252,27 @@ static void recur_compute_centers_ (double R, double a1, double a2,
centers.push_back(std::make_pair(PL,zero));
centers.push_back(std::make_pair(PR,zero));
- //std::map<SPoint2, double> CR;
+ std::vector<SPoint2> centersChild;
+ centersChild.clear();
for (int i=0;i<root->children.size();i++){
multiscaleLaplaceLevel* m = root->children[i];
- if( !m) printf("no child \n");
centers.push_back(std::make_pair(m->center,m));
m->radius = 0.0;
for (std::map<MVertex*,SPoint2>::iterator it = m->coordinates.begin();
it != m->coordinates.end() ; ++it){
- const SPoint2 &p = it->second;
+ SPoint2 p = it->second;
m->radius = std::max(m->radius,sqrt ((m->center.x() - p.x())*(m->center.x() - p.x())+
(m->center.y() - p.y())*(m->center.y() - p.y())));
}
- //CR.insert(std::make_pair(m->center, m->radius));
+ centersChild.push_back(m->center);
}
+
//add the center of real holes ...
std::vector<std::vector<MEdge> > boundaries;
connected_bounds(root->elements, boundaries);
- if (root->children.size()==0 ){ // || boundaries.size()-1 != root->children.size()){
+ int toadd = 0;
+ if (root->children.size()==0 || boundaries.size()-1 != root->children.size() ){
for (unsigned int i = 0; i < boundaries.size(); i++){
std::vector<MEdge> me = boundaries[i];
SPoint2 c(0.0,0.0);
@@ -287,41 +289,29 @@ static void recur_compute_centers_ (double R, double a1, double a2,
(c.y() - p.y())*(c.y() - p.y())));
}
- //check if the center has not been added
-// bool newCenter = false;
-// for (std::map<SPoint2,double>::iterator it = CR.begin(); it != CR.end(); it++){
-// SPoint2 p = it->first;
-// double radius = it->second;
-// double dist = sqrt ((c.x() - p.x())*(c.x() - p.x())+ (c.y() - p.y())*(c.y() - p.y()));
-// if (dist > radius) newCenter = true;
-// }
- if (std::abs(rad/root->radius) < 0.9 && std::abs(rad) < 0.99){
+ //check if the center has not been added
+ bool newCenter = true;
+ for (std::vector<SPoint2>::iterator it2 = centersChild.begin(); it2 != centersChild.end(); it2++){
+ SPoint2 p = *it2;
+ double dist = sqrt ((c.x() - p.x())*(c.x() - p.x())+
+ (c.y() - p.y())*(c.y() - p.y()));
+ if (dist < 0.7*rad) newCenter = false;
+ }
+
+ if (std::abs(rad/root->radius) < 0.6 && std::abs(rad) < 0.95 && newCenter){
+ toadd++;
centers.push_back(std::make_pair(c,zero));
}
}
}
+ if (toadd != boundaries.size()-1-root->children.size() ) {
+ printf("!!!!!!!! ARG added =%d != %d\n", toadd, boundaries.size()-1- root->children.size());
+ //exit(1);
+ }
//sort centers
std::sort(centers.begin(),centers.end(), sort_pred(PL,PR));
- int nbrecur = 0;
- for (int i=1;i<centers.size()-1;i++){
- multiscaleLaplaceLevel* m1 = centers[i-1].second;
- multiscaleLaplaceLevel* m2 = centers[i].second;
- multiscaleLaplaceLevel* m3 = centers[i+1].second;
- if (m2) {
- nbrecur ++;
- }
- }
- if (root->children.size()!= nbrecur) {
- printf("children =%d nbrecur =%d \n", root->children.size(), nbrecur);
- for (int i=0;i<centers.size();i++){
- printf("center i=%d (%g %g)\n", i, centers[i].first.x(), centers[i].first.x());
- if(centers[i].second) printf(" level %d recur %d elems =%d\n", centers[i].second->recur, centers[i].second->region, centers[i].second->elements.size() );
- }
- }
-
-
for (int i=1;i<centers.size()-1;i++){
multiscaleLaplaceLevel* m1 = centers[i-1].second;
multiscaleLaplaceLevel* m2 = centers[i].second;
@@ -530,9 +520,7 @@ static void connectedRegions (std::vector<MElement*> &elements,
static void recur_cut_ (double R, double a1, double a2,
multiscaleLaplaceLevel * root,
std::vector<MElement *> &left,
- std::vector<MElement *> &right , int &totNbElems){
-
- totNbElems += root->elements.size();
+ std::vector<MElement *> &right){
SPoint2 PL (R*cos(a1),R*sin(a1));
SPoint2 PR (R*cos(a2),R*sin(a2));
@@ -567,7 +555,7 @@ static void recur_cut_ (double R, double a1, double a2,
if (m2){
a1 = myatan2 (centers[i-1].first.y() - m2->center.y() , centers[i-1].first.x() - m2->center.x() );
a2 = myatan2 (centers[i+1].first.y() - m2->center.y() , centers[i+1].first.x() - m2->center.x() );
- recur_cut_ (m2->radius, a1, a2, m2, left, right, totNbElems);
+ recur_cut_ (m2->radius, a1, a2, m2, left, right);
}
}
}
@@ -822,10 +810,8 @@ multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements,
root->scale = 1.0;
root->_name = "Root";
- int totNbElems = 0;
- parametrize(*root, totNbElems);
- printf("PARAM: elements =%d, totNbElems = %d \n", elements.size(), totNbElems);
-
+ parametrize(*root);
+
//fill the coordinates
std::vector<double> iScale;
std::vector<SPoint2> iCenter;
@@ -922,7 +908,7 @@ void multiscaleLaplace::fillCoordinates (multiscaleLaplaceLevel & level,
}
-void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level, int &totNbElems){
+void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
//Compute all nodes for the level
@@ -1022,8 +1008,6 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level, int &totNbE
}
}
- totNbElems += level.elements.size();
-
//Save multiscale meshes
std::string name1(level._name+"real.msh");
std::string name2(level._name+"param.msh");
@@ -1068,7 +1052,7 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level, int &totNbE
if (!tooSmallv.empty()){
Msg::Info("Level (%d-%d) Multiscale Laplace (reg[%d] = %d too small)",level.recur,level.region, i, tooSmallv.size());
level.children.push_back(nextLevel);
- parametrize (*nextLevel, totNbElems);
+ parametrize (*nextLevel);
}
}
@@ -1133,25 +1117,21 @@ void multiscaleLaplace::cut (std::vector<MElement *> &elements)
std::vector<MElement*> left,right;
int totNbElems = 0;
- recur_cut_ (1.0, M_PI, 0.0, root,left,right, totNbElems);
+ recur_cut_ (1.0, M_PI, 0.0, root,left,right);
+ connected_left_right(left, right);
+
if ( elements.size() != left.size()+right.size()) {
- Msg::Error("Cutting laplace wrong nb elements (%d) != left + right (%d) TOTRECUR=%d", elements.size(), left.size()+right.size(), totNbElems);
+ Msg::Error("Cutting laplace wrong nb elements (%d) != left + right (%d)", elements.size(), left.size()+right.size());
exit(1);
}
- connected_left_right(left, right);
-
-
- printf("left.size() = %d, right.size() = %d, nbElem = %d elements.size() = %d\n",
- left.size(),right.size(),totNbElems,elements.size());
-
- printLevel ("cut-left.msh",left,0,2.2);
- printLevel ("cut-right.msh",right,0,2.2);
elements.clear();
elements.insert(elements.end(),left.begin(),left.end());
elements.insert(elements.end(),right.begin(),right.end());
- printLevel ("cut-all.msh",elements, 0,2.2);
+ printLevel ("Rootcut-left.msh",left,0,2.2);
+ printLevel ("Rootcut-right.msh",right,0,2.2);
+ printLevel ("Rootcut-all.msh",elements, 0,2.2);
}
diff --git a/Solver/multiscaleLaplace.h b/Solver/multiscaleLaplace.h
index d587dfd9ec30c5c0eb97ae634358fa46ea02de30..5e205190019714b61dd1930036fdeb588913b5b7 100644
--- a/Solver/multiscaleLaplace.h
+++ b/Solver/multiscaleLaplace.h
@@ -36,7 +36,7 @@ public:
std::map<MVertex*, SPoint3> &allCoordinates,
std::vector<double> &iScale,
std::vector<SPoint2> &iCenter);
- void parametrize (multiscaleLaplaceLevel &, int &totNbElems);
+ void parametrize (multiscaleLaplaceLevel &);
void parametrize_method (multiscaleLaplaceLevel & level,
std::set<MVertex*> &allNodes,
std::map<MVertex*,SPoint2> &solution,