diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp
index 4257e08d..de2cd864 100644
--- a/libsrc/meshing/basegeom.cpp
+++ b/libsrc/meshing/basegeom.cpp
@@ -22,7 +22,7 @@ namespace netgen
           ArrayMem<int, 1> points;
           tree.GetIntersecting(p, p, points);
           if(points.Size()==0)
-              throw Exception("cannot find mapped point");
+              throw Exception("cannot find mapped point " + ToString(p));
           return points[0];
       }
 
@@ -935,11 +935,12 @@ namespace netgen
     }
 
     bool have_identifications = false;
+    std::map<tuple<PointIndex, int>, PointIndex> mapto;
     for(auto & face : faces)
         if(face->primary != face.get())
         {
             have_identifications = true;
-            MapSurfaceMesh(mesh, *face);
+            MapSurfaceMesh(mesh, *face, mapto);
         }
 
     // identify points on faces
@@ -983,7 +984,8 @@ namespace netgen
                 if(ident.from == face.get())
                     for(auto pi : pi_of_face[face->nr])
                     {
-                        auto pi_other = tree.Find(ident.trafo(mesh[pi]));
+                        auto pi_primary = ident.from->primary->nr == ident.from->nr ? pi : mapto[{pi, ident.to->primary->nr}];
+                        auto pi_other = ident.to->primary->nr == ident.to->nr ? pi_primary : mapto[{pi_primary, ident.to->nr}];
                         mesh_ident.Add(pi, pi_other, ident.name, ident.type);
                     }
             }
@@ -993,7 +995,7 @@ namespace netgen
     multithread.task = savetask;
   }
 
-  void NetgenGeometry :: MapSurfaceMesh( Mesh & mesh, const GeometryFace & dst ) const
+  void NetgenGeometry :: MapSurfaceMesh( Mesh & mesh, const GeometryFace & dst, std::map<tuple<PointIndex, int>, PointIndex> & mapto ) const
   {
     static Timer timer("MapSurfaceMesh");
     RegionTimer rt(timer);
@@ -1074,6 +1076,8 @@ namespace netgen
                 pmap[pi] = mesh.AddPoint(trafo(mesh[pi]), 1, SURFACEPOINT);
               }
               sel_new[i] = pmap[pi];
+              mapto[{pi, dst.nr}] = pmap[pi];
+              mapto[{pmap[pi], src.nr}] = pi;
           }
           if(do_invert.IsTrue())
               sel_new.Invert();
diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp
index 7dbc1e17..67badd6e 100644
--- a/libsrc/meshing/basegeom.hpp
+++ b/libsrc/meshing/basegeom.hpp
@@ -243,7 +243,7 @@ namespace netgen
     virtual void MeshSurface(Mesh& mesh, const MeshingParameters& mparam) const;
     virtual bool MeshFace(Mesh& mesh, const MeshingParameters& mparam,
                      int nr, FlatArray<int, PointIndex> glob2loc) const;
-    virtual void MapSurfaceMesh( Mesh & mesh, const GeometryFace & dst ) const;
+    virtual void MapSurfaceMesh( Mesh & mesh, const GeometryFace & dst, std::map<tuple<PointIndex, int>, PointIndex> & mapto) const;
     virtual void OptimizeSurface(Mesh& mesh, const MeshingParameters& mparam) const;
 
     virtual void FinalizeMesh(Mesh& mesh) const;