diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py
index ddc7535288cd82c1bbc56cb8decd9b0644cea73a..00e4f0c8496dced20a21e808e68986ecec1d5be2 100644
--- a/examples/advection/var-velocity.py
+++ b/examples/advection/var-velocity.py
@@ -1,17 +1,26 @@
-# 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/>.
+from __future__ import division, absolute_import
+
+__copyright__ = "Copyright (C) 2017 Bogdan Enache"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
 
 
 import numpy as np
@@ -37,20 +46,19 @@ def main(write_output=True, order=4):
 
     dim = 2
 
+    resolution = 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=(resolution, resolution), order=order)
 
-    dt_factor = 4
-    h = 1/20
-
-    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
+    dt_factor = 5
+    h = 1/resolution
 
     sym_x = sym.nodes(2)
 
     advec_v = join_fields(-1*sym_x[1], sym_x[0]) / 2
 
-    flux_type = "central"
+    flux_type = "upwind"
 
     source_center = np.array([0.1, 0.1])
     source_width = 0.05
@@ -58,30 +66,46 @@ def main(write_output=True, order=4):
     sym_x = sym.nodes(2)
     sym_source_center_dist = sym_x - source_center
 
-    def f(x):
+    def f_gaussian(x):
         return sym.exp(
                     -np.dot(sym_source_center_dist, sym_source_center_dist)
                     / source_width**2)
 
+    def f_step(x):
+        return sym.If(
+                sym.Comparison(
+                    np.dot(sym_source_center_dist, sym_source_center_dist),
+                    "<",
+                    (4*source_width)**2),
+                1, 0)
+
     def u_analytic(x):
         return 0
 
     from grudge.models.advection import VariableCoefficientAdvectionOperator
+    from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory  # noqa
+
+    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order,
+            quad_tag_to_group_factory={
+                #"product": None,
+                "product": QuadratureSimplexGroupFactory(order=4*order)
+                })
 
-    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
     op = VariableCoefficientAdvectionOperator(2, advec_v,
-        inflow_u=u_analytic(sym.nodes(dim, sym.BTAG_ALL)),
+        u_analytic(sym.nodes(dim, sym.BTAG_ALL)), quad_tag="product",
         flux_type=flux_type)
 
-    bound_op = bind(discr, op.sym_operator())
+    bound_op = bind(
+            discr, op.sym_operator(),
+            #debug_flags=["dump_sym_operator_stages"]
+            )
 
-    u = bind(discr, f(sym.nodes(dim)))(queue, t=0)
+    u = bind(discr, f_gaussian(sym.nodes(dim)))(queue, t=0)
 
     def rhs(t, u):
         return bound_op(queue, t=t, u=u)
 
