diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp
index ff994594..63b47bc7 100644
--- a/libsrc/meshing/basegeom.cpp
+++ b/libsrc/meshing/basegeom.cpp
@@ -454,7 +454,7 @@ namespace netgen
     mesh.LoadLocalMeshSize(mparam.meshsizefilename);
   }
 
-  void DivideEdge(GeometryEdge * edge, const MeshingParameters & mparam, const Mesh & mesh, Array<Point<3>> & points, Array<double> & params)
+  void DivideEdge(GeometryEdge * edge, const MeshingParameters & mparam, const Mesh & mesh, Array<Point<3>> & points, Array<double> & params, int layer)
   {
       static Timer tdivedgesections("Divide edge sections");
       static Timer tdivide("Divide Edges");
@@ -470,7 +470,7 @@ namespace netgen
       for(auto i : Range(divide_edge_sections))
       {
           auto pt = edge->GetPoint(double(i+1)/divide_edge_sections);
-          hvalue[i+1] = hvalue[i] + 1./mesh.GetH(pt) * (pt-old_pt).Length();
+          hvalue[i+1] = hvalue[i] + 1./mesh.GetH(pt, layer) * (pt-old_pt).Length();
           old_pt = pt;
       }
       int nsubedges = max2(1, int(floor(hvalue[divide_edge_sections]+0.5)));
@@ -525,7 +525,7 @@ namespace netgen
     std::map<size_t, PointIndex> vert2meshpt;
     for(auto & vert : vertices)
       {
-        auto pi = mesh.AddPoint(vert->GetPoint());
+        auto pi = mesh.AddPoint(vert->GetPoint(), vert->properties.layer);
         tree.Insert(mesh[pi], pi);
         vert2meshpt[vert->GetHash()] = pi;
         mesh[pi].Singularity(vert->properties.hpref);
@@ -591,7 +591,7 @@ namespace netgen
             }
             else
             {
-                DivideEdge(edge, mparam, mesh, edge_points, params);
+                DivideEdge(edge, mparam, mesh, edge_points, params, edge->properties.layer);
             }
         }
         else
@@ -635,7 +635,7 @@ namespace netgen
 
         for(auto i : Range(edge_points))
         {
-            auto pi = mesh.AddPoint(edge_points[i]);
+            auto pi = mesh.AddPoint(edge_points[i], edge->properties.layer);
             tree.Insert(mesh[pi], pi);
             pnums[i+1] = pi;
         }
@@ -724,7 +724,7 @@ namespace netgen
 
 
     static Timer t("GenerateMesh"); RegionTimer reg(t);
-    MESHING2_RESULT res = meshing.GenerateMesh(mesh, mparam, mparam.maxh, k+1);
+    MESHING2_RESULT res = meshing.GenerateMesh(mesh, mparam, mparam.maxh, k+1, face.properties.layer);
 
     for(auto i : Range(noldsurfels, mesh.GetNSE()))
       {
diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp
index 31eeabba..e0bf18b6 100644
--- a/libsrc/meshing/basegeom.hpp
+++ b/libsrc/meshing/basegeom.hpp
@@ -17,6 +17,7 @@ namespace netgen
     optional<Vec<4>> col;
     double maxh = 1e99;
     double hpref = 0;  // number of hp refinement levels (will be multiplied by factor later)
+    int layer = 1;
     optional<bool> quad_dominated;
     void Merge(const ShapeProperties & prop2)
     {
@@ -25,6 +26,7 @@ namespace netgen
       maxh = min2(maxh, prop2.maxh);
       hpref = max2(hpref, prop2.hpref);
       if(!quad_dominated.has_value()) quad_dominated = prop2.quad_dominated;
+      layer = max(layer, prop2.layer);
     }
 
     string GetName() const { return name ? *name : "default"; }
@@ -32,7 +34,7 @@ namespace netgen
 
     void DoArchive(Archive& ar)
     {
-        ar & name & col & maxh & hpref;
+        ar & name & col & maxh & hpref & layer;
     }
   };
 
