diff --git a/grudge/execution.py b/grudge/execution.py
index 04bc71d911525c43916e0dd00ced7d14e91b001c..0b5b21e545a6ceec57f353ab57792e7233cdab11 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -239,8 +239,68 @@ class ExecutionMapper(mappers.Evaluator,
         return self.discr.opposite_face_connection(qtag)(
                 self.queue, self.rec(field_expr)).with_queue(self.queue)
 
-    def map_face_mass_operator(self, op, field_expr):
-        raise NotImplementedError
+    # {{{ face mass operator
+
+    def map_ref_face_mass_operator(self, op, field_expr):
+        field = self.rec(field_expr)
+
+        from grudge.tools import is_zero
+        if is_zero(field):
+            return 0
+
+        @memoize_in(self.bound_op, "face_mass_knl")
+        def knl():
+            knl = lp.make_kernel(
+                """{[k,i,f,j]:
+                    0<=k<nelements and
+                    0<=f<nfaces and
+                    0<=i<nvol_nodes and
+                    0<=j<nface_nodes}""",
+                "result[k,i] = sum(f, sum(j, mat[i, f, j] * vec[f, k, j]))",
+                default_offset=lp.auto, name="face_mass")
+
+            knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
+            return lp.tag_inames(knl, dict(k="g.0"))
+
+        qtag = op.dd_in.quadrature_tag
+        if qtag is None:
+            # FIXME: Remove once proper quadrature support arrives
+            qtag = sym.QTAG_NONE
+
+        all_faces_conn = self.discr.all_faces_volume_connection(qtag)
+        all_faces_discr = all_faces_conn.to_discr
+        vol_discr = all_faces_conn.from_discr
+
+        result = vol_discr.empty(
+                queue=self.queue,
+                dtype=field.dtype, allocator=self.bound_op.allocator)
+
+        assert len(all_faces_discr.groups) == len(vol_discr.groups)
+
+        for cgrp, afgrp, volgrp in zip(all_faces_conn.groups,
+                all_faces_discr.groups, vol_discr.groups):
+            cache_key = "face_mass", afgrp, op, field.dtype
+
+            nfaces = volgrp.mesh_el_group.nfaces
+
+            try:
+                matrix = self.bound_op.operator_data_cache[cache_key]
+            except KeyError:
+                matrix = op.matrix(afgrp, volgrp, field.dtype)
+                matrix = (
+                        cl.array.to_device(self.queue, matrix)
+                        .with_queue(None))
+
+                self.bound_op.operator_data_cache[cache_key] = matrix
+
+            input_view = afgrp.view(field).reshape(
+                    nfaces, volgrp.nelements, afgrp.nunit_nodes)
+            knl()(self.queue, mat=matrix, result=volgrp.view(result),
+                    vec=input_view)
+
+        return result
+
+    # }}}
 
     # }}}
 
diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py
index 06518a3be1cbb2ffa4be3004ccd8d11a0611f456..d7f3e71c567821f164574d33cc671203526ddcf0 100644
--- a/grudge/symbolic/mappers/__init__.py
+++ b/grudge/symbolic/mappers/__init__.py
@@ -149,6 +149,7 @@ class OperatorReducerMixin(LocalOpReducerMixin, FluxOpReducerMixin):
 
     map_opposite_interior_face_swap = _map_op_base
     map_face_mass_operator = _map_op_base
+    map_ref_face_mass_operator = _map_op_base
 
 
 class CombineMapperMixin(object):
@@ -196,6 +197,7 @@ class IdentityMapperMixin(LocalOpReducerMixin, FluxOpReducerMixin):
 
     map_opposite_interior_face_swap = map_elementwise_linear
     map_face_mass_operator = map_elementwise_linear
+    map_ref_face_mass_operator = map_elementwise_linear
 
     # }}}
 
@@ -521,6 +523,12 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper):
             return op.RefInverseMassOperator(expr.op.dd_in, expr.op.dd_out)(
                 1/jac_in * self.rec(expr.field))
 
+        elif isinstance(expr.op, op.FaceMassOperator):
+            jac_in_surf = sym.area_element(self.ambient_dim, self.dim - 1,
+                    dd=expr.op.dd_in)
+            return op.RefFaceMassOperator(expr.op.dd_in, expr.op.dd_out)(
+                    jac_in_surf * self.rec(expr.field))
+
         elif isinstance(expr.op, op.StiffnessOperator):
             return op.RefMassOperator()(jac_noquad *
                     self.rec(
@@ -652,6 +660,9 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
     def map_face_mass_operator(self, expr, enclosing_prec):
         return "FaceM" + self._format_op_dd(expr)
 
+    def map_ref_face_mass_operator(self, expr, enclosing_prec):
+        return "RefFaceM" + self._format_op_dd(expr)
+
     def map_opposite_interior_face_swap(self, expr, enclosing_prec):
         return "OppSwap" + self._format_op_dd(expr)
 
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index 631fac36f31c1ed39a95947f219e19da10cfb862..32d82a2770fdd208bb16d10790ca0178f9d1477f 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -383,7 +383,7 @@ class OppositeInteriorFaceSwap(Operator):
     mapper_method = intern("map_opposite_interior_face_swap")
 
 
-class FaceMassOperator(ElementwiseLinearOperator):
+class FaceMassOperatorBase(ElementwiseLinearOperator):
     def __init__(self, dd_in=None, dd_out=None):
         sym = _sym()
 
@@ -391,17 +391,53 @@ class FaceMassOperator(ElementwiseLinearOperator):
             dd_in = sym.DOFDesc(sym.FRESTR_ALL_FACES, None)
 
         if dd_out is None:
-            dd_out = sym.DOFDesc("vol", None)
+            dd_out = sym.DOFDesc("vol", sym.QTAG_NONE)
+        if dd_out.uses_quadrature():
+            raise ValueError("face mass operator outputs are not on "
+                    "quadrature grids")
 
         if not dd_out.is_volume():
             raise ValueError("dd_out must be a volume domain")
         if dd_in.domain_tag is not sym.FRESTR_ALL_FACES:
             raise ValueError("dd_in must be an interior faces domain")
 
-        super(FaceMassOperator, self).__init__(dd_in, dd_out)
+        super(FaceMassOperatorBase, self).__init__(dd_in, dd_out)
 
+
+class FaceMassOperator(FaceMassOperatorBase):
     mapper_method = intern("map_face_mass_operator")
 
+
+class RefFaceMassOperator(ElementwiseLinearOperator):
+    def matrix(self, afgrp, volgrp, dtype):
+        nfaces = volgrp.mesh_el_group.nfaces
+        assert afgrp.nelements == nfaces * volgrp.nelements
+
+        matrix = np.empty(
+                (volgrp.nunit_nodes,
+                    nfaces,
+                    afgrp.nunit_nodes),
+                dtype=dtype)
+
+        from modepy.tools import UNIT_VERTICES
+        import modepy as mp
+        for iface, fvi in enumerate(
+                volgrp.mesh_el_group.face_vertex_indices()):
+            face_vertices = UNIT_VERTICES[volgrp.dim][np.array(fvi)].T
+            matrix[:, iface, :] = mp.nodal_face_mass_matrix(
+                    volgrp.basis(), volgrp.unit_nodes, afgrp.unit_nodes,
+                    volgrp.order,
+                    face_vertices)
+
+        # np.set_printoptions(linewidth=200, precision=3)
+        # matrix[np.abs(matrix) < 1e-13] = 0
+        # print(matrix)
+        # 1/0
+
+        return matrix
+
+    mapper_method = intern("map_ref_face_mass_operator")
+
 # }}}