-    final_time = 2
-
+    final_time = 50
     dt = dt_factor * h/order**2
     nsteps = (final_time // dt) + 1
     dt = final_time/nsteps + 1e-15
@@ -90,7 +114,7 @@ def main(write_output=True, order=4):
     dt_stepper = set_up_rk4("u", dt, u, rhs)
 
     from grudge.shortcuts import make_visualizer
-    vis = make_visualizer(discr, vis_order=order)
+    vis = make_visualizer(discr, vis_order=2*order)
 
     step = 0
 
@@ -98,10 +122,11 @@ def main(write_output=True, order=4):
         if isinstance(event, dt_stepper.StateComputed):
 
             step += 1
+            if step % 30 == 0:
+                print(step)
 
-            #print(step, event.t, norm(queue, u=event.state_component[0]))
-            vis.write_vtk_file("fld-%04d.vtu" % step,
-                    [("u", event.state_component)])
+                vis.write_vtk_file("fld-%04d.vtu" % step,
+                        [("u", event.state_component)])
 
 
 if __name__ == "__main__":
diff --git a/grudge/discretization.py b/grudge/discretization.py
index 6f39762c00be545aa387253b9f39ea02f5dbeb8a..04cf43302a0da235a66e6329d0040ae86ceda761 100644
--- a/grudge/discretization.py
+++ b/grudge/discretization.py
@@ -1,6 +1,6 @@
 from __future__ import division, absolute_import
 
-__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
+__copyright__ = "Copyright (C) 2015-2017 Andreas Kloeckner, Bogdan Enache"
 
 __license__ = """
 Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -47,17 +47,21 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
     .. automethod :: zeros
     """
 
-    def __init__(self, cl_ctx, mesh, order, quad_min_degrees=None):
+    def __init__(self, cl_ctx, mesh, order, quad_tag_to_group_factory=None):
         """
-        :param quad_min_degrees: A mapping from quadrature tags to the degrees
-            to which the desired quadrature is supposed to be exact.
+        :param quad_min_degrees: A mapping from quadrature tags (typically
+            strings--but may be any hashable/comparable object) to a
+            :class:`meshmode.discretization.ElementGroupFactory` indicating with
+            which quadrature discretization the operations are to be carried out,
+            or *None* to indicate that operations with this quadrature tag should
+            be carried out with the standard volume discretization.
         """
 
-        if quad_min_degrees is None:
-            quad_min_degrees = {}
+        if quad_tag_to_group_factory is None:
+            quad_tag_to_group_factory = {}
 
         self.order = order
-        self.quad_min_degrees = quad_min_degrees
+        self.quad_tag_to_group_factory = quad_tag_to_group_factory
 
         from meshmode.discretization import Discretization
 
@@ -82,16 +86,26 @@ 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:
-            return self._all_faces_discr(qtag)
+        if qtag is not sym.QTAG_NONE:
+            no_quad_discr = self.discr_from_dd(sym.DOFDesc(dd.domain_tag))
+
+            from meshmode.discretization import Discretization
+            return Discretization(
+                    self._volume_discr.cl_context,
+                    no_quad_discr.mesh,
+                    self.group_factory_for_quadrature_tag(qtag))
+
+        assert qtag is sym.QTAG_NONE
+
+        if dd.domain_tag is sym.FACE_RESTR_ALL:
+            return self._all_faces_volume_connection().to_discr
         elif dd.domain_tag is sym.FACE_RESTR_INTERIOR:
-            return self._interior_faces_discr(qtag)
+            return self._interior_faces_connection().to_discr
         elif dd.is_boundary():
-            return self._boundary_discr(dd.domain_tag, qtag)
+            return self._boundary_connection(dd.domain_tag).to_discr
         else:
             raise ValueError("DOF desc tag not understood: " + str(dd))
 
@@ -100,18 +114,32 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
         from_dd = sym.as_dofdesc(from_dd)
         to_dd = sym.as_dofdesc(to_dd)
 
-        if from_dd.quadrature_tag is not sym.QTAG_NONE:
-            raise ValueError("cannot interpolate *from* a "
-                    "(non-interpolatory) quadrature grid")
-
         to_qtag = to_dd.quadrature_tag
 
+        if (
+                not from_dd.is_volume()
+                and from_dd.quadrature_tag == to_dd.quadrature_tag
+                and to_dd.domain_tag is sym.FACE_RESTR_ALL):
+            faces_conn = self.connection_from_dds(
+                    sym.DOFDesc("vol"),
+                    sym.DOFDesc(from_dd.domain_tag))
+
+            from meshmode.discretization.connection import \
+                    make_face_to_all_faces_embedding
+
+            return make_face_to_all_faces_embedding(
+                    faces_conn, self.discr_from_dd(to_dd),
+                    self.discr_from_dd(from_dd))
+
         # {{{ simplify domain + qtag change into chained
 
-        if (from_dd.domain_tag != to_dd.domain_tag
-                and from_dd.quadrature_tag != to_dd.quadrature_tag):
+        if (
+                from_dd.domain_tag != to_dd.domain_tag
+                and from_dd.quadrature_tag is sym.QTAG_NONE
+                and to_dd.quadrature_tag is not sym.QTAG_NONE):
 
-            from meshmode.connection import ChainedDiscretizationConnection
+            from meshmode.discretization.connection import \
+                    ChainedDiscretizationConnection
             intermediate_dd = sym.DOFDesc(to_dd.domain_tag)
             return ChainedDiscretizationConnection(
                     [
@@ -130,8 +158,10 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
 
         # {{{ generic to-quad
 
-        if (from_dd.domain_tag == to_dd.domain_tag
-                and from_dd.quadrature_tag != to_dd.quadrature_tag):
+        if (
+                from_dd.domain_tag == to_dd.domain_tag
+                and from_dd.quadrature_tag is sym.QTAG_NONE
+                and to_dd.quadrature_tag is not sym.QTAG_NONE):
             from meshmode.discretization.connection.same_mesh import \
                     make_same_mesh_connection
 
@@ -141,38 +171,30 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
 
         # }}}
 
+        if from_dd.quadrature_tag is not sym.QTAG_NONE:
+            raise ValueError("cannot interpolate *from* a "
+                    "(non-interpolatory) quadrature grid")
+
+        assert to_qtag is sym.QTAG_NONE
+
         if from_dd.is_volume():
             if to_dd.domain_tag is sym.FACE_RESTR_ALL:
-                return self._all_faces_volume_connection(to_qtag)
+                return self._all_faces_volume_connection()
             if to_dd.domain_tag is sym.FACE_RESTR_INTERIOR:
-                return self._interior_faces_connection(to_qtag)
+                return self._interior_faces_connection()
             elif to_dd.is_boundary():
                 assert from_dd.quadrature_tag is sym.QTAG_NONE
                 return self._boundary_connection(to_dd.domain_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))
 
-        elif from_dd.domain_tag is sym.FACE_RESTR_INTERIOR:
-            if to_dd.domain_tag is sym.FACE_RESTR_ALL and to_qtag is sym.QTAG_NONE:
-                return self._all_faces_connection(None)
-            else:
-                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 to_qtag is sym.QTAG_NONE:
-                return self._all_faces_connection(None)
-
-        elif from_dd.is_boundary():
-            if to_dd.domain_tag is sym.FACE_RESTR_ALL and to_qtag is sym.QTAG_NONE:
-                return self._all_faces_connection(from_dd.domain_tag)
-            else:
-                raise ValueError(
-                        "cannot interpolate from interior faces to: "
-                        + str(to_dd))
-
         else:
             raise ValueError("cannot interpolate from: " + str(from_dd))
 
@@ -185,11 +207,10 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
             quadrature_tag = sym.QTAG_NONE
 
         from meshmode.discretization.poly_element import \
-                PolynomialWarpAndBlendGroupFactory, \
-                QuadratureSimplexGroupFactory
+                PolynomialWarpAndBlendGroupFactory
 
         if quadrature_tag is not sym.QTAG_NONE:
-            return QuadratureSimplexGroupFactory(order=self.order)
+            return self.quad_tag_to_group_factory[quadrature_tag]
         else:
             return PolynomialWarpAndBlendGroupFactory(order=self.order)
 
@@ -210,33 +231,17 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
                         self.group_factory_for_quadrature_tag(sym.QTAG_NONE),
                         boundary_tag=boundary_tag)
 
-    @memoize_method
-    def _boundary_discr(self, boundary_tag, quadrature_tag=None):
-        if quadrature_tag is None:
-            quadrature_tag = sym.QTAG_NONE
-
-        if quadrature_tag is sym.QTAG_NONE:
-            return self._boundary_connection(boundary_tag).to_discr
-        else:
-            no_quad_bdry_discr = self.boundary_discr(boundary_tag, sym.QTAG_NONE)
-
-            from meshmode.discretization import Discretization
-            return Discretization(
-                    self._volume_discr.cl_context,
-                    no_quad_bdry_discr.mesh,
-                    self.group_factory_for_quadrature_tag(quadrature_tag))
-
     # }}}
 
     # {{{ interior faces
 
     @memoize_method
-    def _interior_faces_connection(self, quadrature_tag=None):
+    def _interior_faces_connection(self):
         from meshmode.discretization.connection import (
                 make_face_restriction, FACE_RESTR_INTERIOR)
         return make_face_restriction(
                         self._volume_discr,
-                        self.group_factory_for_quadrature_tag(quadrature_tag),
+                        self.group_factory_for_quadrature_tag(sym.QTAG_NONE),
                         FACE_RESTR_INTERIOR,
 
                         # FIXME: This will need to change as soon as we support
@@ -244,32 +249,24 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
                         # types.
                         per_face_groups=False)
 
-    def _interior_faces_discr(self, quadrature_tag=None):
-        return self._interior_faces_connection(quadrature_tag).to_discr
-
     @memoize_method
-    def opposite_face_connection(self, quadrature_tag):
-        if quadrature_tag is not sym.QTAG_NONE:
-            # FIXME
-            raise NotImplementedError("quadrature")
-
+    def opposite_face_connection(self):
         from meshmode.discretization.connection import \
                 make_opposite_face_connection
 
-        return make_opposite_face_connection(
-                self._interior_faces_connection(quadrature_tag))
+        return make_opposite_face_connection(self._interior_faces_connection())
 
     # }}}
 
     # {{{ all-faces
 
     @memoize_method
-    def _all_faces_volume_connection(self, quadrature_tag=None):
+    def _all_faces_volume_connection(self):
         from meshmode.discretization.connection import (
                 make_face_restriction, FACE_RESTR_ALL)
         return make_face_restriction(
                         self._volume_discr,
-                        self.group_factory_for_quadrature_tag(quadrature_tag),
+                        self.group_factory_for_quadrature_tag(sym.QTAG_NONE),
                         FACE_RESTR_ALL,
 
                         # FIXME: This will need to change as soon as we support
@@ -277,30 +274,6 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
                         # types.
                         per_face_groups=False)
 
-    def _all_faces_discr(self, quadrature_tag=None):
-        return self._all_faces_volume_connection(quadrature_tag).to_discr
-
-    @memoize_method
-    def _all_faces_connection(self, boundary_tag):
-        """Return a
-        :class:`meshmode.discretization.connection.DiscretizationConnection`
-        that goes from either
-        :meth:`_interior_faces_discr` (if *boundary_tag* is None)
-        or
-        :meth:`_boundary_discr` (if *boundary_tag* is not None)
-        to a discretization containing all the faces of the volume
-        discretization.
-        """
-        from meshmode.discretization.connection import \
-                make_face_to_all_faces_embedding
-
-        if boundary_tag is None:
-            faces_conn = self._interior_faces_connection()
-        else:
-            faces_conn = self._boundary_connection(boundary_tag)
-
-        return make_face_to_all_faces_embedding(faces_conn, self._all_faces_discr())
-
     # }}}
 
     @property
diff --git a/grudge/execution.py b/grudge/execution.py
index e4aebe5ca18b495acf4e608ff9f4cdfd218f57b1..2da2e508cc8ddf116eb4ffb30ecfd3f8835f957b 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -1,6 +1,6 @@
 from __future__ import division, absolute_import
 
-__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
+__copyright__ = "Copyright (C) 2015-2017 Andreas Kloeckner, Bogdan Enache"
 
 __license__ = """
 Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -198,36 +198,37 @@ class ExecutionMapper(mappers.Evaluator,
             knl = lp.make_kernel(
                 """{[k,i,j]:
                     0<=k<nelements and
-                    0<=i,j<ndiscr_nodes}""",
+                    0<=i<ndiscr_nodes_out and
+                    0<=j<ndiscr_nodes_in}""",
                 "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)
+        in_discr = self.discrwb.discr_from_dd(op.dd_in)
+        out_discr = self.discrwb.discr_from_dd(op.dd_out)
 
-        # FIXME: This shouldn't really assume that it's dealing with a volume
-        # input. What about quadrature? What about boundaries?
-        result = discr.empty(
+        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
+        for in_grp, out_grp in zip(in_discr.groups, out_discr.groups):
+
+            cache_key = "elwise_linear", in_grp, out_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))
-                        .with_queue(None))
+                    cl.array.to_device(
+                        self.queue,
+                        np.asarray(op.matrix(out_grp, in_grp), dtype=field.dtype))
+                    .with_queue(None))
 
                 self.bound_op.operator_data_cache[cache_key] = matrix
 