@@ -51,6 +53,7 @@ namespace netgen
   {
   public:
     int nr = -1;
+    int layer = 1;
     ShapeProperties properties;
     Array<ShapeIdentification> identifications;
     GeometryShape * primary;
diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index 729554a5..855fe169 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -229,7 +229,7 @@ namespace netgen
     surfelementht = nullptr; 
     segmentht = nullptr;
 
-    lochfunc = nullptr;
+    lochfunc = {nullptr};
     // mglevels = 1;
     elementsearchtree = nullptr;
     elementsearchtreets = NextTimeStamp();
@@ -3264,7 +3264,7 @@ namespace netgen
 
 
 
-  void Mesh :: SetLocalH (netgen::Point<3> pmin, netgen::Point<3> pmax, double grading)
+  void Mesh :: SetLocalH (netgen::Point<3> pmin, netgen::Point<3> pmax, double grading, int layer)
   {
     using netgen::Point;
     Point<3> c = Center (pmin, pmax);
@@ -3275,30 +3275,30 @@ namespace netgen
     Point<3> pmin2 = c - Vec<3> (d, d, d);
     Point<3> pmax2 = c + Vec<3> (d, d, d);
 
-    lochfunc = make_unique<LocalH> (pmin2, pmax2, grading, dimension);
+    SetLocalH(make_unique<LocalH> (pmin2, pmax2, grading, dimension), layer);
   }
 
-  void Mesh :: RestrictLocalH (const Point3d & p, double hloc)
+  void Mesh :: RestrictLocalH (const Point3d & p, double hloc, int layer)
   {
     if(hloc < hmin)
       hloc = hmin;
 
     //cout << "restrict h in " << p << " to " << hloc << endl;
-    if (!lochfunc)
+    if (!lochfunc[layer-1])
       {
         PrintWarning("RestrictLocalH called, creating mesh-size tree");
 
         Point3d boxmin, boxmax;
         GetBox (boxmin, boxmax);
-        SetLocalH (boxmin, boxmax, 0.8);
+        SetLocalH (boxmin, boxmax, 0.8, layer);
       }
 
-    lochfunc -> SetH (p, hloc);
+    lochfunc[layer-1] -> SetH (p, hloc);
   }
 
   void Mesh :: RestrictLocalHLine (const Point3d & p1, 
                                    const Point3d & p2,
-                                   double hloc)
+                                   double hloc, int layer)
   {
     if(hloc < hmin)
       hloc = hmin;
@@ -3311,7 +3311,7 @@ namespace netgen
     for (i = 0; i <= steps; i++)
       {
         Point3d p = p1 + (double(i)/double(steps) * v);
-        RestrictLocalH (p, hloc);
+        RestrictLocalH (p, hloc, layer);
       }
   }
 
@@ -3343,24 +3343,26 @@ namespace netgen
   }
 
 
-  double Mesh :: GetH (const Point3d & p) const
+  double Mesh :: GetH (const Point3d & p, int layer) const
   {
+    const auto& lh = GetLocalH(layer);
     double hmin = hglob;
-    if (lochfunc)
+    if (lh)
       {
-        double hl = lochfunc->GetH (p);
+        double hl = lh->GetH (p);
         if (hl < hglob)
           hmin = hl;
       }
     return hmin;
   }
 
