diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py
index ddc7535288cd82c1bbc56cb8decd9b0644cea73a..0a69967e8613adc76a9ff1d93ba02db90fb2dc9a 100644
--- a/examples/advection/var-velocity.py
+++ b/examples/advection/var-velocity.py
@@ -37,12 +37,13 @@ def main(write_output=True, order=4):
 
     dim = 2
 
+    n_divs = 15
     from meshmode.mesh.generation import generate_regular_rect_mesh
     mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5),
-            n=(15, 15), order=order)
+            n=(n_divs, n_divs), order=order)
 
-    dt_factor = 4
-    h = 1/20
+    dt_factor = 5
+    h = 1/n_divs
 
     discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
 
diff --git a/examples/quad/flux_quad.py b/examples/quad/flux_quad.py
new file mode 100644
index 0000000000000000000000000000000000000000..3d7f9a322d39176626f39bd1ea19267e0d42a43a
--- /dev/null
+++ b/examples/quad/flux_quad.py
@@ -0,0 +1,91 @@
+# Copyright (C) 2007 Andreas Kloeckner
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+
+import numpy as np
+import pyopencl as cl  # noqa
+import pyopencl.array  # noqa
+import pyopencl.clmath  # noqa
+
+import pytest  # noqa
+
+from pyopencl.tools import (  # noqa
+                pytest_generate_tests_for_pyopencl as pytest_generate_tests)
+
+import logging
+logger = logging.getLogger(__name__)
+
+from grudge import sym, bind, DGDiscretizationWithBoundaries
+
+import numpy.linalg as la
+
+
+
+
+def main(order=1):
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+
+    dim = 2
+
+
+    from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh
+    #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9]))
+    mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5),
+            n=(3, 3), order=7)
+
+
+
+
+    pquad_dd = sym.DOFDesc("vol", "product")
+    
+    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} )
+
+
+    magic_thing = sym.DOFDesc(sym.FACE_RESTR_ALL, "product")
+    #magic_thing = sym.DOFDesc(sym.FACE_RESTR_ALL)
+    def fluxy_stuff(u):
+        return sym.interp(u.dd, magic_thing)(u.int)
+
+
+    ones = bind(discr, sym.nodes(dim)[0] / sym.nodes(dim)[0])(queue, t=0)
+    nds = sym.nodes(dim)
+    op = sym.FaceMassOperator(magic_thing, "vol")(fluxy_stuff(sym.int_tpair(nds, "product")))
+    #op = (sym.int_tpair(nds, "product").int)
+    #op = sym.FaceMassOperator()(fluxy_stuff(sym.int_tpair(nds)))
+    print(op)
+
+    print(sym.pretty(op))
+    u = bind(discr, op)(queue, t=0)
+
+
+    print(u)
+
+
+
+
+    #vis.write_vtk_file("fld-000o.vtu", [  ("u", o) ])
+    #vis.write_vtk_file("fld-000u.vtu", [  ("u", u) ])
+
+
+
+
+
+
+
+if __name__ == "__main__":
+    main(4)
+
+
diff --git a/examples/quad/mass_quad.py b/examples/quad/mass_quad.py
new file mode 100644
index 0000000000000000000000000000000000000000..05c306bd9af5bd1838a2b09e325c663163f425b5
--- /dev/null
+++ b/examples/quad/mass_quad.py
@@ -0,0 +1,90 @@
+# Copyright (C) 2007 Andreas Kloeckner
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+
+import numpy as np
+import pyopencl as cl  # noqa
+import pyopencl.array  # noqa
+import pyopencl.clmath  # noqa
+
+import pytest  # noqa
+
+from pyopencl.tools import (  # noqa
+                pytest_generate_tests_for_pyopencl as pytest_generate_tests)
+
+import logging
+logger = logging.getLogger(__name__)
+
+from grudge import sym, bind, DGDiscretizationWithBoundaries
+
+import numpy.linalg as la
+
+
+
+
+def main(order=1):
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+
+    dim = 2
+
+
+    from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh
+    #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9]))
+    mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5),
+            n=(3, 3), order=7)
+
+
+
+
+    pquad_dd = sym.DOFDesc("vol", "product")
+    to_pquad = sym.interp("vol", pquad_dd)
+
+    def f(x):
+        x = x[0] / x[0]
+        #return sym.MassOperator()(x)
+        return sym.MassOperator(pquad_dd, "vol")(to_pquad(x))
+        #return np.dot(from_pquad_stiffness_t, to_pquad(x))
+        #return sym.interp("vol", pquad_dd)(sym.sin(3*x[0]))
+    
+    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} )
+
+
+
+
+
+
+
+    ones = bind(discr, sym.nodes(dim)[0] / sym.nodes(dim)[0])(queue, t=0)
+    u = bind(discr, f(sym.nodes(dim)))(queue, t=0)
+
+    print(np.dot(ones.T, u))
+
+
+
+
+    #vis.write_vtk_file("fld-000o.vtu", [  ("u", o) ])
+    #vis.write_vtk_file("fld-000u.vtu", [  ("u", u) ])
+
+
+
+
+
+
+
+if __name__ == "__main__":
+    main(4)
+
+
diff --git a/examples/quad/quad.py b/examples/quad/quad.py
new file mode 100644
index 0000000000000000000000000000000000000000..67c41fa234fe939af00af64c022d3979035f319d
--- /dev/null
+++ b/examples/quad/quad.py
@@ -0,0 +1,80 @@
+# Copyright (C) 2007 Andreas Kloeckner
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+
+import numpy as np
+import pyopencl as cl  # noqa
+import pyopencl.array  # noqa
+import pyopencl.clmath  # noqa
+
+import pytest  # noqa
+
+from pyopencl.tools import (  # noqa
+                pytest_generate_tests_for_pyopencl as pytest_generate_tests)
+
+import logging
+logger = logging.getLogger(__name__)
+
+from grudge import sym, bind, DGDiscretizationWithBoundaries
+
+import numpy.linalg as la
+
+
+
+
+def main(order=1):
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+
+    dim = 2
+
+
+    from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh
+    #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9]))
+    mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5),
+            n=(3, 3), order=7)
+
+
+
+
+    pquad_dd = sym.DOFDesc("vol", "product")
+
+    def f(x):
+        #return sym.sin(3*x[0])
+        return sym.interp("vol", pquad_dd)(sym.sin(3*x[0]))
+    
+    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} )
+
+    u = bind(discr, f(sym.nodes(dim)))(queue, t=0)
+
+
+
+    from meshmode.discretization.visualization import make_visualizer
+    #vis = make_visualizer(queue, discr.quad_volume_discr['product'], 2*order)
+    #vis = make_visualizer(queue, discr.volume_discr, order)
+
+    #vis.write_vtk_file("fld-000.vtu", [  ("u", u) ])
+    print(u)
+
+
+
+
+
+
+
+if __name__ == "__main__":
+    main(4)
+
+
diff --git a/examples/quad/stiffness_quad.py b/examples/quad/stiffness_quad.py
new file mode 100644
index 0000000000000000000000000000000000000000..ff334c4ec0073f1640ad936475ae96088797ecc6
--- /dev/null
+++ b/examples/quad/stiffness_quad.py
@@ -0,0 +1,99 @@
+# Copyright (C) 2007 Andreas Kloeckner
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+
+import numpy as np
+import pyopencl as cl  # noqa
+import pyopencl.array  # noqa
+import pyopencl.clmath  # noqa
+
+import pytest  # noqa
+
+from pyopencl.tools import (  # noqa
+                pytest_generate_tests_for_pyopencl as pytest_generate_tests)
+
+import logging
+logger = logging.getLogger(__name__)
+
+from grudge import sym, bind, DGDiscretizationWithBoundaries
+
+import numpy.linalg as la
+
+
+
+
+def main(order=1):
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+
+    dim = 2
+
+
+    from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh
+    #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9]))
+    mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5),
+            n=(3, 3), order=7)
+
+
+
+
+    pquad_dd = sym.DOFDesc("vol", "product")
+    to_pquad = sym.interp("vol", pquad_dd)
+    from_pquad_stiffness_t = sym.stiffness_t(dim, pquad_dd, "vol")
+    from_pquad_mass = sym.MassOperator()
+
+    def lin(x):
+        return x[0]
+
+    def non_lin(x):
+        return x[1] * x[0] * x[0] + 3* x[0] * x[0] - 4
+
+    def f(x, fun):
+        #return sym.MassOperator()(x)
+        #return np.dot(sym.stiffness_t(2), to_pquad(fun(x)))
+        return np.dot(from_pquad_stiffness_t, fun(to_pquad(x)))
+
+    def reg(x, fun):
+        #return sym.MassOperator()(x)
+        #return np.dot(sym.stiffness_t(2), to_pquad(fun(x)))
+        return np.dot(sym.stiffness_t(2), fun(x))
+    
+    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} )
+
+
+    #ones = bind(discr, sym.nodes(dim)[0] / sym.nodes(dim)[0])(queue, t=0)
+    nds = sym.nodes(dim)
+    #u = bind(discr, to_pquad(nds))(queue, t=0)
+    #u = bind(discr, reg(nds, lin))(queue, t=0)
+    u = bind(discr, f(nds, lin))(queue, t=0)
+
+    print(u)
+
+
+
+
+    #vis.write_vtk_file("fld-000o.vtu", [  ("u", o) ])
+    #vis.write_vtk_file("fld-000u.vtu", [  ("u", u) ])
+
+
+
+
+
+
+
+if __name__ == "__main__":
+    main(4)
+
+
diff --git a/grudge/discretization.py b/grudge/discretization.py
index a1ca6303ce40f787ddce458438c88e983b8b94d3..05360ecdfbe5a53d254b29c077eda7f07e75f211 100644
--- a/grudge/discretization.py
+++ b/grudge/discretization.py
@@ -82,8 +82,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
 
         if dd.is_volume():
             if qtag is not sym.QTAG_NONE:
