diff --git a/contrib/maxima/maxwellrectcav3d.mac b/contrib/maxima/maxwellrectcav3d.mac
new file mode 100644
index 0000000000000000000000000000000000000000..abda9687d244b25b3cf546bea9f83383e222bc34
--- /dev/null
+++ b/contrib/maxima/maxwellrectcav3d.mac
@@ -0,0 +1,73 @@
+kill(all);
+load("eigen");
+load("itensor");
+load("diag");
+
+/* -------------------------------------------------------------------------- */
+/* Utilities */
+/* -------------------------------------------------------------------------- */
+mu: 1;
+epsilon: 1;
+
+coords:[x,y,z];
+
+curl(x):=crossfunc(lambda([i,j], diff(x[j], coords[i])));
+div(f):=sum(diff(f[j], coords[j]), j, 1, length(coords));
+
+faraday(max_E, max_H):= curl(max_E) - %i * omega * max_H;
+ampere(max_E, max_H):= curl(max_H) + %i * omega * mu * epsilon * max_E;
+/*
+ampere(max_E, max_H):= curl(max_H) + %i * omega * max_E;
+*/
+
+
+crossfunc(f):=makelist(
+  sum(sum(
+      levi_civita([i,j,k])*f(j,k),
+   j,1,3),k,1,3),i,1,3)$
+
+crossprod(a,b):=crossfunc(lambda([j,k], a[j]*b[k]));
+
+/* -------------------------------------------------------------------------- */
+/* Attempts */
+/* -------------------------------------------------------------------------- */
+nabla_t_squared(f):=diff(f, x, 2) + diff(f, y, 2);
+nabla_t(f):= [diff(f, x, 1) , diff(f, y, 1), 0];
+
+/*
+gamma : sqrt((k_y^2 + k_x^2)/ (mu *epsilon)) ;
+*/
+omega : sqrt(k_x^2 + k_y^2 + k_z^2);
+gamma : sqrt((k_y^2 + k_x^2)) ;
+
+psi_cand:E_0*sin(k_x*x)*sin(k_y*y);
+
+
+
+wave_eqn(f, gam_s):=nabla_t_squared(f) + mu*epsilon*gam_s*f;
+gam_s : gamma^2;
+
+/* The _t indicates transverse components (i.e. x and y components only) */
+/*
+E_t(psi):=(-k_z/gam_s)*sin(k_z * z)* nabla_t(psi);
+H_t(psi):=crossprod((%i * omega /gam_s)*cos(k_z * z)* [0,0,1], nabla_t(psi));
+*/
+E_t(psi):=(-k_z/gam_s)*sin(k_z * z)* nabla_t(psi);
+H_t(psi):=crossprod((%i * omega * epsilon/gam_s)*cos(k_z * z)* [0,0,1], nabla_t(psi));
+
+
+/* These are used as the analytic solution for a rectangular cavity resonator 
+    with travsverse magnetic waves */
+
+E : E_t(psi_cand) + [0,0,psi_cand * cos(k_z * z)];
+H :H_t(psi_cand);
+
+Etime : E * %e ^ (- %i * omega * t);
+Htime : H * %e ^ (- %i * omega * t);
+
+
+trigrat(div(E));
+trigrat(div(H));
+trigrat(faraday(E,H));
+trigrat(ampere(E,H));
+
diff --git a/examples/maxwell/analytic_solutions.py b/examples/maxwell/analytic_solutions.py
deleted file mode 100644
index 7ff8884f1d19e46d320e0cc7a8c5301f74d14551..0000000000000000000000000000000000000000
--- a/examples/maxwell/analytic_solutions.py
+++ /dev/null
@@ -1,460 +0,0 @@
-# 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
-from __future__ import absolute_import
-from __future__ import print_function
-from grudge.tools import \
-        cyl_bessel_j, \
-        cyl_bessel_j_prime
-from math import sqrt, pi, sin, cos, atan2, acos
-import numpy
-import numpy.linalg as la
-import cmath
-from six.moves import range
-from six.moves import zip
-
-
-
-
-# solution adapters -----------------------------------------------------------
-class RealPartAdapter:
-    def __init__(self, adaptee):
-        self.adaptee = adaptee
-
-    @property
-    def shape(self):
-        return self.adaptee.shape
-
-    def __call__(self, x, el):
-        return [xi.real for xi in self.adaptee(x, el)]
-
-class SplitComplexAdapter:
-    def __init__(self, adaptee):
-        self.adaptee = adaptee
-
-    @property
-    def shape(self):
-        (n,) = self.adaptee.shape
-        return (n*2,)
-
-    def __call__(self, x, el):
-        ad_x = self.adaptee(x, el)
-        return [xi.real for xi in ad_x] + [xi.imag for xi in ad_x]
-
-
-
-
-class CylindricalFieldAdapter:
-    """Turns a field returned as (r, phi, z) components into
-    (x,y,z) components.
-    """
-    def __init__(self, adaptee):
-        self.adaptee = adaptee
-
-    @property
-    def shape(self):
-        return self.adaptee.shape
-
-    def __call__(self, x, el):
-        xy = x[:2]
-        r = la.norm(xy)
-        phi = atan2(x[1], x[0])
-
-        prev_result = self.adaptee(x, el)
-        result = []
-        i = 0
-        while i < len(prev_result):
-            fr, fphi, fz = prev_result[i:i+3]
-            result.extend([
-                    cos(phi)*fr - sin(phi)*fphi, # ex
-                    sin(phi)*fr + cos(phi)*fphi, # ey
-                    fz,
-                    ])
-            i += 3
-
-        return result
-
-
-
-
-class SphericalFieldAdapter:
-    """Turns the adaptee's field C{(Er, Etheta, Ephi)} components into
-    C{(Ex,Ey,Ez)} components.
-
-    Conventions are those of
-    http://mathworld.wolfram.com/SphericalCoordinates.html
-    """
-    def __init__(self, adaptee):
-        self.adaptee = adaptee
-
-    @property
-    def shape(self):
-        return self.adaptee.shape
-
-    @staticmethod
-    def r_theta_phi_from_xyz(x):
-        r = la.norm(x)
-        return r, atan2(x[1], x[0]), acos(x[2]/r)
-
-    def __call__(self, x, el):
-        r, theta, phi = self.r_theta_phi_from_xyz(x)
-
-        er = numpy.array([
-            cos(theta)*sin(phi),
-            sin(theta)*sin(phi),
-            cos(phi)])
-        etheta = numpy.array([
-            -sin(theta),
-            cos(theta),
-            0])
-        ephi = numpy.array([
-            cos(theta)*cos(phi),
-            sin(theta)*cos(phi),
-            -sin(phi)])
-
-        adaptee_result = self.adaptee(x, el)
-        result = numpy.empty_like(adaptee_result)
-        i = 0
-        while i < len(adaptee_result):
-            fr, ftheta, fphi = adaptee_result[i:i+3]
-            result[i:i+3] = fr*er + ftheta*etheta + fphi*ephi
-            i += 3
-
-        return result
-
-
-
-
-# actual solutions ------------------------------------------------------------
-class CylindricalCavityMode:
-    """A cylindrical TM cavity mode.
-
-    Taken from:
-    J.D. Jackson, Classical Electrodynamics, Wiley.
-    3rd edition, 2001.
-    ch 8.7, p. 368f.
-    """
-    
-    def __init__(self,  m, n, p, radius, height, epsilon, mu):
-        try:
-            from bessel_zeros import bessel_zeros
-        except ImportError:
-            print("*** You need to generate the bessel root data file.")
-            print("*** Execute generate-bessel-zeros.py at the command line.")
-            raise
-
-        assert m >= 0 and m == int(m)
-        assert n >= 1 and n == int(n)
-        assert p >= 0 and p == int(p)
-        self.m = m
-        self.n = n
-        self.p = p
-        self.phi_sign = 1
-
-        R = self.radius = radius
-        d = self.height = height
-
-        self.epsilon = epsilon
-        self.mu = mu
-
-        self.t = 0
-
-        x_mn = bessel_zeros[m][n-1]
-
-        self.omega = 1 / sqrt(mu*epsilon) * sqrt(
-                x_mn**2 / R**2
-                + p**2 * pi**2 / d**2)
-
-        self.gamma_mn = x_mn/R
-
-    def set_time(self, t):
-        self.t = t
-
-    shape = (6,)
-
-    def __call__(self, x, el):
-        # coordinates -----------------------------------------------------
-        xy = x[:2]
-        r = sqrt(xy*xy)
-        phi = atan2(x[1], x[0])
-        z = x[2]
-
-        # copy instance variables for easier access -----------------------
-        m = self.m
-        p = self.p
-        gamma = self.gamma_mn
-        phi_sign = self.phi_sign
-        omega = self.omega
-        d = self.height
-        epsilon = self.epsilon
-
-        # common subexpressions -------------------------------------------
-        tdep = cmath.exp(-1j * omega * self.t)
-        phi_factor = cmath.exp(phi_sign * 1j * m * phi)
-
-        # psi and derivatives ---------------------------------------------
-        psi = cyl_bessel_j(m, gamma * r) * phi_factor
-        psi_dr = gamma*cyl_bessel_j_prime(m, gamma*r) * phi_factor
-        psi_dphi = (cyl_bessel_j(m, gamma * r) 
-                * 1/r * phi_sign*1j*m * phi_factor)
-
-        # field components in polar coordinates ---------------------------
-        ez   = tdep * cos(p * pi * z / d) * psi
-
-        e_transverse_factor = (tdep
-                * (-p*pi/(d*gamma**2))
-                * sin(p * pi * z / d))
-
-        er   = e_transverse_factor * psi_dr
-        ephi = e_transverse_factor * psi_dphi
-
-        hz   = 0j
-
-        # z x grad psi = z x (psi_x, psi_y)   = (-psi_y,   psi_x)
-        # z x grad psi = z x (psi_r, psi_phi) = (-psi_phi, psi_r)
-        h_transverse_factor = (tdep
-                * 1j*epsilon*omega/gamma**2
-                * cos(p * pi * z / d))
-
-        hr   = h_transverse_factor * (-psi_dphi)
-        hphi = h_transverse_factor * psi_dr
-
-        return [er, ephi, ez, hr, hphi, hz]
-
-
-
-
-class RectangularWaveguideMode:
-    """A rectangular TM cavity mode."""
-    
-    def __init__(self, epsilon, mu, mode_indices, 
-            dimensions=(1,1,1), coefficients=(1,0,0),
-            forward_coeff=1, backward_coeff=0):
-        for n in mode_indices:
-            assert n >= 0 and n == int(n)
-        self.mode_indices = mode_indices
-        self.dimensions = dimensions
-        self.coefficients = coefficients
-        self.forward_coeff = forward_coeff
-        self.backward_coeff = backward_coeff
-
-        self.epsilon = epsilon
-        self.mu = mu
-
-        self.t = 0
-
-        self.factors = [n*pi/a for n,  a in zip(self.mode_indices, self.dimensions)]
-
-        c = 1/sqrt(mu*epsilon)
-        self.k = sqrt(sum(f**2 for f in self.factors))
-        self.omega = self.k*c
-
-    def set_time(self, t):
-        self.t = t
-
-    def __call__(self, discr):
-        f,g,h = self.factors
-        omega = self.omega
-        k = self.k
-
-        x = discr.nodes[:,0]
-        y = discr.nodes[:,1]
-        if discr.dimensions > 2:
-            z = discr.nodes[:,2]
-        else:
-            z = discr.volume_zeros()
-
-        sx = numpy.sin(f*x)
-        sy = numpy.sin(g*y)
-        cx = numpy.cos(f*x)
-        cy = numpy.cos(g*y)
-
-        zdep_add = numpy.exp(1j*h*z)+numpy.exp(-1j*h*z)
-        zdep_sub = numpy.exp(1j*h*z)-numpy.exp(-1j*h*z)
-
-        tdep = numpy.exp(-1j * omega * self.t)
-
-        C = 1j/(f**2+g**2)
-
-        result = numpy.zeros((6,len(x)), dtype=zdep_add.dtype)
-
-        result[0] = C*f*h*cx*sy*zdep_sub*tdep # ex
-        result[1] = C*g*h*sx*cy*zdep_sub*tdep # ey
-        result[2] = sx*sy*zdep_add*tdep # ez
-        result[3] = -C*g*self.epsilon*omega*sx*cy*zdep_add*tdep # hx
-        result[4] = C*f*self.epsilon*omega*cx*sy*zdep_add*tdep
-
-        return result
-
-
-
-
-class RectangularCavityMode(RectangularWaveguideMode):
-    """A rectangular TM cavity mode."""
-    
-    def __init__(self, *args, **kwargs):
-        if "scale" in kwargs:
-            kwargs["forward_coeff"] = scale
-            kwargs["backward_coeff"] = scale
-        else:
-            kwargs["forward_coeff"] = 1
-            kwargs["backward_coeff"] = 1
-        RectangularWaveguideMode.__init__(self, *args, **kwargs)
-
-
-
-
-
-# dipole solution -------------------------------------------------------------
-class DipoleFarField:
-    """Dipole EM field.
-
-    See U{http://piki.tiker.net/wiki/Dipole}.
-    """
-    def __init__(self, q, d, omega, epsilon, mu):
-        self.q = q
-        self.d = d
-
-        self.epsilon = epsilon
-        self.mu = mu
-        self.c = c = 1/sqrt(mu*epsilon)
-
-        self.omega = omega
-
-        freq = omega/(2*pi)
-        self.wavelength = c/freq
-
-    def far_field_mask(self, x, el):
-        if la.norm(x) > 0.5*self.wavelength:
-            return 1
-        else:
-            return 0
-
-    def set_time(self, t):
-        self.t = t
-
-    shape = (6,)
-
-    def __call__(self, x, el):
-        # coordinates ---------------------------------------------------------
-        r, theta, phi = SphericalFieldAdapter.r_theta_phi_from_xyz(x)
-
-        # copy instance variables for easier access ---------------------------
-        q = self.q
-        d = self.d
-
-        epsilon = self.epsilon
-        mu = self.mu
-        c = self.c
-
-        omega = self.omega
-        t = self.t
-
-        # field components ----------------------------------------------------
-        tdep = omega*t-omega*r/c
-        e_r  = (2*cos(phi)*q*d) / (4*pi*epsilon) * (
-                1/r**3 * sin(tdep)
-                +omega/(c*r**2) * cos(tdep))
-        e_phi = (sin(phi)*q*d) / (4*pi*epsilon) * (
-                (1/r**3 - omega**2/(c**2*r)*sin(tdep))
-                + omega/(c*r**2)*cos(tdep))
-        h_theta = (omega*sin(phi)*q*d)/(4*pi) * (
-                - omega/(c*r)*sin(tdep)
-                + 1/r**2*cos(tdep))
-
-        return [e_r, 0, e_phi, 0, h_theta, 0]
-
-    def source_modulation(self, t):
-        j0 = self.q*self.d*self.omega/self.epsilon
-        return j0*cos(self.omega*t)
-
-
-
-
-# analytic solution tools -----------------------------------------------------
-def check_time_harmonic_solution(discr, mode, c_sol):
-    from grudge.discretization import bind_nabla, bind_mass_matrix
-    from grudge.visualization import SiloVisualizer
-    from grudge.silo import SiloFile
-    from grudge.tools import dot, cross
-    from grudge.silo import DB_VARTYPE_VECTOR
-
-    def curl(field):
-        return cross(nabla, field)
-
-    vis = SiloVisualizer(discr)
-
-    nabla = bind_nabla(discr)
-    mass = bind_mass_matrix(discr)
-
-    def rel_l2_error(err, base):
-        def l2_norm(field):
-            return sqrt(dot(field, mass*field))
-
-        base_l2 = l2_norm(base)
-        err_l2 = l2_norm(err)
-        if base_l2 == 0:
-            if err_l2 == 0:
-                return 0.
-            else:
-                return float("inf")
-        else:
-            return err_l2/base_l2
-
-    dt = 0.1
-
-    for step in range(10):
-        t = step*dt
-        mode.set_time(t)
-        fields = discr.interpolate_volume_function(c_sol)
-
-        er = fields[0:3]
-        hr = fields[3:6]
-        ei = fields[6:9]
-        hi = fields[9:12]
-
-        silo = SiloFile("em-complex-%04d.silo" % step)
-        vis.add_to_silo(silo,
-                vectors=[
-                    ("curl_er", curl(er)), 
-                    ("om_hi", -mode.mu*mode.omega*hi), 
-                    ("curl_hr", curl(hr)), 
-                    ("om_ei", mode.epsilon*mode.omega*hi), 
-                    ],
-                expressions=[
-                ("diff_er", "curl_er-om_hi", DB_VARTYPE_VECTOR),
-                ("diff_hr", "curl_hr-om_ei", DB_VARTYPE_VECTOR),
-                ],
-                write_coarse_mesh=True,
-                time=t, step=step
-                )
-
-        er_res = curl(er) + mode.mu     *mode.omega*hi
-        ei_res = curl(ei) - mode.mu     *mode.omega*hr
-        hr_res = curl(hr) - mode.epsilon*mode.omega*ei
-        hi_res = curl(hi) + mode.epsilon*mode.omega*er
-
-        print("time=%f, rel l2 residual in Re[E]=%g\tIm[E]=%g\tRe[H]=%g\tIm[H]=%g" % (
-                t,
-                rel_l2_error(er_res, er),
-                rel_l2_error(ei_res, ei),
-                rel_l2_error(hr_res, hr),
-                rel_l2_error(hi_res, hi),
-                ))
-
diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py
index 6fb7d54ac4eb69cf87705fce00268f68debfd4ba..85a7c3dec0cef99b4a72fe9c23439c4212c0c32d 100644
--- a/examples/maxwell/cavities.py
+++ b/examples/maxwell/cavities.py
@@ -1,260 +1,131 @@
-# 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
-from __future__ import absolute_import
-from __future__ import print_function
-import numpy as np
-
-import logging
-logger = logging.getLogger(__name__)
-
-
-def main(write_output=True, allow_features=None, flux_type_arg=1,
-        bdry_flux_type_arg=None, extra_discr_args={}):
-    from grudge.mesh.generator import make_cylinder_mesh, make_box_mesh
-    from grudge.tools import EOCRecorder, to_obj_array
-    from math import sqrt, pi  # noqa
-    from analytic_solutions import (  # noqa
-            RealPartAdapter,
-            SplitComplexAdapter,
-            CylindricalFieldAdapter,
-            CylindricalCavityMode,
-            RectangularWaveguideMode,
-            RectangularCavityMode)
-    from grudge.models.em import MaxwellOperator
-
-    logging.basicConfig(level=logging.DEBUG)
-
-    from grudge.backends import guess_run_context
-    rcon = guess_run_context(allow_features)
-
-    epsilon0 = 8.8541878176e-12  # C**2 / (N m**2)
-    mu0 = 4*pi*1e-7  # N/A**2.
-    epsilon = 1*epsilon0
-    mu = 1*mu0
-
-    eoc_rec = EOCRecorder()
-
-    cylindrical = False
-    periodic = False
-
-    if cylindrical:
-        R = 1
-        d = 2
-        mode = CylindricalCavityMode(m=1, n=1, p=1,
-                radius=R, height=d,
-                epsilon=epsilon, mu=mu)
-        # r_sol = CylindricalFieldAdapter(RealPartAdapter(mode))
-        # c_sol = SplitComplexAdapter(CylindricalFieldAdapter(mode))
-
-        if rcon.is_head_rank:
-            mesh = make_cylinder_mesh(radius=R, height=d, max_volume=0.01)
-    else:
-        if periodic:
-            mode = RectangularWaveguideMode(epsilon, mu, (3, 2, 1))
-            periodicity = (False, False, True)
-        else:
-            periodicity = None
-        mode = RectangularCavityMode(epsilon, mu, (1, 2, 2))
-
-        if rcon.is_head_rank:
-            mesh = make_box_mesh(max_volume=0.001, periodicity=periodicity)
-
-    if rcon.is_head_rank:
-        mesh_data = rcon.distribute_mesh(mesh)
-    else:
-        mesh_data = rcon.receive_mesh()
-
-    for order in [4, 5, 6]:
-    #for order in [1,2,3,4,5,6]:
-        extra_discr_args.setdefault("debug", []).extend([
-            "cuda_no_plan", "cuda_dump_kernels"])
-
-        op = MaxwellOperator(epsilon, mu,
-                flux_type=flux_type_arg,
-                bdry_flux_type=bdry_flux_type_arg)
-
-        discr = rcon.make_discretization(mesh_data, order=order,
-                tune_for=op.sym_operator(),
-                **extra_discr_args)
-
-        from grudge.visualization import VtkVisualizer
-        if write_output:
-            vis = VtkVisualizer(discr, rcon, "em-%d" % order)
-
-        mode.set_time(0)
-
-        def get_true_field():
-            return discr.convert_volume(
-                to_obj_array(mode(discr)
-                    .real.astype(discr.default_scalar_type).copy()),
-                kind=discr.compute_kind)
-
-        fields = get_true_field()
+"""Minimal example of a grudge driver."""
 