-  double Mesh :: GetMinH (const Point3d & pmin, const Point3d & pmax)
+  double Mesh :: GetMinH (const Point3d & pmin, const Point3d & pmax, int layer)
   {
+    const auto& lh = GetLocalH(layer);
     double hmin = hglob;
-    if (lochfunc)
+    if (lh)
       {
-        double hl = lochfunc->GetMinH (pmin, pmax);
+        double hl = lh->GetMinH (pmin, pmax);
         if (hl < hmin)
           hmin = hl;
       }
@@ -3404,16 +3406,16 @@ namespace netgen
 
 
 
-  void Mesh :: CalcLocalH (double grading) 
+  void Mesh :: CalcLocalH (double grading, int layer)
   {
     static Timer t("Mesh::CalcLocalH"); RegionTimer reg(t);
     
-    if (!lochfunc)
+    if (!lochfunc[layer-1])
       {
         Point3d pmin, pmax;
         GetBox (pmin, pmax);
         // SetLocalH (pmin, pmax, mparam.grading);
-	SetLocalH (pmin, pmax, grading);
+	SetLocalH (pmin, pmax, grading, layer);
       }
 
     PrintMessage (3,
@@ -3459,7 +3461,7 @@ namespace netgen
                 const Point3d & p1 = points[el.PNum(1)];
                 const Point3d & p2 = points[el.PNum(2)];
                 const Point3d & p3 = points[el.PNum(3)];
-                lochfunc->SetH (Center (p1, p2, p3), hel);
+                lochfunc[layer-1]->SetH (Center (p1, p2, p3), hel);
               }
           }
         else
@@ -3467,12 +3469,12 @@ namespace netgen
             {
               const Point3d & p1 = points[el.PNum(1)];
               const Point3d & p2 = points[el.PNum(2)];
-              lochfunc->SetH (Center (p1, p2), 2 * Dist (p1, p2));
+              lochfunc[layer-1]->SetH (Center (p1, p2), 2 * Dist (p1, p2));
             }
             {
               const Point3d & p1 = points[el.PNum(3)];
               const Point3d & p2 = points[el.PNum(4)];
-              lochfunc->SetH (Center (p1, p2), 2 * Dist (p1, p2));
+              lochfunc[layer-1]->SetH (Center (p1, p2), 2 * Dist (p1, p2));
             }
           }
       }
@@ -3490,7 +3492,7 @@ namespace netgen
         */
         if (!ident -> UsedSymmetric (seg[0], seg[1]))
           {
-            lochfunc->SetH (Center (p1, p2), Dist (p1, p2));
+            lochfunc[layer-1]->SetH (Center (p1, p2), Dist (p1, p2));
           }
       }
     /*
@@ -3540,17 +3542,17 @@ namespace netgen
   }
 
 
-  void Mesh :: CalcLocalHFromPointDistances(double grading)
+  void Mesh :: CalcLocalHFromPointDistances(double grading, int layer)
   {
     PrintMessage (3, "Calculating local h from point distances");
 
-    if (!lochfunc)
+    if (!lochfunc[layer-1])
       {
         Point3d pmin, pmax;
         GetBox (pmin, pmax);
 
         // SetLocalH (pmin, pmax, mparam.grading);
-	SetLocalH (pmin, pmax, grading);
+	SetLocalH (pmin, pmax, grading, layer);
       }
 
     PointIndex i,j;
@@ -3575,17 +3577,17 @@ namespace netgen
   }
 
 
-  void Mesh :: CalcLocalHFromSurfaceCurvature (double grading, double elperr) 
+  void Mesh :: CalcLocalHFromSurfaceCurvature (double grading, double elperr, int layer) 
   {
     PrintMessage (3, "Calculating local h from surface curvature");
 
-    if (!lochfunc)
+    if (!lochfunc[layer-1])
       {
         Point3d pmin, pmax;
         GetBox (pmin, pmax);
 
         // SetLocalH (pmin, pmax, mparam.grading);
-	SetLocalH (pmin, pmax, grading);
+	SetLocalH (pmin, pmax, grading, layer);
       }
 
 
@@ -3839,6 +3841,18 @@ namespace netgen
 
 
 
+  void Mesh :: SetLocalH(shared_ptr<LocalH> loch, int layer)
+  {
+      if(layer>lochfunc.Size())
+      {
+          auto pre_size = lochfunc.Size();
+          lochfunc.SetSize(layer);
+          for(auto & func : lochfunc.Range(pre_size, layer-1))
+              func = lochfunc[0];
+      }
+      lochfunc[layer-1] = loch;
+  }
+
   void Mesh :: GetBox (Point3d & pmin, Point3d & pmax, int dom) const
   {
     if (points.Size() == 0)
diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp
index 9346ffba..01197895 100644
--- a/libsrc/meshing/meshclass.hpp
+++ b/libsrc/meshing/meshclass.hpp
@@ -69,9 +69,9 @@ namespace netgen
 
 
     /**
-       Representation of local mesh-size h
+       Representation of local mesh-size h (one function per mesh layer)
     */
-    shared_ptr<LocalH> lochfunc;
+    Array<shared_ptr<LocalH>> lochfunc;
     ///
     double hglob;
     ///
@@ -435,18 +435,18 @@ namespace netgen
     */
     DLL_HEADER double AverageH (int surfnr = 0) const;
     /// Calculates localh 