-            knl()(self.queue, mat=matrix, result=grp.view(result),
-                    vec=grp.view(field))
+            knl()(self.queue, mat=matrix, result=out_grp.view(result),
+                    vec=in_grp.view(field))
 
         return result
 
@@ -246,14 +247,7 @@ class ExecutionMapper(mappers.Evaluator,
         return conn(self.queue, self.rec(field_expr)).with_queue(self.queue)
 
     def map_opposite_interior_face_swap(self, op, field_expr):
-        dd = op.dd_in
-
-        qtag = dd.quadrature_tag
-        if qtag is None:
-            # FIXME: Remove once proper quadrature support arrives
-            qtag = sym.QTAG_NONE
-
-        return self.discrwb.opposite_face_connection(qtag)(
+        return self.discrwb.opposite_face_connection()(
                 self.queue, self.rec(field_expr)).with_queue(self.queue)
 
     # {{{ face mass operator
@@ -289,8 +283,7 @@ class ExecutionMapper(mappers.Evaluator,
 
         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):
+        for afgrp, volgrp in zip(all_faces_discr.groups, vol_discr.groups):
             cache_key = "face_mass", afgrp, op, field.dtype
 
             nfaces = volgrp.mesh_el_group.nfaces
@@ -365,10 +358,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
         assert repr_op.dd_in.domain_tag == repr_op.dd_out.domain_tag
 
         @memoize_in(self.discrwb, "reference_derivative_knl")
@@ -377,7 +366,8 @@ class ExecutionMapper(mappers.Evaluator,
                 """{[imatrix,k,i,j]:
                     0<=imatrix<nmatrices and
                     0<=k<nelements and
-                    0<=i,j<nunit_nodes}""",
+                    0<=i<nunit_nodes_out and
+                    0<=j<nunit_nodes_in}""",
                 """
                 result[imatrix, k, i] = sum(
                         j, diff_mat[imatrix, i, j] * vec[k, j])
@@ -388,26 +378,31 @@ 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 in_grp, out_grp in zip(in_discr.groups, out_discr.groups):
+
+            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)], []
 
