diff --git a/include/geometrycentral/surface/mutation_manager.h b/include/geometrycentral/surface/mutation_manager.h
index 569b5b7c..b912b79b 100644
--- a/include/geometrycentral/surface/mutation_manager.h
+++ b/include/geometrycentral/surface/mutation_manager.h
@@ -128,20 +128,23 @@ class MutationManager {
   // to `tSplit`.
   // In general, both the `tSplit` and the `newVertexPosition` are used; `tSplit` is necessary to allow callbacks to
   // interpolate data; if called with either the other will be inferred.
-  void splitEdge(Edge e, double tSplit);
-  void splitEdge(Edge e, Vector3 newVertexPosition);
-  void splitEdge(Edge e, double tSplit, Vector3 newVertexPosition);
+  // Returns the new halfedge which points away from the new vertex (so he.vertex() is the new vertex), and is the same
+  // direction as e.halfedge() on the original edge. The halfedge direction of the other part of the new split edge is
+  // also preserved.
+  Halfedge splitEdge(Edge e, double tSplit);
+  Halfedge splitEdge(Edge e, Vector3 newVertexPosition);
+  Halfedge splitEdge(Edge e, double tSplit, Vector3 newVertexPosition);
 
   // Collapse an edge.
-  // Returns true if the edge could actually be collapsed.
-  bool collapseEdge(Edge e, double tCollapse);
-  bool collapseEdge(Edge e, Vector3 newVertexPosition);
-  bool collapseEdge(Edge e, double tCollapse, Vector3 newVertexPosition);
-
-  // Split a face (i.e. insert a vertex into the face)
-  void splitFace(Face f, const std::vector<double>& bSplit);
-  void splitFace(Face f, Vector3 newVertexPosition);
-  void splitFace(Face f, const std::vector<double>& bSplit, Vector3 newVertexPosition);
+  // Returns the new vertex if the edge could be collapsed, and Vertex() otherwise
+  Vertex collapseEdge(Edge e, double tCollapse);
+  Vertex collapseEdge(Edge e, Vector3 newVertexPosition);
+  Vertex collapseEdge(Edge e, double tCollapse, Vector3 newVertexPosition);
+
+  // Split a face (i.e. insert a vertex into the face), and return the new vertex
+  Vertex splitFace(Face f, const std::vector<double>& bSplit);
+  Vertex splitFace(Face f, Vector3 newVertexPosition);
+  Vertex splitFace(Face f, const std::vector<double>& bSplit, Vector3 newVertexPosition);
 
   // ======================================================
   // ======== High-level mutations
diff --git a/include/geometrycentral/surface/quadric_error_simplification.h b/include/geometrycentral/surface/quadric_error_simplification.h
new file mode 100644
index 00000000..4fda48c4
--- /dev/null
+++ b/include/geometrycentral/surface/quadric_error_simplification.h
@@ -0,0 +1,37 @@
+#pragma once
+
+#include "geometrycentral/surface/manifold_surface_mesh.h"
+#include "geometrycentral/surface/mutation_manager.h"
+#include "geometrycentral/surface/vertex_position_geometry.h"
+
+#include "geometrycentral/utilities/elementary_geometry.h"
+
+#include <queue>
+
+namespace geometrycentral {
+namespace surface {
+class Quadric {
+public:
+  Quadric();
+  Quadric(const Eigen::Matrix3d& A_, const Eigen::Vector3d& b_, double c_);
+  Quadric(const Quadric& Q1, const Quadric& Q2);
+  double cost(const Eigen::Vector3d& v);
+  Eigen::Vector3d optimalPoint();
+
+  Quadric operator+=(const Quadric& Q);
+
+protected:
+  Eigen::Matrix3d A;
+  Eigen::Vector3d b;
+  double c;
+
+  friend Quadric operator+(const Quadric& Q1, const Quadric& Q2);
+};
+
+Quadric operator+(const Quadric& Q1, const Quadric& Q2);
+
+void quadricErrorSimplify(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geo, double tol = 0.05);
+void quadricErrorSimplify(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geo, double tol, MutationManager& mm);
+
+} // namespace surface
+} // namespace geometrycentral
diff --git a/include/geometrycentral/surface/subdivide.h b/include/geometrycentral/surface/subdivide.h
new file mode 100644
index 00000000..0febb06d
--- /dev/null
+++ b/include/geometrycentral/surface/subdivide.h
@@ -0,0 +1,19 @@
+#pragma once
+
+#include "geometrycentral/surface/manifold_surface_mesh.h"
+#include "geometrycentral/surface/mutation_manager.h"
+#include "geometrycentral/surface/vertex_position_geometry.h"
+
+namespace geometrycentral {
+namespace surface {
+
+// Subdivide to form quad mesh
+// Note: these do not take MutationManager arguments since MutationManager only supports triangular meshes at the moment
+void linearSubdivide(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geo);
+void catmullClarkSubdivide(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geo);
+
+void loopSubdivide(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geo);
+void loopSubdivide(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geo, MutationManager& mm);
+
+} // namespace surface
+} // namespace geometrycentral
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 990d41d3..18efb7a2 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -46,6 +46,8 @@ SET(SRCS
   surface/tufted_laplacian.cpp
   surface/flip_geodesics.cpp
   surface/transfer_functions.cpp
+  surface/quadric_error_simplification.cpp
+  surface/subdivide.cpp
   #surface/detect_symmetry.cpp
   #surface/mesh_ray_tracer.cpp
   
@@ -105,6 +107,7 @@ SET(HEADERS
   ${INCLUDE_ROOT}/surface/mesh_graph_algorithms.h
   ${INCLUDE_ROOT}/surface/mesh_ray_tracer.h
   ${INCLUDE_ROOT}/surface/parameterize.h
+  ${INCLUDE_ROOT}/surface/quadric_error_simplification.h
   ${INCLUDE_ROOT}/surface/rich_surface_mesh_data.h
   ${INCLUDE_ROOT}/surface/rich_surface_mesh_data.ipp
   ${INCLUDE_ROOT}/surface/polygon_soup_mesh.h
@@ -112,6 +115,7 @@ SET(HEADERS
   ${INCLUDE_ROOT}/surface/signpost_intrinsic_triangulation.ipp
   ${INCLUDE_ROOT}/surface/simple_idt.h
   ${INCLUDE_ROOT}/surface/simple_polygon_mesh.h
+  ${INCLUDE_ROOT}/surface/subdivide.h
   ${INCLUDE_ROOT}/surface/surface_centers.h
   ${INCLUDE_ROOT}/surface/surface_mesh.h
   ${INCLUDE_ROOT}/surface/surface_mesh.ipp
diff --git a/src/surface/manifold_surface_mesh.cpp b/src/surface/manifold_surface_mesh.cpp
index 35a741e4..02f65ee6 100644
--- a/src/surface/manifold_surface_mesh.cpp
+++ b/src/surface/manifold_surface_mesh.cpp
@@ -1263,6 +1263,64 @@ Vertex ManifoldSurfaceMesh::collapseEdgeTriangular(Edge e) {
   modificationTick++;
 }
 
+Face ManifoldSurfaceMesh::removeEdge(Edge e) {
+  if (e.isBoundary()) {
+    throw std::runtime_error("not implemented");
+  }
+
+  // Halfedges/edges/faces that will be removed
+  // (except first face)
+  std::vector<Halfedge> toRemove{e.halfedge(), e.halfedge().twin()};
+  std::vector<Halfedge> ringHalfedges;
+  for (Halfedge heStart : toRemove) {
+    Halfedge he = heStart.next();
+    while (he != heStart) {
+      // The one-ring must not contain any other copies of e, or we cannot remove the edge
+      if (he.edge() == e) {
+        return Face();
+      }
+      ringHalfedges.push_back(he);
+      he = he.next();
+    }
+  }
+
+  // If both faces are the same, we cannot remove the edge. This should have been caught above
+  if (toRemove[0].face() == toRemove[1].face()) {
+    return Face();
+  }
+  Face keepFace = toRemove[0].face();
+
+
+  // Record these before we break pointers
+  Vertex src = e.halfedge().vertex();
+  Vertex dst = e.halfedge().twin().vertex();
+  Halfedge altSrcHedge = e.halfedge().twin().next();
+  Halfedge altDstHedge = e.halfedge().next();
+
+  // Hook up next and face refs for the halfedges along the ring
+  size_t N = ringHalfedges.size();
+  for (size_t i = 0; i < N; i++) {
+    heNextArr[ringHalfedges[i].getIndex()] = ringHalfedges[(i + 1) % N].getIndex();
+    heFaceArr[ringHalfedges[i].getIndex()] = keepFace.getIndex();
+  }
+
+  // only update vHalfedgeArr if needed to avoid disturbing boundary halfedges
+  if (src.halfedge().edge() == e) {
+    vHalfedgeArr[src.getIndex()] = altSrcHedge.getIndex();
+  }
+  if (dst.halfedge().edge() == e) {
+    vHalfedgeArr[dst.getIndex()] = altDstHedge.getIndex();
+  }
+
+  fHalfedgeArr[keepFace.getIndex()] = ringHalfedges[0].getIndex();
+
+  // Actually delete all of the elements
+  deleteElement(toRemove[1].face());
+  deleteEdgeBundle(e);
+
+  modificationTick++;
+  return keepFace;
+}
 
 
 bool ManifoldSurfaceMesh::removeFaceAlongBoundary(Face f) {
diff --git a/src/surface/mutation_manager.cpp b/src/surface/mutation_manager.cpp
index eed1b273..11f4ccd8 100644
--- a/src/surface/mutation_manager.cpp
+++ b/src/surface/mutation_manager.cpp
@@ -59,16 +59,16 @@ bool MutationManager::flipEdge(Edge e) {
   return true;
 }
 
-void MutationManager::splitEdge(Edge e, double tSplit) {
+Halfedge MutationManager::splitEdge(Edge e, double tSplit) {
   Vector3 newPos{0., 0., 0.};
   if (geometry) {
     VertexData<Vector3>& pos = geometry->vertexPositions;
     newPos = (1. - tSplit) * pos[e.halfedge().tailVertex()] + tSplit * pos[e.halfedge().tipVertex()];
   }
-  splitEdge(e, tSplit, newPos);
+  return splitEdge(e, tSplit, newPos);
 }
 
-void MutationManager::splitEdge(Edge e, Vector3 newVertexPosition) {
+Halfedge MutationManager::splitEdge(Edge e, Vector3 newVertexPosition) {
 
   double tSplit = -1;
   GC_SAFETY_ASSERT(geometry, "must have geometry to split by position");
@@ -80,10 +80,10 @@ void MutationManager::splitEdge(Edge e, Vector3 newVertexPosition) {
     tSplit = pointLineSegmentNeaestLocation(newVertexPosition, posTail, posTip);
   }
 
-  splitEdge(e, tSplit, newVertexPosition);
+  return splitEdge(e, tSplit, newVertexPosition);
 }
 
-void MutationManager::splitEdge(Edge e, double tSplit, Vector3 newVertexPosition) {
+Halfedge MutationManager::splitEdge(Edge e, double tSplit, Vector3 newVertexPosition) {
 
   // Invoke before callbacks
   for (EdgeSplitPolicy* policy : edgeSplitPolicies) {
@@ -102,11 +102,13 @@ void MutationManager::splitEdge(Edge e, double tSplit, Vector3 newVertexPosition
   for (EdgeSplitPolicy* policy : edgeSplitPolicies) {
     policy->afterEdgeSplit(newHeFront, newHeBack, tSplit);
   }
+
+  return newHeFront;
 }
 
 // Collapse an edge.
-// Returns true if the edge could actually be collapsed.
-bool MutationManager::collapseEdge(Edge e, double tCollapse) {
+// Returns the new vertex if the edge could be collapsed, and Vertex() otherwise
+Vertex MutationManager::collapseEdge(Edge e, double tCollapse) {
   Vector3 newPos{0., 0., 0.};
   if (geometry) {
     // Find the nearest tCoord
@@ -116,7 +118,7 @@ bool MutationManager::collapseEdge(Edge e, double tCollapse) {
   return collapseEdge(e, tCollapse, newPos);
 }
 
-bool MutationManager::collapseEdge(Edge e, Vector3 newVertexPosition) {
+Vertex MutationManager::collapseEdge(Edge e, Vector3 newVertexPosition) {
 
   double tCollapse = -1;
   GC_SAFETY_ASSERT(geometry, "must have geometry to split by position");
@@ -130,8 +132,7 @@ bool MutationManager::collapseEdge(Edge e, Vector3 newVertexPosition) {
 
   return collapseEdge(e, tCollapse, newVertexPosition);
 }
-
-bool MutationManager::collapseEdge(Edge e, double tCollapse, Vector3 newVertexPosition) {
+Vertex MutationManager::collapseEdge(Edge e, double tCollapse, Vector3 newVertexPosition) {
 
   // Invoke before callbacks
   // TODO need to handle possiblity that collapse fails -- check before calling
@@ -140,7 +141,7 @@ bool MutationManager::collapseEdge(Edge e, double tCollapse, Vector3 newVertexPo
   }
 
   Vertex newV = mesh.collapseEdgeTriangular(e);
-  if (newV == Vertex()) return false;
+  if (newV == Vertex()) return Vertex();
 
   if (geometry) {
     VertexData<Vector3>& pos = geometry->vertexPositions;
@@ -152,11 +153,11 @@ bool MutationManager::collapseEdge(Edge e, double tCollapse, Vector3 newVertexPo
     policy->afterEdgeCollapse(newV, tCollapse);
   }
 
-  return true;
+  return newV;
 }
 
 // Split a face (i.e. insert a vertex into the face)
-void MutationManager::splitFace(Face f, const std::vector<double>& bSplit) {
+Vertex MutationManager::splitFace(Face f, const std::vector<double>& bSplit) {
   Vector3 newPos = Vector3::zero();
   if (geometry) {
     size_t iV = 0;
@@ -167,16 +168,17 @@ void MutationManager::splitFace(Face f, const std::vector<double>& bSplit) {
     }
   }
 
-  splitFace(f, bSplit, newPos);
+  return splitFace(f, bSplit, newPos);
 }
 
-void MutationManager::splitFace(Face f, Vector3 newVertexPosition) {
+Vertex MutationManager::splitFace(Face f, Vector3 newVertexPosition) {
   // TODO
   throw std::runtime_error("Face split based on vertex position not implemented yet");
+  return Vertex();
 }
 
 
-void MutationManager::splitFace(Face f, const std::vector<double>& bSplit, Vector3 newVertexPosition) {
+Vertex MutationManager::splitFace(Face f, const std::vector<double>& bSplit, Vector3 newVertexPosition) {
   // Invoke before callbacks
   for (FaceSplitPolicy* policy : faceSplitPolicies) {
     policy->beforeFaceSplit(f, bSplit);
@@ -192,6 +194,8 @@ void MutationManager::splitFace(Face f, const std::vector<double>& bSplit, Vecto
   for (FaceSplitPolicy* policy : faceSplitPolicies) {
     policy->afterFaceSplit(newV, bSplit);
   }
+
+  return newV;
 }
 
 
@@ -290,7 +294,7 @@ MutationPolicyHandle MutationManager::registerPolicy(MutationPolicy* policyObjec
 }
 
 void MutationManager::removePolicy(const MutationPolicyHandle& toRemove) {
-  // Remove from all lists 
+  // Remove from all lists
   removeFromVector(vertexRepositionPolicies, toRemove.policy);
   removeFromVector(edgeFlipPolicies, toRemove.policy);
   removeFromVector(edgeSplitPolicies, toRemove.policy);
diff --git a/src/surface/quadric_error_simplification.cpp b/src/surface/quadric_error_simplification.cpp
new file mode 100644
index 00000000..7ef9d3a4
--- /dev/null
+++ b/src/surface/quadric_error_simplification.cpp
@@ -0,0 +1,114 @@
+#include "geometrycentral/surface/quadric_error_simplification.h"
+
+namespace geometrycentral {
+namespace surface {
+
+Quadric::Quadric() {
+  A = Eigen::Matrix3d::Zero();
+  b = Eigen::Vector3d::Zero();
+  c = 0;
+}
+
+Quadric::Quadric(const Quadric& Q1, const Quadric& Q2) {
+  A = Q1.A + Q2.A;
+  b = Q1.b + Q2.b;
+  c = Q1.c + Q2.c;
+}
+
+Quadric::Quadric(const Eigen::Matrix3d& A_, const Eigen::Vector3d& b_, double c_) : A(A_), b(b_), c(c_) {}
+
+double Quadric::cost(const Eigen::Vector3d& v) { return v.dot(A * v) + 2 * b.dot(v) + c; }
+
+Eigen::Vector3d Quadric::optimalPoint() { return -A.inverse() * b; }
+
+Quadric Quadric::operator+=(const Quadric& Q) {
+  A += Q.A;
+  b += Q.b;
+  c += Q.c;
+
+  return *this;
+}
+
+Quadric operator+(const Quadric& Q1, const Quadric& Q2) { return Quadric(Q1.A + Q2.A, Q1.b + Q2.b, Q1.c + Q2.c); }
+
+void quadricErrorSimplify(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geo, double tol) {
+  MutationManager mm(mesh, geo);
+  quadricErrorSimplify(mesh, geo, tol, mm);
+}
+
+void quadricErrorSimplify(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geo, double tol, MutationManager& mm) {
+  auto toEigen = [](Vector3 v) -> Eigen::Vector3d {
+    Eigen::Vector3d ret;
+    ret << v.x, v.y, v.z;
+    return ret;
+  };
+  auto fromEigen = [](Eigen::Vector3d v) -> Vector3 { return Vector3{v(0), v(1), v(2)}; };
+
+  VertexData<Quadric> Q(mesh, Quadric());
+
+  geo.requireFaceNormals();
+
+  for (Face f : mesh.faces()) {
+    Eigen::Vector3d n = toEigen(geo.faceNormals[f]);
+    Eigen::Matrix3d M = n * n.transpose();
+    for (Vertex v : f.adjacentVertices()) {
+      Eigen::Vector3d q = toEigen(geo.inputVertexPositions[v]);
+      double d = -n.dot(q);
+
+      Q[v] += Quadric(M, d * n, d * d);
+    }
+  }
+
+  using PotentialEdge = std::tuple<double, Edge>;
+
+  auto cmp = [](const PotentialEdge& a, const PotentialEdge& b) -> bool { return std::get<0>(a) > std::get<0>(b); };
+
+  std::priority_queue<PotentialEdge, std::vector<PotentialEdge>, decltype(cmp)> edgesToCheck(cmp);
+
+  for (Edge e : mesh.edges()) {
+    Quadric Qe = Q[e.halfedge().tailVertex()] + Q[e.halfedge().tipVertex()];
+    Eigen::Vector3d q = Qe.optimalPoint();
+    double cost = Qe.cost(q);
+    edgesToCheck.push(std::make_tuple(cost, e));
+  }
+
+  while (!edgesToCheck.empty()) {
+    PotentialEdge best = edgesToCheck.top();
+    edgesToCheck.pop();
+
+    // Stop when collapse becomes too expensive
+    double cost = std::get<0>(best);
+    if (cost > tol) break;
+
+    Edge e = std::get<1>(best);
+    if (!e.isDead()) {
+
+      Vertex v1 = e.halfedge().tailVertex();
+      Vertex v2 = e.halfedge().tipVertex();
+
+      // Get edge quadric
+      Quadric Qe(Q[v1], Q[v2]);
+      Eigen::Vector3d q = Qe.optimalPoint();
+
+      // If either vertex has been collapsed since the edge was pushed
+      // onto the queue, the old cost was wrong. In that case, give up
+      if (abs(cost - Qe.cost(q)) > 1e-8) continue;
+
+      Vertex v = mm.collapseEdge(e, fromEigen(q));
+      if (v == Vertex()) continue;
+      Q[v] = Qe;
+
+      for (Edge f : v.adjacentEdges()) {
+        Quadric Qf(Q[f.halfedge().tailVertex()], Q[f.halfedge().tipVertex()]);
+        Eigen::Vector3d q = Qf.optimalPoint();
+        double cost = Qf.cost(q);
+        edgesToCheck.push(std::make_tuple(cost, f));
+      }
+    }
+  }
+
+  mesh.compress();
+  return;
+}
+} // namespace surface
+} // namespace geometrycentral
diff --git a/src/surface/subdivide.cpp b/src/surface/subdivide.cpp
new file mode 100644
index 00000000..adfe0f0c
--- /dev/null
+++ b/src/surface/subdivide.cpp
@@ -0,0 +1,255 @@
+#include "geometrycentral/surface/subdivide.h"
+
+namespace geometrycentral {
+namespace surface {
+
+void linearSubdivide(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geo) {
+
+  VertexData<Vector3>& pos = geo.inputVertexPositions;
+
+  // Compute new positions for original vertices
+  VertexData<Vector3> newPositions = geo.inputVertexPositions;
+
+  FaceData<Vector3> splitFacePositions(mesh);
+  for (Face f : mesh.faces()) {
+    double D = (double)f.degree();
+    splitFacePositions[f] = Vector3::zero();
+    for (Vertex v : f.adjacentVertices()) {
+      splitFacePositions[f] += geo.inputVertexPositions[v] / D;
+    }
+  }
+
+  EdgeData<Vector3> splitEdgePositions(mesh);
+  for (Edge e : mesh.edges()) {
+    splitEdgePositions[e] = (pos[e.halfedge().tailVertex()] + pos[e.halfedge().tipVertex()]) / 2.;
+  }
+
+  // Subdivide edges
+  VertexData<bool> isOrigVert(mesh, true);
+  EdgeData<bool> isOrigEdge(mesh, true);
+  for (Edge e : mesh.edges()) {
+    if (!isOrigEdge[e]) continue;
+
+    Vector3 newPos = splitEdgePositions[e];
+
+    // split the edge
+    Halfedge newHe = mesh.insertVertexAlongEdge(e);
+    Vertex newV = newHe.vertex();
+
+    isOrigVert[newV] = false;
+    isOrigEdge[newHe.edge()] = false;
+    isOrigEdge[newHe.twin().next().edge()] = false;
+
+    GC_SAFETY_ASSERT(isOrigVert[newHe.twin().vertex()] && isOrigVert[newHe.twin().next().twin().vertex()], "???");
+
+    newPositions[newV] = newPos;
+  }
+
+  // Subdivide faces
+  FaceData<bool> isOrigFace(mesh, true);
+  for (Face f : mesh.faces()) {
+    if (!isOrigFace[f]) continue;
+
+    Vector3 newPos = splitFacePositions[f];
+
+    // split the face
+    Vertex newV = mesh.insertVertex(f);
+    isOrigVert[newV] = false;
+    newPositions[newV] = newPos;
+
+    for (Face f : newV.adjacentFaces()) {
+      isOrigFace[f] = false;
+    }
+
+    // get incoming halfedges before we start deleting edges
+    std::vector<Halfedge> incomingHalfedges;
+    for (Halfedge he : newV.incomingHalfedges()) {
+      incomingHalfedges.push_back(he);
+    }
+
+    for (Halfedge he : incomingHalfedges) {
+      if (isOrigVert[he.vertex()]) {
+        Face mergedF = mesh.removeEdge(he.edge());
+        if (mergedF == Face()) {
+          std::cout << "merge was impossible?" << std::endl;
+        } else {
+          isOrigFace[mergedF] = false;
+        }
+      }
+    }
+  }
+
+  mesh.compress();
+  geo.inputVertexPositions = newPositions;
+  geo.refreshQuantities();
+}
+
+void catmullClarkSubdivide(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geo) {
+  VertexData<Vector3>& pos = geo.inputVertexPositions;
+
+  // Compute new positions for original vertices
+  VertexData<Vector3> newPositions(mesh);
+
+  FaceData<Vector3> splitFacePositions(mesh);
+  for (Face f : mesh.faces()) {
+    double D = (double)f.degree();
+    splitFacePositions[f] = Vector3::zero();
+    for (Vertex v : f.adjacentVertices()) {
+      splitFacePositions[f] += geo.inputVertexPositions[v] / D;
+    }
+  }
+
+  EdgeData<Vector3> splitEdgePositions(mesh);
+  for (Edge e : mesh.edges()) {
+    std::array<Face, 2> neigh{e.halfedge().face(), e.halfedge().twin().face()};
+    splitEdgePositions[e] = (splitFacePositions[neigh[0]] + splitFacePositions[neigh[1]]) / 2.;
+  }
+
+  for (Vertex v : mesh.vertices()) {
+    double D = (double)v.degree();
+    Vector3 S = geo.inputVertexPositions[v];
+    Vector3 R = Vector3::zero();
+    for (Edge e : v.adjacentEdges()) {
+      R += splitEdgePositions[e] / D;
+    }
+    Vector3 Q = Vector3::zero();
+    for (Face f : v.adjacentFaces()) {
+      Q += splitFacePositions[f] / D;
+    }
+    newPositions[v] = (Q + 2 * R + (D - 3) * S) / D;
+  }
+
+  // Subdivide edges
+  VertexData<bool> isOrigVert(mesh, true);
+  EdgeData<bool> isOrigEdge(mesh, true);
+  for (Edge e : mesh.edges()) {
+    if (!isOrigEdge[e]) continue;
+
+    Vector3 newPos = splitEdgePositions[e];
+
+    // split the edge
+    Halfedge newHe = mesh.insertVertexAlongEdge(e);
+    Vertex newV = newHe.vertex();
+
+    isOrigVert[newV] = false;
+    isOrigEdge[newHe.edge()] = false;
+    isOrigEdge[newHe.twin().next().edge()] = false;
+
+    GC_SAFETY_ASSERT(isOrigVert[newHe.twin().vertex()] && isOrigVert[newHe.twin().next().twin().vertex()], "???");
+
+    newPositions[newV] = newPos;
+  }
+
+  // Subdivide faces
+  FaceData<bool> isOrigFace(mesh, true);
+  for (Face f : mesh.faces()) {
+    if (!isOrigFace[f]) continue;
+
+    Vector3 newPos = splitFacePositions[f];
+
+    // split the face
+    Vertex newV = mesh.insertVertex(f);
+    isOrigVert[newV] = false;
+    newPositions[newV] = newPos;
+
+    for (Face f : newV.adjacentFaces()) {
+      isOrigFace[f] = false;
+    }
+
+    // get incoming halfedges before we start deleting edges
+    std::vector<Halfedge> incomingHalfedges;
+    for (Halfedge he : newV.incomingHalfedges()) {
+      incomingHalfedges.push_back(he);
+    }
+
+    for (Halfedge he : incomingHalfedges) {
+      if (isOrigVert[he.vertex()]) {
+        Face mergedF = mesh.removeEdge(he.edge());
+        if (mergedF == Face()) {
+          std::cout << "merge was impossible?" << std::endl;
+        } else {
+          isOrigFace[mergedF] = false;
+        }
+      }
+    }
+  }
+
+  mesh.compress();
+  geo.inputVertexPositions = newPositions;
+  geo.refreshQuantities();
+}
+
+void loopSubdivide(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geo) {
+  MutationManager mm(mesh, geo);
+  return loopSubdivide(mesh, geo, mm);
+}
+
+void loopSubdivide(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geo, MutationManager& mm) {
+  GC_SAFETY_ASSERT(mesh.isTriangular(), "Cannot run loop subdivision on a mesh with non-triangular faces");
+
+  VertexData<Vector3>& pos = geo.inputVertexPositions;
+  VertexData<bool> isOrigVert(mesh, true);
+  EdgeData<bool> isOrigEdge(mesh, true);
+
+  // Compute new positions for original vertices
+  VertexData<Vector3> newPositions(mesh);
+  for (Vertex v : mesh.vertices()) {
+    double n = (double)v.degree();
+    double u = (n == 3) ? 3. / 16. : 3. / (8. * n);
+    newPositions[v] = (1 - n * u) * pos[v];
+    for (Vertex w : v.adjacentVertices()) {
+      newPositions[v] += u * pos[w];
+    }
+  }
+
+  EdgeData<Vector3> splitVertexPositions(mesh);
+  for (Edge e : mesh.edges()) {
+    if (e.isBoundary()) {
+      splitVertexPositions[e] = (pos[e.halfedge().tailVertex()] + pos[e.halfedge().tipVertex()]) / 2.;
+    } else {
+      std::array<Vertex, 2> ends{e.halfedge().tailVertex(), e.halfedge().tipVertex()};
+      std::array<Vertex, 2> sides{e.halfedge().next().next().vertex(), e.halfedge().twin().next().next().vertex()};
+      splitVertexPositions[e] = 3. / 8. * (pos[ends[0]] + pos[ends[1]]) + 1. / 8. * (pos[sides[0]] + pos[sides[1]]);
+    }
+  }
+
+  // Subdivide edges
+  std::vector<Edge> toFlip;
+  for (Edge e : mesh.edges()) {
+    if (!isOrigEdge[e]) continue;
+
+    // gather both vertices incident on the edge
+    Vertex oldA = e.halfedge().tipVertex();
+    Vertex oldB = e.halfedge().tailVertex();
+
+    Vector3 newPos = splitVertexPositions[e];
+
+    // split the edge
+    Vertex newV = mm.splitEdge(e, newPos).vertex();
+    isOrigVert[newV] = false;
+    newPositions[newV] = newPos;
+
+    // iterate through the edges incident on the new vertex
+    for (Edge e : newV.adjacentEdges()) {
+      isOrigEdge[e] = false;               // mark the new edges
+      Vertex otherV = e.otherVertex(newV); // other side of edge
+
+      // if this is a new edge between an old an new vertex, save for
+      // flipping
+      if (isOrigVert[otherV] && otherV != oldA && otherV != oldB) {
+        toFlip.push_back(e);
+      }
+    }
+  }
+
+  for (Edge e : toFlip) {
+    mesh.flip(e);
+  }
+
+  mesh.compress();
+  geo.inputVertexPositions = newPositions;
+  geo.refreshQuantities();
+}
+
+} // namespace surface
+} // namespace geometrycentral