diff --git a/CMakeLists.txt b/CMakeLists.txt
index d536474135c0bbdae1523d1a34ebe39577502ed9..d3be9d1a18d77da079083c2c4fb1d18bd3be7317 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -573,7 +573,7 @@ if(ENABLE_TAUCS)
   endif(NOT HAVE_METIS)
   find_library(TAUCS_LIB taucs PATH_SUFFIXES lib)
   if(TAUCS_LIB)
-    find_path(TAUCS_INC "taucs.h" PATH_SUFFIXES src include)
+    find_path(TAUCS_INC "taucs.h" PATH_SUFFIXES src include taucs)
     if(TAUCS_INC)
       set_config_option(HAVE_TAUCS "Taucs")
       add_definitions(-DTAUCS_CILK)
diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 81baf9f855638d886e416afcf41bcf789c5562bf..11a8a33eb8a134c03e7e4c512885275e3bf4f9f8 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -39,6 +39,7 @@ void GEdge::deleteMesh()
   mesh_vertices.clear();
   for(unsigned int i = 0; i < lines.size(); i++) delete lines[i];
   lines.clear();
+  _normals.clear();
   deleteVertexArrays();
   model()->destroyMeshCaches();
 }
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index 9efe5be4bdab72c9992da8f198870b0f0dec7cc1..f27407a8c3fcb59b18e7952e2a87397edeb1c80e 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -31,9 +31,11 @@ class GEdge : public GEntity {
 
  protected:
   GVertex *v0, *v1;
+  // FIXME: normals need to be mutable at the moment, because thay can
+  // be created in const member functions
+  mutable std::map<MVertex*, SVector3, std::less<MVertex*> > _normals;
   GEdgeCompound *compound; // this model edge belongs to a compound 
   std::list<GFace *> l_faces;
-  std::set<GFace *> bl_faces;
   // for specific solid modelers that need to re-do the internal curve
   // if a topological change ending points is done (gluing)
   virtual void replaceEndingPointsInternals(GVertex *, GVertex *) {}
@@ -153,13 +155,7 @@ class GEdge : public GEntity {
   // true if entity is periodic in the "dim" direction.
   virtual bool periodic(int dim) const { return v0 == v1; }
 
-  // true if edge is used in hyperbolic layer on face gf
-  virtual bool inHyperbolicLayer(GFace *gf)
-  {
-    return bl_faces.find(gf) != bl_faces.end();
-  }
-
-  virtual void flagInHyperbolicLayer(GFace *gf) { bl_faces.insert(gf); }
+  std::map<MVertex*, SVector3, std::less<MVertex*> > &getNormals() { return _normals; }
 
   // get bounds of parametric coordinate 
   virtual Range<double> parBounds(int i) const = 0;
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 35b9a7131df97c77c069dd2247c795664a99dff7..eb3f8959d652fe03c3ed46d40b333337288ccc80 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -387,6 +387,8 @@ class GModel
                    double angle);
   GEntity *extrude(GEntity *e, std::vector<double> p1, std::vector<double> p2);
   GEntity *addPipe(GEntity *e, std::vector<GEdge *>  edges);
+  void createBoundaryLayer(std::vector<GEntity *> e, double h);
+
   void addRuledFaces (std::vector<std::vector<GEdge *> > edges);
   GFace* addFace (std::vector<GEdge *> edges, std::vector< std::vector<double > > points);
   GFace* addPlanarFace (std::vector<std::vector<GEdge *> > edges);
diff --git a/Geo/GModelFactory.h b/Geo/GModelFactory.h
index 79172fedcab76e5e7a300af961e9660a2254d86d..b21647ccc4b9646bfcbedbac445d5d462802d0ef 100644
--- a/Geo/GModelFactory.h
+++ b/Geo/GModelFactory.h
@@ -15,6 +15,7 @@
 class GVertex;
 class GEdge;
 class GModel;
+class Field;
 
 // Abstract CAD creation factory.
 class GModelFactory {
diff --git a/Geo/boundaryLayerEdge.h b/Geo/boundaryLayerEdge.h
new file mode 100644
index 0000000000000000000000000000000000000000..6c1a42064e655e4149cf47aa03fe414ddcf484cc
--- /dev/null
+++ b/Geo/boundaryLayerEdge.h
@@ -0,0 +1,27 @@
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _BOUNDARY_LAYER_EDGE_H_
+#define _BOUNDARY_LAYER_EDGE_H_
+
+#include "discreteEdge.h"
+
+class Field;
+
+class boundaryLayerEdge : public discreteEdge {
+ private:
+  // original entity from which this part of the boundary layer is
+  // extruded
+  GEntity *_orig;
+ public:
+  boundaryLayerEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1, 
+                    GEntity *orig) 
+    : discreteEdge(model, num, _v0, _v1), _orig(orig)  {}
+  virtual ~boundaryLayerEdge() {}
+  virtual GeomType geomType() const { return BoundaryLayerCurve; }
+  GEntity *getOriginalEntity(){ return _orig; }
+};
+
+#endif
diff --git a/Geo/boundaryLayerFace.h b/Geo/boundaryLayerFace.h
new file mode 100644
index 0000000000000000000000000000000000000000..32ddd71e09093226923b6951734e27db58341fbe
--- /dev/null
+++ b/Geo/boundaryLayerFace.h
@@ -0,0 +1,33 @@
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _BOUNDARY_LAYER_FACE_H_
+#define _BOUNDARY_LAYER_FACE_H_
+
+#include "discreteFace.h"
+
+class Field;
+
+class boundaryLayerFace : public discreteFace {
+ private:
+  // list of entities from which the whole boundary layer is created
+  std::vector<GEntity*> _source;
+  // field providing extrude information
+  Field *_field;
+  // original entity from which this part of the boundary layer is
+  // extruded
+  GEntity *_orig;
+ public:
+  boundaryLayerFace(GModel *model, int num, std::vector<GEntity*> source, 
+                    Field *field, GEntity *orig) 
+    : discreteFace(model, num), _source(source), _field(field), _orig(orig){}
+  virtual ~boundaryLayerFace() {}
+  GEntity::GeomType geomType() const { return BoundaryLayerSurface; }
+  std::vector<GEntity*> &getSourceEntities(){ return _source; }
+  Field * getField(){ return _field; }
+  GEntity *getSourceEntity(){ return _orig; }
+};
+
+#endif
diff --git a/Geo/boundaryLayerVertex.h b/Geo/boundaryLayerVertex.h
new file mode 100644
index 0000000000000000000000000000000000000000..892ba3611f15a90519ca6d60edb0b1b30a2f35db
--- /dev/null
+++ b/Geo/boundaryLayerVertex.h
@@ -0,0 +1,29 @@
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _BOUNDARY_LAYER_VERTEX_H_
+#define _BOUNDARY_LAYER_VERTEX_H_
+
+#include "discreteVertex.h"
+
+class boundaryLayerVertex : public discreteVertex {
+ private:
+  // original entity from which this part of the boundary layer is
+  // extruded
+  GVertex *_orig;
+ public:
+  boundaryLayerVertex(GModel *m, int num, GVertex *orig) 
+    : discreteVertex(m, num), _orig(orig)
+  {
+    // don't wait for Mesh0D to create the mesh vertex, as it is used
+    // to remove duplicates during model creation (before any mesh
+    // generation is performed)
+    mesh_vertices.push_back(new MVertex(orig->x(), orig->y(), orig->z(), this));
+  }
+  virtual ~boundaryLayerVertex(){}
+  virtual GeomType geomType() const { return BoundaryLayerPoint; }
+};
+
+#endif
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 05108fbac30ad055b265ff1cf80e7588e7097030..c112f4df7d47da90ab358fec1f6b76985c99794b 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -391,7 +391,7 @@ void discreteEdge::computeNormals () const
 //   delete normals;
 }
 