@@ -495,6 +490,10 @@ def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None,
     dumper("before-cfold", sym_operator)
     sym_operator = mappers.CommutativeConstantFoldingMapper()(sym_operator)
 
+    dumper("before-qcheck", sym_operator)
+    sym_operator = mappers.QuadratureCheckerAndRemover(
+            discrwb.quad_tag_to_group_factory)(sym_operator)
+
     # Work around https://github.com/numpy/numpy/issues/9438
     #
     # The idea is that we need 1j as an expression to survive
@@ -510,35 +509,9 @@ def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None,
             complex_type=discrwb.complex_dtype.type,
             )(sym_operator)
 
-    # Ordering restriction:
-    #
-    # - Must run constant fold before first type inference pass, because zeros,
-    # while allowed, violate typing constraints (because they can't be assigned
-    # a unique type), and need to be killed before the type inferrer sees them.
-
-    # FIXME: Reenable type inference
-
-    # from grudge.symbolic.mappers.type_inference import TypeInferrer
-    # dumper("before-specializer", sym_operator)
-    # sym_operator = mappers.OperatorSpecializer(
-    #         TypeInferrer()(sym_operator)
-    #         )(sym_operator)
-
-    # Ordering restriction:
-    #
-    # - Must run OperatorSpecializer before performing the GlobalToReferenceMapper,
-    # because otherwise it won't differentiate the type of grids (node or quadrature
-    # grids) that the operators will apply on.
-
     dumper("before-global-to-reference", sym_operator)
     sym_operator = mappers.GlobalToReferenceMapper(discrwb.ambient_dim)(sym_operator)
 
