diff --git a/grudge/geometry/metrics.py b/grudge/geometry/metrics.py
index 1ce02ce2e64797c2b3d08e7436de0c41a520d5ea..24c5f8a28fe9320993a82b0d876b2a1bf542abde 100644
--- a/grudge/geometry/metrics.py
+++ b/grudge/geometry/metrics.py
@@ -414,7 +414,8 @@ def inverse_surface_metric_derivative(
 
 
 def inverse_surface_metric_derivative_mat(
-        actx: ArrayContext, dcoll: DiscretizationCollection, dd=None):
+        actx: ArrayContext, dcoll: DiscretizationCollection, dd=None,
+        *, times_area_element=False):
     r"""Computes the matrix of inverse surface metric derivatives, indexed by
     ``(xyz_axis, rst_axis)``. It returns all values of
     :func:`inverse_surface_metric_derivative_mat` in cached matrix form.
@@ -423,16 +424,25 @@ def inverse_surface_metric_derivative_mat(
 
     :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one.
         Defaults to the base volume discretization.
+    :arg times_area_element: If *True*, each entry of the matrix is premultiplied
+        with the value of :func:`area_element`, reflecting the typical use
+        of the matrix in integrals evaluating weak derivatives.
     :returns: a :class:`~meshmode.dof_array.DOFArray` containing the
         inverse metric derivatives in per-group arrays of shape
         ``(xyz_dimension, rst_dimension, nelements, ndof)``.
     """
 
-    @memoize_in(dcoll, (inverse_surface_metric_derivative_mat, dd))
+    @memoize_in(dcoll, (inverse_surface_metric_derivative_mat, dd,
+        times_area_element))
     def _inv_surf_metric_deriv():
+        if times_area_element:
+            multiplier = area_element(actx, dcoll, dd=dd)
+        else:
+            multiplier = 1
+
         mat = actx.np.stack([
                 actx.np.stack(
-                    [inverse_surface_metric_derivative(actx, dcoll,
+                    [multiplier * inverse_surface_metric_derivative(actx, dcoll,
                         rst_axis, xyz_axis, dd=dd)
                         for rst_axis in range(dcoll.dim)])
                 for xyz_axis in range(dcoll.ambient_dim)])
diff --git a/grudge/op.py b/grudge/op.py
index 5176a710eb8ff52a1874f06cac97ace9ba438b33..4aac98815a89eaf63dc138d2279d8f0017e4b575 100644
--- a/grudge/op.py
+++ b/grudge/op.py
@@ -306,32 +306,30 @@ def weak_local_grad(dcoll: DiscretizationCollection, *args, nested=False):
         else:
             return np.stack(grad, axis=0)
 
-    from grudge.geometry import \
-        inverse_surface_metric_derivative_mat, area_element
+    from grudge.geometry import inverse_surface_metric_derivative_mat
 
     in_discr = dcoll.discr_from_dd(dd_in)
     out_discr = dcoll.discr_from_dd(dof_desc.DD_VOLUME)
 
     actx = vec.array_context
-    area_elements = area_element(actx, dcoll, dd=dd_in)
-    inverse_jac_mat = inverse_surface_metric_derivative_mat(actx, dcoll, dd=dd_in)
+    inverse_jac_mat = inverse_surface_metric_derivative_mat(actx, dcoll, dd=dd_in,
+            times_area_element=True)
 
     per_group_grads = [
         # r for rst axis
         # x for xyz axis
-        actx.einsum("rij,ej,ej,xrej->xei",
+        actx.einsum("rij,ej,xrej->xei",
                     reference_stiffness_transpose_matrix(
                         actx,
                         out_element_group=out_grp,
                         in_element_group=in_grp
                     ),
-                    ae_i,
                     vec_i,
                     ijm_i,
-                    arg_names=("ref_stiffT_mat", "jac", "vec", "inv_jac_t"),
+                    arg_names=("ref_stiffT_mat", "vec", "inv_jac_t"),
                     tagged=(FirstAxisIsElementsTag(),))
-        for out_grp, in_grp, vec_i, ae_i, ijm_i in zip(
-            out_discr.groups, in_discr.groups, vec, area_elements,
+        for out_grp, in_grp, vec_i, ijm_i in zip(
+            out_discr.groups, in_discr.groups, vec,
             inverse_jac_mat)]
 
     return make_obj_array([
@@ -374,33 +372,31 @@ def weak_local_d_dx(dcoll: DiscretizationCollection, *args):
     else:
         raise TypeError("invalid number of arguments")
 
-    from grudge.geometry import \
-        inverse_surface_metric_derivative_mat, area_element
+    from grudge.geometry import inverse_surface_metric_derivative_mat
 
     in_discr = dcoll.discr_from_dd(dd_in)
     out_discr = dcoll.discr_from_dd(dof_desc.DD_VOLUME)
 
     actx = vec.array_context
-    area_elements = area_element(actx, dcoll, dd=dd_in)
-    inverse_jac_mat = inverse_surface_metric_derivative_mat(actx, dcoll, dd=dd_in)
+    inverse_jac_mat = inverse_surface_metric_derivative_mat(actx, dcoll, dd=dd_in,
+            times_area_element=True)
 
     return DOFArray(
         actx,
         data=tuple(
-            actx.einsum("dij,ej,ej,dej->ei",
+            actx.einsum("dij,ej,dej->ei",
                         reference_stiffness_transpose_matrix(
                             actx,
                             out_element_group=out_grp,
                             in_element_group=in_grp
                         ),
-                        ae_i,
                         vec_i,
                         ijm_i[xyz_axis],
-                        arg_names=("ref_stiffT_mat", "jac", "vec", "inv_jac_t"),
+                        arg_names=("ref_stiffT_mat", "vec", "inv_jac_t"),
                         tagged=(FirstAxisIsElementsTag(),))
 
-            for out_grp, in_grp, vec_i, ae_i, ijm_i in zip(
-                out_discr.groups, in_discr.groups, vec, area_elements,
+            for out_grp, in_grp, vec_i, ijm_i in zip(
+                out_discr.groups, in_discr.groups, vec,
                 inverse_jac_mat)))