-    DLL_HEADER void CalcLocalH (double grading);
+    DLL_HEADER void CalcLocalH (double grading, int layer=1);
     ///
-    DLL_HEADER void SetLocalH (netgen::Point<3> pmin, netgen::Point<3> pmax, double grading);
+    DLL_HEADER void SetLocalH (netgen::Point<3> pmin, netgen::Point<3> pmax, double grading, int layer=1);
     ///
-    DLL_HEADER void RestrictLocalH (const Point3d & p, double hloc);
+    DLL_HEADER void RestrictLocalH (const Point3d & p, double hloc, int layer=1);
     ///
     DLL_HEADER void RestrictLocalHLine (const Point3d & p1, const Point3d & p2, 
-			     double hloc);
+			     double hloc, int layer=1);
     /// number of elements per radius
-    DLL_HEADER void CalcLocalHFromSurfaceCurvature(double grading, double elperr);
+    DLL_HEADER void CalcLocalHFromSurfaceCurvature(double grading, double elperr, int layer=1);
     ///
-    DLL_HEADER void CalcLocalHFromPointDistances(double grading);
+    DLL_HEADER void CalcLocalHFromPointDistances(double grading, int layer=1);
     ///
     DLL_HEADER void RestrictLocalH (resthtype rht, int nr, double loch);
     ///
@@ -460,19 +460,25 @@ namespace netgen
     ///
 	DLL_HEADER void SetMaxHDomain (const NgArray<double> & mhd);
     ///
-    DLL_HEADER double GetH (const Point3d & p) const;
+    DLL_HEADER double GetH (const Point3d & p, int layer=1) const;
+    DLL_HEADER double GetH (PointIndex pi) const { return GetH(points[pi], points[pi].GetLayer()); }
     ///
-    double GetMinH (const Point3d & pmin, const Point3d & pmax);
+    double GetMinH (const Point3d & pmin, const Point3d & pmax, int layer=1);
     ///
-    bool HasLocalHFunction () { return lochfunc != nullptr; }
+    bool HasLocalHFunction (int layer=1) { return lochfunc[layer-1] != nullptr; }
     ///
-    LocalH & LocalHFunction () { return * lochfunc; }
+    LocalH & LocalHFunction (int layer=1) { return * lochfunc[layer-1]; }
 
-    shared_ptr<LocalH> GetLocalH() const { return lochfunc; }
-    void SetLocalH(shared_ptr<LocalH> loch) { lochfunc = loch; }
+    shared_ptr<LocalH> GetLocalH(int layer=1) const
+    {
+      if(lochfunc.Size() == 1)
+        return lochfunc[0];
+      return lochfunc[layer-1];
+    }
+    void SetLocalH(shared_ptr<LocalH> loch, int layer=1);
 
     ///
-    bool LocalHFunctionGenerated(void) const { return (lochfunc != NULL); }
+    bool LocalHFunctionGenerated(int layer=1) const { return (lochfunc[layer-1] != NULL); }
 
     /// Find bounding box
     DLL_HEADER void GetBox (Point3d & pmin, Point3d & pmax, int dom = -1) const;
diff --git a/libsrc/meshing/meshing2.cpp b/libsrc/meshing/meshing2.cpp
index c7d60e8d..6700b70c 100644
--- a/libsrc/meshing/meshing2.cpp
+++ b/libsrc/meshing/meshing2.cpp
@@ -241,7 +241,7 @@ namespace netgen
 
 
 
