diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp
index dd108592..4dcbedbc 100644
--- a/libsrc/meshing/basegeom.cpp
+++ b/libsrc/meshing/basegeom.cpp
@@ -56,18 +56,18 @@ namespace netgen
           return false;
       auto & e = *other_ptr;
 
-      if(tol < Dist(GetCenter(), e.GetCenter()))
+      if(tol < Dist(trafo(GetCenter()), e.GetCenter()))
           return false;
 
-      auto v0 = GetStartVertex().GetPoint();
-      auto v1 = GetEndVertex().GetPoint();
+      auto v0 = trafo(GetStartVertex().GetPoint());
+      auto v1 = trafo(GetEndVertex().GetPoint());
       auto w0 = e.GetStartVertex().GetPoint();
       auto w1 = e.GetEndVertex().GetPoint();
 
       // have two closed edges, use midpoints to compare
       if(Dist(v0,v1) < tol && Dist(w0,w1) < tol)
       {
-          v1 = GetPoint(0.5);
+          v1 = trafo(GetPoint(0.5));
           w1 = other_ptr->GetPoint(0.5);
       }
 
@@ -75,6 +75,38 @@ namespace netgen
               (Dist(v0, w1) < tol && Dist(v1, w0) < tol) );
   }
 
+  bool GeometryFace :: IsMappedShape( const GeometryShape & other_, const Transformation<3> & trafo, double tol ) const
+  {
+      const auto other_ptr = dynamic_cast<const GeometryFace*>(&other_);
+      if(!other_ptr)
+          return false;
+      auto & f = *other_ptr;
+
+      if(tol < Dist(GetCenter(), f.GetCenter()))
+          return false;
+
+      // simple check: check if there is a bijective mapping of mapped edges
+      auto & other_edges = f.edges;
+      if(edges.Size() != other_edges.Size())
+          return false;
+
+      auto nedges = edges.Size();
+      Array<bool> is_mapped(nedges);
+      is_mapped = false;
+
+      for(auto e : edges)
+      {
+          int found_mapping = 0;
+          for(auto other_e : other_edges)
+              if(e->IsMappedShape(*other_e, trafo, tol))
+                  found_mapping++;
+          if(found_mapping != 1)
+              return false;
+      }
+
+      return true;
+  }
+
   void GeometryFace :: RestrictHTrig(Mesh& mesh,
                                      const PointGeomInfo& gi0,
                                      const PointGeomInfo& gi1,
@@ -169,6 +201,15 @@ namespace netgen
 
   void NetgenGeometry :: ProcessIdentifications()
   {
+      for(auto i : Range(vertices))
+          vertices[i]->nr = i;
+      for(auto i : Range(edges))
+          edges[i]->nr = i;
+      for(auto i : Range(faces))
+          faces[i]->nr = i;
+      for(auto i : Range(solids))
+          solids[i]->nr = i;
+
       auto mirror_identifications = [&] ( auto & shapes )
       {
           for(auto i : Range(shapes))
@@ -181,11 +222,39 @@ namespace netgen
           }
       };
 
+      auto tol = 1e-8 * bounding_box.Diam();
+      for(auto & f : faces)
+        for(auto & ident: f->identifications)
+          for(auto e : static_cast<GeometryFace*>(ident.from)->edges)
+            for(auto e_other : static_cast<GeometryFace*>(ident.to)->edges)
+              if(e->IsMappedShape(*e_other, ident.trafo, tol))
+                e->identifications.Append( {e, e_other, ident.trafo, ident.type, ident.name} );
+
+      for(auto & e : edges)
+        for(auto & ident: e->identifications)
+          {
+              auto & from = static_cast<GeometryEdge&>(*ident.from);
+              auto & to = static_cast<GeometryEdge&>(*ident.to);
+
+              GeometryVertex * pfrom[] = { &from.GetStartVertex(), &from.GetEndVertex() };
+              GeometryVertex * pto[] = { &to.GetStartVertex(), &to.GetEndVertex() };
+
+              // swap points of other edge if necessary
+              Point<3> p_from0 = ident.trafo(from.GetStartVertex().GetPoint());
+              Point<3> p_from1 = ident.trafo(from.GetEndVertex().GetPoint());
+              Point<3> p_to0 = ident.trafo(to.GetStartVertex().GetPoint());
+
+              if(Dist(p_from1, p_to0) < Dist(p_from1, p_to0))
+                  swap(pto[0], pto[1]);
+
+              for(auto i : Range(2))
+                  pfrom[i]->identifications.Append( {pfrom[i], pto[i], ident.trafo, ident.type, ident.name} );
+          }
+
       mirror_identifications(vertices);
       mirror_identifications(edges);
       mirror_identifications(faces);
 
-      // todo: propagate identifications faces -> edges -> vertices
 
       auto find_primary = [&] (auto & shapes)
       {
@@ -729,6 +798,9 @@ namespace netgen
                     sel[1] = s[1];
                     sel[2] = tree.Find(trafo(mesh[s[1]]));
                     sel[3] = tree.Find(trafo(mesh[s[0]]));
+                    for(auto i : Range(4))
+                        sel.GeomInfo()[i] = face.Project(mesh[sel[i]]);
+
                     sel.SetIndex(face.nr+1);
                     mesh.AddSurfaceElement(sel);
                 }
diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp
index e1e3016b..d7371e2b 100644
--- a/libsrc/meshing/basegeom.hpp
+++ b/libsrc/meshing/basegeom.hpp
@@ -70,11 +70,19 @@ namespace netgen
 
   class DLL_HEADER GeometryEdge : public GeometryShape
   {
+  protected:
+      GeometryVertex *start, *end;
   public:
     int domin=-1, domout=-1;
 
-    virtual const GeometryVertex& GetStartVertex() const = 0;
-    virtual const GeometryVertex& GetEndVertex() const = 0;
+    GeometryEdge( GeometryVertex &start_, GeometryVertex &end_ )
+        : start(&start_), end(&end_)
+    {}
+
+    virtual const GeometryVertex& GetStartVertex() const { return *start; }
+    virtual const GeometryVertex& GetEndVertex() const { return *end; }
+    virtual GeometryVertex& GetStartVertex() { return *start; }
+    virtual GeometryVertex& GetEndVertex() { return *end; }
     virtual double GetLength() const = 0;
     virtual Point<3> GetCenter() const = 0;
     virtual Point<3> GetPoint(double t) const = 0;
@@ -100,11 +108,13 @@ namespace netgen
   class DLL_HEADER GeometryFace : public GeometryShape
   {
   public:
+    Array<GeometryEdge*> edges;
     int domin=-1, domout=-1;
 
     virtual Point<3> GetCenter() const = 0;
     virtual size_t GetNBoundaries() const = 0;
     virtual Array<Segment> GetBoundary(const Mesh& mesh) const = 0;
+
     virtual PointGeomInfo Project(Point<3>& p) const = 0;
     // Project point using geo info. Fast if point is close to
     // parametrization in geo info.
@@ -143,6 +153,8 @@ namespace netgen
         newgi = Project(newp);
     }
 
+    virtual bool IsMappedShape( const GeometryShape & other, const Transformation<3> & trafo, double tolerance ) const override;
+
   protected:
     void RestrictHTrig(Mesh& mesh,
                        const PointGeomInfo& gi0,
diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp
index 597a61fa..68dcd347 100644
--- a/libsrc/meshing/meshfunc.cpp
+++ b/libsrc/meshing/meshfunc.cpp
@@ -156,8 +156,6 @@ namespace netgen
                   if(!pi0.IsValid() || !pi1.IsValid())
                       continue;
 
-                  if(pi1<pi0)
-                      Swap(pi0,pi1);
                   m_ident.Add(pi0, pi1, n);
               }
               m_ident.SetType( n, identifications.GetType(n) );
@@ -166,6 +164,59 @@ namespace netgen
       return ret;
   }
 
+  // Add between identified surface elements (only consider closesurface identifications)
+  void FillCloseSurface( MeshingData & md)
+  {
+      static Timer timer("FillCloseSurface"); RegionTimer rtimer(timer);
+
+      auto & mesh = md.mesh;
+      auto & identifications = mesh->GetIdentifications();
+      auto nmax = identifications.GetMaxNr();
+
+      bool have_closesurfaces = false;
+      for(auto i : Range(1,nmax+1))
+          if(identifications.GetType(i) == Identifications::CLOSESURFACES)
+              have_closesurfaces = true;
+      if(!have_closesurfaces)
+          return;
+
+      NgArray<int, PointIndex::BASE> map;
+      for(auto identnr : Range(1,nmax+1))
+      {
+          if(identifications.GetType(identnr) != Identifications::CLOSESURFACES)
+              continue;
+
+          identifications.GetMap(identnr, map);
+
+          for(auto & sel : mesh->SurfaceElements())
+          {
+              bool is_mapped = true;
+              for(auto pi : sel.PNums())
+                  if(!PointIndex(map[pi]).IsValid())
+                      is_mapped = false;
+
+              if(!is_mapped)
+                  continue;
+              
+              // in case we have symmetric mapping (used in csg), only map in one direction
+              if(map[map[sel[0]]] == sel[0] && map[sel[0]] < sel[0])
+                  continue;
+
+              // insert prism
+              auto np = sel.GetNP();
+              Element el(2*np);
+              for(auto i : Range(np))
+              {
+                  el[i] = sel[i];
+                  el[i+np] = map[sel[i]];
+              }
+
+              el.SetIndex(md.domain);
+              mesh->AddVolumeElement(el);
+          }
+      }
+  }
+
   void CloseOpenQuads( MeshingData & md)
   {
     auto & mesh = *md.mesh;
@@ -488,6 +539,9 @@ namespace netgen
            if (md[i].mesh->CheckOverlappingBoundary())
              throw NgException ("Stop meshing since boundary mesh is overlapping");
          
+         // TODO: FillCloseSurface is not working with CSG closesurfaces
+         if(md[i].mesh->GetGeometry()->GetGeomType() == Mesh::GEOM_OCC)
+           FillCloseSurface( md[i] );
          CloseOpenQuads( md[i] );
          MeshDomain(md[i]);
        });
