diff --git a/Geo/GModel.h b/Geo/GModel.h
index 77e53530e3a02795e7df361900da8d98ce9f6c50..ccd36d4a601da672e281ccb89a0cf83b433cb27e 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -158,10 +158,6 @@ class GModel
   std::set<int> meshPartitions;
   int partitionSize[2];
 
-  // boundary layer columns i.e. list of vertices that form columns
-  // in boundary layers
-  BoundaryLayerColumns _columns;
-
  public:
   GModel(std::string name="");
   virtual ~GModel();
diff --git a/Geo/GRegion.h b/Geo/GRegion.h
index 110c232942091fef97f86be4c1e0050a1da4af9a..2e1ccfb222909b00870e296b105a394650877695 100644
--- a/Geo/GRegion.h
+++ b/Geo/GRegion.h
@@ -11,6 +11,7 @@
 #include <vector>
 #include <stdio.h>
 #include "GEntity.h"
+#include "boundaryLayersData.h"
 
 class MElement;
 class MTetrahedron;
@@ -32,6 +33,7 @@ class GRegion : public GEntity {
   // replace faces (for gluing) for specific modelers, we have to
   // re-create internal data ...
   virtual void replaceFacesInternal (std::list<GFace*> &) {}
+  BoundaryLayerColumns _columns;
 
  public:
   GRegion(GModel *model, int tag);
@@ -135,6 +137,9 @@ class GRegion : public GEntity {
   void addPrism(MPrism *p){ prisms.push_back(p); }
   void addPyramid(MPyramid *p){ pyramids.push_back(p); }
   void addPolyhedron(MPolyhedron *p){ polyhedra.push_back(p); }
+
+  // get the boundary layer columns
+  BoundaryLayerColumns *getColumns () {return &_columns;}
 };
 
 #endif
diff --git a/Geo/boundaryLayersData.cpp b/Geo/boundaryLayersData.cpp
index 7bbed0a3712ea8e32803955a83ae9651325b8e16..0fb8ba97c74924747190d4bdc7a0bf1eb1e8c0e6 100644
--- a/Geo/boundaryLayersData.cpp
+++ b/Geo/boundaryLayersData.cpp
@@ -1161,15 +1161,15 @@ static bool preprocessVertex (MVertex *v,
   return onSymmetryPlane;
 }
 
-BoundaryLayerColumns *buildAdditionalPoints3D(GRegion *gr)
+bool buildAdditionalPoints3D(GRegion *gr)
 {
   BoundaryLayerField *blf = getBLField (gr->model());
 
-  if (!blf)return 0;
+  if (!blf) return false;
 
   blf->setupFor3d();
 
-  BoundaryLayerColumns * _columns = new BoundaryLayerColumns;
+  BoundaryLayerColumns *_columns = gr->getColumns();
 
   std::list<GFace*> faces = gr->faces();
   std::list<GFace*>::iterator itf = faces.begin();
@@ -1299,7 +1299,7 @@ BoundaryLayerColumns *buildAdditionalPoints3D(GRegion *gr)
 
   // END OF DEBUG STUFF
 
-  return _columns;
+  return true;
 }
 
 #endif
diff --git a/Geo/boundaryLayersData.h b/Geo/boundaryLayersData.h
index baf7529c9758b347fa61aa2c659c089634cb79e8..057361b1c97fb9b9aad69a7ee88243e5ca027c2e 100644
--- a/Geo/boundaryLayersData.h
+++ b/Geo/boundaryLayersData.h
@@ -148,7 +148,7 @@ public:
 
 BoundaryLayerField* getBLField(GModel *gm);
 bool buildAdditionalPoints2D (GFace *gf ) ;
-BoundaryLayerColumns * buildAdditionalPoints3D (GRegion *gr) ;
+bool buildAdditionalPoints3D (GRegion *gr) ;
 void buildMeshMetric(GFace *gf, double *uv, SMetric3 &m, double metric[3]);
 
 #endif
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 2b5ea1191811700c57d6386ec36610c1e4cb1e44..267c28d463752b56a9ebf7f2c6a865f85df99158 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -776,8 +776,8 @@ static int getWedge(BoundaryLayerColumns* _columns, MVertex *v1, MVertex *v2,
 static bool modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr)
 {
   if (getBLField(gr->model())) insertVerticesInRegion(gr,-1);
-  BoundaryLayerColumns* _columns = buildAdditionalPoints3D (gr);
-  if (!_columns)return false;
+  if (!buildAdditionalPoints3D (gr)) return false;
+  BoundaryLayerColumns* _columns = gr->getColumns();
   std::map<MFace,MElement*,Less_Face> bfaces;
 
   std::vector<MPrism*> blPrisms;