-        if rcon.is_head_rank:
-            print("---------------------------------------------")
-            print("order %d" % order)
-            print("---------------------------------------------")
-            print("#elements=", len(mesh.elements))
+from __future__ import division, print_function
 
-        from grudge.timestep.runge_kutta import LSRK4TimeStepper
-        stepper = LSRK4TimeStepper(dtype=discr.default_scalar_type, rcon=rcon)
-        #from grudge.timestep.dumka3 import Dumka3TimeStepper
-        #stepper = Dumka3TimeStepper(3, dtype=discr.default_scalar_type, rcon=rcon)
+__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
 
-        # {{{ diagnostics setup
+__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:
 
-        from pytools.log import LogManager, add_general_quantities, \
-                add_simulation_quantities, add_run_info
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
 
-        if write_output:
-            log_file_name = "maxwell-%d.dat" % order
-        else:
-            log_file_name = None
+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.
+"""
 
-        logmgr = LogManager(log_file_name, "w", rcon.communicator)
 
-        add_run_info(logmgr)
-        add_general_quantities(logmgr)
-        add_simulation_quantities(logmgr)
-        discr.add_instrumentation(logmgr)
-        stepper.add_instrumentation(logmgr)
-
-        from pytools.log import IntervalTimer
-        vis_timer = IntervalTimer("t_vis", "Time spent visualizing")
-        logmgr.add_quantity(vis_timer)
-
-        from grudge.log import EMFieldGetter, add_em_quantities
-        field_getter = EMFieldGetter(discr, op, lambda: fields)
-        add_em_quantities(logmgr, op, field_getter)
-
-        logmgr.add_watches(
-                ["step.max", "t_sim.max",
-                    ("W_field", "W_el+W_mag"),
-                    "t_step.max"]
-                )
-
-        # }}}
-
-        # {{{ timestep loop
-
-        rhs = op.bind(discr)
-        final_time = 0.5e-9
-
-        try:
-            from grudge.timestep import times_and_steps
-            step_it = times_and_steps(
-                    final_time=final_time, logmgr=logmgr,
-                    max_dt_getter=lambda t: op.estimate_timestep(discr,
-                        stepper=stepper, t=t, fields=fields))
+import numpy as np
+import pyopencl as cl
+from grudge.shortcuts import set_up_rk4
+from grudge import sym, bind, Discretization
 
-            for step, t, dt in step_it:
-                if step % 50 == 0 and write_output:
-                    sub_timer = vis_timer.start_sub_timer()
-                    e, h = op.split_eh(fields)
-                    visf = vis.make_file("em-%d-%04d" % (order, step))
-                    vis.add_data(
-                            visf,
-                            [
-                                ("e",
-                                    discr.convert_volume(e, kind="numpy")),
-                                ("h",
-                                    discr.convert_volume(h, kind="numpy")),
-                                ],
-                            time=t, step=step)
-                    visf.close()
-                    sub_timer.stop().submit()
+from grudge.models.em import get_rectangular_cavity_mode
 
-                fields = stepper(fields, t, dt, rhs)
 
-            mode.set_time(final_time)
+def main(dims, write_output=True, order=4):
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
 
-            eoc_rec.add_data_point(order,
-                    discr.norm(fields-get_true_field()))
+    from meshmode.mesh.generation import generate_regular_rect_mesh
+    mesh = generate_regular_rect_mesh(
+            a=(0.0,)*dims,
+            b=(1.0,)*dims,
+            n=(5,)*dims)
 
-        finally:
-            if write_output:
-                vis.close()
+    discr = Discretization(cl_ctx, mesh, order=order)
 
-            logmgr.close()
-            discr.close()
+    if 0:
+        epsilon0 = 8.8541878176e-12  # C**2 / (N m**2)
+        mu0 = 4*np.pi*1e-7  # N/A**2.
+        epsilon = 1*epsilon0
+        mu = 1*mu0
+    else:
+        epsilon = 1
+        mu = 1
 
-        if rcon.is_head_rank:
-            print()
-            print(eoc_rec.pretty_print("P.Deg.", "L2 Error"))
+    from grudge.models.em import MaxwellOperator
+    op = MaxwellOperator(epsilon, mu, flux_type=0.5, dimensions=dims)
 
-        # }}}
+    queue = cl.CommandQueue(discr.cl_context)
 
-    assert eoc_rec.estimate_order_of_convergence()[0, 1] > 6
+    if dims == 3:
+        sym_mode = get_rectangular_cavity_mode(1, (1, 2, 2))
+        fields = bind(discr, sym_mode)(queue, t=0, epsilon=epsilon, mu=mu)
+    else:
+        sym_mode = get_rectangular_cavity_mode(1, (2, 3))
+        fields = bind(discr, sym_mode)(queue, t=0)
 
+    # FIXME
+    #dt = op.estimate_rk4_timestep(discr, fields=fields)
 
-# {{{ entry points for py.test
+    op.check_bc_coverage(mesh)
 
-from pytools.test import mark_test
+    # print(sym.pretty(op.sym_operator()))
+    bound_op = bind(discr, op.sym_operator())
 
+    def rhs(t, w):
+        return bound_op(queue, t=t, w=w)
 
-@mark_test.long
-def test_maxwell_cavities():
-    main(write_output=False)
+    if mesh.dim == 2:
+        dt = 0.004
+    elif mesh.dim == 3:
+        dt = 0.002
 
+    dt_stepper = set_up_rk4("w", dt, fields, rhs)
 
-@mark_test.long
-def test_maxwell_cavities_lf():
-    main(write_output=False, flux_type_arg="lf", bdry_flux_type_arg=1)
+    final_t = dt * 600
+    nsteps = int(final_t/dt)
 
+    print("dt=%g nsteps=%d" % (dt, nsteps))
 
-@mark_test.mpi
-@mark_test.long
-def test_maxwell_cavities_mpi():
-    from pytools.mpi import run_with_mpi_ranks
-    run_with_mpi_ranks(__file__, 2, main,
-            write_output=False, allow_features=["mpi"])
+    from grudge.shortcuts import make_visualizer
+    vis = make_visualizer(discr, vis_order=order)
 
+    step = 0
 
-def test_cuda():
-    try:
-        from pycuda.tools import mark_cuda_test
-    except ImportError:
-        pass
+    norm = bind(discr, sym.norm(2, sym.var("u")))
 
-    yield "SP CUDA Maxwell", mark_cuda_test(
-            do_test_maxwell_cavities_cuda), np.float32
-    yield "DP CUDA Maxwell", mark_cuda_test(
-            do_test_maxwell_cavities_cuda), np.float64
+    from time import time
+    t_last_step = time()
 
+    e, h = op.split_eh(fields)
 
-def do_test_maxwell_cavities_cuda(dtype):
-    import pytest  # noqa
+    if 1:
+        vis.write_vtk_file("fld-%04d.vtu" % step,
+                [
+                    ("e", e),
+                    ("h", h),
+                    ])
 
-    main(write_output=False, allow_features=["cuda"],
-            extra_discr_args=dict(
-                debug=["cuda_no_plan"],
-                default_scalar_type=dtype,
-                ))
+    for event in dt_stepper.run(t_end=final_t):
+        if isinstance(event, dt_stepper.StateComputed):
+            assert event.component_id == "w"
 
-# }}}
+            step += 1
 
+            print(step, event.t, norm(queue, u=e[0]), norm(queue, u=e[1]),
+                    norm(queue, u=h[0]), norm(queue, u=h[1]),
+                    time()-t_last_step)
+            if step % 10 == 0:
+                e, h = op.split_eh(event.state_component)
+                vis.write_vtk_file("fld-%04d.vtu" % step,
+                        [
+                            ("e", e),
+                            ("h", h),
+                            ])
+            t_last_step = time()
 
-# entry point -----------------------------------------------------------------
 
 if __name__ == "__main__":
-    from pytools.mpi import in_mpi_relaunch
-    if in_mpi_relaunch():
-        test_maxwell_cavities_mpi()
-    else:
-        do_test_maxwell_cavities_cuda(np.float32)
+    main(3)
diff --git a/grudge/discretization.py b/grudge/discretization.py
index cf640c2c4a4338463e41625530fbbf98d7dddb71..a1ca6303ce40f787ddce458438c88e983b8b94d3 100644
--- a/grudge/discretization.py
+++ b/grudge/discretization.py
@@ -25,6 +25,7 @@ THE SOFTWARE.
 
 from pytools import memoize_method
 from grudge import sym
+import numpy as np
 
 
 # FIXME Naming not ideal
@@ -292,6 +293,14 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
     def ambient_dim(self):
         return self._volume_discr.ambient_dim
 
+    @property
+    def real_dtype(self):
+        return self._volume_discr.real_dtype
+
+    @property
+    def complex_dtype(self):
+        return self._volume_discr.complex_dtype
+
     @property
     def mesh(self):
         return self._volume_discr.mesh
@@ -318,6 +327,11 @@ class PointsDiscretization(DiscretizationBase):
 
     def __init__(self, nodes):
         self._nodes = nodes
+        self.real_dtype = np.dtype(np.float64)
+        self.complex_dtype = np.dtype({
+                np.float32: np.complex64,
+                np.float64: np.complex128
+        }[self.real_dtype.type])
 
     def ambient_dim(self):
         return self._nodes.shape[0]
diff --git a/grudge/execution.py b/grudge/execution.py
index 99a3d4e5aa3bb971b83489583247e9625f3638f1..e4aebe5ca18b495acf4e608ff9f4cdfd218f57b1 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -95,14 +95,40 @@ class ExecutionMapper(mappers.Evaluator,
 
         # FIXME: Make a way to register functions
 
-        pars = [self.rec(p) for p in expr.parameters]
-        if any(isinstance(par, cl.array.Array) for par in pars):
-            import pyopencl.clmath as clmath
-            func = getattr(clmath, expr.function.name)
-        else:
+        args = [self.rec(p) for p in expr.parameters]
+        from numbers import Number
+        representative_arg = args[0]
+        if (
+                isinstance(representative_arg,  Number)
+                or (isinstance(representative_arg, np.ndarray)
+                    and representative_arg.shape == ())):
             func = getattr(np, expr.function.name)
+            return func(*args)
+
+        cached_name = "map_call_knl_"
+
+        i = Variable("i")
+        func = Variable(expr.function.name)
+        if expr.function.name == "fabs":  # FIXME
+            func = Variable("abs")
+            cached_name += "abs"
+        else:
+            cached_name += expr.function.name
+
+        @memoize_in(self.bound_op, cached_name)
+        def knl():
+            knl = lp.make_kernel(
+                "{[i]: 0<=i<n}",
+                [
+                    lp.Assignment(Variable("out")[i],
+                        func(Variable("a")[i]))
+                ], default_offset=lp.auto)
+            return lp.split_iname(knl, "i", 128, outer_tag="g.0", inner_tag="l.0")
+
+        assert len(args) == 1
+        evt, (out,) = knl()(self.queue, a=args[0])
 
-        return func(*pars)
+        return out
 
     def map_nodal_sum(self, op, field_expr):
         # FIXME: Could allow array scalars
@@ -291,7 +317,6 @@ class ExecutionMapper(mappers.Evaluator,
     # }}}
 
     # {{{ code execution functions
-
     def map_insn_loopy_kernel(self, insn):
         kwargs = {}
         kdescr = insn.kernel_descriptor
@@ -450,27 +475,41 @@ class BoundOperator(object):
 
 # {{{ process_sym_operator function
 
-def process_sym_operator(sym_operator, post_bind_mapper=None,
-        dumper=lambda name, sym_operator: None, mesh=None):
+def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None,
+        dumper=lambda name, sym_operator: None):
 
     import grudge.symbolic.mappers as mappers
 
     dumper("before-bind", sym_operator)
     sym_operator = mappers.OperatorBinder()(sym_operator)
 
-    mappers.ErrorChecker(mesh)(sym_operator)
+    mappers.ErrorChecker(discrwb.mesh)(sym_operator)
 
     if post_bind_mapper is not None:
         dumper("before-postbind", sym_operator)
         sym_operator = post_bind_mapper(sym_operator)
 
-    if mesh is not None:
-        dumper("before-empty-flux-killer", sym_operator)
-        sym_operator = mappers.EmptyFluxKiller(mesh)(sym_operator)
+    dumper("before-empty-flux-killer", sym_operator)
+    sym_operator = mappers.EmptyFluxKiller(discrwb.mesh)(sym_operator)
 
     dumper("before-cfold", sym_operator)
     sym_operator = mappers.CommutativeConstantFoldingMapper()(sym_operator)
 
+    # Work around https://github.com/numpy/numpy/issues/9438
+    #
+    # The idea is that we need 1j as an expression to survive
+    # until code generation time. If it is evaluated and combined
+    # with other constants, we will need to determine its size
+    # (as np.complex64/128) within the expression. But because
+    # of the above numpy bug, sized numbers are not likely to survive
+    # expression building--so that's why we step in here to fix that.
+
+    dumper("before-csize", sym_operator)
+    sym_operator = mappers.ConstantToNumpyConversionMapper(
+            real_type=discrwb.real_dtype.type,
+            complex_type=discrwb.complex_dtype.type,
+            )(sym_operator)
+
     # Ordering restriction:
     #
     # - Must run constant fold before first type inference pass, because zeros,
@@ -491,9 +530,8 @@ def process_sym_operator(sym_operator, post_bind_mapper=None,
     # because otherwise it won't differentiate the type of grids (node or quadrature
     # grids) that the operators will apply on.
 
-    assert mesh is not None
     dumper("before-global-to-reference", sym_operator)
-    sym_operator = mappers.GlobalToReferenceMapper(mesh.ambient_dim)(sym_operator)
+    sym_operator = mappers.GlobalToReferenceMapper(discrwb.ambient_dim)(sym_operator)
 
     # Ordering restriction:
     #
@@ -534,10 +572,11 @@ def bind(discr, sym_operator, post_bind_mapper=lambda x: x,
                     pretty(sym_operator))
             stage[0] += 1
 
-    sym_operator = process_sym_operator(sym_operator,
+    sym_operator = process_sym_operator(
+            discr,
+            sym_operator,
             post_bind_mapper=post_bind_mapper,
-            dumper=dump_optemplate,
-            mesh=discr.mesh)
+            dumper=dump_optemplate)
 
     from grudge.symbolic.compiler import OperatorCompiler
     discr_code, eval_code = OperatorCompiler(discr)(sym_operator)
diff --git a/grudge/models/em.py b/grudge/models/em.py
index 3a2b70cbac3d00dab20420e060ceb073f6c1447a..f0e44f903d226e28710709651e948293a534d1b8 100644
--- a/grudge/models/em.py
+++ b/grudge/models/em.py
@@ -1,11 +1,13 @@
 # -*- coding: utf8 -*-
 """grudge operators modelling electromagnetic phenomena."""
 
-from __future__ import division
-from __future__ import absolute_import
-from six.moves import range
+from __future__ import division, absolute_import
 
-__copyright__ = "Copyright (C) 2007, 2010 Andreas Kloeckner, David Powell"
+__copyright__ = """
+Copyright (C) 2007-2017 Andreas Kloeckner
+Copyright (C) 2010 David Powell
+Copyright (C) 2017 Bogdan Enache
+"""
 
 __license__ = """
 Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -27,11 +29,13 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