-void discreteEdge::getLocalParameter(const double &t, int &iLine,
+bool discreteEdge::getLocalParameter(const double &t, int &iLine,
                                      double &tLoc) const
 {
   for (iLine = 0; iLine < (int)lines.size(); iLine++){
@@ -400,16 +400,17 @@ void discreteEdge::getLocalParameter(const double &t, int &iLine,
     if (t >= tmin && t <= tmax){
       tLoc = _orientation[iLine] ? (t-tmin)/(tmax-tmin) :
         1 - (t-tmin)/(tmax-tmin);
-      return;
+      return true;
     }
   }
+  return false;
 }
 
 GPoint discreteEdge::point(double par) const
 {
   double tLoc;
   int iEdge;
-  getLocalParameter(par,iEdge,tLoc);
+  if(!getLocalParameter(par,iEdge,tLoc)) return GPoint();
 
   double x, y, z;
   MVertex *vB = lines[iEdge]->getVertex(0);
diff --git a/Geo/discreteEdge.h b/Geo/discreteEdge.h
index 3177df7e404548ce3cf0334b1fe68c25715c7fad..5c475bc28bc051d712cdaeaf608d19e647952135 100644
--- a/Geo/discreteEdge.h
+++ b/Geo/discreteEdge.h
@@ -15,12 +15,11 @@ class discreteEdge : public GEdge {
   std::vector<double> _pars;
   std::vector<int> _orientation;
   std::map<MVertex*, MLine*> boundv;
-  mutable std::map<MVertex*, SVector3, std::less<MVertex*> > _normals;
   bool createdTopo;
  public:
   discreteEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1);
   virtual ~discreteEdge() {}
-  void getLocalParameter(const double &t, int &iEdge, double &tLoc) const;
+  bool getLocalParameter(const double &t, int &iEdge, double &tLoc) const;
   virtual GeomType geomType() const { return DiscreteCurve; }
   virtual GPoint point(double p) const;
   virtual SVector3 firstDer(double par) const;
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index f85b47ce99f075bec6bd02883870f50490ba7f06..5a56ea02efb1c4f85cb3ad0fa3117acc976266ae 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -53,12 +53,12 @@ void discreteFace::findEdges(std::map<MEdge, std::vector<int>, Less_Edge> &map_e
 
 void discreteFace::setBoundEdges(std::vector<int> tagEdges)
 {
- for (std::vector<int>::iterator it = tagEdges.begin();it != tagEdges.end(); it++){
-   GEdge *ge = GModel::current()->getEdgeByTag(*it);
-   l_edges.push_back(ge);
-   l_dirs.push_back(1);
-   ge->addFace(this);
- }
+  for (std::vector<int>::iterator it = tagEdges.begin(); it != tagEdges.end(); it++){
+    GEdge *ge = GModel::current()->getEdgeByTag(*it);
+    l_edges.push_back(ge);
+    l_dirs.push_back(1);
+    ge->addFace(this);
+  }
 }
 
 GPoint discreteFace::point(double par1, double par2) const
diff --git a/Mesh/BoundaryLayers.cpp b/Mesh/BoundaryLayers.cpp
index f5165b0c30a8fac35c1cbf20cac66f0de0bcb6f0..83fe34ea8e3097bf78b9213f79955a302f887a90 100644
--- a/Mesh/BoundaryLayers.cpp
+++ b/Mesh/BoundaryLayers.cpp
@@ -40,7 +40,8 @@ int Mesh2DWithBoundaryLayers(GModel *m)
   std::map<int, bool> sourceFaceInvert;
   for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++){
     GFace *gf = *it;
-    if(gf->geomType() == GEntity::BoundaryLayerSurface){
+    if(gf->getNativeType() == GEntity::GmshModel && 
+       gf->geomType() == GEntity::BoundaryLayerSurface){
       ExtrudeParams *ep = gf->meshAttributes.extrude;
       if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == COPIED_ENTITY){
         GFace *from = m->getFaceByTag(std::abs(ep->geo.Source));
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index ab0fc79380535c0c15db93b74067d21656836461..e3c98321b90c179876b7e1d082713c8eb606d881 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -7,22 +7,16 @@ set(SRC
   Generator.cpp
   Field.cpp
     meshGEdge.cpp 
-    meshGEdgeExtruded.cpp 
+      meshGEdgeExtruded.cpp
     meshGFace.cpp 
-    meshGFaceTransfinite.cpp 
-    meshGFaceExtruded.cpp 
-    meshGFaceBamg.cpp 
-    meshGFaceBDS.cpp 
-    meshGFaceDelaunayInsertion.cpp 
-    meshGFaceLloyd.cpp 
-    meshGFaceOptimize.cpp 
-    meshGFaceQuadrilateralize.cpp 
+      meshGFaceTransfinite.cpp meshGFaceExtruded.cpp 
+      meshGFaceBamg.cpp meshGFaceBDS.cpp meshGFaceDelaunayInsertion.cpp 
+      meshGFaceLloyd.cpp meshGFaceOptimize.cpp 
+      meshGFaceQuadrilateralize.cpp meshGFaceBoundaryLayer.cpp 
     meshGRegion.cpp 
-    meshGRegionDelaunayInsertion.cpp 
-    meshGRegionTransfinite.cpp 
-    meshGRegionExtruded.cpp 
-    meshGRegionCarveHole.cpp 
-    meshGRegionLocalMeshMod.cpp
+      meshGRegionDelaunayInsertion.cpp meshGRegionTransfinite.cpp 
+      meshGRegionExtruded.cpp  meshGRegionCarveHole.cpp 
+      meshGRegionLocalMeshMod.cpp
     BackgroundMesh.cpp 
     qualityMeasures.cpp 
     BoundaryLayers.cpp 
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 637e4182a0c3aa3142fcb96a94e5c54c9473eeb4..40eb85be8e9d85b6a05321b04fc8090feba8c1da 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1495,7 +1495,6 @@ void meshGFace::operator() (GFace *gf)
   }
 
   if(gf->geomType() == GEntity::DiscreteSurface) return;
-  if(gf->geomType() == GEntity::BoundaryLayerSurface) return;
   if(gf->geomType() == GEntity::ProjectionFace) return;
   if(gf->meshAttributes.Method == MESH_NONE) return;
   if(CTX::instance()->mesh.meshOnlyVisible && !gf->getVisibility()) return;
@@ -1504,6 +1503,7 @@ void meshGFace::operator() (GFace *gf)
   deMeshGFace dem;
   dem(gf);
 
+  if(MeshBoundaryLayerSurface(gf)) return;
   if(MeshTransfiniteSurface(gf)) return;
   if(MeshExtrudedSurface(gf)) return;
   if(gf->meshMaster() != gf->tag()){
diff --git a/Mesh/meshGFace.h b/Mesh/meshGFace.h
index 3ea931e09a526eb3932ad3f9a08d2c9491d4a7b2..09158465c33c48443bb02581039abd3d254a0e4c 100644
--- a/Mesh/meshGFace.h
+++ b/Mesh/meshGFace.h
@@ -45,6 +45,7 @@ void findTransfiniteCorners(GFace *gf, std::vector<MVertex*> &corners);
 int MeshTransfiniteSurface(GFace *gf);
 int MeshExtrudedSurface(GFace *gf, std::set<std::pair<MVertex*, MVertex*> > 
                         *constrainedEdges=0);
+int MeshBoundaryLayerSurface(GFace *gf);
 void partitionAndRemesh(GFaceCompound *gf);
 bool checkMeshCompound(GFaceCompound *gf, std::list<GEdge*> &edges);