diff --git a/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h b/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h
index 8ac55c3a3115..67c23660d047 100644
--- a/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h
+++ b/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h
@@ -1112,7 +1112,7 @@ class Surface_mesh
   
   bool join(const Surface_mesh& other)
   {
-    size_type nv = num_vertices(), nh = num_halfedges(), nf = num_faces();
+    const size_type nv = num_vertices(), nh = num_halfedges(), nf = num_faces();
     resize(num_vertices()+  other.num_vertices(),
             num_edges()+  other.num_edges(),
             num_faces()+  other.num_faces());
@@ -1155,7 +1155,9 @@ class Surface_mesh
         Vertex_index vi(nv+other.vertices_freelist_);
         Halfedge_index inf((std::numeric_limits<size_type>::max)());
         while(vconn_[vi].halfedge_ != inf){
-          vi = Vertex_index(size_type(vconn_[vi].halfedge_));
+          Vertex_index corrected_vi = Vertex_index(size_type(vconn_[vi].halfedge_)+nv-nh);
+          vconn_[vi].halfedge_ = Halfedge_index(corrected_vi);
+          vi = corrected_vi;
         }
         vconn_[vi].halfedge_ = Halfedge_index(vertices_freelist_);
       }
@@ -1166,7 +1168,9 @@ class Surface_mesh
         Face_index fi(nf+other.faces_freelist_);
         Halfedge_index inf((std::numeric_limits<size_type>::max)());
         while(fconn_[fi].halfedge_ != inf){
-          fi = Face_index(size_type(fconn_[fi].halfedge_));
+          Face_index corrected_fi = Face_index(size_type(fconn_[fi].halfedge_)+nf-nh);
+          fconn_[fi].halfedge_ = Halfedge_index(corrected_fi);
+          fi = corrected_fi;
         }
         fconn_[fi].halfedge_ = Halfedge_index(faces_freelist_);
       }
@@ -1174,14 +1178,14 @@ class Surface_mesh
     }
     if(other.edges_freelist_ != inf_value){
       if(edges_freelist_ != inf_value){
-        Halfedge_index hi((nh>>1)+other.edges_freelist_);
+        Halfedge_index hi(nh+other.edges_freelist_);
         Halfedge_index inf((std::numeric_limits<size_type>::max)());
         while(hconn_[hi].next_halfedge_ != inf){
           hi = hconn_[hi].next_halfedge_;
         }
         hconn_[hi].next_halfedge_ = Halfedge_index(edges_freelist_);
       }
-      edges_freelist_ = (nh>>1) + other.edges_freelist_; 
+      edges_freelist_ = nh + other.edges_freelist_; 
     }
     garbage_ = garbage_ || other.garbage_;
     removed_vertices_ += other.removed_vertices_;
@@ -1962,6 +1966,22 @@ class Surface_mesh
   }
   /// @}
 
+#if defined(CGAL_SURFACE_MESH_TEST_SUITE)
+  Vertex_index vertex_freelist() const
+  {
+    return Vertex_index(vertices_freelist_);
+  }
+  
+  Face_index face_freelist() const
+  {
+    return Face_index(faces_freelist_);
+  }
+  
+  Edge_index edge_freelist() const
+  {
+    return Edge_index(edges_freelist_>>1);
+  }
+#endif
   
 private: //--------------------------------------------------- helper functions
 
diff --git a/Surface_mesh/test/Surface_mesh/sm_join_test.cpp b/Surface_mesh/test/Surface_mesh/sm_join_test.cpp
new file mode 100644
index 000000000000..937f7b05e5b8
--- /dev/null
+++ b/Surface_mesh/test/Surface_mesh/sm_join_test.cpp
@@ -0,0 +1,77 @@
+#define CGAL_SURFACE_MESH_TEST_SUITE 1  // so that we can access the freelists
+
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/boost/graph/Euler_operations.h>
+
+#include <iostream>
+#include <fstream>
+
+typedef CGAL::Simple_cartesian<double> K;
+typedef K::Point_3 Point_3;
+typedef CGAL::Surface_mesh<Point_3> Sm;
+
+typedef boost::graph_traits<Sm>::vertex_descriptor vertex_descriptor;
+typedef boost::graph_traits<Sm>::halfedge_descriptor halfedge_descriptor;
+typedef boost::graph_traits<Sm>::edge_descriptor edge_descriptor;
+typedef boost::graph_traits<Sm>::face_descriptor face_descriptor;
+
+void
+freelist(const Sm& sm, int vc, int fc, int ec)
+{
+  std::cout << "vertex freelist" << std::endl;
+  vertex_descriptor vd = sm.vertex_freelist();
+  while(vd != sm.null_vertex()){
+    --vc;
+    std::cout << vd << std::endl;
+    halfedge_descriptor hd = halfedge(vd,sm);
+    vd = vertex_descriptor((Sm::size_type)hd);
+  }
+  assert(vc == 0);
+  
+  std::cout << "face freelist" << std::endl;
+  face_descriptor fd = sm.face_freelist();
+  while(fd != sm.null_face()){
+    --fc;
+    std::cout << fd << std::endl;
+    halfedge_descriptor hd = halfedge(fd,sm);
+    fd = face_descriptor((Sm::size_type)hd);
+  }
+  assert(fc == 0);
+
+  std::cout << "edge freelist" << std::endl;
+  edge_descriptor ed = sm.edge_freelist();
+  while(ed != sm.null_edge()){
+    --ec;
+    std::cout << ed << std::endl;
+    halfedge_descriptor hd = next(halfedge(ed,sm),sm);
+    ed = edge(hd,sm);
+  }
+  assert(ec == 0);
+}
+
+
+int main()
+{
+  Sm sm1, sm2;
+  {
+    std::ifstream in("cube.off");
+    in >> sm1;
+    CGAL::Euler::remove_center_vertex(*(halfedges(sm1).first),sm1);
+  }
+  freelist(sm1,1,5,6);
+  
+  {
+    std::ifstream in("cube.off");
+    in >> sm2;
+    CGAL::Euler::remove_center_vertex(*(halfedges(sm2).first),sm2);
+  }
+
+  freelist(sm1,1,5,6);
+
+  sm1.join(sm2);
+  freelist(sm1,2,10,12);
+  
+  return 0;
+}
+