diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index e1335bbfe2b382b088b708b27e327748608a1a1f..b6b5b73ac60fbb6d4a67a15f3e31b58b7fede1af 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -659,22 +659,21 @@ class RefFaceMassOperator(ElementwiseLinearOperator):
         vol_basis = volgrp.basis_obj()
         faces = mp.faces_for_shape(volgrp.shape)
 
-        # NOTE: Assumes each face has the same quadrature/basis
-        face_quadrature = mp.quadrature_for_space(
-            # FIXME: How do we want to handle this in general?
-            # Do we always want a quadrature rule that exactly
-            # integrates 2N polynomials on the faces?
-            # For tensor-product and simplices?
-            mp.space_for_shape(faces[-1], 2*afgrp.order),
-            faces[-1]
-        )
-        face_basis = mp.basis_for_space(
-            mp.space_for_shape(faces[-1], afgrp.order),
-            faces[-1]
-        )
-
-        for iface, fvi in enumerate(
-                volgrp.mesh_el_group.face_vertex_indices()):
+        for iface, face in enumerate(faces):
+            # NOTE: Quadrature degree determined such that
+            # we either guarantee exact integration of degree 2N
+            # polynomials, or use a higher-order quadrature rule
+            # provided by the user (for example, when doing overintegration)
+            face_quad_degree = max(2*volgrp.order, afgrp.order)
+            face_quadrature = mp.quadrature_for_space(
+                mp.space_for_shape(face, face_quad_degree),
+                face
+            )
+            assert face_quadrature.exact_to >= afgrp.order
+            face_basis = mp.basis_for_space(
+                mp.space_for_shape(face, afgrp.order),
+                face
+            )
             matrix[:, iface, :] = mp.nodal_mass_matrix_for_face(
                 faces[iface], face_quadrature,
                 face_basis.functions, vol_basis.functions,