-    # Ordering restriction:
-    #
-    # - Must specialize quadrature operators before performing inverse mass
-    # contraction, because there are no inverse-mass-contracted variants of the
-    # quadrature operators.
-
     dumper("before-imass", sym_operator)
     sym_operator = mappers.InverseMassContractor()(sym_operator)
 
@@ -564,19 +537,20 @@ def bind(discr, sym_operator, post_bind_mapper=lambda x: x,
 
     stage = [0]
 
-    def dump_optemplate(name, sym_operator):
-        if "dump_optemplate_stages" in debug_flags:
-            from grudge.tools import open_unique_debug_file
-            from grudge.optemplate import pretty
-            open_unique_debug_file("%02d-%s" % (stage[0], name), ".txt").write(
-                    pretty(sym_operator))
+    def dump_sym_operator(name, sym_operator):
+        if "dump_sym_operator_stages" in debug_flags:
+            from pytools.debug import open_unique_debug_file
+            outf, name = open_unique_debug_file("%02d-%s" % (stage[0], name), ".txt")
+            with outf:
+                outf.write(sym.pretty(sym_operator))
+
             stage[0] += 1
 
     sym_operator = process_sym_operator(
             discr,
             sym_operator,
             post_bind_mapper=post_bind_mapper,
-            dumper=dump_optemplate)
+            dumper=dump_sym_operator)
 
     from grudge.symbolic.compiler import OperatorCompiler
     discr_code, eval_code = OperatorCompiler(discr)(sym_operator)
diff --git a/grudge/models/advection.py b/grudge/models/advection.py
index 1c83ed9675d5a45c75b22e3b4d5ca7926acf9133..091b9870c4b82a3fde2fc21b84bf92e5ec62ab17 100644
--- a/grudge/models/advection.py
+++ b/grudge/models/advection.py
@@ -3,7 +3,7 @@
 
 from __future__ import division, absolute_import
 
