only force even num of points for Blossom if we actually recombine something!

Mesh/meshGEdge.cpp

@@ -27,20 +27,16 @@ typedef struct {
 } IntPoint;
 
 
-static double smoothPrimitive (GEdge *ge,
-			       double alpha,
+static double smoothPrimitive(GEdge *ge, double alpha,
                               std::vector<IntPoint> &Points)
 {
-  //  printf("alpha = %g\n",alpha);
-
   int ITER = 0;
   while (1){
     int count1 = 0;
     int count2 = 0;
 
-    // use a gauss-seidel iteration
-    // iterate forward and then backward
-    // convergence is usually very fast.
+    // use a gauss-seidel iteration; iterate forward and then backward;
+    // convergence is usually very fast
     for (int i=1; i< Points.size(); i++){
       double dh = (Points[i].xp/Points[i].lc - Points[i-1].xp/Points[i-1].lc);
       double dt = Points[i].t - Points[i-1].t;
@@ -60,18 +56,14 @@ static double smoothPrimitive (GEdge *ge,
       	double hnew = Points[i].xp / Points[i].lc + dt * (alpha-1) * Points[i].xp;
       	Points[i-1].lc = Points[i-1].xp / hnew;
       	count2 ++;
-	//	dh = (Points[i-1].xp/Points[i-1].lc - Points[i].xp/Points[i].lc);
-	//	dt = fabs(Points[i-1].t - Points[i].t);
-	//	dhdt =  dh/dt;      
-	//	printf("dhdt / Points[i-1].xp
       }
     }
 
-    
     ++ITER;
-    if (ITER > 2000){break;}
+    if (ITER > 2000) break;
     if (!(count1 + count2)) break;
   }
+
   // recompute the primitive
   for (int i = 1; i < Points.size(); i++){
     IntPoint &pt2 = Points[i];
@@ -100,11 +92,6 @@ static double F_Lc(GEdge *ge, double t)
   SVector3 der = ge->firstDer(t);
   const double d = norm(der);
   return d / lc_here;
-
-  //  SMetric3 metric(1. / (lc_here * lc_here));
-  //  double lSquared = dot(der, metric, der);
-
-  //  return sqrt(lSquared);
 }
 
 static double F_Lc_aniso(GEdge *ge, double t)
@@ -124,17 +111,7 @@ static double F_Lc_aniso(GEdge *ge, double t)
     lc_here = BGM_MeshMetric(ge, t, 0, p.x(), p.y(), p.z());
 
   SVector3 der = ge->firstDer(t);
-
   double lSquared = dot(der, lc_here, der);
-
-  /*
-  if (ge->tag() == 3){
-    printf("%12.5E %12.5E\n",p.x(),1./sqrt(lSquared));
-  }
-  */
-  //  der.normalize();
-  //  printf("in the function %g n %g %g\n", sqrt(lSquared),der.x(),der.y());
-
   return sqrt(lSquared);
 }
 
@@ -389,34 +366,39 @@ void meshGEdge::operator() (GEdge *ge)
     N = ge->meshAttributes.nbPointsTransfinite;
   }
   else{
-    if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG || blf
-        /* FIXME || CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL_QUAD */) {
+    if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG || blf)
       a = Integration(ge, t_begin, t_end, F_Lc_aniso, Points,
                       CTX::instance()->mesh.lcIntegrationPrecision);
-      //      printFandPrimitive(ge->tag()+10000, Points);
-    }
     else
       a = Integration(ge, t_begin, t_end, F_Lc, Points,
                       CTX::instance()->mesh.lcIntegrationPrecision);
 
-    // FIXME JF : MAYBE WE SHOULD NOT ALWAYS SMOOTH THE 1D MESH SIZE FIELD ??
+    // we should maybe provide an option to disable the smoothing
     for (int i = 0; i < Points.size(); i++){
       IntPoint &pt = Points[i];
       SVector3 der = ge->firstDer(pt.t);
       pt.xp = der.norm();
     }
-    //printFandPrimitive(ge->tag(), Points);
-    a =  smoothPrimitive (ge,/*CTX::instance()->mesh.smoothRatio*/CTX::instance()->mesh.smoothRatio,Points);
-    //    printf(" coucou %g\n",a);
-    //printFandPrimitive(ge->tag()+10000, Points);
-    
+    a = smoothPrimitive(ge, CTX::instance()->mesh.smoothRatio, Points);
     N = std::max(ge->minimumMeshSegments() + 1, (int)(a + 1.));
   }
 
-  // force odd number of points for if blossom is used for
-  // recombination
+  // force odd number of points for if blossom is used for recombination
   if(ge->meshAttributes.Method != MESH_TRANSFINITE &&
-     CTX::instance()->mesh.algoRecombine == 1 && N % 2 == 0) N++;
+     CTX::instance()->mesh.algoRecombine == 1 && N % 2 == 0){
+    if(CTX::instance()->mesh.recombineAll){
+      N++;
+    }
+    else{
+      std::list<GFace*> faces = ge->faces();
+      for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end(); it++){
+        if((*it)->meshAttributes.recombine){
+          N++;
+          break;
+        }
+      }
+    }
+  }
 
   // printFandPrimitive(ge->tag(),Points);