+from six.moves import range
+
 from pytools import memoize_method
 
 from grudge.models import HyperbolicOperator
-from grudge.symbolic.primitives import make_common_subexpression as cse
 from meshmode.mesh import BTAG_ALL, BTAG_NONE
+from grudge import sym
 from pytools.obj_array import join_fields, make_obj_array
 
 # TODO: Check PML
@@ -106,17 +110,8 @@ class MaxwellOperator(HyperbolicOperator):
         self.current = current
         self.incident_bc_data = incident_bc
 
-    @property
-    def c(self):
-        from warnings import warn
-        warn("MaxwellOperator.c is deprecated", DeprecationWarning)
-        if not self.fixed_material:
-            raise RuntimeError("Cannot compute speed of light "
-                    "for non-constant material")
-
-        return 1/(self.mu*self.epsilon)**0.5
 
-    def flux(self, flux_type):
+    def flux(self, w):
         """The template for the numerical flux for variable coefficients.
 
         :param flux_type: can be in [0,1] for anything between central and upwind,
@@ -124,88 +119,72 @@ class MaxwellOperator(HyperbolicOperator):
 
         As per Hesthaven and Warburton page 433.
         """
-        from grudge.flux import (make_normal, FluxVectorPlaceholder,
-                FluxConstantPlaceholder)
 