-__copyright__ = "Copyright (C) 2009 Andreas Kloeckner"
+__copyright__ = "Copyright (C) 2009-2017 Andreas Kloeckner, Bogdan Enache"
 
 __license__ = """
 Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -132,13 +132,16 @@ class WeakAdvectionOperator(AdvectionOperatorBase):
 # {{{ variable-coefficient advection
 
 class VariableCoefficientAdvectionOperator(HyperbolicOperator):
-    def __init__(self, dim, v, inflow_u, flux_type="central"):
+    def __init__(self, dim, v, inflow_u, flux_type="central", quad_tag="product"):
         self.ambient_dim = dim
         self.v = v
         self.inflow_u = inflow_u
         self.flux_type = flux_type
 
+        self.quad_tag = quad_tag
+
     def flux(self, u):
+
         normal = sym.normal(u. dd, self.ambient_dim)
 
         surf_v = sym.interp("vol", u.dd)(self.v)
@@ -167,15 +170,26 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator):
         # boundary conditions -------------------------------------------------
         bc_in = self.inflow_u
 
+        all_faces_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag)
+        boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag)
+
         def flux(pair):
-            return sym.interp(pair.dd, "all_faces")(
+            return sym.interp(pair.dd, all_faces_dd)(
                     self.flux(pair))
 
+        quad_dd = sym.DOFDesc("vol", self.quad_tag)
+
+        to_quad = sym.interp("vol", quad_dd)
+
+        stiff_t = sym.stiffness_t(self.ambient_dim, quad_dd, "vol")
+        quad_v = to_quad(self.v)
+        quad_u = to_quad(u)
+
         return sym.InverseMassOperator()(
-                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))
+                (stiff_t[0](quad_u * quad_v[0]) + stiff_t[1](quad_u * quad_v[1]))
+                - sym.FaceMassOperator(all_faces_dd, "vol")(
+                    flux(sym.int_tpair(u, self.quad_tag))
+                    + flux(sym.bv_tpair(boundary_dd, 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/compiler.py b/grudge/symbolic/compiler.py
index f6e238083266d4ba0c510243b116b2e49b024631..c555cea0a56e947b4a94e0cd1736e8fa2b97ce63 100644
--- a/grudge/symbolic/compiler.py
+++ b/grudge/symbolic/compiler.py
@@ -338,7 +338,7 @@ class Code(object):
         self.static_schedule_attempts = 5
 
     def dump_dataflow_graph(self):
-        from grudge.tools import open_unique_debug_file
+        from pytools.debug import open_unique_debug_file
 
         open_unique_debug_file("dataflow", ".dot")\
                 .write(dot_dataflow_graph(self, max_node_label_length=None))
@@ -792,7 +792,7 @@ def aggregate_assignments(inf_mapper, instructions, result,
 # }}}
 
 
-# {{{
+# {{{ to-loopy mapper
 
 def set_once(d, k, v):
     try:
diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py
index 405cb5f06706116266322c3d9ab0ab1c206a6c22..ca11cbf0d92c96b6c9d00a6e035b7b840c2a85cd 100644
--- a/grudge/symbolic/mappers/__init__.py
+++ b/grudge/symbolic/mappers/__init__.py
@@ -717,96 +717,67 @@ class NoCSEStringifyMapper(StringifyMapper):
 
 # {{{ quadrature support
 
-class QuadratureUpsamplerRemover(CSECachingMapperMixin, IdentityMapper):
-    def __init__(self, quad_min_degrees, do_warn=True):
+class QuadratureCheckerAndRemover(CSECachingMapperMixin, IdentityMapper):
+    """Checks whether all quadratu
+    """
+
+    def __init__(self, quad_tag_to_group_factory):
         IdentityMapper.__init__(self)
         CSECachingMapperMixin.__init__(self)
-        self.quad_min_degrees = quad_min_degrees
-        self.do_warn = do_warn
+        self.quad_tag_to_group_factory = quad_tag_to_group_factory
 
     map_common_subexpression_uncached = \
             IdentityMapper.map_common_subexpression
 
-    def map_operator_binding(self, expr):
-        if isinstance(expr.op, (op.QuadratureGridUpsampler,
-                op.QuadratureInteriorFacesGridUpsampler,
-                op.QuadratureBoundaryGridUpsampler)):
-            try:
-                min_degree = self.quad_min_degrees[expr.op.quadrature_tag]
-            except KeyError:
-                if self.do_warn:
-                    from warnings import warn
-                    warn("No minimum degree for quadrature tag '%s' specified--"
-                            "falling back to nodal evaluation"
-                            % expr.op.quadrature_tag)
-                return self.rec(expr.field)
-            else:
-                if min_degree is None:
-                    return self.rec(expr.field)
-                else:
-                    return expr.op(self.rec(expr.field))
-        else:
-            return IdentityMapper.map_operator_binding(self, expr)
+    def _process_dd(self, dd, location_descr):
+        from grudge.symbolic.primitives import DOFDesc, QTAG_NONE
+        if dd.quadrature_tag is not QTAG_NONE:
+            if dd.quadrature_tag not in self.quad_tag_to_group_factory:
+                raise ValueError("found unknown quadrature tag '%s' in '%s'"
+                        % (dd.quadrature_tag, location_descr))
 
+            grp_factory = self.quad_tag_to_group_factory[dd.quadrature_tag]
+            if grp_factory is None:
+                dd = DOFDesc(dd.domain_tag, QTAG_NONE)
 
-class QuadratureDetector(CSECachingMapperMixin, CombineMapper):
-    """For a given expression, this mapper returns the upsampling
-    operator in effect at the root of the expression, or *None*
-    if there isn't one.
-    """
-    class QuadStatusNotKnown:
-        pass
-
-    map_common_subexpression_uncached = \
-            CombineMapper.map_common_subexpression
+        return dd
 
-    def combine(self, values):
-        from pytools import single_valued
-        return single_valued([
-            v for v in values if v is not self.QuadStatusNotKnown])
+    def map_operator_binding(self, expr):
+        dd_in_orig = dd_in = expr.op.dd_in
+        dd_out_orig = dd_out = expr.op.dd_out
 
-    def map_variable(self, expr):
-        return None
+        dd_in = self._process_dd(dd_in, "dd_in of %s" % type(expr.op).__name__)
+        dd_out = self._process_dd(dd_out, "dd_out of %s" % type(expr.op).__name__)
 
-    def map_constant(self, expr):
-        return self.QuadStatusNotKnown
+        if dd_in_orig == dd_in and dd_out_orig == dd_out:
+            # unchanged
+            return IdentityMapper.map_operator_binding(self, expr)
 
-    def map_operator_binding(self, expr):
-        if isinstance(expr.op, (
-                op.DiffOperatorBase, op.FluxOperatorBase,
-                op.MassOperatorBase)):
-            return None
-        elif isinstance(expr.op, (op.QuadratureGridUpsampler,
-                op.QuadratureInteriorFacesGridUpsampler)):
-            return expr.op
-        else:
-            return CombineMapper.map_operator_binding(self, expr)
+        import grudge.symbolic.operators as op
+        # changed
 
+        if dd_in == dd_out and isinstance(expr.op, op.InterpolationOperator):
+            # This used to be to-quad interpolation and has become a no-op.
+            # Remove it.
+            return self.rec(expr.field)
 
-class QuadratureUpsamplerChanger(CSECachingMapperMixin, IdentityMapper):
-    """This mapper descends the expression tree, down to each
-    quadrature-consuming operator (diff, mass) along each branch.
-    It then change
-    """
-    def __init__(self, desired_quad_op):
-        IdentityMapper.__init__(self)
-        CSECachingMapperMixin.__init__(self)
+        if isinstance(expr.op, op.StiffnessTOperator):
+            new_op = type(expr.op)(expr.op.xyz_axis, dd_in, dd_out)
+        elif isinstance(expr.op, (op.FaceMassOperator, op.InterpolationOperator)):
+            new_op = type(expr.op)(dd_in, dd_out)
+        else:
+            raise NotImplementedError("do not know how to modify dd_in and dd_out "
+                    "in %s" % type(expr.op).__name__)
 
-        self.desired_quad_op = desired_quad_op
+        return new_op(self.rec(expr.field))
 
-    map_common_subexpression_uncached = \
-            IdentityMapper.map_common_subexpression
+    def map_ones(self, expr):
+        dd = self._process_dd(expr.dd, location_descr=type(expr).__name__)
+        return type(expr)(dd)
 
-    def map_operator_binding(self, expr):
-        if isinstance(expr.op, (
-                op.DiffOperatorBase, op.FluxOperatorBase,
-                op.MassOperatorBase)):
-            return expr
-        elif isinstance(expr.op, (op.QuadratureGridUpsampler,
-                op.QuadratureInteriorFacesGridUpsampler)):
-            return self.desired_quad_op(expr.field)
-        else:
-            return IdentityMapper.map_operator_binding(self, expr)
+    def map_node_coordinate_component(self, expr):
+        dd = self._process_dd(expr.dd, location_descr=type(expr).__name__)
+        return type(expr)(expr.axis, dd)
 
 # }}}
 
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index 84dd4679bad9ba642ba00ffc4fcf391df4185dde..194b071614dba86a7ff150c075e5a3c987f23a4b 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -2,7 +2,7 @@
 
 from __future__ import division, absolute_import
 
-__copyright__ = "Copyright (C) 2008 Andreas Kloeckner"
+__copyright__ = "Copyright (C) 2008-2017 Andreas Kloeckner, Bogdan Enache"
 
 __license__ = """
 Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -97,6 +97,14 @@ class ElementwiseLinearOperator(Operator):
 
 
 class InterpolationOperator(Operator):