-                # FIXME
-                raise NotImplementedError("quadrature")
+                return self._quad_volume_discr(qtag)
             return self._volume_discr
 
         elif dd.domain_tag is sym.FACE_RESTR_ALL:
@@ -120,6 +119,12 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
                                 self._boundary_to_quad_connection(
                                     to_dd.domain_tag.tag, to_dd.quadrature_tag)])
 
+            elif to_dd.is_volume():
+                from meshmode.discretization.connection import 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)
+
             else:
                 raise ValueError("cannot interpolate from volume to: " + str(to_dd))
 
@@ -130,6 +135,24 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
                 raise ValueError(
                         "cannot interpolate from interior faces to: "
                         + str(to_dd))
+        elif from_dd.domain_tag is sym.FACE_RESTR_ALL:
+            if to_dd.domain_tag is sym.FACE_RESTR_ALL and qtag is not sym.QTAG_NONE:
+                from meshmode.discretization.connection import make_face_restriction, FACE_RESTR_ALL
+                return make_face_restriction(
+                        self._all_faces_discr(),
+                        self.group_factory_for_quadrature_tag(qtag),
+                        FACE_RESTR_ALL,
+
+                        # FIXME: This will need to change as soon as we support
+                        # pyramids or other elements with non-identical face
+                        # types.
+                        per_face_groups=False)
+
+                #return self._all_faces_connection(None, qtag)
+            else:
+                raise ValueError(
+                        "cannot interpolate from interior faces to: "
+                        + str(to_dd))
 
         elif from_dd.is_boundary():
             if to_dd.domain_tag is sym.FACE_RESTR_ALL and qtag is sym.QTAG_NONE:
@@ -155,7 +178,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
                 QuadratureSimplexGroupFactory
 
         if quadrature_tag is not sym.QTAG_NONE:
-            return QuadratureSimplexGroupFactory(order=self.order)
+            return QuadratureSimplexGroupFactory(order=self.quad_min_degrees[quadrature_tag])
         else:
             return PolynomialWarpAndBlendGroupFactory(order=self.order)
 
@@ -259,7 +282,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
         return self._all_faces_volume_connection(quadrature_tag).to_discr
 
     @memoize_method
-    def _all_faces_connection(self, boundary_tag):
+    def _all_faces_connection(self, boundary_tag, to_qtag=None):
         """Return a
         :class:`meshmode.discretization.connection.DiscretizationConnection`
         that goes from either
@@ -277,7 +300,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
         else:
             faces_conn = self._boundary_connection(boundary_tag)
 
-        return make_face_to_all_faces_embedding(faces_conn, self._all_faces_discr())
+        return make_face_to_all_faces_embedding(faces_conn, self._all_faces_discr(to_qtag))
 
     # }}}
 
diff --git a/grudge/execution.py b/grudge/execution.py
index e4aebe5ca18b495acf4e608ff9f4cdfd218f57b1..40dcb19a1178239b6b617dceb0b676d3a85fd468 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -198,36 +198,42 @@ class ExecutionMapper(mappers.Evaluator,
             knl = lp.make_kernel(
                 """{[k,i,j]:
                     0<=k<nelements and
-                    0<=i,j<ndiscr_nodes}""",
+                    0<=i<out_ndiscr_nodes and
+                    0<=j<in_ndiscr_nodes}""",
                 "result[k,i] = sum(j, mat[i, j] * vec[k, j])",
                 default_offset=lp.auto, name="diff")
 
             knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
             return lp.tag_inames(knl, dict(k="g.0"))
 