-  MESHING2_RESULT Meshing2 :: GenerateMesh (Mesh & mesh, const MeshingParameters & mp, double gh, int facenr)
+  MESHING2_RESULT Meshing2 :: GenerateMesh (Mesh & mesh, const MeshingParameters & mp, double gh, int facenr, int layer)
   {
     static Timer timer("surface meshing"); RegionTimer reg(timer);
 
@@ -477,7 +477,7 @@ namespace netgen
 	double his = Dist (p1, p2);
 
 	Point<3> pmid = Center (p1, p2);
-	double hshould = CalcLocalH (pmid, mesh.GetH (pmid));
+	double hshould = CalcLocalH (pmid, mesh.GetH (pmid, layer));
 	if (gh < hshould) hshould = gh;
 
 	mesh.RestrictLocalH (pmid, hshould);
diff --git a/libsrc/meshing/meshing2.hpp b/libsrc/meshing/meshing2.hpp
index 8eb233a4..a05cd3dd 100644
--- a/libsrc/meshing/meshing2.hpp
+++ b/libsrc/meshing/meshing2.hpp
@@ -59,7 +59,7 @@ public:
   void LoadRules (const char * filename, bool quad);
 
   /// 
-  DLL_HEADER MESHING2_RESULT GenerateMesh (Mesh & mesh, const MeshingParameters & mp, double gh, int facenr);
+  DLL_HEADER MESHING2_RESULT GenerateMesh (Mesh & mesh, const MeshingParameters & mp, double gh, int facenr, int layer=1);
 
   DLL_HEADER void Delaunay (Mesh & mesh, int domainnr, const MeshingParameters & mp);
   DLL_HEADER void BlockFillLocalH (Mesh & mesh, const MeshingParameters & mp);
diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp
index 9a1ac204..0fb605d0 100644
--- a/libsrc/meshing/meshtype.hpp
+++ b/libsrc/meshing/meshtype.hpp
@@ -1350,6 +1350,7 @@ namespace netgen
     public:
       Point<3> pnt;
       double h;
+      int layer = 1;
       MeshSizePoint (Point<3> _pnt, double _h) : pnt(_pnt), h(_h) { ; }
       MeshSizePoint () = default;
       MeshSizePoint (const MeshSizePoint &) = default;
diff --git a/libsrc/meshing/smoothing3.cpp b/libsrc/meshing/smoothing3.cpp
index 165f37f1..b30b8156 100644
--- a/libsrc/meshing/smoothing3.cpp
+++ b/libsrc/meshing/smoothing3.cpp
@@ -1243,8 +1243,8 @@ void Mesh :: ImproveMesh (const CSG eometry & geometry, OPTIMIZEGOAL goal)
 
       //      if (uselocalh)
       {
-	double lh = GetH (points.Get(i));
-	pf->SetLocalH (GetH (points.Get(i)));
+	double lh = GetH(points.Get(i), points.Get(i).GetLayer());
+	pf->SetLocalH (lh);
 	par.typx = lh / 10;
 	//	  (*testout) << "lh(" << points.Get(i) << ") = " << lh << "\n";
       }
@@ -1366,10 +1366,10 @@ void Mesh :: ImproveMeshSequential (const MeshingParameters & mp, OPTIMIZEGOAL g
 
   NgArray<double, PointIndex::BASE> pointh (points.Size());
 
-  if(lochfunc)
+  if(HasLocalHFunction())
     {
       for (PointIndex pi : points.Range())
-	pointh[pi] = GetH(points[pi]);
+	pointh[pi] = GetH(pi);
     }
   else
     {
@@ -1505,13 +1505,13 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
 
   Array<double, PointIndex> pointh (points.Size());
 
-  if(lochfunc)
+  if(HasLocalHFunction())
     {
       RegionTimer rt(tloch);
       ParallelForRange(points.Range(), [&] (auto myrange)
          {
            for(auto pi : myrange)
-             pointh[pi] = GetH(points[pi]);
+             pointh[pi] = GetH(pi);
          });
     }
   else
@@ -1640,11 +1640,11 @@ void Mesh :: ImproveMeshJacobian (const MeshingParameters & mp,
 
   NgArray<double, PointIndex::BASE, PointIndex> pointh (points.Size());
 
-  if(lochfunc)
+  if(HasLocalHFunction())
     {
       // for(i = 1; i<=points.Size(); i++)
       for (PointIndex pi : points.Range())
-	pointh[pi] = GetH(points[pi]);
+	pointh[pi] = GetH(pi);
     }
   else
     {
@@ -1797,11 +1797,11 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp,
 
   NgArray<double, PointIndex::BASE> pointh (points.Size());
  
-  if(lochfunc)
+  if(HasLocalHFunction())
     {
       // for(i=1; i<=points.Size(); i++)
       for (PointIndex pi : points.Range())
-	pointh[pi] = GetH(points[pi]);
+	pointh[pi] = GetH(pi);
     }
   else
     {
diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp
index 14da5c9a..d36b5d88 100644
--- a/libsrc/occ/occgenmesh.cpp
+++ b/libsrc/occ/occgenmesh.cpp
@@ -78,8 +78,7 @@ namespace netgen
 
 
   void RestrictHTriangle (gp_Pnt2d & par0, gp_Pnt2d & par1, gp_Pnt2d & par2,
-                          BRepLProp_SLProps * prop, BRepLProp_SLProps * prop2, Mesh & mesh, int depth, double h,
-                          const MeshingParameters & mparam)
+                          BRepLProp_SLProps * prop, BRepLProp_SLProps * prop2, Mesh & mesh, int depth, double h, int layer, const MeshingParameters & mparam)
   {
     int ls = -1;
 
@@ -190,20 +189,20 @@ namespace netgen
         if(ls == 0)
           {
             pm.SetX(0.5*(par1.X()+par2.X())); pm.SetY(0.5*(par1.Y()+par2.Y()));
-            RestrictHTriangle(pm, par2, par0, prop, prop2, mesh, depth+1, h, mparam);
-            RestrictHTriangle(pm, par0, par1, prop, prop2, mesh, depth+1, h, mparam);
+            RestrictHTriangle(pm, par2, par0, prop, prop2, mesh, depth+1, h, layer, mparam);
+            RestrictHTriangle(pm, par0, par1, prop, prop2, mesh, depth+1, h, layer, mparam);
           }
         else if(ls == 1)
           {
             pm.SetX(0.5*(par0.X()+par2.X())); pm.SetY(0.5*(par0.Y()+par2.Y()));
-            RestrictHTriangle(pm, par1, par2, prop, prop2, mesh, depth+1, h, mparam);
-            RestrictHTriangle(pm, par0, par1, prop, prop2, mesh, depth+1, h, mparam);
+            RestrictHTriangle(pm, par1, par2, prop, prop2, mesh, depth+1, h, layer, mparam);
+            RestrictHTriangle(pm, par0, par1, prop, prop2, mesh, depth+1, h, layer, mparam);
           }
         else if(ls == 2)
           {
             pm.SetX(0.5*(par0.X()+par1.X())); pm.SetY(0.5*(par0.Y()+par1.Y()));
-            RestrictHTriangle(pm, par1, par2, prop, prop2, mesh, depth+1, h, mparam);
-            RestrictHTriangle(pm, par2, par0, prop, prop2, mesh, depth+1, h, mparam);
+            RestrictHTriangle(pm, par1, par2, prop, prop2, mesh, depth+1, h, layer, mparam);
+            RestrictHTriangle(pm, par2, par0, prop, prop2, mesh, depth+1, h, layer, mparam);
           }
 
       }
@@ -215,16 +214,16 @@ namespace netgen
         prop->SetParameters (parmid.X(), parmid.Y());
         pnt = prop->Value();
         p3d = Point3d(pnt.X(), pnt.Y(), pnt.Z());
-        mesh.RestrictLocalH (p3d, h);
+        mesh.RestrictLocalH (p3d, h, layer);
 
         p3d = Point3d(pnt0.X(), pnt0.Y(), pnt0.Z());
-        mesh.RestrictLocalH (p3d, h);
+        mesh.RestrictLocalH (p3d, h, layer);
 
         p3d = Point3d(pnt1.X(), pnt1.Y(), pnt1.Z());
-        mesh.RestrictLocalH (p3d, h);
+        mesh.RestrictLocalH (p3d, h, layer);
 
         p3d = Point3d(pnt2.X(), pnt2.Y(), pnt2.Z());
-        mesh.RestrictLocalH (p3d, h);
+        mesh.RestrictLocalH (p3d, h, layer);
 
         //(*testout) << "p = " << p3d << ", h = " << h << ", maxside = " << maxside << endl;
 
@@ -375,6 +374,7 @@ namespace netgen
     //double maxh = mparam.maxh;
     //      int noldpoints = mesh->GetNP();
     int noldsurfel = mesh.GetNSE();
+    int layer = OCCGeometry::global_shape_properties[TopoDS::Face(geom.fmap(k)).TShape()].layer;
 
     static Timer tsurfprop("surfprop");
     tsurfprop.Start();
@@ -391,7 +391,7 @@ namespace netgen
     
     try {
       static Timer t("GenerateMesh"); RegionTimer reg(t);
-      res = meshing.GenerateMesh (mesh, mparam_without_overlap, maxh, k);
+      res = meshing.GenerateMesh (mesh, mparam_without_overlap, maxh, k, layer);
     }
 
     catch (SingularMatrixException)
@@ -437,24 +437,28 @@ namespace netgen
     NgArray<double> maxhdom;
     maxhdom.SetSize (geom.NrSolids());
     maxhdom = mparam.maxh;
+    int maxlayer = 1;
 
     int dom = 0;
     for (TopExp_Explorer e(geom.GetShape(), TopAbs_SOLID); e.More(); e.Next(), dom++)
+    {
       maxhdom[dom] = min2(maxhdom[dom], OCCGeometry::global_shape_properties[e.Current().TShape()].maxh);
+      maxlayer = max2(maxlayer, OCCGeometry::global_shape_properties[e.Current().TShape()].layer);
+    }
+
 
     mesh.SetMaxHDomain (maxhdom);
 
     Box<3> bb = geom.GetBoundingBox();
     bb.Increase (bb.Diam()/10);
 
-    mesh.SetLocalH (bb.PMin(), bb.PMax(), 0.5);
-
     if (mparam.uselocalh)
       {
         const char * savetask = multithread.task;
         multithread.percent = 0;
 
-        mesh.SetLocalH (bb.PMin(), bb.PMax(), mparam.grading);
+        for(auto layer : Range(1, maxlayer+1))
+            mesh.SetLocalH (bb.PMin(), bb.PMax(), mparam.grading, layer);
 
         int nedges = geom.emap.Extent();
 
@@ -482,6 +486,7 @@ namespace netgen
         for (int i = 1; i <= nedges && !multithread.terminate; i++)
           {
             TopoDS_Edge e = TopoDS::Edge (geom.emap(i));
+            int layer = OCCGeometry::global_shape_properties[e.TShape()].layer;
             multithread.percent = 100 * (i-1)/double(nedges);
             if (BRep_Tool::Degenerated(e)) continue;
 
@@ -540,7 +545,7 @@ namespace netgen
             for (int j = 0; j <= maxj; j++)
               {
                 gp_Pnt pnt = c->Value (s0+double(j)/maxj*(s1-s0));
-                mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), localh);
+                mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), localh, layer);
               }
           }
 
@@ -555,6 +560,7 @@ namespace netgen
             double maxcur = 0;
             multithread.percent = 100 * (i-1)/double(nedges);
             TopoDS_Edge edge = TopoDS::Edge (geom.emap(i));
+            int layer = OCCGeometry::global_shape_properties[edge.TShape()].layer;
             if (BRep_Tool::Degenerated(edge)) continue;
             double s0, s1;
             Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1);
@@ -575,7 +581,7 @@ namespace netgen
 
                 gp_Pnt pnt = c->Value (s);
 
-                mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), ComputeH (fabs(curvature), mparam));
+                mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), ComputeH (fabs(curvature), mparam), layer);
               }
           }
 