diff --git a/libsrc/occ/occ_edge.cpp b/libsrc/occ/occ_edge.cpp
index 901417bc..7f2c3d66 100644
--- a/libsrc/occ/occ_edge.cpp
+++ b/libsrc/occ/occ_edge.cpp
@@ -7,8 +7,9 @@
 
 namespace netgen
 {
-    OCCEdge::OCCEdge(TopoDS_Shape edge_)
-        : tedge(edge_.TShape()),
+    OCCEdge::OCCEdge(TopoDS_Shape edge_, GeometryVertex & start_, GeometryVertex & end_)
+        : GeometryEdge(start_, end_),
+          tedge(edge_.TShape()),
           edge(TopoDS::Edge(edge_))
     {
         curve = BRep_Tool::Curve(edge, s0, s1);
@@ -18,26 +19,13 @@ namespace netgen
         if(verts.size() != 2)
             throw Exception("OCC edge does not have 2 vertices");
 
-        start = OCCVertex(verts[0]);
-        end = OCCVertex(verts[1]);
-
         // swap start/end if necessary
-        double d00 = Dist(GetPoint(0), start.GetPoint());
-        double d01 = Dist(GetPoint(0), end.GetPoint());
+        double d00 = Dist(GetPoint(0), start->GetPoint());
+        double d01 = Dist(GetPoint(0), end->GetPoint());
         if(d01 < d00)
             swap(start, end);
     }
 
-    const GeometryVertex& OCCEdge::GetStartVertex() const 
-    {
-        return start;
-    }
-
-    const GeometryVertex& OCCEdge::GetEndVertex() const
-    {
-        return end;
-    }
-
     double OCCEdge::GetLength() const
     {
         return props.Mass();
diff --git a/libsrc/occ/occ_edge.hpp b/libsrc/occ/occ_edge.hpp
index c1d95428..f34cee56 100644
--- a/libsrc/occ/occ_edge.hpp
+++ b/libsrc/occ/occ_edge.hpp
@@ -14,6 +14,7 @@ namespace netgen
 {
     class OCCEdge : public GeometryEdge 
     {
+        public:
         T_Shape tedge;
         TopoDS_Edge edge;
         Handle(Geom_Curve) curve;
@@ -21,16 +22,11 @@ namespace netgen
         GProp_GProps props;
 
         public:
-        OCCVertex start;
-        OCCVertex end;
-
-        OCCEdge(TopoDS_Shape edge_);
+        OCCEdge(TopoDS_Shape edge_, GeometryVertex & start_, GeometryVertex & end_);
 
         auto Shape() const { return edge; }
         T_Shape TShape() const { return tedge; }
 
-        const GeometryVertex& GetStartVertex() const override;
-        const GeometryVertex& GetEndVertex() const override;
         double GetLength() const override;
         Point<3> GetCenter() const override;
         Point<3> GetPoint(double t) const override;
diff --git a/libsrc/occ/occ_face.cpp b/libsrc/occ/occ_face.cpp
index d4d1a7c4..7b985f06 100644
--- a/libsrc/occ/occ_face.cpp
+++ b/libsrc/occ/occ_face.cpp
@@ -98,7 +98,6 @@ namespace netgen
 
                 // auto cof = curve_on_face[ORIENTATION][edgenr];
                 auto edge = edge_on_face[ORIENTATION][edgenr];
-                OCCEdge gedge(edge);
                 double s0, s1;
                 auto cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1);
 
diff --git a/libsrc/occ/occ_utils.cpp b/libsrc/occ/occ_utils.cpp
index 8b747137..210e358e 100644
--- a/libsrc/occ/occ_utils.cpp
+++ b/libsrc/occ/occ_utils.cpp
@@ -11,6 +11,22 @@ namespace netgen
         return occ2ng( Handle(BRep_TVertex)::DownCast(shape)->Pnt() );
     }
 
+    Transformation<3> occ2ng (const gp_Trsf & occ_trafo)
+    {
+        Transformation<3> trafo;
+        auto v = occ_trafo.TranslationPart();
+        auto m = occ_trafo.VectorialPart();
+        auto & tv = trafo.GetVector();
+        auto & tm = trafo.GetMatrix();
+        for(auto i : Range(3))
+        {
+            tv[i] = v.Coord(i+1);
+            for(auto k : Range(3))
+                tm(i,k) = m(i+1,k+1);
+        }
+        return trafo;
+    }
+
     Box<3> GetBoundingBox( const TopoDS_Shape & shape )
     {
         Bnd_Box bb;
diff --git a/libsrc/occ/occ_utils.hpp b/libsrc/occ/occ_utils.hpp
index 1e9fa9a1..e51a2b51 100644
--- a/libsrc/occ/occ_utils.hpp
+++ b/libsrc/occ/occ_utils.hpp
@@ -9,6 +9,7 @@
 #include <TopTools_IndexedMapOfShape.hxx>
 #include <TopoDS.hxx>
 #include <TopoDS_Vertex.hxx>
+#include <gp_Trsf.hxx>
 
 #include "meshing.hpp"
 
@@ -47,6 +48,8 @@ namespace netgen
         return occ2ng (BRep_Tool::Pnt (v));
     }
 
+    DLL_HEADER Transformation<3> occ2ng (const gp_Trsf & t);
+
     inline gp_Pnt ng2occ (const Point<3> & p)
     {
         return gp_Pnt(p(0), p(1), p(2));
diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp
index 290c4ece..1b84b5ac 100644
--- a/libsrc/occ/occgeom.cpp
+++ b/libsrc/occ/occgeom.cpp
@@ -1108,20 +1108,24 @@ namespace netgen
       fsingular = esingular = vsingular = false;
 
       // Add shapes
-      for(auto v : GetVertices(shape))
+      for(auto i1 : Range(1, vmap.Extent()+1))
       {
+          auto v = vmap(i1);
           auto tshape = v.TShape();
           if(vertex_map.count(tshape)!=0)
               continue;
-          vertex_map[tshape] = vertices.Size();
           auto occ_vertex = make_unique<OCCVertex>(TopoDS::Vertex(v));
+          occ_vertex->nr = vertices.Size();
+          vertex_map[tshape] = occ_vertex->nr;
+
           if(global_shape_properties.count(tshape)>0)
               occ_vertex->properties = global_shape_properties[tshape];
           vertices.Append(std::move(occ_vertex));
       }
 
-      for(auto e : GetEdges(shape))
+      for(auto i1 : Range(1, emap.Extent()+1))
       {
+          auto e = emap(i1);
           auto tshape = e.TShape();
           auto edge = TopoDS::Edge(e);
           if(edge_map.count(tshape)!=0)
@@ -1129,15 +1133,15 @@ namespace netgen
           if(BRep_Tool::Degenerated(edge))
               continue;
           edge_map[tshape] = edges.Size();
-          auto occ_edge = make_unique<OCCEdge>(edge);
+          auto verts = GetVertices(e);
+          auto occ_edge = make_unique<OCCEdge>(edge, *vertices[vertex_map[verts[0].TShape()]], *vertices[vertex_map[verts[1].TShape()]] );
           occ_edge->properties = global_shape_properties[tshape];
-          occ_edge->start.nr = vertex_map[occ_edge->start.TShape()];
-          occ_edge->end.nr = vertex_map[occ_edge->end.TShape()];
           edges.Append(std::move(occ_edge));
       }
 
-      for(auto f : GetFaces(shape))
+      for(auto i1 : Range(1, fmap.Extent()+1))
       {
+          auto f = fmap(i1);
           auto tshape = f.TShape();
           if(face_map.count(tshape)==0)
           {
@@ -1145,6 +1149,10 @@ namespace netgen
               auto k = faces.Size();
               face_map[tshape] = k;
               auto occ_face = make_unique<OCCFace>(f);
+
+              for(auto e : GetEdges(f))
+                  occ_face->edges.Append( edges[edge_map[e.TShape()]].get() );
+
               if(global_shape_properties.count(tshape)>0)
                   occ_face->properties = global_shape_properties[tshape];
               faces.Append(std::move(occ_face));
@@ -1162,8 +1170,9 @@ namespace netgen
       }
 
 
-      for(auto s : GetSolids(shape))
+      for(auto i1 : Range(1, somap.Extent()+1))
       {
+          auto s = somap(i1);
           auto tshape = s.TShape();
           int k;
           if(solid_map.count(tshape)==0)
@@ -1194,6 +1203,9 @@ namespace netgen
             if(identifications.count(tshape))
                 for(auto & ident : identifications[tshape])
                 {
+                    if(shape_map.count(ident.from)==0 || shape_map.count(ident.to)==0)
+                        continue;
+
                     ShapeIdentification si{
                         shapes[shape_map[ident.from]].get(),
                         shapes[shape_map[ident.to]].get(),
@@ -1208,9 +1220,8 @@ namespace netgen
       add_identifications( edges, edge_map );
       add_identifications( faces, face_map );
 
-      ProcessIdentifications();
-
       bounding_box = ::netgen::GetBoundingBox( shape );
+      ProcessIdentifications();
    }
 
 
@@ -1940,25 +1951,6 @@ namespace netgen
       return occ2ng( props.CentreOfMass() );
   }
 
-  void OCCGeometry :: IdentifyEdges(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type) 
-  {
-    auto cme = GetCenter(me);
-    auto cyou = GetCenter(you);
-    Transformation<3> trafo{cyou-cme};
-    identifications[me.TShape()].push_back( {me.TShape(), you.TShape(), Transformation<3>(cyou - cme), name, type} );
-
-    auto vme = GetVertices(me);
-    auto vyou = GetVertices(you);
-    Point<3> pme0 = trafo(occ2ng(vme[0]));
-    Point<3> pme1 = trafo(occ2ng(vme[1]));
-    Point<3> pyou = occ2ng(vyou[0]);
-
-    bool do_swap = Dist(pme1, pyou) < Dist(pme0, pyou);
-
-    for(auto i : Range(2))
-        identifications[vme[i].TShape()].push_back( {vme[i].TShape(), vyou[do_swap ? 1-i : i].TShape(), trafo, name, type} );
-  }
-
   bool IsMappedShape(const Transformation<3> & trafo, const TopoDS_Shape & me, const TopoDS_Shape & you)
   {
       if(me.ShapeType() != you.ShapeType()) return false;
@@ -1976,6 +1968,11 @@ namespace netgen
       std::map<T_Shape, T_Shape> vmap;
 
       auto verts_me = GetVertices(me);
+      auto verts_you = GetVertices(you);
+
+      if(verts_me.size() != verts_you.size())
+          return false;
+
       for (auto i : Range(verts_me.size()))
       {
           auto s = verts_me[i].TShape();
@@ -1986,8 +1983,7 @@ namespace netgen
           vmap[s] = nullptr;
       }
           
-      bool all_verts_mapped = true;
-      for (auto vert : GetVertices(you))
+      for (auto vert : verts_you)
       {
           auto s = vert.TShape();
           auto p = occ2ng(s);
@@ -1998,30 +1994,62 @@ namespace netgen
                   return true;
           });
           if(!vert_mapped)
-          {
-              all_verts_mapped = false;
-              break;
-          }
+              return false;
       }
-      return all_verts_mapped;
+      return true;
   }
 
-  void OCCGeometry :: IdentifyFaces(const TopoDS_Shape & solid, const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type) 
+  void Identify(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional<gp_Trsf> opt_trafo) 
   {
-    auto cme = GetCenter(me);
-    auto cyou = GetCenter(you);
-    Transformation<3> trafo(cyou-cme);
+    gp_Trsf trafo;
+    if(opt_trafo)
+    {
+        trafo = *opt_trafo;
+    }
+    else
+    {
+        auto v = GetCenter(you) - GetCenter(me);
+        trafo.SetTranslation(gp_Vec(v[0], v[1], v[2]));
+    }
 
-    identifications[me.TShape()].push_back
-        (OCCIdentification { me.TShape(), you.TShape(), trafo, name, type });
+    ListOfShapes list_me, list_you;
+    list_me.push_back(me);
+    list_you.push_back(you);
+    Identify(list_me, list_you, name, type, trafo);
+  }
 
-    auto edges_me = GetEdges(me);
-    auto edges_you = GetEdges(you);
+  void Identify(const ListOfShapes & me, const ListOfShapes & you, string name, Identifications::ID_TYPE type, gp_Trsf occ_trafo) 
+  {
+    Transformation<3> trafo = occ2ng(occ_trafo);
 
-    for (auto e_me : edges_me)
-      for (auto e_you : edges_you)
-        if(IsMappedShape(trafo, e_me, e_you))
-            IdentifyEdges(e_me, e_you, name, type);
+    ListOfShapes id_me;
+    ListOfShapes id_you;
+
+    if(auto faces_me = me.Faces(); faces_me.size()>0)
+    {
+        id_me = faces_me;
+        id_you = you.Faces();
+    }
+    else if(auto edges_me = me.Edges(); edges_me.size()>0)
+    {
+        id_me = edges_me;
+        id_you = you.Edges();
+    }
+    else
+    {
+        id_me = me.Vertices();
+        id_you = you.Vertices();
+    }
+
+    for(auto shape_me : id_me)
+        for(auto shape_you : id_you)
+        {
+            if(!IsMappedShape(trafo, shape_me, shape_you))
+                continue;
+
+            OCCGeometry::identifications[shape_me.TShape()].push_back
+                (OCCIdentification { shape_me.TShape(), shape_you.TShape(), trafo, name, type });
+        }
   }
 
   void OCCParameters :: Print(ostream & ost) const
diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp
index 27eae5b7..0619cb97 100644
--- a/libsrc/occ/occgeom.hpp
+++ b/libsrc/occ/occgeom.hpp
@@ -44,6 +44,7 @@ namespace netgen
 #define OCCGEOMETRYVISUALIZATIONFULLCHANGE 1  // Compute transformation matrices and redraw
 #define OCCGEOMETRYVISUALIZATIONHALFCHANGE 2  // Redraw
 
+  bool IsMappedShape(const Transformation<3> & trafo, const TopoDS_Shape & me, const TopoDS_Shape & you);
 
   class EntityVisualizationCode
   {
@@ -341,12 +342,13 @@ namespace netgen
     bool ErrorInSurfaceMeshing ();
 
     //      void WriteOCC_STL(char * filename);
-    static void IdentifyEdges(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type);
-    static void IdentifyFaces(const TopoDS_Shape & solid,const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type);
 
   private:
     //bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const;
   };
+
+  void Identify(const ListOfShapes & me, const ListOfShapes & you, string name, Identifications::ID_TYPE type, gp_Trsf occ_trafo);
+  void Identify(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional<gp_Trsf> opt_trafo);
    
 
   void PrintContents (OCCGeometry * geom);
diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp
index e52d270b..178dfbc5 100644
--- a/libsrc/occ/python_occ_shapes.cpp
+++ b/libsrc/occ/python_occ_shapes.cpp
@@ -199,14 +199,14 @@ py::object CastShape(const TopoDS_Shape & s)
 template <class TBuilder>
 void PropagateIdentifications (TBuilder & builder, TopoDS_Shape shape)
 {
-  std::map<T_Shape, T_Shape> mod_map;
+  std::map<T_Shape, std::set<T_Shape>> mod_map;
   std::map<T_Shape, bool> tshape_handled;
 
   for (auto typ : { TopAbs_SOLID, TopAbs_FACE,  TopAbs_EDGE, TopAbs_VERTEX })
     for (TopExp_Explorer e(shape, typ); e.More(); e.Next())
     {
         auto tshape = e.Current().TShape();
-        mod_map[tshape] = tshape;
+        mod_map[tshape].insert(tshape);
         tshape_handled[tshape] = false;
     }
 
@@ -215,11 +215,8 @@ void PropagateIdentifications (TBuilder & builder, TopoDS_Shape shape)
       {
         auto tshape = e.Current().TShape();
 
-        auto & modified = builder.Modified(e.Current());
-        if(modified.Size()!=1)
-            continue;
-
-        mod_map[tshape] = modified.First().TShape();
+        for (auto mods : builder.Modified(e.Current()))
+            mod_map[tshape].insert(mods.TShape());
       }
 
   for (auto typ : { TopAbs_SOLID, TopAbs_FACE,  TopAbs_EDGE, TopAbs_VERTEX })
@@ -236,21 +233,33 @@ void PropagateIdentifications (TBuilder & builder, TopoDS_Shape shape)
 
         auto tshape_mapped = mod_map[tshape];
 
-        for(auto & ident : OCCGeometry::identifications[tshape])
+        for(auto ident : OCCGeometry::identifications[tshape])
         {
-            // update existing identification
-            if(tshape == tshape_mapped)
-            {
-                ident.to = mod_map[ident.to];
-                ident.from = mod_map[ident.from];
-            }
-            else
-            {
-                OCCIdentification id_new = ident;
-                id_new.to = mod_map[id_new.to];
-                id_new.from = mod_map[id_new.from];
-                OCCGeometry::identifications[mod_map[tshape]].push_back(id_new);
-            }
+            // nothing happened
+            if(mod_map[ident.to].size()==1 && mod_map[ident.from].size() ==1)
+                continue;
+
+            auto from = ident.from;
+            auto to = ident.to;
+
+            for(auto from_mapped : mod_map[from])
+                for(auto to_mapped : mod_map[to])
+                {
+                    if(from==from_mapped && to==to_mapped)
+                        continue;
+
+                    TopoDS_Shape s_from; s_from.TShape(from_mapped);
+                    TopoDS_Shape s_to; s_to.TShape(to_mapped);
+
+                    if(!IsMappedShape(ident.trafo, s_from, s_to))
+                        continue;
+
+                    OCCIdentification id_new = ident;
+                    id_new.to = to_mapped;
+                    id_new.from = from_mapped;
+                    auto id_owner = from == tshape ? from_mapped : to_mapped;
+                    OCCGeometry::identifications[id_owner].push_back(id_new);
+                }
         }
     }
 }
@@ -1071,27 +1080,11 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
            BRepMesh_IncrementalMesh (shape, deflection, true);
          })
 