-        discr = self.discrwb.discr_from_dd(op.dd_in)
 
-        # FIXME: This shouldn't really assume that it's dealing with a volume
-        # input. What about quadrature? What about boundaries?
-        result = discr.empty(
+        in_discr = self.discrwb.discr_from_dd(op.dd_in)
+        out_discr = self.discrwb.discr_from_dd(op.dd_out)
+
+        result = out_discr.empty(
                 queue=self.queue,
                 dtype=field.dtype, allocator=self.bound_op.allocator)
 
-        for grp in discr.groups:
-            cache_key = "elwise_linear", grp, op, field.dtype
-            try:
-                matrix = self.bound_op.operator_data_cache[cache_key]
-            except KeyError:
-                matrix = (
+        for i in range(len(out_discr.groups)):
+            in_grp = in_discr.groups[i]
+            out_grp = out_discr.groups[i]
+	    # FIXME: This cache thing
+            #cache_key = "elwise_linear", grp, op, field.dtype
+            #try:
+                #matrix = self.bound_op.operator_data_cache[cache_key]
+            #except KeyError:
+            matrix = (
                         cl.array.to_device(
                             self.queue,
-                            np.asarray(op.matrix(grp), dtype=field.dtype))
+                            np.asarray(op.matrix(out_grp, in_grp), dtype=field.dtype))
                         .with_queue(None))
 
-                self.bound_op.operator_data_cache[cache_key] = matrix
+                #self.bound_op.operator_data_cache[cache_key] = matrix
+
+            knl()(self.queue, mat=matrix, result=out_grp.view(result),
+                    vec=in_grp.view(field))
+
 
-            knl()(self.queue, mat=matrix, result=grp.view(result),
-                    vec=grp.view(field))
 
         return result
 
@@ -365,7 +371,6 @@ class ExecutionMapper(mappers.Evaluator,
         # This should be unified with map_elementwise_linear, which should
         # be extended to support batching.
 
-        discr = self.discrwb.discr_from_dd(repr_op.dd_in)
 
         # FIXME: Enable
         # assert repr_op.dd_in == repr_op.dd_out
@@ -377,7 +382,8 @@ class ExecutionMapper(mappers.Evaluator,
                 """{[imatrix,k,i,j]:
                     0<=imatrix<nmatrices and
                     0<=k<nelements and
-                    0<=i,j<nunit_nodes}""",
+                     0<=i<o_nunit_nodes and
+                    0<=j<i_nunit_nodes}""",
                 """
                 result[imatrix, k, i] = sum(
                         j, diff_mat[imatrix, i, j] * vec[k, j])
@@ -388,26 +394,33 @@ class ExecutionMapper(mappers.Evaluator,
             return lp.tag_inames(knl, dict(k="g.0"))
 
         noperators = len(insn.operators)
-        result = discr.empty(
+
+        in_discr = self.discrwb.discr_from_dd(repr_op.dd_in)
+        out_discr = self.discrwb.discr_from_dd(repr_op.dd_out)
+
+        result = out_discr.empty(
                 queue=self.queue,
                 dtype=field.dtype, extra_dims=(noperators,),
                 allocator=self.bound_op.allocator)
 
-        for grp in discr.groups:
-            if grp.nelements == 0:
+        for i in range(len(out_discr.groups)):
+            in_grp = in_discr.groups[i]
+            out_grp = out_discr.groups[i]
+
+            if in_grp.nelements == 0:
                 continue
 
-            matrices = repr_op.matrices(grp)
+            matrices = repr_op.matrices(out_grp, in_grp)
 
             # FIXME: Should transfer matrices to device and cache them
             matrices_ary = np.empty((
-                noperators, grp.nunit_nodes, grp.nunit_nodes))
+                noperators, out_grp.nunit_nodes, in_grp.nunit_nodes))
             for i, op in enumerate(insn.operators):
                 matrices_ary[i] = matrices[op.rst_axis]
 
             knl()(self.queue,
                     diff_mat=matrices_ary,
-                    result=grp.view(result), vec=grp.view(field))
+                    result=out_grp.view(result), vec=in_grp.view(field))
 
         return [(name, result[i]) for i, name in enumerate(insn.names)], []
 
diff --git a/grudge/models/advection.py b/grudge/models/advection.py
index 1c83ed9675d5a45c75b22e3b4d5ca7926acf9133..e0203378214983cceb11ecb46d6230a2dfc7574e 100644
--- a/grudge/models/advection.py
+++ b/grudge/models/advection.py
@@ -175,7 +175,7 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator):
                 np.dot(sym.stiffness_t(self.ambient_dim), sym.cse(self.v*u))
                 - sym.FaceMassOperator()(
                     flux(sym.int_tpair(u))
-                    + flux(sym.bv_tpair(sym.BTAG_ALL, u, bc_in))
+                    #+ flux(sym.bv_tpair(sym.BTAG_ALL, u, bc_in))
 
                     # FIXME: Add back support for inflow/outflow tags
                     #+ flux(sym.bv_tpair(self.inflow_tag, u, bc_in))
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index 34ec47aaa556369ecea8e688f81bdcbaa0a6db3f..9115360453fffc498bc7e39de515be7f6d009cf9 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -97,6 +97,12 @@ class ElementwiseLinearOperator(Operator):
 
 
 class InterpolationOperator(Operator):
+    #def __init__(self, dd_in, dd_out):
+        #print("ayy")
+        #print(dd_in)
+        #print(dd_out)
+
+        #super().__init__(dd_in, dd_out)
     mapper_method = intern("map_interpolation")
 
 
@@ -143,6 +149,8 @@ class DiffOperatorBase(Operator):
             dd_in = _sym().DD_VOLUME
         if dd_out is None:
             dd_out = dd_in.with_qtag(_sym().QTAG_NONE)
+        if dd_out == "vol": #FIXME I'm pretty sure doing it this way is wrong
+            dd_out = _sym().DD_VOLUME
         if dd_out.uses_quadrature():
             raise ValueError("differentiation outputs are not on "
                     "quadrature grids")
@@ -215,18 +223,41 @@ class RefDiffOperator(RefDiffOperatorBase):
     mapper_method = intern("map_ref_diff")
 
     @staticmethod
-    def matrices(element_group):
-        return element_group.diff_matrices()
+    def matrices(out_element_group, in_element_group):
+        assert in_element_group == out_element_group
+        return in_element_group.diff_matrices()
 
 
 class RefStiffnessTOperator(RefDiffOperatorBase):
     mapper_method = intern("map_ref_stiffness_t")
 
     @staticmethod
-    def matrices(element_group):
-        assert element_group.is_orthogonal_basis()
-        mmat = element_group.mass_matrix()
-        return [dmat.T.dot(mmat.T) for dmat in element_group.diff_matrices()]
+    def matrices(out_element_group, in_element_group):
+        #FIXME: Is this right? Do we expose weird bugs if we delete this?
+        if in_element_group == out_element_group:
+            assert in_element_group.is_orthogonal_basis()
+            mmat = in_element_group.mass_matrix()
+            return [dmat.T.dot(mmat.T) for dmat in in_element_group.diff_matrices()]
+
+        n_quad_nodes = in_element_group.nunit_nodes
+        n_reg_nodes = out_element_group.nunit_nodes
+
+        from modepy import vandermonde
+        vand = vandermonde(out_element_group.basis(), out_element_group.unit_nodes)
+        grad_vand = vandermonde(out_element_group.grad_basis(), in_element_group.unit_nodes)
+        vand_inv_T = np.linalg.inv(vand).T
+
+        grad_basis = out_element_group.grad_basis()
+        weights = in_element_group.weights
+        nodes = in_element_group.unit_nodes
+
+        f_res = [np.zeros((n_reg_nodes, n_quad_nodes)) for d in range(out_element_group.dim)]
+        for d in range(out_element_group.dim):
+            for i in range(n_reg_nodes):
+                for j in range(n_quad_nodes):
+                    f_res[d][i,j] = weights[j] * np.dot(vand_inv_T[i,:], grad_vand[d][j,:])
+        return f_res
+
 
 # }}}
 
@@ -338,9 +369,6 @@ class MassOperatorBase(Operator):
         if dd_out is None:
             dd_out = dd_in
 
-        if dd_out != dd_in:
-            raise ValueError("dd_out and dd_in must be identical")
-
         super(MassOperatorBase, self).__init__(dd_in, dd_out)
 
 
@@ -358,19 +386,42 @@ class RefMassOperatorBase(ElementwiseLinearOperator):
 
 class RefMassOperator(RefMassOperatorBase):
     @staticmethod
-    def matrix(element_group):
-        return element_group.mass_matrix()
+    def matrix(out_element_group, in_element_group):
+        #FIXME: Is this normal? Will bugs surface if I remove this?
+        if out_element_group == in_element_group:
+            return element_group.mass_matrix()
+
+        n_quad_nodes = in_element_group.nunit_nodes
+        n_reg_nodes = out_element_group.nunit_nodes
+
+        from modepy import vandermonde
+        vand = vandermonde(out_element_group.basis(), out_element_group.unit_nodes)
+        o_vand = vandermonde(out_element_group.basis(), in_element_group.unit_nodes)
+        vand_inv_T = np.linalg.inv(vand).T
+
+        grad_basis = out_element_group.grad_basis()
+        weights = in_element_group.weights
+        nodes = in_element_group.unit_nodes
+
+        f_res = np.zeros((n_reg_nodes, n_quad_nodes))
+        for i in range(n_reg_nodes):
+            for j in range(n_quad_nodes):
+                f_res[i,j] = weights[j] * np.dot(vand_inv_T[i,:], o_vand[j,:])
+
+        return f_res
+
 
     mapper_method = intern("map_ref_mass")
 
 
 class RefInverseMassOperator(RefMassOperatorBase):
     @staticmethod
-    def matrix(element_group):
+    def matrix(in_element_group, out_element_group):
+        assert in_element_group == out_element_group
         import modepy as mp
         return mp.inverse_mass_matrix(
-                element_group.basis(),
-                element_group.unit_nodes)
+                in_element_group.basis(),
+                in_element_group.unit_nodes)
 
     mapper_method = intern("map_ref_inverse_mass")
 
