diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp
index bc005da4..d3bf1ccc 100644
--- a/libsrc/meshing/boundarylayer.cpp
+++ b/libsrc/meshing/boundarylayer.cpp
@@ -584,8 +584,6 @@ struct GrowthVectorLimiter {
 
     for(auto i : Range(growthvectors))
       growthvectors[i] *= limits[i];
-  cout << "limits " << limits << endl;
-  cout << "old np " << np << endl;
 
 }
 
@@ -1541,7 +1539,6 @@ struct GrowthVectorLimiter {
 
     auto & identifications = mesh.GetIdentifications();
     const int identnr = identifications.GetNr("boundarylayer");
-    identifications.SetType(identnr, Identifications::CLOSESURFACES);
 
     auto add_points = [&](PointIndex pi, Vec<3> & growth_vector, Array<PointIndex> & new_points)
     {
@@ -1553,13 +1550,8 @@ struct GrowthVectorLimiter {
           auto pi_new = mesh.AddPoint(p);
           new_points.Append(pi_new);
           growth_vector_map[pi_new] = { &growth_vector, params.heights[i] };
-          if(special_boundary_points.count(pi) == 0)
-            mesh.GetIdentifications().Add(pi_last, pi_new, identnr);
-          else
-          {
-            cout << "add locked point " << pi_new << endl;
+          if(special_boundary_points.count(pi) > 0)
             mesh.AddLockedPoint(pi_new);
-          }
           pi_last = pi_new;
         }
     };
@@ -1793,6 +1785,14 @@ struct GrowthVectorLimiter {
                 el.SetIndex(new_mat_nrs[sel.GetIndex()]);
                 if(add_volume_element)
                   mesh.AddVolumeElement(el);
+                else
+                {
+                  // Let the volume mesher fill the hole with pyramids/tets
+                  // To insert pyramids, we need close surface identifications on open quads
+                  for(auto i : Range(points))
+                    if(numGroups(sel[i]) == 1)
+                      identifications.Add(el[i], el[i+points.Size()], identnr);
+                }
               }
             Element2d newel = sel;
             for(auto i: Range(points))
diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index e3891cda..5685a17b 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -4081,7 +4081,6 @@ namespace netgen
           }
         else
           {
-            cout << "unused point " << pi << endl;
             op2np[pi].Invalidate(); 
           }
       
@@ -4141,11 +4140,6 @@ namespace netgen
     for (int i = 0; i < lockedpoints.Size(); i++)
       lockedpoints[i] = op2np[lockedpoints[i]];
 
-  cout << "op2np " << endl;
-  for(auto i : Range(op2np)) {
-    if(i!=op2np[i])
-    cout << i << '\t'  << op2np[i] << endl;
-  }
     GetIdentifications().MapPoints(op2np);
     /*
     for (int i = 0; i < facedecoding.Size(); i++)
diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp
index 4803e863..2663996f 100644
--- a/libsrc/meshing/meshfunc.cpp
+++ b/libsrc/meshing/meshfunc.cpp
@@ -548,10 +548,6 @@ namespace netgen
      auto & ident_points = mesh.GetIdentifications().GetIdentifiedPoints();
      ident_points.DeleteData();
      
-     for(auto & m_ : md) {
-     m_.mesh->GetIdentifications().Delete();
-      m_.mesh->Save("mesh_domain_"+ToString(m_.domain)+".vol.gz");
-  }
 
      for(auto & m_ : md)
      {
@@ -615,9 +611,7 @@ namespace netgen
   {
     static Timer t("MeshVolume"); RegionTimer reg(t);
 
-  cout << "initial compress " << endl;
      mesh3d.Compress();
-  cout << "initial compress done" << endl;
 
 
 
diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp
index 04d50338..bec964d5 100644
--- a/libsrc/meshing/meshtype.cpp
+++ b/libsrc/meshing/meshtype.cpp
@@ -2787,6 +2787,20 @@ namespace netgen
   }
 
 
+  Array<INDEX_3> Identifications :: GetPairs () const
+  {
+    Array<INDEX_3> pairs;
+    for (auto i : IntRange(1, identifiedpoints.GetNBags()+1))
+      for (auto j : IntRange(1, identifiedpoints.GetBagSize(i)+1))
+        {
+          INDEX_2 i2;
+          int nr;
+          identifiedpoints.GetData (i, j, i2, nr);
+          pairs.Append ({i2.I1(), i2.I2(), nr});
+        }
+    return pairs;
+  }
+
   void Identifications :: GetPairs (int identnr, 
                                     NgArray<INDEX_2> & identpairs) const
   {
@@ -2832,32 +2846,20 @@ namespace netgen
         }
   }
 
+  // Map points in the identifications to new point numbers
+  // deletes identifications with invalid (mapped) points
   void Identifications :: MapPoints(FlatArray<PointIndex, PointIndex> op2np)
   {
-    Array<NgArray<INDEX_2>> pairs;
-    pairs.SetSize(maxidentnr+1);
-    for(auto i : Range(maxidentnr+1))
-      GetPairs(i, pairs[i]);
-    
+    auto pairs = GetPairs();
     Delete();
-    for (auto i : Range(maxidentnr+1))
-      for(auto pair : pairs[i]) {
-        pair.I1() = op2np[pair.I1()];
-        pair.I2() = op2np[pair.I2()];
-        auto invalid = PointIndex(PointIndex::INVALID);
-        if(pair.I1() != invalid && pair.I2() != invalid)
-          Add(pair.I1(), pair.I2(), i);
+
+    for(auto pair : pairs)
+      {
+        auto p1 = op2np[pair.I1()];
+        auto p2 = op2np[pair.I2()];
+        if(p1.IsValid() && p2.IsValid())
+          Add(p1, p2, pair.I3());
       }
-    // for (auto i : IntRange(1, identifiedpoints.GetNBags()+1))
-    //   for (auto j : IntRange(1, identifiedpoints.GetBagSize(i)+1))
-    //     {
-    //       INDEX_2 i2;
-    //       int nr;
-    //       identifiedpoints.GetData (i, j, i2, nr);
-    //       i2.I1() = op2np[i2.I1()];
-    //       i2.I2() = op2np[i2.I2()];
-    //       identifiedpoints.SetData (i, j, i2, nr);	    
-    //     }
   }
 
   void Identifications :: Print (ostream & ost) const
diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp
index 8e24aed7..f452c1b3 100644
--- a/libsrc/meshing/meshtype.hpp
+++ b/libsrc/meshing/meshtype.hpp
@@ -1592,6 +1592,7 @@ namespace netgen
     
     ///
     DLL_HEADER void GetPairs (int identnr, NgArray<INDEX_2> & identpairs) const;
+    DLL_HEADER Array<INDEX_3> GetPairs () const;
     ///
     int GetMaxNr () const { return maxidentnr; }