@@ -589,6 +595,7 @@ namespace netgen
           {
             multithread.percent = 100 * (i-1)/double(nfaces);
             TopoDS_Face face = TopoDS::Face(geom.fmap(i));
+            int layer = OCCGeometry::global_shape_properties[face.TShape()].layer;
             TopLoc_Location loc;
             Handle(Geom_Surface) surf = BRep_Tool::Surface (face);
             Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc);
@@ -628,7 +635,7 @@ namespace netgen
                 //maxside = max (maxside, p[1].Distance(p[2]));
                 //cout << "\rFace " << i << " pos11 ntriangles " << ntriangles << " maxside " << maxside << flush;
 
-                RestrictHTriangle (par[0], par[1], par[2], &prop, &prop2, mesh, 0, 0, mparam);
+                RestrictHTriangle (par[0], par[1], par[2], &prop, &prop2, mesh, 0, 0, layer, mparam);
                 //cout << "\rFace " << i << " pos12 ntriangles " << ntriangles << flush;
               }
           }
@@ -654,6 +661,7 @@ namespace netgen
             for (int i = 1; i <= nedges && !multithread.terminate; i++)
               {
                 TopoDS_Edge edge = TopoDS::Edge (geom.emap(i));
+                int layer = OCCGeometry::global_shape_properties[edge.TShape()].layer;
                 if (BRep_Tool::Degenerated(edge)) continue;
 
                 double s0, s1;
@@ -678,6 +686,7 @@ namespace netgen
                         gp_Pnt p1 = c->Value (s);
                         lines[nlines].p0 = Point<3> (p0.X(), p0.Y(), p0.Z());
                         lines[nlines].p1 = Point<3> (p1.X(), p1.Y(), p1.Z());
+                        lines[nlines].layer = layer;
 
                         Box3d box;
                         box.SetPoint (Point3d(lines[nlines].p0));
@@ -719,8 +728,8 @@ namespace netgen
                 Box3d box;
                 box.SetPoint (Point3d(line.p0));
                 box.AddPoint (Point3d(line.p1));
-                double maxhline = max (mesh.GetH(box.PMin()),
-                                       mesh.GetH(box.PMax()));
+                double maxhline = max (mesh.GetH(box.PMin(), line.layer),
+                                       mesh.GetH(box.PMax(), line.layer));
                 box.Increase(maxhline);
 
                 double mindist = 1e99;
@@ -731,6 +740,7 @@ namespace netgen
                   {
                     int num = linenums[j]-1;
                     if (i == num) continue;
+                    if (line.layer != lines[num].layer) continue;
                     if( is_identified_edge(edgenumber[i], edgenumber[num]) ) continue;
                     if ((line.p0-lines[num].p0).Length2() < 1e-15) continue;
                     if ((line.p0-lines[num].p1).Length2() < 1e-15) continue;
@@ -749,12 +759,12 @@ namespace netgen
                     mindist = 1e-3 * bb.Diam();
                   }
 
-                mesh.RestrictLocalHLine(line.p0, line.p1, mindist);
+                mesh.RestrictLocalHLine(line.p0, line.p1, mindist, line.layer);
               }
           }
 
         for (auto mspnt : mparam.meshsize_points)
