diff --git a/include/geometrycentral/surface/embedded_geometry_interface.h b/include/geometrycentral/surface/embedded_geometry_interface.h
index 64d2710d..fed68671 100644
--- a/include/geometrycentral/surface/embedded_geometry_interface.h
+++ b/include/geometrycentral/surface/embedded_geometry_interface.h
@@ -28,6 +28,11 @@ class EmbeddedGeometryInterface : public ExtrinsicGeometryInterface {
   void requireVertexPositions();
   void unrequireVertexPositions();
 
+  // Face centroids
+  FaceData<Vector3> faceCentroids;
+  void requireFaceCentroids();
+  void unrequireFaceCentroids();
+
   // Face normal
   FaceData<Vector3> faceNormals;
   void requireFaceNormals();
@@ -116,6 +121,9 @@ class EmbeddedGeometryInterface : public ExtrinsicGeometryInterface {
   DependentQuantityD<VertexData<Vector3>> vertexPositionsQ;
   virtual void computeVertexPositions() = 0;
 
+  DependentQuantityD<FaceData<Vector3>> faceCentroidsQ;
+  virtual void computeFaceCentroids();
+
   DependentQuantityD<FaceData<Vector3>> faceNormalsQ;
   virtual void computeFaceNormals();
 
diff --git a/src/surface/embedded_geometry_interface.cpp b/src/surface/embedded_geometry_interface.cpp
index a32035b3..9b0dbd48 100644
--- a/src/surface/embedded_geometry_interface.cpp
+++ b/src/surface/embedded_geometry_interface.cpp
@@ -13,6 +13,7 @@ EmbeddedGeometryInterface::EmbeddedGeometryInterface(SurfaceMesh& mesh_) :
   ExtrinsicGeometryInterface(mesh_),
 
   vertexPositionsQ                (&vertexPositions,                std::bind(&EmbeddedGeometryInterface::computeVertexPositions, this),                quantities),
+  faceCentroidsQ                  (&faceCentroids,                  std::bind(&EmbeddedGeometryInterface::computeFaceCentroids, this),                  quantities),
   faceNormalsQ                    (&faceNormals,                    std::bind(&EmbeddedGeometryInterface::computeFaceNormals, this),                    quantities),
   vertexNormalsQ                  (&vertexNormals,                  std::bind(&EmbeddedGeometryInterface::computeVertexNormals, this),                  quantities),
   faceTangentBasisQ               (&faceTangentBasis,               std::bind(&EmbeddedGeometryInterface::computeFaceTangentBasis, this),               quantities),
@@ -76,6 +77,25 @@ void EmbeddedGeometryInterface::computeEdgeDihedralAngles() {
 void EmbeddedGeometryInterface::requireVertexPositions() { vertexPositionsQ.require(); }
 void EmbeddedGeometryInterface::unrequireVertexPositions() { vertexPositionsQ.unrequire(); }
 
+// Face centroids
+void EmbeddedGeometryInterface::computeFaceCentroids() {
+  vertexPositionsQ.ensureHave();
+
+  faceCentroids = FaceData<Vector3>(mesh);
+
+  for (Face f : mesh.faces()) {
+    Vector3 centroid = Vector3::zero();
+    size_t n = 0;
+    for (Vertex v : f.adjacentVertices()) {
+      centroid += vertexPositions[v];
+      n++;
+    }
+    centroid /= static_cast<double>(n);
+    faceCentroids[f] = centroid;
+  }
+}
+void EmbeddedGeometryInterface::requireFaceCentroids() { faceCentroidsQ.require(); }
+void EmbeddedGeometryInterface::unrequireFaceCentroids() { faceCentroidsQ.unrequire(); }
 
 void EmbeddedGeometryInterface::computeFaceNormals() {
   vertexPositionsQ.ensureHave();