+    def __init__(self, dd_in, dd_out):
+        official_dd_in = _sym().as_dofdesc(dd_in)
+        official_dd_out = _sym().as_dofdesc(dd_out)
+        if official_dd_in == official_dd_out:
+            raise ValueError("Interpolating from {} to {}"
+            " does not do anything.".format(official_dd_in, official_dd_out))
+
+        super(InterpolationOperator, self).__init__(dd_in, dd_out)
     mapper_method = intern("map_interpolation")
 
 
@@ -141,8 +149,12 @@ class DiffOperatorBase(Operator):
     def __init__(self, xyz_axis, dd_in=None, dd_out=None):
         if dd_in is None:
             dd_in = _sym().DD_VOLUME
+
         if dd_out is None:
             dd_out = dd_in.with_qtag(_sym().QTAG_NONE)
+        else:
+            dd_out = _sym().as_dofdesc(dd_out)
+
         if dd_out.uses_quadrature():
             raise ValueError("differentiation outputs are not on "
                     "quadrature grids")
@@ -215,18 +227,30 @@ 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_elem_grp, in_elem_grp):
+        if in_elem_grp == out_elem_grp:
+            assert in_elem_grp.is_orthogonal_basis()
+            mmat = in_elem_grp.mass_matrix()
+            return [dmat.T.dot(mmat.T) for dmat in in_elem_grp.diff_matrices()]
+
+        from modepy import vandermonde
+        vand = vandermonde(out_elem_grp.basis(), out_elem_grp.unit_nodes)
+        grad_vand = vandermonde(out_elem_grp.grad_basis(), in_elem_grp.unit_nodes)
+        vand_inv_t = np.linalg.inv(vand).T
+
+        weights = in_elem_grp.weights
+
+        return np.einsum('c,bz,acz->abc', weights, vand_inv_t, grad_vand)
+
 
 # }}}
 
@@ -322,8 +346,6 @@ class VandermondeOperator(ElementwiseLinearOperator):
 
 # }}}
 
-# }}}
-
 
 # {{{ mass operators
 
@@ -338,9 +360,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 +377,30 @@ class RefMassOperatorBase(ElementwiseLinearOperator):
 
 class RefMassOperator(RefMassOperatorBase):
     @staticmethod