-        normal = make_normal(self.dimensions)
+        normal = sym.normal(w.dd, self.dimensions)
 
         if self.fixed_material:
-            from grudge.tools import count_subset
-            w = FluxVectorPlaceholder(count_subset(self.get_eh_subset()))
-
             e, h = self.split_eh(w)
-            epsilon = FluxConstantPlaceholder(self.epsilon)
-            mu = FluxConstantPlaceholder(self.mu)
-
-        else:
-            from grudge.tools import count_subset
-            w = FluxVectorPlaceholder(count_subset(self.get_eh_subset())+2)
+            epsilon = self.epsilon
+            mu = self.mu
 
-            epsilon, mu, e, h = self.split_eps_mu_eh(w)
 
-        Z_int = (mu.int/epsilon.int)**0.5
+        Z_int = (mu/epsilon)**0.5
         Y_int = 1/Z_int
-        Z_ext = (mu.ext/epsilon.ext)**0.5
+        Z_ext = (mu/epsilon)**0.5
         Y_ext = 1/Z_ext
 
-        if flux_type == "lf":
+        if self.flux_type == "lf":
             if self.fixed_material:
                 max_c = (self.epsilon*self.mu)**(-0.5)
-            else:
-                from grudge.flux import Max
-                c_int = (epsilon.int*mu.int)**(-0.5)
-                c_ext = (epsilon.ext*mu.ext)**(-0.5)
-                max_c = Max(c_int, c_ext)  # noqa
 
             return join_fields(
                     # flux e,
                     1/2*(
                         -self.space_cross_h(normal, h.int-h.ext)
                         # multiplication by epsilon undoes material divisor below
-                        #-max_c*(epsilon.int*e.int - epsilon.ext*e.ext)
+                        #-max_c*(epsilon*e.int - epsilon*e.ext)
                     ),
                     # flux h
                     1/2*(
                         self.space_cross_e(normal, e.int-e.ext)
                         # multiplication by mu undoes material divisor below
-                        #-max_c*(mu.int*h.int - mu.ext*h.ext)
+                        #-max_c*(mu*h.int - mu*h.ext)
                     ))
