basic 2D boundary layers (using old code)

Mesh/BoundaryLayers.cpp

@@ -4,6 +4,7 @@
 // bugs and problems to <gmsh@geuz.org>.
 
 #include "GModel.h"
+#include "MLine.h"
 #include "MTriangle.h"
 #include "MQuadrangle.h"
 #include "BoundaryLayers.h"
@@ -38,16 +39,17 @@ static void addExtrudeNormals(std::vector<T*> &elements, int invert,
   else{ // get extrusion data from Gouraud-shaded element normals
     for(unsigned int i = 0; i < elements.size(); i++){
       MElement *ele = elements[i];
-      for(int j = 0; j < ele->getNumFaces(); j++){
-        MFace fac = ele->getFace(j);
-        SVector3 n = fac.normal();
-        if(invert) n *= -1.;
-        if(n[0] || n[1] || n[2]){
-          double nn[3] = {n[0], n[1], n[2]};
-          for(int k = 0; k < fac.getNumVertices(); k++){
-            MVertex *v = fac.getVertex(k);
-            ExtrudeParams::normals->add(v->x(), v->y(), v->z(), 3, nn);
-          }
+      SVector3 n(0, 0, 0);
+      if(ele->getDim() == 2)
+        n = ele->getFace(0).normal();
+      else if(ele->getDim() == 1) // FIXME only valid in XY-plane
+        n = crossprod(ele->getEdge(0).tangent(), SVector3(0, 0, 1));
+      if(invert) n *= -1.;
+      if(n[0] || n[1] || n[2]){
+        double nn[3] = {n[0], n[1], n[2]};
+        for(int k = 0; k < ele->getNumVertices(); k++){
+          MVertex *v = ele->getVertex(k);
+          ExtrudeParams::normals->add(v->x(), v->y(), v->z(), 3, nn);
         }
       }
     }
@@ -58,8 +60,29 @@ int Mesh2DWithBoundaryLayers(GModel *m)
 {
   std::set<GFace*> sourceFaces, otherFaces;
   std::set<GEdge*> sourceEdges, otherEdges;
-  std::map<int, bool> sourceFaceInvert;
-  std::map<int, int> sourceUseView;
+  std::map<int, bool> sourceFaceInvert, sourceEdgeInvert;
+  std::map<int, int> sourceFaceUseView, sourceEdgeUseView;
+
+  // 2D boundary layers
+  for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++){
+    GEdge *ge = *it;
+    if(ge->getNativeType() == GEntity::GmshModel && 
+       ge->geomType() == GEntity::BoundaryLayerCurve){
+      ExtrudeParams *ep = ge->meshAttributes.extrude;
+      if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == COPIED_ENTITY){
+        GEdge *from = m->getEdgeByTag(std::abs(ep->geo.Source));
+        if(!from){
+          Msg::Error("Unknown source curve %d for boundary layer", ep->geo.Source);
+          return 0;
+        }
+        if(ep->geo.Source < 0) sourceEdgeInvert[from->tag()] = true;
+        if(ep->mesh.ViewIndex >= 0) sourceEdgeUseView[from->tag()] = ep->mesh.ViewIndex;
+        sourceEdges.insert(from);
+      }
+    }
+  }
+
+  // 3D boundary layers
   for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++){
     GFace *gf = *it;
     if(gf->getNativeType() == GEntity::GmshModel && 
@@ -72,14 +95,16 @@ int Mesh2DWithBoundaryLayers(GModel *m)
           return 0;
         }
         if(ep->geo.Source < 0) sourceFaceInvert[from->tag()] = true;
-        if(ep->mesh.ViewIndex >= 0) sourceUseView[from->tag()] = ep->mesh.ViewIndex;
+        if(ep->mesh.ViewIndex >= 0) sourceFaceUseView[from->tag()] = ep->mesh.ViewIndex;
         sourceFaces.insert(from);
         std::list<GEdge*> e = from->edges();
         sourceEdges.insert(e.begin(), e.end());
       }
     }
   }