-    .def("Identify", [](const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE idtype) {
-        // only edges supported, by now
-        auto type = me.ShapeType();
-        auto tyou = you.ShapeType();
-        if(type != tyou)
-            throw NgException ("Identify: cannot identify different shape types");
 
-        switch(type)
-        {
-          case TopAbs_VERTEX:
-          case TopAbs_EDGE:
-            OCCGeometry::IdentifyEdges(me, you, name, idtype);
-            break;
-          default:
-            throw NgException ("Identify: unsupported shape type");
-            break;
-        }
-      }, py::arg("other"), py::arg("name"), py::arg("type")=Identifications::PERIODIC, "Identify shapes for periodic meshing")
-
-    .def("Identify", OCCGeometry::IdentifyFaces, "Identify faces",
-            py::arg("from"), py::arg("to"), py::arg("name"), py::arg("type")=Identifications::PERIODIC)
+    .def("Identify", py::overload_cast<const TopoDS_Shape &, const TopoDS_Shape &, string, Identifications::ID_TYPE, std::optional<gp_Trsf>>(&Identify),
+            py::arg("other"), py::arg("name"),
+            py::arg("type")=Identifications::PERIODIC, py::arg("trafo")=nullopt,
+            "Identify shapes for periodic meshing")
 
     .def("Distance", [](const TopoDS_Shape& self,
                         const TopoDS_Shape& other)
@@ -1634,6 +1627,11 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
           }
       }, "set hpref for all elements of list")
     