-        elif isinstance(flux_type, (int, float)):
+        elif isinstance(self.flux_type, (int, float)):
             # see doc/maxima/maxwell.mac
             return join_fields(
                     # flux e,
                     (
                         -1/(Z_int+Z_ext)*self.space_cross_h(normal,
                             Z_ext*(h.int-h.ext)
-                            - flux_type*self.space_cross_e(normal, e.int-e.ext))
+                            - self.flux_type*self.space_cross_e(normal, e.int-e.ext))
                         ),
                     # flux h
                     (
                         1/(Y_int + Y_ext)*self.space_cross_e(normal,
                             Y_ext*(e.int-e.ext)
-                            + flux_type*self.space_cross_h(normal, h.int-h.ext))
+                            + self.flux_type*self.space_cross_h(normal, h.int-h.ext))
                         ),
                     )
         else:
             raise ValueError("maxwell: invalid flux_type (%s)"
                     % self.flux_type)
 
-    def local_derivatives(self, w=None):
+    def local_derivatives(self, w):
         """Template for the spatial derivatives of the relevant components of
         :math:`E` and :math:`H`
         """
 
-        e, h = self.split_eh(self.field_placeholder(w))
+        e, h = self.split_eh(w)
 
+        nabla = sym.nabla(self.dimensions)
         def e_curl(field):
             return self.space_cross_e(nabla, field)
 
         def h_curl(field):
             return self.space_cross_h(nabla, field)
 
