diff --git a/grudge/discretization.py b/grudge/discretization.py
index ce130aef8dbad0cb165653774960ae96b722ea8c..e700e1812f705ae6fd5aedfd93df41b5c3ebf3b6 100644
--- a/grudge/discretization.py
+++ b/grudge/discretization.py
@@ -266,7 +266,8 @@ class DGDiscretizationWithBoundaries:
                         make_same_mesh_connection
                 to_discr = self._quad_volume_discr(to_dd.quadrature_tag)
                 from_discr = self._volume_discr
-                return make_same_mesh_connection(to_discr, from_discr)
+                return make_same_mesh_connection(self._setup_actx, to_discr,
+                            from_discr)
 
             else:
                 raise ValueError("cannot interpolate from volume to: " + str(to_dd))
diff --git a/grudge/eager.py b/grudge/eager.py
index 4079cdd0f8be6f5f8b4275623ca39bff0a694850..76c366ebf57a5c25f81ddba73527357a7c965ca6 100644
--- a/grudge/eager.py
+++ b/grudge/eager.py
@@ -58,6 +58,7 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries):
     .. automethod:: weak_div
 
     .. automethod:: normal
+    .. automethod:: mass
     .. automethod:: inverse_mass
     .. automethod:: face_mass
 
@@ -83,6 +84,11 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries):
         :arg tgt: a :class:`~grudge.sym.DOFDesc`, or a value convertible to one
         :arg vec: a :class:`~meshmode.dof_array.DOFArray`
         """
+        src = sym.as_dofdesc(src)
+        tgt = sym.as_dofdesc(tgt)
+        if src == tgt:
+            return vec
+
         if isinstance(vec, np.ndarray):
             return obj_array_vectorize(
                     lambda el: self.project(src, tgt, el), vec)
@@ -262,6 +268,26 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries):
                     local_only=True)
                 (array_context=actx))
 
+    @memoize_method
+    def _bound_mass(self, dd):
+        return bind(self, sym.MassOperator(dd_in=dd)(sym.Variable("u", dd)),
+                local_only=True)
+
+    def mass(self, *args):
+        if len(args) == 1:
+            vec, = args
+            dd = sym.DOFDesc("vol", sym.QTAG_NONE)
+        elif len(args) == 2:
+            dd, vec = args
+        else:
+            raise TypeError("invalid number of arguments")
+
+        if isinstance(vec, np.ndarray):
+            return obj_array_vectorize(
+                    lambda el: self.mass(dd, el), vec)
+
+        return self._bound_mass(dd)(u=vec)
+
     @memoize_method
     def _bound_inverse_mass(self):
         return bind(self, sym.InverseMassOperator()(sym.Variable("u")),
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index 5155cbb6ad213691f56054001b697e505db40878..b6db10eb62869217a201bdbd58bc0da9bfe72956 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -486,7 +486,7 @@ class MassOperatorBase(Operator):
         if dd_in is None:
             dd_in = prim.DD_VOLUME
         if dd_out is None:
-            dd_out = dd_in
+            dd_out = prim.DD_VOLUME
 
         super().__init__(dd_in, dd_out)
 
@@ -705,7 +705,7 @@ def integral(arg, dd=None):
 
     return NodalSum(dd)(
             arg * prim.cse(
-                MassOperator(dd_in=dd)(prim.Ones(dd)),
+                MassOperator(dd_in=dd, dd_out=dd)(prim.Ones(dd)),
                 "mass_quad_weights",
                 prim.cse_scope.DISCRETIZATION))
 
@@ -764,7 +764,7 @@ def h_max_from_volume(ambient_dim, dim=None, dd=None):
 
     return NodalMax(dd_in=dd)(
             ElementwiseSumOperator(dd)(
-                MassOperator(dd_in=dd)(prim.Ones(dd))
+                MassOperator(dd_in=dd, dd_out=dd)(prim.Ones(dd))
                 )
             )**(1.0/dim)
 
@@ -785,7 +785,7 @@ def h_min_from_volume(ambient_dim, dim=None, dd=None):
 
     return NodalMin(dd_in=dd)(
             ElementwiseSumOperator(dd)(
-                MassOperator(dd_in=dd)(prim.Ones(dd))
+                MassOperator(dd_in=dd, dd_out=dd)(prim.Ones(dd))
                 )
             )**(1.0/dim)