+    .def("Identify", py::overload_cast<const ListOfShapes&, const ListOfShapes&, string, Identifications::ID_TYPE, gp_Trsf>(&Identify),
+            py::arg("other"), py::arg("name"),
+            py::arg("type")=Identifications::PERIODIC, py::arg("trafo"),
+            "Identify shapes for periodic meshing")
+
     ;
          
 
diff --git a/tests/pytest/test_occ_identifications.py b/tests/pytest/test_occ_identifications.py
new file mode 100644
index 00000000..d15abc5e
--- /dev/null
+++ b/tests/pytest/test_occ_identifications.py
@@ -0,0 +1,76 @@
+import pytest
+
+from netgen.meshing import IdentificationType
+idtype = IdentificationType.CLOSESURFACES
+
+def test_two_boxes():
+    occ = pytest.importorskip("netgen.occ")
+    inner = occ.Box((0,0,0), (1,1,1))
+    trafo = occ.gp_Trsf().Scale(inner.center, 1.1)
+    outer = trafo(inner)
+
+    inner.Identify(outer, "", idtype, trafo)
+    shape = occ.Glue([outer-inner, inner])
+
+    geo = occ.OCCGeometry(shape)
+    mesh = geo.GenerateMesh(maxh=0.3)
+    have_prisms = False
+
+    for el in mesh.Elements3D():
+        if len(el.vertices)==6:
+            have_prisms = True
+            break
+
+    assert have_prisms
+
+def test_two_circles():
+    occ = pytest.importorskip("netgen.occ")
+    circ1 = occ.WorkPlane().Circle(1).Face()
+    trafo = occ.gp_Trsf().Scale(circ1.center, 1.1)
+
+    circ2 = trafo(circ1)
+    circ1.edges[0].Identify(circ2.edges[0], "", idtype, trafo)
+    circ2 -= circ1
+    shape = occ.Glue([circ1, circ2])
+
+    geo = occ.OCCGeometry(shape, 2)
+    mesh = geo.GenerateMesh(maxh=0.2)
+    have_quads = False
+
+    for el in mesh.Elements2D():
+        if len(el.vertices)==4:
+            have_quads = True
+            break
+
+    assert have_quads
+
+def test_cut_identified_face():
+    occ = pytest.importorskip("netgen.occ")
+    from netgen.occ import Z, Box, Cylinder, Glue, OCCGeometry
+    box = Box((-1,-1,0), (1,1,1))
+    cyl = Cylinder( (0,0,0), Z, 0.5, 1 )
+
+    box.faces.Min(Z).Identify(box.faces.Max(Z), "", idtype)
+    shape = Glue([cyl, box])
+    geo = OCCGeometry(shape)
+    mesh = geo.GenerateMesh(maxh=0.5)
+
+    for el in mesh.Elements3D():
+        assert len(el.vertices)==6
+
+def test_identify_multiple_faces():
+    occ = pytest.importorskip("netgen.occ")
+    from netgen.occ import Z, Box, Cylinder, Glue, OCCGeometry, gp_Trsf
+    box = Box((-1,-1,0), (1,1,1))
+    cyl = Cylinder( (0,0,0), Z, 0.5, 1 )
+
+    shape = Glue([box, cyl])
+    bot_faces = shape.faces[Z < 0.1]
+    top_faces = shape.faces[Z > 0.1]
+    bot_faces.Identify(top_faces, "", idtype, gp_Trsf.Translation((0,0,1)))
+
+    geo = OCCGeometry(shape)
+    mesh = geo.GenerateMesh(maxh=0.3)
+
+    for el in mesh.Elements3D():
+        assert len(el.vertices)==6