-        from grudge.symbolic import make_nabla
 
-        nabla = make_nabla(self.dimensions)
 
         # in conservation form: u_t + A u_x = 0
         return join_fields(
@@ -213,62 +192,44 @@ class MaxwellOperator(HyperbolicOperator):
                 e_curl(e)
                 )
 
-    def field_placeholder(self, w=None):
-        "A placeholder for E and H."
-        from grudge.tools import count_subset
-        fld_cnt = count_subset(self.get_eh_subset())
-        if w is None:
-            from grudge.symbolic import make_sym_vector
-            w = make_sym_vector("w", fld_cnt)
-
-        return w
-
-    def pec_bc(self, w=None):
+    def pec_bc(self, w):
         "Construct part of the flux operator template for PEC boundary conditions"
-        e, h = self.split_eh(self.field_placeholder(w))
+        e, h = self.split_eh(w)
+
+        pec_e = sym.cse(sym.interp("vol", self.pec_tag)(e))
+        pec_h = sym.cse(sym.interp("vol", self.pec_tag)(h))
 
-        from grudge.symbolic import RestrictToBoundary
-        pec_e = RestrictToBoundary(self.pec_tag)(e)
-        pec_h = RestrictToBoundary(self.pec_tag)(h)
 
         return join_fields(-pec_e, pec_h)
 
-    def pmc_bc(self, w=None):
+    def pmc_bc(self, w):
         "Construct part of the flux operator template for PMC boundary conditions"
-        e, h = self.split_eh(self.field_placeholder(w))
+        e, h = self.split_eh(w)
 
-        from grudge.symbolic import RestrictToBoundary
-        pmc_e = RestrictToBoundary(self.pmc_tag)(e)
-        pmc_h = RestrictToBoundary(self.pmc_tag)(h)
+        pmc_e = sym.cse(sym.interp("vol", self.pmc_tag)(e))
+        pmc_h = sym.cse(sym.interp("vol", self.pmc_tag)(h))
 
         return join_fields(pmc_e, -pmc_h)
 
-    def absorbing_bc(self, w=None):
+    def absorbing_bc(self, w):
         """Construct part of the flux operator template for 1st order
         absorbing boundary conditions.
         """
 
-        from grudge.symbolic import normal
-        absorb_normal = normal(self.absorb_tag, self.dimensions)
+        absorb_normal = sym.normal(self.absorb_tag, self.dimensions)
 
-        from grudge.symbolic import RestrictToBoundary, Field
 
-        e, h = self.split_eh(self.field_placeholder(w))
+        e, h = self.split_eh(w)
 
         if self.fixed_material:
             epsilon = self.epsilon
             mu = self.mu
-        else:
-            epsilon = cse(
-                    RestrictToBoundary(self.absorb_tag)(Field("epsilon")))
-            mu = cse(
-                    RestrictToBoundary(self.absorb_tag)(Field("mu")))
 
         absorb_Z = (mu/epsilon)**0.5
         absorb_Y = 1/absorb_Z
 