@@ -405,7 +456,7 @@ class FaceMassOperatorBase(ElementwiseLinearOperator):
         if dd_in is None:
             dd_in = sym.DOFDesc(sym.FACE_RESTR_ALL, None)
 
-        if dd_out is None:
+        if dd_out is None or dd_out == "vol":
             dd_out = sym.DOFDesc("vol", sym.QTAG_NONE)
         if dd_out.uses_quadrature():
             raise ValueError("face mass operator outputs are not on "
@@ -476,10 +527,10 @@ def stiffness(dim):
         [StiffnessOperator(i) for i in range(dim)])
 
 
-def stiffness_t(dim):
+def stiffness_t(dim, dd_in=None, dd_out=None):
     from pytools.obj_array import make_obj_array
     return make_obj_array(
-        [StiffnessTOperator(i) for i in range(dim)])
+        [StiffnessTOperator(i, dd_in, dd_out) for i in range(dim)])
 
 
 def integral(arg, dd=None):
diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
index d64aef3a39f7927c80ccc85b48b5850d83c03bfb..91122a89e215ba26c6ea22b39f7efd21ea49f47e 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -653,10 +653,15 @@ class TracePair:
         return 0.5*(self.int + self.ext)
 
 
-def int_tpair(expression):
-    i = cse(_sym().interp("vol", "int_faces")(expression))
+def int_tpair(expression, qtag=None):
+    dd = _sym().DOFDesc("int_faces", qtag)
+    i = cse(_sym().interp("vol", dd)(expression))
     e = cse(_sym().OppositeInteriorFaceSwap()(i))
-    return TracePair("int_faces", i, e)
+    return TracePair(dd, i, e)
+
+    #i = cse(_sym().interp("vol", "int_faces")(expression))
+    #e = cse(_sym().OppositeInteriorFaceSwap()(i))
+    #return TracePair("int_faces", i, e)
 
 
 def bdry_tpair(dd, interior, exterior):