-  if(sourceFaces.empty()) return 0;
+
+  if(sourceEdges.empty() && sourceFaces.empty())
+    return 0;
 
   // compute set of non-source edges and faces
   for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++)
@@ -96,27 +121,48 @@ int Mesh2DWithBoundaryLayers(GModel *m)
   // extrusion)
   std::for_each(sourceFaces.begin(), sourceFaces.end(), orientMeshGFace());
 
-  // compute a normal field for all the source faces
+  // compute a normal field for all the source edges or faces
   if(ExtrudeParams::normals) delete ExtrudeParams::normals;
   ExtrudeParams::normals = new smooth_data();
-  for(std::set<GFace*>::iterator it = sourceFaces.begin(); 
-      it != sourceFaces.end(); it++){
-    GFace *gf = *it;
-    int invert = sourceFaceInvert.count(gf->tag());
-    OctreePost *octree = 0;
+  
+  if(sourceFaces.empty()){
+    for(std::set<GEdge*>::iterator it = sourceEdges.begin(); 
+        it != sourceEdges.end(); it++){
+      GEdge *ge = *it;
+      int invert = sourceEdgeInvert.count(ge->tag());
+      OctreePost *octree = 0;
 #if defined(HAVE_POST)
-    if(sourceUseView.count(gf->tag())){
-      int index = sourceUseView[gf->tag()];
-      if(index >= 0 && index < PView::list.size())
-        octree = new OctreePost(PView::list[index]);
-      else
-        Msg::Error("Unknown View[%d]: using normals instead", index);
+      if(sourceEdgeUseView.count(ge->tag())){
+        int index = sourceEdgeUseView[ge->tag()];
+        if(index >= 0 && index < PView::list.size())
+          octree = new OctreePost(PView::list[index]);
+        else
+          Msg::Error("Unknown View[%d]: using normals instead", index);
+      }
+#endif
+      addExtrudeNormals(ge->lines, invert, octree);
     }
+  }
+  else{
+    for(std::set<GFace*>::iterator it = sourceFaces.begin(); 
+        it != sourceFaces.end(); it++){
+      GFace *gf = *it;
+      int invert = sourceFaceInvert.count(gf->tag());
+      OctreePost *octree = 0;
+#if defined(HAVE_POST)
+      if(sourceFaceUseView.count(gf->tag())){
+        int index = sourceFaceUseView[gf->tag()];
+        if(index >= 0 && index < PView::list.size())
+          octree = new OctreePost(PView::list[index]);
+        else
+          Msg::Error("Unknown View[%d]: using normals instead", index);
+      }
 #endif
-    addExtrudeNormals(gf->triangles, invert, octree);
-    addExtrudeNormals(gf->quadrangles, invert, octree);
+      addExtrudeNormals(gf->triangles, invert, octree);
+      addExtrudeNormals(gf->quadrangles, invert, octree);
+    }
   }
-  if(sourceUseView.empty())
+  if(sourceEdgeUseView.empty() && sourceFaceUseView.empty())
     ExtrudeParams::normals->normalize();
 
   // set the position of boundary layer points using the smooth normal

Mesh/meshGEdge.cpp

@@ -20,66 +20,6 @@ typedef struct {
   double t, lc, p;
 } IntPoint;
 
-struct xi2lc {
-  double xi, lc;
-  xi2lc(const double &_xi, const double _lc)
-    : xi(_xi), lc(_lc)
-  { 
-  }
-  bool operator < (const xi2lc &other) const
-  {
-    return xi < other.xi; 
-  }
-};
-
-static std::vector<xi2lc> interpLc;
-
-static void buildInterpLc(const std::vector<IntPoint> &lcPoints)
-{
-  IntPoint p;
-  interpLc.clear();
-  for(unsigned int i = 0; i < lcPoints.size(); i++){
-    p = lcPoints[i];
-    interpLc.push_back(xi2lc(p.t, p.lc));
-  }
-}
-
-static double F_Lc_usingInterpLc(GEdge *ge, double t)
-{
-  std::vector<xi2lc>::iterator it = std::lower_bound(interpLc.begin(),
-                                                     interpLc.end(), xi2lc(t, 0));
-  double t1 = it->xi;
-  double l1 = it->lc;
-  it++;
-  if(it == interpLc.end()) return l1;
-  double t2 = it->xi;
-  double l2 = it->lc;
-  double l = l1 + ((t - t1) / (t2 - t1)) * (l2 - l1);
-  return l;
-}
-
-static double F_Lc_usingInterpLcBis(GEdge *ge, double t)
-{
-  GPoint p = ge->point(t);
-  double lc_here;
-
-  Range<double> bounds = ge->parBounds(0);
-  double t_begin = bounds.low();
-  double t_end = bounds.high();
-
-  SVector3 der = ge->firstDer(t);
-  const double d = norm(der);
-
-  if(t == t_begin)
-    lc_here = BGM_MeshSize(ge->getBeginVertex(), t, 0, p.x(), p.y(), p.z());
-  else if(t == t_end)
-    lc_here = BGM_MeshSize(ge->getEndVertex(), t, 0, p.x(), p.y(), p.z());
-  else
-    lc_here = BGM_MeshSize(ge, t, 0, p.x(), p.y(), p.z());
-
-  return d / lc_here;
-}
-
 static double F_Lc(GEdge *ge, double t)
 {
   GPoint p = ge->point(t);
@@ -122,7 +62,7 @@ static double F_Lc_aniso(GEdge *ge, double t)
 
   SVector3 der = ge->firstDer(t);
 
-  double lSquared = dot (der,lc_here,der);
+  double lSquared = dot(der, lc_here, der);
 
   //  der.normalize();
   //  printf("in the function %g n %g %g\n", sqrt(lSquared),der.x(),der.y());
@@ -367,25 +307,18 @@ void meshGEdge::operator() (GEdge *ge)
     N = ge->meshAttributes.nbPointsTransfinite;
   }
   else{
-    /*if(CTX::instance()->mesh.lcIntegrationPrecision > 1.e-2){ JF says this code was a test but it is not useful
-      std::vector<IntPoint> lcPoints;
-      Integration(ge, t_begin, t_end, F_Lc_usingInterpLcBis, lcPoints, 
-                  CTX::instance()->mesh.lcIntegrationPrecision);
-      buildInterpLc(lcPoints);
-      a = Integration(ge, t_begin, t_end, F_Lc_usingInterpLc, Points, 1.e-8);
-    }
-    else*/{
-      if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG) 
-	a = Integration(ge, t_begin, t_end, F_Lc_aniso, Points,
-			CTX::instance()->mesh.lcIntegrationPrecision);
-      else
-	a = Integration(ge, t_begin, t_end, F_Lc, Points,
-			CTX::instance()->mesh.lcIntegrationPrecision);
-    }
+    if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG) 
+      a = Integration(ge, t_begin, t_end, F_Lc_aniso, Points,
+                      CTX::instance()->mesh.lcIntegrationPrecision);
+    else
+      a = Integration(ge, t_begin, t_end, F_Lc, Points,
+                      CTX::instance()->mesh.lcIntegrationPrecision);
     N = std::max(ge->minimumMeshSegments() + 1, (int)(a + 1.));
   }
-  
-  if(CTX::instance()->mesh.algoRecombine == 1 && N%2 == 0) N++;
+
+  // force odd number of points for if blossom is used for
+  // recombination
+  if(CTX::instance()->mesh.algoRecombine == 1 && N % 2 == 0) N++;
 
   //  printFandPrimitive(ge->tag(),Points);