-        absorb_e = RestrictToBoundary(self.absorb_tag)(e)
-        absorb_h = RestrictToBoundary(self.absorb_tag)(h)
+        absorb_e = sym.cse(sym.interp("vol", self.absorb_tag)(e))
+        absorb_h = sym.cse(sym.interp("vol", self.absorb_tag)(h))
 
         bc = join_fields(
                 absorb_e + 1/2*(self.space_cross_h(absorb_normal, self.space_cross_e(
@@ -281,17 +242,12 @@ class MaxwellOperator(HyperbolicOperator):
 
         return bc
 
-    def incident_bc(self, w=None):
+    def incident_bc(self, w):
         "Flux terms for incident boundary conditions"
         # NOTE: Untested for inhomogeneous materials, but would usually be
         # physically meaningless anyway (are there exceptions to this?)
 
-        e, h = self.split_eh(self.field_placeholder(w))
-        if not self.fixed_material:
-            from warnings import warn
-            if self.incident_tag != BTAG_NONE:
-                warn("Incident boundary conditions assume homogeneous"
-                     " background material, results may be unphysical")
+        e, h = self.split_eh(w)
 
         from grudge.tools import count_subset
         fld_cnt = count_subset(self.get_eh_subset())
@@ -301,7 +257,7 @@ class MaxwellOperator(HyperbolicOperator):
         if is_zero(incident_bc_data):
             return make_obj_array([0]*fld_cnt)
         else:
-            return cse(-incident_bc_data)
+            return sym.cse(-incident_bc_data)
 
     def sym_operator(self, w=None):
         """The full operator template - the high level description of
@@ -310,23 +266,9 @@ class MaxwellOperator(HyperbolicOperator):
         Combines the relevant operator templates for spatial
         derivatives, flux, boundary conditions etc.
         """
-        w = self.field_placeholder(w)
-
-        if self.fixed_material:
-            flux_w = w
-        else:
-            epsilon = self.epsilon
-            mu = self.mu
-
-            flux_w = join_fields(epsilon, mu, w)
-
-        from grudge.symbolic import BoundaryPair, \
-                InverseMassOperator, get_flux_operator
-
-        flux_op = get_flux_operator(self.flux(self.flux_type))
-        bdry_flux_op = get_flux_operator(self.flux(self.bdry_flux_type))
+        from grudge.tools import count_subset
+        w = sym.make_sym_array("w", count_subset(self.get_eh_subset()))
 
-        from grudge.tools.indexing import count_subset
         elec_components = count_subset(self.get_eh_subset()[0:3])
         mag_components = count_subset(self.get_eh_subset()[3:6])
 
@@ -334,10 +276,6 @@ class MaxwellOperator(HyperbolicOperator):
             # need to check this
             material_divisor = (
                     [self.epsilon]*elec_components+[self.mu]*mag_components)
-        else:
-            material_divisor = join_fields(
-                    [epsilon]*elec_components,
-                    [mu]*mag_components)
 
         tags_and_bcs = [
                 (self.pec_tag, self.pec_bc(w)),
@@ -346,67 +284,20 @@ class MaxwellOperator(HyperbolicOperator):
                 (self.incident_tag, self.incident_bc(w)),
                 ]
 
-        def make_flux_bc_vector(tag, bc):
-            if self.fixed_material:
-                return bc
-            else:
-                from grudge.symbolic import RestrictToBoundary
-                return join_fields(
-                        cse(RestrictToBoundary(tag)(epsilon)),
-                        cse(RestrictToBoundary(tag)(mu)),
-                        bc)
+
+        def flux(pair):
+            return sym.interp(pair.dd, "all_faces")(self.flux(pair))
 
         return (
                 - self.local_derivatives(w)
-                + InverseMassOperator()(
-                    flux_op(flux_w)
+                - sym.InverseMassOperator()(sym.FaceMassOperator()(
+                    flux(sym.int_tpair(w))
                     + sum(
-                        bdry_flux_op(BoundaryPair(
-                            flux_w, make_flux_bc_vector(tag, bc), tag))
-                        for tag, bc in tags_and_bcs))
-                    ) / material_divisor
-
-    def bind(self, discr):
-        "Convert the abstract operator template into compiled code."
-        from grudge.mesh import check_bc_coverage
-        check_bc_coverage(discr.mesh, [
-            self.pec_tag, self.absorb_tag, self.incident_tag])
-
-        compiled_sym_operator = discr.compile(self.op_template())
-
-        def rhs(t, w, **extra_context):
-            kwargs = {}
-            kwargs.update(extra_context)
-
-            return compiled_sym_operator(w=w, t=t, **kwargs)
-
-        return rhs
-
-    def assemble_eh(self, e=None, h=None, discr=None):
-        "Combines separate E and H vectors into a single array."
-        if discr is None:
-            def zero():
-                return 0
-        else:
-            def zero():
-                return discr.volume_zeros()
-
-        from grudge.tools import count_subset
-        e_components = count_subset(self.get_eh_subset()[0:3])
-        h_components = count_subset(self.get_eh_subset()[3:6])
-
-        def default_fld(fld, comp):
-            if fld is None:
-                return [zero() for i in range(comp)]
-            else:
-                return fld
-
-        e = default_fld(e, e_components)
-        h = default_fld(h, h_components)
+                        flux(sym.bv_tpair(tag, w, bc))
+                            for tag, bc in tags_and_bcs)
+                    ) )) / material_divisor
 
-        return join_fields(e, h)
 
-    assemble_fields = assemble_eh
 
     @memoize_method
     def partial_to_eh_subsets(self):
@@ -445,11 +336,8 @@ class MaxwellOperator(HyperbolicOperator):
         e_idx, h_idx = self.partial_to_eh_subsets()
         e, h = w[e_idx], w[h_idx]
 
-        from grudge.flux import FluxVectorPlaceholder as FVP
-        if isinstance(w, FVP):
-            return FVP(scalars=e), FVP(scalars=h)
-        else:
-            return make_obj_array(e), make_obj_array(h)
+        return e,h
+        #return make_obj_array(e), make_obj_array(h)
 
     def get_eh_subset(self):
         """Return a 6-tuple of :class:`bool` objects indicating whether field
@@ -476,6 +364,14 @@ class MaxwellOperator(HyperbolicOperator):
             raise ValueError("max_eigenvalue is no longer supported for "
                     "variable-coefficient problems--use max_eigenvalue_expr")
 
+    def check_bc_coverage(self, mesh):
+        from meshmode.mesh import check_bc_coverage
+        check_bc_coverage(mesh, [
+            self.pec_tag,
+            self.pmc_tag,
+            self.absorb_tag,
+            self.incident_tag])
+
 
 class TMMaxwellOperator(MaxwellOperator):
     """A 2D TM Maxwell operator with PEC boundaries.
@@ -539,3 +435,61 @@ class SourceFree1DMaxwellOperator(MaxwellOperator):
                 +
                 (False, False, True)
                 )
+
+def get_rectangular_cavity_mode(E_0, mode_indices):
+    """A rectangular TM cavity mode for a rectangle / cube
+    with one corner at the origin and the other at (1,1[,1])."""
+    dims = len(mode_indices)
+    if dims != 2 and dims != 3:
+        raise ValueError("Improper mode_indices dimensions")
+    import numpy
+
+    factors = [n*numpy.pi for n in mode_indices]
+
+    kx, ky = factors[0:2]
+    if dims == 3:
+        kz = factors[2]
+
+    omega = numpy.sqrt(sum(f**2 for f in factors))
+
+
+    nodes = sym.nodes(dims)
+    x = nodes[0]
+    y = nodes[1]
+    if dims == 3:
+        z = nodes[2]
+
+    sx = sym.sin(kx*x)
+    cx = sym.cos(kx*x)
+    sy = sym.sin(ky*y)
+    cy = sym.cos(ky*y)
+    if dims == 3:
+        sz = sym.sin(kz*z)
+        cz = sym.cos(kz*z)
+
+    if dims == 2:
+        tfac = sym.ScalarVariable("t") * omega
+
+        result = sym.join_fields(
+                0,
+                0,
+                sym.sin(kx * x) * sym.sin(ky * y) * sym.cos(tfac),  # ez
+                -ky * sym.sin(kx * x) * sym.cos(ky * y) * sym.sin(tfac) / omega,  # hx
+                kx * sym.cos(kx * x) * sym.sin(ky * y) * sym.sin(tfac) / omega,  # hy
+                0,
+                )
+    else:
+        tdep = sym.exp(-1j * omega * sym.ScalarVariable("t"))
+
+        gamma_squared = ky**2 + kx**2
+        result = sym.join_fields(
+            -kx * kz * E_0*cx*sy*sz*tdep / gamma_squared,  # ex
+            -ky * kz * E_0*sx*cy*sz*tdep / gamma_squared,  # ey
+            E_0 * sx*sy*cz*tdep,  # ez
+
+            -1j * omega * ky*E_0*sx*cy*cz*tdep / gamma_squared,  # hx
+            1j * omega * kx*E_0*cx*sy*cz*tdep / gamma_squared,
+            0,
+            )
+
+    return result
diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py
index ac7765ed919f87ed7b8a17099a302263b94ce347..f6e238083266d4ba0c510243b116b2e49b024631 100644
--- a/grudge/symbolic/compiler.py
+++ b/grudge/symbolic/compiler.py
@@ -852,7 +852,11 @@ class ToLoopyExpressionMapper(mappers.IdentityMapper):
     def map_call(self, expr):
         if isinstance(expr.function, sym.CFunction):
             from pymbolic import var
-            return var(expr.function.name)(
+            func_name = expr.function.name
+            if func_name == "fabs":
+                func_name = "abs"
+
+            return var(func_name)(
                     *[self.rec(par) for par in expr.parameters])
         else:
             raise NotImplementedError(
diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py
index 0fd8fbd1e22e4c3ef2f9a93f758c46abf3f46c3a..405cb5f06706116266322c3d9ab0ab1c206a6c22 100644
--- a/grudge/symbolic/mappers/__init__.py
+++ b/grudge/symbolic/mappers/__init__.py
@@ -34,6 +34,7 @@ import pymbolic.mapper.evaluator
 import pymbolic.mapper.dependency
 import pymbolic.mapper.substitutor
 import pymbolic.mapper.constant_folder
+import pymbolic.mapper.constant_converter
 import pymbolic.mapper.flop_counter
 from pymbolic.mapper import CSECachingMapperMixin
 
@@ -812,6 +813,12 @@ class QuadratureUpsamplerChanger(CSECachingMapperMixin, IdentityMapper):
 
 # {{{ simplification / optimization
 
+class ConstantToNumpyConversionMapper(
+        pymbolic.mapper.constant_converter.ConstantToNumpyConversionMapper,
+        IdentityMapperMixin):
+    pass
+
+
 class CommutativeConstantFoldingMapper(
         pymbolic.mapper.constant_folder.CommutativeConstantFoldingMapper,
         IdentityMapperMixin):
diff --git a/grudge/tools.py b/grudge/tools.py
index cd1cc49f1828ba82f95c0e805262315e5f8e2be6..88d6714364e250d137ee7988558c91138acc1eb7 100644
--- a/grudge/tools.py
+++ b/grudge/tools.py
@@ -25,6 +25,7 @@ THE SOFTWARE.
 """
 
 import numpy as np
+from pytools import levi_civita
 
 
 def is_zero(x):
@@ -36,3 +37,90 @@ def is_zero(x):
         return x == 0
     else:
         return False
+
+
+class SubsettableCrossProduct:
+    """A cross product that can operate on an arbitrary subsets of its
+    two operands and return an arbitrary subset of its result.
+    """
+
+    full_subset = (True, True, True)
+
+    def __init__(self, op1_subset=full_subset, op2_subset=full_subset,
+            result_subset=full_subset):
+        """Construct a subset-able cross product.
+        :param op1_subset: The subset of indices of operand 1 to be taken into
+            account.  Given as a 3-sequence of bools.
+        :param op2_subset: The subset of indices of operand 2 to be taken into
+            account.  Given as a 3-sequence of bools.
+        :param result_subset: The subset of indices of the result that are
+            calculated.  Given as a 3-sequence of bools.
+        """
+        def subset_indices(subset):
+            return [i for i, use_component in enumerate(subset)
+                    if use_component]
+
+        self.op1_subset = op1_subset
+        self.op2_subset = op2_subset
+        self.result_subset = result_subset
+
+        import pymbolic
+        op1 = pymbolic.var("x")
+        op2 = pymbolic.var("y")
+
+        self.functions = []
+        self.component_lcjk = []
+        for i, use_component in enumerate(result_subset):
+            if use_component:
+                this_expr = 0
+                this_component = []
+                for j, j_real in enumerate(subset_indices(op1_subset)):
+                    for k, k_real in enumerate(subset_indices(op2_subset)):
+                        lc = levi_civita((i, j_real, k_real))
+                        if lc != 0:
+                            this_expr += lc*op1.index(j)*op2.index(k)
+                            this_component.append((lc, j, k))
+                self.functions.append(pymbolic.compile(this_expr,
+                    variables=[op1, op2]))
+                self.component_lcjk.append(this_component)
+
+    def __call__(self, x, y, three_mult=None):
+        """Compute the subsetted cross product on the indexables *x* and *y*.
+        :param three_mult: a function of three arguments *sign, xj, yk*
+          used in place of the product *sign*xj*yk*. Defaults to just this
+          product if not given.
+        """
+        from pytools.obj_array import join_fields
+        if three_mult is None:
+            return join_fields(*[f(x, y) for f in self.functions])
+        else:
+            return join_fields(
+                    *[sum(three_mult(lc, x[j], y[k]) for lc, j, k in lcjk)
+                    for lcjk in self.component_lcjk])
+
+
+cross = SubsettableCrossProduct()
+
+
+def count_subset(subset):
+    from pytools import len_iterable
+    return len_iterable(uc for uc in subset if uc)
+
+
+def partial_to_all_subset_indices(subsets, base=0):
+    """Takes a sequence of bools and generates it into an array of indices
+    to be used to insert the subset into the full set.
+    Example:
+    >>> list(partial_to_all_subset_indices([[False, True, True], [True,False,True]]))
+    [array([0 1]), array([2 3]
+    """
+
+    idx = base
+    for subset in subsets:
+        result = []
+        for is_in in subset:
+            if is_in:
+                result.append(idx)
+                idx += 1
+
+        yield np.array(result, dtype=np.intp)
diff --git a/test/test_grudge.py b/test/test_grudge.py
index 6602699039c357fadac3fc78cc7262285ce596f2..84c46168ede26efe1a42935501d69fba99b7700a 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -333,6 +333,76 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type
     assert eoc_rec.order_estimate() > order
 
 
+@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):
+    """Test whether 3D maxwells actually converges"""
+
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+
+    from pytools.convergence import EOCRecorder
+    eoc_rec = EOCRecorder()
+
+    dims = 3
+    ns = [4, 6, 8]
+    for n in ns:
+        from meshmode.mesh.generation import generate_regular_rect_mesh
+        mesh = generate_regular_rect_mesh(
+                a=(0.0,)*dims,
+                b=(1.0,)*dims,
+                n=(n,)*dims)
+
+        discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
+
+        epsilon = 1
+        mu = 1
+
+        from grudge.models.em import get_rectangular_cavity_mode
+        sym_mode = get_rectangular_cavity_mode(1, (1, 2, 2))
+
+        analytic_sol = bind(discr, sym_mode)
+        fields = analytic_sol(queue, t=0, epsilon=epsilon, mu=mu)
+
+        from grudge.models.em import MaxwellOperator
+        op = MaxwellOperator(epsilon, mu, flux_type=0.5, dimensions=dims)
+        op.check_bc_coverage(mesh)
+        bound_op = bind(discr, op.sym_operator())
+
+        def rhs(t, w):
+            return bound_op(queue, t=t, w=w)
+
+        dt = 0.002
+        final_t = dt * 5
+        nsteps = int(final_t/dt)
+
+        from grudge.shortcuts import set_up_rk4
+        dt_stepper = set_up_rk4("w", dt, fields, rhs)
+
+        print("dt=%g nsteps=%d" % (dt, nsteps))
+
+        norm = bind(discr, sym.norm(2, sym.var("u")))
+
+        step = 0
+        for event in dt_stepper.run(t_end=final_t):
+            if isinstance(event, dt_stepper.StateComputed):
+                assert event.component_id == "w"
+                esc = event.state_component
+
+                step += 1
+                print(step)
+
+        sol = analytic_sol(queue, mu=mu, epsilon=epsilon, t=step * dt)
+        vals = [norm(queue, u=(esc[i] - sol[i])) / norm(queue, u=sol[i]) for i in range(5)] # noqa E501
+        total_error = sum(vals)
+        eoc_rec.add_data_point(1.0/n, total_error)
+
+    print(eoc_rec.pretty_print(abscissa_label="h",
+            error_label="L2 Error"))
+
+    assert eoc_rec.order_estimate() > order
+
+
 def test_foreign_points(ctx_factory):
     pytest.importorskip("sumpy")
     import sumpy.point_calculus as pc