diff --git a/grudge/reductions.py b/grudge/reductions.py
index 6dcde1316..6f6108462 100644
--- a/grudge/reductions.py
+++ b/grudge/reductions.py
@@ -295,14 +295,23 @@ def integral(dcoll: DiscretizationCollection, dd, vec) -> Scalar:
         :class:`~arraycontext.ArrayContainer` of them.
     :returns: a device scalar denoting the evaluated integral.
     """
-    from grudge.op import _apply_mass_operator
-
     dd = dof_desc.as_dofdesc(dd)
+    discr = dcoll.discr_from_dd(dd)
 
-    ones = dcoll.discr_from_dd(dd).zeros(vec.array_context) + 1.0
-    return nodal_sum(
-        dcoll, dd, vec * _apply_mass_operator(dcoll, dd, dd, ones)
+    from grudge.geometry import area_element
+    actx = vec.array_context
+    area_elements = area_element(
+        actx, dcoll, dd=dd,
+        _use_geoderiv_connection=actx.supports_nonscalar_broadcasting)
+    qwts = DOFArray(
+        actx,
+        data=tuple(
+            actx.from_numpy(
+                np.tile(grp.quadrature_rule().weights,
+                        (ae_i.shape[0], 1)))*ae_i
+            for grp, ae_i in zip(discr.groups, area_elements))
     )
+    return nodal_sum(dcoll, dd, vec*qwts)
 
 # }}}
 
@@ -509,13 +518,22 @@ def elementwise_integral(
         raise TypeError("invalid number of arguments")
 
     dd = dof_desc.as_dofdesc(dd)
+    discr = dcoll.discr_from_dd(dd)
 
-    from grudge.op import _apply_mass_operator
-
-    ones = dcoll.discr_from_dd(dd).zeros(vec.array_context) + 1.0
-    return elementwise_sum(
-        dcoll, dd, vec * _apply_mass_operator(dcoll, dd, dd, ones)
+    from grudge.geometry import area_element
+    actx = vec.array_context
+    area_elements = area_element(
+        actx, dcoll, dd=dd,
+        _use_geoderiv_connection=actx.supports_nonscalar_broadcasting)
+    qwts = DOFArray(
+        actx,
+        data=tuple(
+            actx.from_numpy(
+                np.tile(grp.quadrature_rule().weights,
+                        (ae_i.shape[0], 1)))*ae_i
+            for grp, ae_i in zip(discr.groups, area_elements))
     )
+    return elementwise_sum(dcoll, dd, vec*qwts)
 
 # }}}
 
diff --git a/test/mesh_data.py b/test/mesh_data.py
index 58c9b80a7..26b222ecc 100644
--- a/test/mesh_data.py
+++ b/test/mesh_data.py
@@ -117,14 +117,20 @@ class _BoxMeshBuilderBase(MeshBuilder):
     a = (-0.5, -0.5, -0.5)
     b = (+0.5, +0.5, +0.5)
 
+    tpe: bool
+
+    def __init__(self, tpe=False):
+        self._tpe = tpe
+
     def get_mesh(self, resolution, mesh_order=4):
         if not isinstance(resolution, (list, tuple)):
             resolution = (resolution,) * self.ambient_dim
-
+        from meshmode.mesh import TensorProductElementGroup
+        group_cls = TensorProductElementGroup if self._tpe else None
         return mgen.generate_regular_rect_mesh(
                 a=self.a, b=self.b,
                 nelements_per_axis=resolution,
-                order=mesh_order)
+                order=mesh_order, group_cls=group_cls)
 
 
 class BoxMeshBuilder1D(_BoxMeshBuilderBase):
diff --git a/test/test_grudge.py b/test/test_grudge.py
index b0849310b..181557276 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -154,27 +154,35 @@ def _spheroid_surface_area(radius, aspect_ratio):
         return 2.0 * np.pi * radius**2 * (1 + (c/a)**2 / e * np.arctanh(e))
 
 
+# This tests uses op.integral which has been changed to directly use quadrature
+# weights instead of the mass matrix.  We either need to reimplement this test
+# to do what it was intended to do, or rename it appropriately.
 @pytest.mark.parametrize("name", [
-    "2-1-ellipse", "spheroid", "box2d", "box3d"
+    "box2d-tpe", "box3d-tpe", "box2d", "box3d", "2-1-ellipse", "spheroid",
     ])
-def test_mass_surface_area(actx_factory, name):
+@pytest.mark.parametrize("oi", [False, True])
+def test_mass_surface_area(actx_factory, name, oi):
+    from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD, as_dofdesc
     actx = actx_factory()
+    qtag = DISCR_TAG_QUAD if oi else DISCR_TAG_BASE
+    vol_dd_base = as_dofdesc(dof_desc.DTAG_VOLUME_ALL)
+    vol_dd_quad = vol_dd_base.with_discr_tag(qtag)
 
     # {{{ cases
 
     order = 4
-
+    tpe = name.endswith("-tpe")
     if name == "2-1-ellipse":
         builder = mesh_data.EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0)
         surface_area = _ellipse_surface_area(builder.radius, builder.aspect_ratio)
     elif name == "spheroid":
         builder = mesh_data.SpheroidMeshBuilder()
         surface_area = _spheroid_surface_area(builder.radius, builder.aspect_ratio)
-    elif name == "box2d":
-        builder = mesh_data.BoxMeshBuilder2D()
+    elif name.startswith("box2d"):
+        builder = mesh_data.BoxMeshBuilder2D(tpe=tpe)
         surface_area = 1.0
-    elif name == "box3d":
-        builder = mesh_data.BoxMeshBuilder3D()
+    elif name.startswith("box3d"):
+        builder = mesh_data.BoxMeshBuilder3D(tpe=tpe)
         surface_area = 1.0
     else:
         raise ValueError(f"unknown geometry name: {name}")
@@ -189,15 +197,20 @@ def test_mass_surface_area(actx_factory, name):
 
     for resolution in builder.resolutions:
         mesh = builder.get_mesh(resolution, order)
-        dcoll = make_discretization_collection(actx, mesh, order=order)
-        volume_discr = dcoll.discr_from_dd(dof_desc.DD_VOLUME_ALL)
+        dcoll = make_discretization_collection(
+                   actx, mesh,
+                   discr_tag_to_group_factory={
+                       DISCR_TAG_BASE: InterpolatoryEdgeClusteredGroupFactory(order),
+                       DISCR_TAG_QUAD: QuadratureGroupFactory(3 * order)
+                   })
+        volume_discr = dcoll.discr_from_dd(vol_dd_quad)
 
         logger.info("ndofs:     %d", volume_discr.ndofs)
         logger.info("nelements: %d", volume_discr.mesh.nelements)
 
         # {{{ compute surface area
 
-        dd = dof_desc.DD_VOLUME_ALL
+        dd = vol_dd_quad
         ones_volm = volume_discr.zeros(actx) + 1
         approx_surface_area = actx.to_numpy(op.integral(dcoll, dd, ones_volm))