diff --git a/src/catmull_clark.cpp b/src/catmull_clark.cpp
index de8ce26..d34502d 100644
--- a/src/catmull_clark.cpp
+++ b/src/catmull_clark.cpp
@@ -2,6 +2,22 @@
 #include <unordered_map>
 #include <utility>
 #include <functional>
+#include <vector>
+
+/**
+ * Compute the centroids of the faces of a mesh.
+ *
+ * Inputs:
+ * @param V  #V by 3 list of vertex positions
+ * @param F  #F by 3 list of triangle mesh indices into V
+ * Outputs:
+ * @param FC  #F by 3 list of face centroids
+ */
+void compute_face_centroids(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  Eigen::MatrixXd & FC
+);
 
 void catmull_clark(
   const Eigen::MatrixXd & V,
@@ -10,9 +26,59 @@ void catmull_clark(
   Eigen::MatrixXd & SV,
   Eigen::MatrixXi & SF)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Replace with your code here:
+  Eigen::MatrixXd VF;
+
+  /* Calculate new vertices. */
+
+  // Calculate face centroids into VF.
+  compute_face_centroids(V, F, VF);
+
+  // Calculate edge points into VE.
+  // TODO
+
+  // Move original vertices to new positions into VV.
+  // TODO
+
+  /* Form edges and faces in the new mesh. */
+
+  // Connect each new face point (VF) to the new edge points (VE) of all original edges defining the original face.
+  // TODO
+
+  // Connect each new vertex point to the new edge points of all original edges incident on the original vertex.
+  // TODO
+
+  // Define the new faces as enclosed by edges.
+  for (int i = 0; )
+
+  // TODO Replace me
   SV = V;
   SF = F;
-  ////////////////////////////////////////////////////////////////////////////
+}
+
+void extract_edges(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  std::unordered_map<std::pair<int, int>, std::pair<int, int>, std::function<size_t(const std::pair<int, int> &)>> & edge_map
+) {
+  for (int i = 0; i < F.rows(); i++) {
+    for (int j = 0; j < F.cols(); j++) {
+      std::pair<int, int> edge(F(i, j), F(i, (j + 1) % F.cols()));
+      if (edge_map.find(edge) == edge_map.end()) {
+        edge_map[edge] = std::make_pair(i, -1);
+      } else {
+        edge_map[edge].second = i;
+      }
+    }
+  }
+}
+
+void compute_face_centroids(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  Eigen::MatrixXd & FC
+) {
+  FC.resize(F.rows(), 3);
+  for (int i = 0; i < F.rows(); i++) {
+    FC.row(i) = (V.row(F(i, 0)) + V.row(F(i, 1)) + V.row(F(i, 2))) / 3;
+  }
 }
diff --git a/src/vertex_triangle_adjacency.cpp b/src/vertex_triangle_adjacency.cpp
index 4b002fb..3f5148c 100644
--- a/src/vertex_triangle_adjacency.cpp
+++ b/src/vertex_triangle_adjacency.cpp
@@ -13,4 +13,3 @@ void vertex_triangle_adjacency(
     }
   }
 }
-