2D boundary layers now are ablt take into account axis of symmetry

Mesh/Field.cpp

@@ -2249,18 +2249,54 @@ void BoundaryLayerField::removeAttractors()
   update_needed = true;
 }
 
+void BoundaryLayerField::setupFor1d(int iE) {
+
+  if (faces_id_saved.empty() &&
+      edges_id_saved.empty() &&
+      faces_id_saved.empty() ){     
+    faces_id_saved = faces_id;
+    edges_id_saved = edges_id;
+    nodes_id_saved = nodes_id;
+  }
+
+  nodes_id.clear();
+  edges_id.clear();
+  faces_id.clear();
+
+  bool found = std::find(edges_id_saved.begin(), edges_id_saved.end(), iE) !=
+    edges_id_saved.end();
+  
+  if (!found) {
+    GEdge *ge = GModel::current()->getEdgeByTag(iE);
+    GVertex *gv0 = ge->getBeginVertex();
+    found = std::find(nodes_id_saved.begin(), nodes_id_saved.end(), gv0->tag()) !=
+      nodes_id_saved.end();
+    if (found)nodes_id.push_back(gv0->tag());
+    GVertex *gv1 = ge->getEndVertex();
+    found = std::find(nodes_id_saved.begin(), nodes_id_saved.end(), gv1->tag()) !=
+      nodes_id_saved.end();
+    if (found)nodes_id.push_back(gv1->tag());
+  }
+  //  printf("edge %d %d nodes added\n",iE,nodes_id.size());
+  //  getchar();
+  removeAttractors();
+}
+
+
 void BoundaryLayerField::setupFor2d(int iF)
 {
-  if (1 || faces_id.size()){
   /* preprocess data in the following way
      remove GFaces from the attarctors (only used in 2D)
      for edges and vertices
   */
-    if ( !faces_id_saved.size()){
+  if (faces_id_saved.empty() &&
+      edges_id_saved.empty() &&
+      faces_id_saved.empty() ){     
     faces_id_saved = faces_id;
     edges_id_saved = edges_id;
     nodes_id_saved = nodes_id;
   }
+
   nodes_id.clear();
   edges_id.clear();
   faces_id.clear();
@@ -2313,7 +2349,6 @@ void BoundaryLayerField::setupFor2d(int iF)
     // printf("face %d %d BL Edges\n", iF, (int)edges_id.size());
   removeAttractors();
 }
-}
 
 void BoundaryLayerField::setupFor3d()
 {
@@ -2339,7 +2374,9 @@ double BoundaryLayerField::operator() (double x, double y, double z, GEntity *ge
     update_needed = false;
   }
 
+
   double dist = 1.e22;
+  if (_att_fields.empty())return dist;
   //  AttractorField *cc;
   for (std::list<AttractorField*>::iterator it = _att_fields.begin();
        it != _att_fields.end(); ++it){
@@ -2349,11 +2386,15 @@ double BoundaryLayerField::operator() (double x, double y, double z, GEntity *ge
       dist = cdist;
     }
   }
+
+  if (dist > thickness*ratio) return 1.e22;
   current_distance = dist;
   //  const double dist = (*field) (x, y, z);
   //  current_distance = dist;
-  const double lc = dist*(ratio-1) + hwall_t;
-  //    double lc = hwall * pow (ratio, dist / hwall);
+  double lc = dist*(ratio-1) + hwall_n;
+
+  //  double lc =  hwall_n;
+  //  double lc = hwall_n * pow (ratio, dist / hwall_t);  
   return std::min (hfar,lc);
 }
 

Mesh/Field.h

@@ -182,6 +182,7 @@ class BoundaryLayerField : public Field {
     return std::find(nodes_id.begin(),nodes_id.end(),iV) != nodes_id.end();
   }
   void computeFor1dMesh(double x, double y, double z, SMetric3 &metr);
+  void setupFor1d(int iE);
   void setupFor2d(int iF);
   void setupFor3d();
   void removeAttractors();

Mesh/meshGEdge.cpp

@@ -78,12 +78,34 @@ static double smoothPrimitive(GEdge *ge, double alpha,
 
 static double F_LcB(GEdge *ge, double t)
 {
+  BoundaryLayerField *blf = 0;
+#if defined(HAVE_ANN)
+  FieldManager *fields = ge->model()->getFields();
+  Field *bl_field = fields->get(fields->getBoundaryLayerField());
+  blf = dynamic_cast<BoundaryLayerField*> (bl_field);
+#endif
+
   GPoint p = ge->point(t);
-  return BGM_MeshSize(ge, t, 0, p.x(), p.y(), p.z());
+  double lc = BGM_MeshSize(ge, t, 0, p.x(), p.y(), p.z());
+
+  if (blf){
+    double lc2 = (*blf)( p.x(), p.y(), p.z() , ge);    
+    //    printf("p %g %g lc %g\n",p.x(),p.y(),lc2);
+    lc = std::min(lc, lc2);
+  }
+
+  return lc;
 }
 
 static double F_Lc(GEdge *ge, double t)
 {
+  BoundaryLayerField *blf = 0;
+#if defined(HAVE_ANN)
+  FieldManager *fields = ge->model()->getFields();
+  Field *bl_field = fields->get(fields->getBoundaryLayerField());
+  blf = dynamic_cast<BoundaryLayerField*> (bl_field);
+#endif
+  
   GPoint p = ge->point(t);
   double lc_here;
 
@@ -98,6 +120,12 @@ static double F_Lc(GEdge *ge, double t)
   else
     lc_here = BGM_MeshSize(ge, t, 0, p.x(), p.y(), p.z());
 
+  if (blf){
+    double lc2 = (*blf)( p.x(), p.y(), p.z() , ge);    
+    //    printf("p %g %g lc %g\n",p.x(),p.y(),lc2);
+    lc_here = std::min(lc_here, lc2);
+  }
+
   SVector3 der = ge->firstDer(t);
   const double d = norm(der);
 
@@ -115,6 +143,8 @@ static double F_Lc_aniso(GEdge *ge, double t)
   bool blf = false;
 #endif
 
+  printf("coucou\n");
+
 
   GPoint p = ge->point(t);
   SMetric3 lc_here;
@@ -397,11 +427,12 @@ static void filterPoints (GEdge*ge) {
 
 void meshGEdge::operator() (GEdge *ge)
 {
+  BoundaryLayerField *blf = 0;
 #if defined(HAVE_ANN)
   FieldManager *fields = ge->model()->getFields();
-  BoundaryLayerField *blf = 0;
   Field *bl_field = fields->get(fields->getBoundaryLayerField());
   blf = dynamic_cast<BoundaryLayerField*> (bl_field);
+  if (blf)blf->setupFor1d(ge->tag());
 #else
   bool blf = false;
 #endif
@@ -476,7 +507,7 @@ void meshGEdge::operator() (GEdge *ge)
       N /= CTX::instance()->mesh.lcFactor;
   }
   else{
-    if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG || blf){
+    if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG/* || blf*/){
       a = Integration(ge, t_begin, t_end, F_Lc_aniso, Points,
                       CTX::instance()->mesh.lcIntegrationPrecision);
     }
@@ -491,7 +522,7 @@ void meshGEdge::operator() (GEdge *ge)
       SVector3 der = ge->firstDer(pt.t);
       pt.xp = der.norm();
     }
-    a = smoothPrimitive(ge, sqrt(CTX::instance()->mesh.smoothRatio), Points);
+    //    a = smoothPrimitive(ge, sqrt(CTX::instance()->mesh.smoothRatio), Points);
     N = std::max(ge->minimumMeshSegments() + 1, (int)(a + 1.99));
   }
 
@@ -517,6 +548,7 @@ void meshGEdge::operator() (GEdge *ge)
     }
   }
 
+
   //printFandPrimitive(ge->tag(),Points);
 
   // if the curve is periodic and if the begin vertex is identical to