diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9fc53163..19f48136 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -60,7 +60,7 @@ option(LIBIGL_CYCODEBASE           "Build igl::cycodebase bindings"        ON)
 FetchContent_Declare(
   libigl
   GIT_REPOSITORY https://github.com/libigl/libigl.git
-  GIT_TAG 477e15a3d566a21f415aa5ee62992b12a836b01b
+  GIT_TAG 72217e00bae6a7e041221e187aafbdc793f905d2
 )
 FetchContent_MakeAvailable(libigl)
 
@@ -208,4 +208,3 @@ if(LIBIGL_CYCODEBASE)
 endif()
 
 
-
diff --git a/src/marching_cubes.cpp b/src/marching_cubes.cpp
index 72c14972..6de13b1c 100644
--- a/src/marching_cubes.cpp
+++ b/src/marching_cubes.cpp
@@ -35,8 +35,9 @@ namespace pyigl
   {
     Eigen::MatrixXN V;
     Eigen::MatrixXI F;
-    igl::marching_cubes(S,GV,GI,isovalue,V,F);
-    return std::make_tuple(V,F);
+    std::unordered_map<std::int64_t,int> E2V;
+    igl::marching_cubes(S,GV,GI,isovalue,V,F,E2V);
+    return std::make_tuple(V,F,E2V);
   }
 }
 
@@ -68,7 +69,7 @@ points, and generates a mesh defined by vertices and faces
 @param[out] E2V  map from edge key to index into rows of V
 
 # unpack keys into (i,j,v) index triplets
-EV = np.array([[k & 0xFFFFFFFF, k >> 32, v] for k, v in E2V.items()], dtype=np.int64)
+E2V_triplets = np.array([[k & 0xFFFFFFFF, k >> 32, v] for k, v in E2V.items()], dtype=np.int64)
 )");
 
   m.def(
@@ -86,5 +87,10 @@ points, and generates a mesh defined by vertices and faces
 @param[in] GI  #GI by 8 list of grid corner indices into rows of GV
 @param[in] isovalue  the isovalue of the surface to reconstruct
 @param[out] V  #V by 3 list of mesh vertex positions
-@param[out] F  #F by 3 list of mesh triangle indices into rows of V)");
+@param[out] F  #F by 3 list of mesh triangle indices into rows of V
+@param[out] E2V  map from edge key to index into rows of V
+
+# unpack keys into (i,j,v) index triplets
+E2V_triplets = np.array([[k & 0xFFFFFFFF, k >> 32, v] for k, v in E2V.items()], dtype=np.int64)
+)");
 }
diff --git a/tests/test_all.py b/tests/test_all.py
index 9728d16d..232946ab 100644
--- a/tests/test_all.py
+++ b/tests/test_all.py
@@ -424,7 +424,8 @@ def test_implicit():
     S = np.sqrt(((GV - np.array([0.5,0.5,0.5],dtype=np.float64))**2).sum(axis=1))-0.25;
     V,F,E2V = igl.marching_cubes(S,GV,res[0],res[1],res[2])
     # unpack keys into (i,j,v) index triplets
-    EV = np.array([[k & 0xFFFFFFFF, k >> 32, v] for k, v in E2V.items()], dtype=np.int64)
+    EV_list = [[k & 0xFFFFFFFF, k >> 32, v] for k, v in E2V.items()]
+    EV = np.array(EV_list, dtype=np.int64) if EV_list else np.empty((0,3), dtype=np.int64)
 
     h = igl.avg_edge_length(V,F)
     m0,m1,m2 = igl.moments(V,F)
@@ -625,7 +626,14 @@ def udf_sphere(Q):
     h = h0 / (2**max_depth)
     unique_ijk, J, unique_corners = igl.unique_sparse_voxel_corners(origin,h0,max_depth,ijk)
     unique_S = sdf_sphere(unique_corners)
-    V,F = igl.marching_cubes(unique_S,unique_corners,J,0.0)
+    V,F,E2V = igl.marching_cubes(unique_S,unique_corners,J,0.0)
+    assert V.shape[0] > 0
+    assert F.shape[0] > 0
+    EV_list = [[k & 0xFFFFFFFF, k >> 32, v] for k, v in E2V.items()]
+    EV = np.array(EV_list, dtype=np.int64) if EV_list else np.empty((0,3), dtype=np.int64)
+    assert len(E2V) == EV.shape[0]
+    assert np.all(EV[:,2] >= 0)
+    assert np.all(EV[:,2] < V.shape[0])
 
 def test_is_intrinsic_delaunay() -> None:
     # vs and fs come from a simple plane from pyvista