-          mesh.RestrictLocalH(mspnt.pnt, mspnt.h);
+          mesh.RestrictLocalH(mspnt.pnt, mspnt.h, mspnt.layer);
 
         multithread.task = savetask;
 
diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp
index c36342aa..d266c69d 100644
--- a/libsrc/occ/occgeom.hpp
+++ b/libsrc/occ/occgeom.hpp
@@ -91,6 +91,7 @@ namespace netgen
   {
   public:
     Point<3> p0, p1;
+    int layer = 1;
     double Dist (Line l);
     double Length () { return (p1-p0).Length(); }
   };
diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp
index 86a5ebc2..4a4e021d 100644
--- a/libsrc/occ/python_occ_shapes.cpp
+++ b/libsrc/occ/python_occ_shapes.cpp
@@ -1456,14 +1456,29 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
   py::class_<TopoDS_Solid, TopoDS_Shape> (m, "Solid");
   
   py::class_<TopoDS_Compound, TopoDS_Shape> (m, "Compound")
-    .def(py::init([](std::vector<TopoDS_Shape> shapes) {
+    .def(py::init([](std::vector<TopoDS_Shape> shapes, bool separate_layers) {
           BRep_Builder builder;
           TopoDS_Compound comp;
           builder.MakeCompound(comp);
-          for(auto& s : shapes)
-            builder.Add(comp, s);
+          for(auto i : Range(shapes.size()))
+          {
+            builder.Add(comp, shapes[i]);
+            if(separate_layers)
+            {
+                auto & props = OCCGeometry::global_shape_properties;
+                for(auto & s : GetSolids(shapes[i]))
+                    props[s.TShape()].layer = i+1;
+                for(auto & s : GetFaces(shapes[i]))
+                    props[s.TShape()].layer = i+1;
+                for(auto & s : GetEdges(shapes[i]))
+                    props[s.TShape()].layer = i+1;
+                for(auto & s : GetVertices(shapes[i]))
+                    props[s.TShape()].layer = i+1;
+            }
+          }
+
           return comp;
-        }))
+        }), py::arg("shapes"), py::arg("separate_layers")=false)
     ;