-    def matrix(element_group):
-        return element_group.mass_matrix()
+    def matrix(out_element_group, in_element_group):
+        if out_element_group == in_element_group:
+            return in_element_group.mass_matrix()
+
+        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
+
+        weights = in_element_group.weights
+
+        return np.einsum('c,bz,acz->abc', weights, vand_inv_t, o_vand)
 
     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 +435,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 +506,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..00f7b8d386af6130810cf900bbaaebb0027f5adc 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -2,7 +2,7 @@
 
 from __future__ import division, absolute_import
 
-__copyright__ = "Copyright (C) 2008 Andreas Kloeckner"
+__copyright__ = "Copyright (C) 2008-2017 Andreas Kloeckner, Bogdan Enache"
 
 __license__ = """
 Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -653,10 +653,18 @@ 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):
+    i = _sym().interp("vol", "int_faces")(expression)
     e = cse(_sym().OppositeInteriorFaceSwap()(i))
-    return TracePair("int_faces", i, e)
+
+    if qtag is not None and qtag != _sym().QTAG_NONE:
+        q_dd = _sym().DOFDesc("int_faces", qtag)
+        i = cse(_sym().interp("int_faces", q_dd)(i))
+        e = cse(_sym().interp("int_faces", q_dd)(e))
+    else:
+        q_dd = "int_faces"
+
+    return TracePair(q_dd, i, e)
 
 
 def bdry_tpair(dd, interior, exterior):
diff --git a/test/test_grudge.py b/test/test_grudge.py
index 84c46168ede26efe1a42935501d69fba99b7700a..34532bfc80371e29bde28ea9674c48e36ab889aa 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -334,8 +334,7 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type
 
 
 @pytest.mark.parametrize("order", [3, 4, 5])
-# test: 'test_convergence_advec(cl._csc, "disk", [0.1, 0.05], "strong", "upwind", 3)'
-def test_convergence_maxwell(ctx_factory,  order, visualize=False):
+def test_convergence_maxwell(ctx_factory,  order):
     """Test whether 3D maxwells actually converges"""
 
     cl_ctx = cl.create_some_context()
@@ -403,6 +402,73 @@ def test_convergence_maxwell(ctx_factory,  order, visualize=False):
     assert eoc_rec.order_estimate() > order
 
 
+@pytest.mark.parametrize("order", [2, 3, 4])
+def test_improvement_quadrature(ctx_factory, order):
+    """Test whether quadrature improves things and converges"""
+    from meshmode.mesh.generation import generate_regular_rect_mesh
+    from grudge.models.advection import VariableCoefficientAdvectionOperator
+    from pytools.convergence import EOCRecorder
+    from pytools.obj_array import join_fields
+    from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory
+
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+
+    dims = 2
+    sym_nds = sym.nodes(dims)
+    advec_v = join_fields(-1*sym_nds[1], sym_nds[0])
+
+    flux = "upwind"
+    op = VariableCoefficientAdvectionOperator(2, advec_v, 0, flux_type=flux)
+
+    def gaussian_mode():
+        source_width = 0.1
+        sym_x = sym.nodes(2)
+        return sym.exp(-np.dot(sym_x, sym_x) / source_width**2)
+
+    def conv_test(descr, use_quad):
+        print("-"*75)
+        print(descr)
+        print("-"*75)
+        eoc_rec = EOCRecorder()
+
+        ns = [20, 25]
+        for n in ns:
+            mesh = generate_regular_rect_mesh(
+                a=(-0.5,)*dims,
+                b=(0.5,)*dims,
+                n=(n,)*dims,
+                order=order)
+
+            if use_quad:
+                quad_tag_to_group_factory = {
+                    "product": QuadratureSimplexGroupFactory(order=4*order)
+                    }
+            else:
+                quad_tag_to_group_factory = {"product": None}
+
+            discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order,
+                    quad_tag_to_group_factory=quad_tag_to_group_factory)
+
+            bound_op = bind(discr, op.sym_operator())
+            fields = bind(discr, gaussian_mode())(queue, t=0)
+            norm = bind(discr, sym.norm(2, sym.var("u")))
+
+            esc = bound_op(queue, u=fields)
+            total_error = norm(queue, u=esc)
+            eoc_rec.add_data_point(1.0/n, total_error)
+
+        print(eoc_rec.pretty_print(abscissa_label="h", error_label="LInf Error"))
+
+        return eoc_rec.order_estimate(), np.array([x[1] for x in eoc_rec.history])
+
+    eoc, errs = conv_test("no quadrature", False)
+    q_eoc, q_errs = conv_test("with quadrature", True)
+    assert q_eoc > eoc
+    assert (q_errs < errs).all()
+    assert q_eoc > order
+
+
 def test_foreign_points(ctx_factory):
     pytest.importorskip("sumpy")
     import sumpy.point_calculus as pc