diff --git a/examples/ellipsoid.step b/examples/ellipsoid.step
new file mode 100644
index 0000000000000000000000000000000000000000..883cd89ff3599e534d3f66fde3f3c71ea9cb2b31
--- /dev/null
+++ b/examples/ellipsoid.step
@@ -0,0 +1,138 @@
+ISO-10303-21;
+HEADER;
+FILE_DESCRIPTION(('FreeCAD Model'),'2;1');
+FILE_NAME('/home/andreas/ellipsoid.step','2017-06-05T10:39:15',('Author'
+    ),(''),'Open CASCADE STEP processor 6.8','FreeCAD','Unknown');
+FILE_SCHEMA(('AUTOMOTIVE_DESIGN_CC2 { 1 2 10303 214 -1 1 5 4 }'));
+ENDSEC;
+DATA;
+#1 = APPLICATION_PROTOCOL_DEFINITION('committee draft',
+  'automotive_design',1997,#2);
+#2 = APPLICATION_CONTEXT(
+  'core data for automotive mechanical design processes');
+#3 = SHAPE_DEFINITION_REPRESENTATION(#4,#10);
+#4 = PRODUCT_DEFINITION_SHAPE('','',#5);
+#5 = PRODUCT_DEFINITION('design','',#6,#9);
+#6 = PRODUCT_DEFINITION_FORMATION('','',#7);
+#7 = PRODUCT('ASSEMBLY','ASSEMBLY','',(#8));
+#8 = MECHANICAL_CONTEXT('',#2,'mechanical');
+#9 = PRODUCT_DEFINITION_CONTEXT('part definition',#2,'design');
+#10 = SHAPE_REPRESENTATION('',(#11,#15),#19);
+#11 = AXIS2_PLACEMENT_3D('',#12,#13,#14);
+#12 = CARTESIAN_POINT('',(0.,0.,0.));
+#13 = DIRECTION('',(0.,0.,1.));
+#14 = DIRECTION('',(1.,0.,-0.));
+#15 = AXIS2_PLACEMENT_3D('',#16,#17,#18);
+#16 = CARTESIAN_POINT('',(0.,0.,0.));
+#17 = DIRECTION('',(0.,0.,1.));
+#18 = DIRECTION('',(1.,0.,0.));
+#19 = ( GEOMETRIC_REPRESENTATION_CONTEXT(3) 
+GLOBAL_UNCERTAINTY_ASSIGNED_CONTEXT((#23)) GLOBAL_UNIT_ASSIGNED_CONTEXT(
+(#20,#21,#22)) REPRESENTATION_CONTEXT('Context #1',
+  '3D Context with UNIT and UNCERTAINTY') );
+#20 = ( LENGTH_UNIT() NAMED_UNIT(*) SI_UNIT(.MILLI.,.METRE.) );
+#21 = ( NAMED_UNIT(*) PLANE_ANGLE_UNIT() SI_UNIT($,.RADIAN.) );
+#22 = ( NAMED_UNIT(*) SI_UNIT($,.STERADIAN.) SOLID_ANGLE_UNIT() );
+#23 = UNCERTAINTY_MEASURE_WITH_UNIT(LENGTH_MEASURE(1.E-07),#20,
+  'distance_accuracy_value','confusion accuracy');
+#24 = PRODUCT_TYPE('part',$,(#7));
+#25 = ADVANCED_BREP_SHAPE_REPRESENTATION('',(#11,#26),#69);
+#26 = MANIFOLD_SOLID_BREP('',#27);
+#27 = CLOSED_SHELL('',(#28));
+#28 = ADVANCED_FACE('',(#29),#33,.T.);
+#29 = FACE_BOUND('',#30,.T.);
+#30 = VERTEX_LOOP('',#31);
+#31 = VERTEX_POINT('',#32);
+#32 = CARTESIAN_POINT('',(2.449293598295E-16,-5.999039130647E-32,-2.));
+#33 = ( BOUNDED_SURFACE() B_SPLINE_SURFACE(2,2,(
+    (#34,#35,#36,#37,#38)
+    ,(#39,#40,#41,#42,#43)
+    ,(#44,#45,#46,#47,#48)
+    ,(#49,#50,#51,#52,#53)
+    ,(#54,#55,#56,#57,#58)
+    ,(#59,#60,#61,#62,#63)
+    ,(#64,#65,#66,#67,#68
+)),.UNSPECIFIED.,.T.,.F.,.F.) B_SPLINE_SURFACE_WITH_KNOTS((1,2,2,2,2,1),
+  (3,2,3),(-2.094395102393,0.,2.094395102393,4.188790204786,
+    6.28318530718,8.377580409573),(-1.570796326795,0.,1.570796326795),
+.UNSPECIFIED.) GEOMETRIC_REPRESENTATION_ITEM() RATIONAL_B_SPLINE_SURFACE
+((
+    (1.,0.707106781187,1.,0.707106781187,1.)
+    ,(0.5,0.353553390593,0.5,0.353553390593,0.5)
+    ,(1.,0.707106781187,1.,0.707106781187,1.)
+    ,(0.5,0.353553390593,0.5,0.353553390593,0.5)
+    ,(1.,0.707106781187,1.,0.707106781187,1.)
+    ,(0.5,0.353553390593,0.5,0.353553390593,0.5)
+,(1.,0.707106781187,1.,0.707106781187,1.
+))) REPRESENTATION_ITEM('') SURFACE() );
+#34 = CARTESIAN_POINT('',(2.449293598295E-16,0.,-2.));
+#35 = CARTESIAN_POINT('',(4.,0.,-2.));
+#36 = CARTESIAN_POINT('',(4.,0.,0.));
+#37 = CARTESIAN_POINT('',(4.,0.,2.));
+#38 = CARTESIAN_POINT('',(2.449293598295E-16,0.,2.));
+#39 = CARTESIAN_POINT('',(2.449293598295E-16,4.2423009549E-16,-2.));
+#40 = CARTESIAN_POINT('',(4.,6.928203230276,-2.));
+#41 = CARTESIAN_POINT('',(4.,6.928203230276,0.));
+#42 = CARTESIAN_POINT('',(4.,6.928203230276,2.));
+#43 = CARTESIAN_POINT('',(2.449293598295E-16,4.2423009549E-16,2.));
+#44 = CARTESIAN_POINT('',(-1.224646799147E-16,2.12115047745E-16,-2.));
+#45 = CARTESIAN_POINT('',(-2.,3.464101615138,-2.));
+#46 = CARTESIAN_POINT('',(-2.,3.464101615138,0.));
+#47 = CARTESIAN_POINT('',(-2.,3.464101615138,2.));
+#48 = CARTESIAN_POINT('',(-1.224646799147E-16,2.12115047745E-16,2.));
+#49 = CARTESIAN_POINT('',(-4.898587196589E-16,5.999039130647E-32,-2.));
+#50 = CARTESIAN_POINT('',(-8.,9.797174393179E-16,-2.));
+#51 = CARTESIAN_POINT('',(-8.,9.797174393179E-16,0.));
+#52 = CARTESIAN_POINT('',(-8.,9.797174393179E-16,2.));
+#53 = CARTESIAN_POINT('',(-4.898587196589E-16,5.999039130647E-32,2.));
+#54 = CARTESIAN_POINT('',(-1.224646799147E-16,-2.12115047745E-16,-2.));
+#55 = CARTESIAN_POINT('',(-2.,-3.464101615138,-2.));
+#56 = CARTESIAN_POINT('',(-2.,-3.464101615138,0.));
+#57 = CARTESIAN_POINT('',(-2.,-3.464101615138,2.));
+#58 = CARTESIAN_POINT('',(-1.224646799147E-16,-2.12115047745E-16,2.));
+#59 = CARTESIAN_POINT('',(2.449293598295E-16,-4.2423009549E-16,-2.));
+#60 = CARTESIAN_POINT('',(4.,-6.928203230276,-2.));
+#61 = CARTESIAN_POINT('',(4.,-6.928203230276,0.));
+#62 = CARTESIAN_POINT('',(4.,-6.928203230276,2.));
+#63 = CARTESIAN_POINT('',(2.449293598295E-16,-4.2423009549E-16,2.));
+#64 = CARTESIAN_POINT('',(2.449293598295E-16,0.,-2.));
+#65 = CARTESIAN_POINT('',(4.,0.,-2.));
+#66 = CARTESIAN_POINT('',(4.,0.,0.));
+#67 = CARTESIAN_POINT('',(4.,0.,2.));
+#68 = CARTESIAN_POINT('',(2.449293598295E-16,0.,2.));
+#69 = ( GEOMETRIC_REPRESENTATION_CONTEXT(3) 
+GLOBAL_UNCERTAINTY_ASSIGNED_CONTEXT((#73)) GLOBAL_UNIT_ASSIGNED_CONTEXT(
+(#70,#71,#72)) REPRESENTATION_CONTEXT('Context #1',
+  '3D Context with UNIT and UNCERTAINTY') );
+#70 = ( LENGTH_UNIT() NAMED_UNIT(*) SI_UNIT(.MILLI.,.METRE.) );
+#71 = ( NAMED_UNIT(*) PLANE_ANGLE_UNIT() SI_UNIT($,.RADIAN.) );
+#72 = ( NAMED_UNIT(*) SI_UNIT($,.STERADIAN.) SOLID_ANGLE_UNIT() );
+#73 = UNCERTAINTY_MEASURE_WITH_UNIT(LENGTH_MEASURE(1.E-07),#70,
+  'distance_accuracy_value','confusion accuracy');
+#74 = SHAPE_DEFINITION_REPRESENTATION(#75,#25);
+#75 = PRODUCT_DEFINITION_SHAPE('','',#76);
+#76 = PRODUCT_DEFINITION('design','',#77,#80);
+#77 = PRODUCT_DEFINITION_FORMATION('','',#78);
+#78 = PRODUCT('Ellipsoid','Ellipsoid','',(#79));
+#79 = MECHANICAL_CONTEXT('',#2,'mechanical');
+#80 = PRODUCT_DEFINITION_CONTEXT('part definition',#2,'design');
+#81 = CONTEXT_DEPENDENT_SHAPE_REPRESENTATION(#82,#84);
+#82 = ( REPRESENTATION_RELATIONSHIP('','',#25,#10) 
+REPRESENTATION_RELATIONSHIP_WITH_TRANSFORMATION(#83) 
+SHAPE_REPRESENTATION_RELATIONSHIP() );
+#83 = ITEM_DEFINED_TRANSFORMATION('','',#11,#15);
+#84 = PRODUCT_DEFINITION_SHAPE('Placement','Placement of an item',#85);
+#85 = NEXT_ASSEMBLY_USAGE_OCCURRENCE('1','=>[0:1:1:2]','',#5,#76,$);
+#86 = PRODUCT_TYPE('part',$,(#78));
+#87 = MECHANICAL_DESIGN_GEOMETRIC_PRESENTATION_REPRESENTATION('',(#88),
+  #69);
+#88 = STYLED_ITEM('color',(#89),#28);
+#89 = PRESENTATION_STYLE_ASSIGNMENT((#90));
+#90 = SURFACE_STYLE_USAGE(.BOTH.,#91);
+#91 = SURFACE_SIDE_STYLE('',(#92));
+#92 = SURFACE_STYLE_FILL_AREA(#93);
+#93 = FILL_AREA_STYLE('',(#94));
+#94 = FILL_AREA_STYLE_COLOUR('',#95);
+#95 = COLOUR_RGB('',0.800000011921,0.800000011921,0.800000011921);
+ENDSEC;
+END-ISO-10303-21;
diff --git a/examples/layerpot-3d.py b/examples/layerpot-3d.py
index da1f50e7df7a620cb5acd84dc089dd1790ce2773..68d9e33e81e8fa6041f7a873e31f824bd2d8177f 100644
--- a/examples/layerpot-3d.py
+++ b/examples/layerpot-3d.py
@@ -90,7 +90,9 @@ if 1:
                 ]
             )
 
-    bdry_normals = bind(density_discr, sym.normal())(queue).as_vector(dtype=object)
+    bdry_normals = bind(
+            density_discr,
+            sym.normal(density_discr.ambient_dim))(queue).as_vector(dtype=object)
 
     from meshmode.discretization.visualization import make_visualizer
     bdry_vis = make_visualizer(queue, density_discr, target_order)
diff --git a/examples/maxwell.py b/examples/maxwell.py
new file mode 100644
index 0000000000000000000000000000000000000000..108a96c3f82cb170292d1f2c77be7ae107609144
--- /dev/null
+++ b/examples/maxwell.py
@@ -0,0 +1,141 @@
+from __future__ import division
+import numpy as np
+import pyopencl as cl
+from sumpy.visualization import FieldPlotter
+#from mayavi import mlab
+from sumpy.kernel import one_kernel_2d, LaplaceKernel, HelmholtzKernel  # noqa
+
+import faulthandler
+from six.moves import range
+faulthandler.enable()
+
+import logging
+logging.basicConfig(level=logging.INFO)
+logger = logging.getLogger(__name__)
+
+cl_ctx = cl.create_some_context()
+queue = cl.CommandQueue(cl_ctx)
+
+target_order = 5
+qbx_order = 3
+nelements = 60
+mode_nr = 0
+
+h = 1
+
+k = 4
+
+
+def main():
+    from meshmode.mesh.io import generate_gmsh, FileSource
+    mesh = generate_gmsh(
+            FileSource("ellipsoid.step"), 2, order=2,
+            other_options=["-string", "Mesh.CharacteristicLengthMax = %g;" % h])
+
+    from meshmode.mesh.processing import perform_flips
+    # Flip elements--gmsh generates inside-out geometry.
+    mesh = perform_flips(mesh, np.ones(mesh.nelements))
+
+    print("%d elements" % mesh.nelements)
+
+    from meshmode.mesh.processing import find_bounding_box
+    bbox_min, bbox_max = find_bounding_box(mesh)
+    bbox_center = 0.5*(bbox_min+bbox_max)
+    bbox_size = max(bbox_max-bbox_min) / 2
+
+    logger.info("%d elements" % mesh.nelements)
+
+    from pytential.qbx import QBXLayerPotentialSource
+    from meshmode.discretization import Discretization
+    from meshmode.discretization.poly_element import \
+            InterpolatoryQuadratureSimplexGroupFactory
+
+    density_discr = Discretization(
+            cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order))
+
+    qbx = QBXLayerPotentialSource(density_discr, 4*target_order, qbx_order,
+            fmm_order=qbx_order + 10, fmm_backend="fmmlib")
+
+    from pytential.symbolic.pde.maxwell import PECAugmentedMFIEOperator
+    pde_op = PECAugmentedMFIEOperator()
+    from pytential import bind, sym
+
+    rho_sym = sym.var("rho")
+    # jt_sym = sym.make_sym_vector("jt", 2)
+    # repr_op = pde_op.scattered_volume_field(jt_sym, rho_sym)
+
+    # {{{ make a density
+
+    nodes = density_discr.nodes().with_queue(queue)
+
+    angle = cl.clmath.atan2(nodes[1], nodes[0])
+
+    sigma = cl.clmath.cos(mode_nr*angle)
+    if 0:
+        sigma = 0*angle
+        from random import randrange
+        for i in range(5):
+            sigma[randrange(len(sigma))] = 1
+
+    sigma = sigma.astype(np.complex128)
+
+    jt = sym.make_obj_array([sigma, sigma])
+    rho = sigma
+
+    # }}}
+
+    #mlab.figure(bgcolor=(1, 1, 1))
+    if 1:
+        from meshmode.discretization.visualization import make_visualizer
+        bdry_vis = make_visualizer(queue, density_discr, target_order)
+
+        bdry_normals = bind(density_discr, sym.normal(3))(queue)\
+                .as_vector(dtype=object)
+
+        bdry_vis.write_vtk_file("source.vtu", [
+            ("sigma", sigma),
+            ("bdry_normals", bdry_normals),
+            ])
+
+        fplot = FieldPlotter(bbox_center, extent=2*bbox_size, npoints=(150, 150, 1))
+
+        qbx_stick_out = qbx.copy(target_stick_out_factor=0.1)
+        from pytential.target import PointsTarget
+        from pytential.qbx import QBXTargetAssociationFailedException
+        try:
+            fld_in_vol = bind(
+                    (qbx_stick_out, PointsTarget(fplot.points)),
+                    sym.make_obj_array([
+                        # sym.S(pde_op.kernel, rho_sym, k=sym.var("k"),
+                        #     qbx_forced_limit=None),
+                        # sym.d_dx(3, sym.S(pde_op.kernel, rho_sym, k=sym.var("k"),
+                        #     qbx_forced_limit=None)),
+                        sym.d_dy(3, sym.S(pde_op.kernel, rho_sym, k=sym.var("k"),
+                            qbx_forced_limit=None)),
+                        # sym.d_dz(3, sym.S(pde_op.kernel, rho_sym, k=sym.var("k"),
+                        #     qbx_forced_limit=None)),
+                        ])
+                    )(queue, jt=jt, rho=rho, k=k)
+        except QBXTargetAssociationFailedException as e:
+            fplot.write_vtk_file(
+                    "failed-targets.vts",
+                    [
+                        ("failed_targets", e.failed_target_flags.get(queue))
+                        ])
+            raise
+
+        fld_in_vol = sym.make_obj_array(
+            [fiv.get() for fiv in fld_in_vol])
+
+        #fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5)
+        fplot.write_vtk_file(
+                "potential.vts",
+                [
+                    ("potential", fld_in_vol[0]),
+                    ("grad", fld_in_vol[1:]),
+                    ]
+                )
+
+
+if __name__ == "__main__":
+    main()
diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py
index ad6578771f7805d9ea17c2beb4bd533e0cf3676e..5292f7bb6685e766ff36861eee0d7dcb3e657bfc 100644
--- a/pytential/qbx/__init__.py
+++ b/pytential/qbx/__init__.py
@@ -78,7 +78,8 @@ class QBXLayerPotentialSource(LayerPotentialSource):
             debug=True,
             refined_for_global_qbx=False,
             expansion_disks_in_tree_have_extent=False,
-            performance_data_file=None):
+            performance_data_file=None,
+            fmm_backend="sumpy"):
         """
         :arg fine_order: The total degree to which the (upsampled)
              underlying quadrature is exact.
@@ -115,6 +116,7 @@ class QBXLayerPotentialSource(LayerPotentialSource):
         self.density_discr = density_discr
         self.fmm_level_to_order = fmm_level_to_order
         self.target_stick_out_factor = target_stick_out_factor
+        self.fmm_backend = fmm_backend
 
         # Default values are lazily provided if these are None
         self._base_resampler = base_resampler
@@ -164,7 +166,8 @@ class QBXLayerPotentialSource(LayerPotentialSource):
                     else self.refined_for_global_qbx),
                 expansion_disks_in_tree_have_extent=(
                     self.expansion_disks_in_tree_have_extent),
-                performance_data_file=self.performance_data_file)
+                performance_data_file=self.performance_data_file,
+                fmm_backend=self.fmm_backend)
 
     # }}}
 
@@ -418,12 +421,24 @@ class QBXLayerPotentialSource(LayerPotentialSource):
         fmm_local_factory = partial(local_expn_class, base_kernel)
         qbx_local_factory = partial(local_expn_class, base_kernel)
 
-        from pytential.qbx.fmm import \
-                QBXExpansionWranglerCodeContainer
-        return QBXExpansionWranglerCodeContainer(
-                self.cl_context,
-                fmm_mpole_factory, fmm_local_factory, qbx_local_factory,
-                out_kernels)
+        if self.fmm_backend == "sumpy":
+            from pytential.qbx.fmm import \
+                    QBXSumpyExpansionWranglerCodeContainer
+            return QBXSumpyExpansionWranglerCodeContainer(
+                    self.cl_context,
+                    fmm_mpole_factory, fmm_local_factory, qbx_local_factory,
+                    out_kernels)
+
+        elif self.fmm_backend == "fmmlib":
+            from pytential.qbx.fmmlib import \
+                    QBXFMMLibExpansionWranglerCodeContainer
+            return QBXFMMLibExpansionWranglerCodeContainer(
+                    self.cl_context,
+                    fmm_mpole_factory, fmm_local_factory, qbx_local_factory,
+                    out_kernels)
+
+        else:
+            raise ValueError("invalid FMM backend: %s" % self.fmm_backend)
 
     def exec_compute_potential_insn_fmm(self, queue, insn, bound_expr, evaluate):
         # {{{ build list of unique target discretizations used
diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py
index 426de45d3789a02fd392fe5fadc92ab286e72b33..c7d6c029448ce9c5c47ebedfc92a91a509a35157 100644
--- a/pytential/qbx/fmm.py
+++ b/pytential/qbx/fmm.py
@@ -37,7 +37,7 @@ logger = logging.getLogger(__name__)
 
 
 __doc__ = """
-.. autoclass:: QBXExpansionWranglerCodeContainer
+.. autoclass:: QBXSumpyExpansionWranglerCodeContainer
 
 .. autoclass:: QBXExpansionWrangler
 
@@ -45,9 +45,9 @@ __doc__ = """
 """
 
 
-# {{{ expansion wrangler
+# {{{ sumpy expansion wrangler
 
-class QBXExpansionWranglerCodeContainer(SumpyExpansionWranglerCodeContainer):
+class QBXSumpyExpansionWranglerCodeContainer(SumpyExpansionWranglerCodeContainer):
     def __init__(self, cl_context,
             multipole_expansion_factory, local_expansion_factory,
             qbx_local_expansion_factory, out_kernels):
@@ -126,10 +126,10 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
 
     # {{{ data vector utilities
 
-    def potential_zeros(self):
-        """This ought to be called ``non_qbx_potential_zeros``, but since
+    def output_zeros(self):
+        """This ought to be called ``non_qbx_output_zeros``, but since
         it has to override the superclass's behavior to integrate seamlessly,
-        it needs to be called just :meth:`potential_zeros`.
+        it needs to be called just :meth:`output_zeros`.
         """
 
         nqbtl = self.geo_data.non_qbx_box_target_lists()
@@ -142,10 +142,10 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
                     dtype=self.dtype)
                 for k in self.code.out_kernels])
 
-    def full_potential_zeros(self):
+    def full_output_zeros(self):
         # The superclass generates a full field of zeros, for all
         # (not just non-QBX) targets.
-        return SumpyExpansionWrangler.potential_zeros(self)
+        return SumpyExpansionWrangler.output_zeros(self)
 
     def qbx_local_expansion_zeros(self):
         order = self.qbx_order
@@ -175,12 +175,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
 
     # {{{ source/target dispatch
 
-    def box_source_list_kwargs(self):
-        return dict(
-                box_source_starts=self.tree.box_source_starts,
-                box_source_counts_nonchild=(
-                    self.tree.box_source_counts_nonchild),
-                sources=self.tree.sources)
+    # box_source_list_kwargs inherited from superclass
 
     def box_target_list_kwargs(self):
         # This only covers the non-QBX targets.
@@ -311,7 +306,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
         return qbx_expansions
 
     def eval_qbx_expansions(self, qbx_expansions):
-        pot = self.full_potential_zeros()
+        pot = self.full_output_zeros()
 
         geo_data = self.geo_data
         if len(geo_data.global_qbx_centers()) == 0:
@@ -342,6 +337,9 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
 
     # }}}
 
+    def finalize_potential(self, potential):
+        return potential
+
 # }}}
 
 
@@ -371,12 +369,12 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     logger.info("start qbx fmm")
 
-    logger.debug("reorder source weights")
+    logger.info("reorder source weights")
     src_weights = wrangler.reorder_sources(src_weights)
 
     # {{{ construct local multipoles
 
-    logger.debug("construct local multipoles")
+    logger.info("construct local multipoles")
     mpole_exps = wrangler.form_multipoles(
             traversal.level_start_source_box_nrs,
             traversal.source_boxes,
@@ -386,7 +384,7 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     # {{{ propagate multipoles upward
 
-    logger.debug("propagate multipoles upward")
+    logger.info("propagate multipoles upward")
     wrangler.coarsen_multipoles(
             traversal.level_start_source_parent_box_nrs,
             traversal.source_parent_boxes,
@@ -396,7 +394,7 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     # {{{ direct evaluation from neighbor source boxes ("list 1")
 
-    logger.debug("direct evaluation from neighbor source boxes ('list 1')")
+    logger.info("direct evaluation from neighbor source boxes ('list 1')")
     non_qbx_potentials = wrangler.eval_direct(
             traversal.target_boxes,
             traversal.neighbor_source_boxes_starts,
@@ -407,7 +405,7 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     # {{{ translate separated siblings' ("list 2") mpoles to local
 
-    logger.debug("translate separated siblings' ('list 2') mpoles to local")
+    logger.info("translate separated siblings' ('list 2') mpoles to local")
     local_exps = wrangler.multipole_to_local(
             traversal.level_start_target_or_target_parent_box_nrs,
             traversal.target_or_target_parent_boxes,
@@ -419,7 +417,7 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     # {{{ evaluate sep. smaller mpoles ("list 3") at particles
 
-    logger.debug("evaluate sep. smaller mpoles at particles ('list 3 far')")
+    logger.info("evaluate sep. smaller mpoles at particles ('list 3 far')")
 
     # (the point of aiming this stage at particles is specifically to keep its
     # contribution *out* of the downward-propagating local expansions)
@@ -437,7 +435,7 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     # {{{ form locals for separated bigger mpoles ("list 4")
 
-    logger.debug("form locals for separated bigger mpoles ('list 4 far')")
+    logger.info("form locals for separated bigger mpoles ('list 4 far')")
 
     local_exps = local_exps + wrangler.form_locals(
             traversal.level_start_target_or_target_parent_box_nrs,
@@ -453,7 +451,7 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     # {{{ propagate local_exps downward
 
-    logger.debug("propagate local_exps downward")
+    logger.info("propagate local_exps downward")
     wrangler.refine_locals(
             traversal.level_start_target_or_target_parent_box_nrs,
             traversal.target_or_target_parent_boxes,
@@ -463,7 +461,7 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     # {{{ evaluate locals
 
-    logger.debug("evaluate locals")
+    logger.info("evaluate locals")
     non_qbx_potentials = non_qbx_potentials + wrangler.eval_locals(
             traversal.level_start_target_box_nrs,
             traversal.target_boxes,
@@ -473,22 +471,22 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     # {{{ wrangle qbx expansions
 
-    logger.debug("form global qbx expansions from list 1")
+    logger.info("form global qbx expansions from list 1")
     qbx_expansions = wrangler.form_global_qbx_locals(
             traversal.neighbor_source_boxes_starts,
             traversal.neighbor_source_boxes_lists,
             src_weights)
 
-    logger.debug("translate from list 3 multipoles to qbx local expansions")
+    logger.info("translate from list 3 multipoles to qbx local expansions")
     qbx_expansions = qbx_expansions + \
             wrangler.translate_box_multipoles_to_qbx_local(mpole_exps)
 
-    logger.debug("translate from box local expansions to contained "
+    logger.info("translate from box local expansions to contained "
             "qbx local expansions")
     qbx_expansions = qbx_expansions + \
             wrangler.translate_box_local_to_qbx_local(local_exps)
 
-    logger.debug("evaluate qbx local expansions")
+    logger.info("evaluate qbx local expansions")
     qbx_potentials = wrangler.eval_qbx_expansions(
             qbx_expansions)
 
@@ -496,17 +494,11 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     # {{{ reorder potentials
 
-    logger.debug("reorder potentials")
+    logger.info("reorder potentials")
 
     nqbtl = geo_data.non_qbx_box_target_lists()
 
-    from pytools.obj_array import make_obj_array
-    all_potentials_in_tree_order = make_obj_array([
-            cl.array.zeros(
-                wrangler.queue,
-                tree.ntargets,
-                dtype=wrangler.dtype)
-            for k in wrangler.code.out_kernels])
+    all_potentials_in_tree_order = wrangler.full_output_zeros()
 
     for ap_i, nqp_i in zip(all_potentials_in_tree_order, non_qbx_potentials):
         ap_i[nqbtl.unfiltered_from_filtered_target_indices] = nqp_i
@@ -514,7 +506,9 @@ def drive_fmm(expansion_wrangler, src_weights):
     all_potentials_in_tree_order += qbx_potentials
 
     def reorder_potentials(x):
-        return x[tree.sorted_target_ids]
+        # "finalize" gives host FMMs (like FMMlib) a chance to turn the
+        # potential back into a CL array.
+        return wrangler.finalize_potential(x[tree.sorted_target_ids])
 
     from pytools.obj_array import with_object_array_or_scalar
     result = with_object_array_or_scalar(
@@ -529,6 +523,8 @@ def drive_fmm(expansion_wrangler, src_weights):
 # }}}
 
 
+# {{{ performance model
+
 def write_performance_model(outf, geo_data):
     from pymbolic import var
     p_fmm = var("p_fmm")
@@ -714,5 +710,6 @@ def write_performance_model(outf, geo_data):
 
     # }}}
 
+# }}}
 
 # vim: foldmethod=marker
diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py
new file mode 100644
index 0000000000000000000000000000000000000000..14d523f11ee39cd84064ff527a028905afa60d1a
--- /dev/null
+++ b/pytential/qbx/fmmlib.py
@@ -0,0 +1,427 @@
+from __future__ import division, absolute_import
+
+__copyright__ = "Copyright (C) 2017 Andreas Kloeckner"
+
+__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
+from pytools import memoize_method
+import pyopencl as cl  # noqa
+import pyopencl.array  # noqa: F401
+from boxtree.pyfmmlib_integration import HelmholtzExpansionWrangler
+from sumpy.kernel import HelmholtzKernel
+
+
+class QBXFMMLibExpansionWranglerCodeContainer(object):
+    def __init__(self, cl_context,
+            multipole_expansion_factory, local_expansion_factory,
+            qbx_local_expansion_factory, out_kernels):
+        self.cl_context = cl_context
+        self.multipole_expansion_factory = multipole_expansion_factory
+        self.local_expansion_factory = local_expansion_factory
+        self.qbx_local_expansion_factory = qbx_local_expansion_factory
+
+        self.out_kernels = out_kernels
+
+    def get_wrangler(self, queue, geo_data, dtype,
+            qbx_order, fmm_level_to_order,
+            source_extra_kwargs={},
+            kernel_extra_kwargs=None):
+
+        return QBXFMMLibHelmholtzExpansionWrangler(self, queue, geo_data, dtype,
+                qbx_order, fmm_level_to_order,
+                source_extra_kwargs,
+                kernel_extra_kwargs)
+
+# }}}
+
+
+# {{{ host geo data wrapper
+
+class ToHostTransferredGeoDataWrapper(object):
+    def __init__(self, queue, geo_data):
+        self.queue = queue
+        self.geo_data = geo_data
+
+    @memoize_method
+    def tree(self):
+        return self.traversal().tree
+
+    @memoize_method
+    def traversal(self):
+        return self.geo_data.traversal().get(queue=self.queue)
+
+    @property
+    def ncenters(self):
+        return self.geo_data.ncenters
+
+    @memoize_method
+    def centers(self):
+        return np.array([
+            ci.get(queue=self.queue)
+            for ci in self.geo_data.centers()])
+
+    @memoize_method
+    def global_qbx_centers(self):
+        return self.geo_data.global_qbx_centers().get(queue=self.queue)
+
+    @memoize_method
+    def qbx_center_to_target_box(self):
+        return self.geo_data.qbx_center_to_target_box().get(queue=self.queue)
+
+    @memoize_method
+    def non_qbx_box_target_lists(self):
+        return self.geo_data.non_qbx_box_target_lists().get(queue=self.queue)
+
+    @memoize_method
+    def center_to_tree_targets(self):
+        return self.geo_data.center_to_tree_targets().get(queue=self.queue)
+
+    @memoize_method
+    def all_targets(self):
+        """All (not just non-QBX) targets packaged into a single array."""
+        return np.array(list(self.tree().targets))
+
+# }}}
+
+
+# {{{ fmmlib expansion wrangler
+
+class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler):
+    def __init__(self, code, queue, geo_data, dtype,
+            qbx_order, fmm_level_to_order,
+            source_extra_kwargs,
+            kernel_extra_kwargs):
+
+        self.code = code
+        self.queue = queue
+
+        # FMMLib is CPU-only. This wrapper gets the geometry out of
+        # OpenCL-land.
+        self.geo_data = ToHostTransferredGeoDataWrapper(queue, geo_data)
+
+        self.qbx_order = qbx_order
+
+        # {{{ digest out_kernels
+
+        from sumpy.kernel import AxisTargetDerivative
+
+        k_names = []
+
+        def is_supported_helmknl(knl):
+            result = isinstance(knl, HelmholtzKernel) and knl.dim == 3
+            if result:
+                k_names.append(knl.helmholtz_k_name)
+            return result
+
+        ifgrad = False
+        outputs = []
+        for out_knl in self.code.out_kernels:
+            if is_supported_helmknl(out_knl):
+                outputs.append(())
+            elif (isinstance(out_knl, AxisTargetDerivative)
+                    and is_supported_helmknl(out_knl.inner_kernel)):
+                outputs.append((out_knl.axis,))
+                ifgrad = True
+            else:
+                raise NotImplementedError(
+                        "only the 3D Helmholtz kernel and its target derivatives "
+                        "are supported for now")
+
+        self.outputs = outputs
+
+        # }}}
+
+        from pytools import single_valued
+        k_name = single_valued(k_names)
+        helmholtz_k = kernel_extra_kwargs[k_name]
+
+        self.level_orders = [
+                fmm_level_to_order(level)
+                for level in range(self.geo_data.tree().nlevels)]
+
+        # FIXME: For now
+        from pytools import single_valued
+        assert single_valued(self.level_orders)
+
+        super(QBXFMMLibHelmholtzExpansionWrangler, self).__init__(
+                self.geo_data.tree(),
+
+                helmholtz_k=helmholtz_k,
+
+                # FIXME
+                nterms=fmm_level_to_order(0),
+
+                ifgrad=ifgrad)
+
+    # {{{ data vector helpers
+
+    def output_zeros(self):
+        """This ought to be called ``non_qbx_output_zeros``, but since
+        it has to override the superclass's behavior to integrate seamlessly,
+        it needs to be called just :meth:`output_zeros`.
+        """
+
+        nqbtl = self.geo_data.non_qbx_box_target_lists()
+
+        from pytools.obj_array import make_obj_array
+        return make_obj_array([
+                np.zeros(nqbtl.nfiltered_targets, self.dtype)
+                for k in self.outputs])
+
+    def full_output_zeros(self):
+        """This includes QBX and non-QBX targets."""
+
+        from pytools.obj_array import make_obj_array
+        return make_obj_array([
+                np.zeros(self.tree.ntargets, self.dtype)
+                for k in self.outputs])
+
+    def reorder_sources(self, source_array):
+        source_array = source_array.get(queue=self.queue)
+        return (super(QBXFMMLibHelmholtzExpansionWrangler, self)
+                .reorder_sources(source_array))
+
+    def reorder_potentials(self, potentials):
+        raise NotImplementedError("reorder_potentials should not "
+            "be called on a QBXFMMLibHelmholtzExpansionWrangler")
+
+        # Because this is a multi-stage, more complicated process that combines
+        # potentials from non-QBX targets and QBX targets.
+
+    def add_potgrad_onto_output(self, output, output_slice, pot, grad):
+        for i_out, out in enumerate(self.outputs):
+            if len(out) == 0:
+                output[i_out][output_slice] += pot
+            elif len(out) == 1:
+                axis, = out
+                if isinstance(grad, np.ndarray):
+                    output[i_out][output_slice] += grad[axis]
+                else:
+                    assert grad == 0
+            else:
+                raise ValueError("element '%s' of outputs array not "
+                        "understood" % out)
+
+    # }}}
+
+    # {{{ override target lists to only hit non-QBX targets
+
+    def box_target_starts(self):
+        nqbtl = self.geo_data.non_qbx_box_target_lists()
+        return nqbtl.box_target_starts
+
+    def box_target_counts_nonchild(self):
+        nqbtl = self.geo_data.non_qbx_box_target_lists()
+        return nqbtl.box_target_counts_nonchild
+
+    def targets(self):
+        nqbtl = self.geo_data.non_qbx_box_target_lists()
+        return nqbtl.targets
+
+    # }}}
+
+    # {{{ qbx-related
+
+    def qbx_local_expansion_zeros(self):
+        return np.zeros(
+                    (self.geo_data.ncenters,) + self.expansion_shape(self.qbx_order),
+                    dtype=self.dtype)
+
+    def form_global_qbx_locals(self, starts, lists, src_weights):
+        local_exps = self.qbx_local_expansion_zeros()
+
+        rscale = 1  # FIXME
+        geo_data = self.geo_data
+
+        if len(geo_data.global_qbx_centers()) == 0:
+            return local_exps
+
+        qbx_center_to_target_box = geo_data.qbx_center_to_target_box()
+        qbx_centers = geo_data.centers()
+
+        formta = self.get_routine("%ddformta")
+
+        for itgt_center, tgt_icenter in enumerate(geo_data.global_qbx_centers()):
+            itgt_box = qbx_center_to_target_box[tgt_icenter]
+
+            isrc_box_start = starts[itgt_box]
+            isrc_box_stop = starts[itgt_box+1]
+
+            tgt_center = qbx_centers[:, tgt_icenter]
+
+            ctr_coeffs = 0
+
+            for isrc_box in range(isrc_box_start, isrc_box_stop):
+                src_ibox = lists[isrc_box]
+
+                src_pslice = self._get_source_slice(src_ibox)
+
+                ier, coeffs = formta(
+                        self.helmholtz_k, rscale,
+                        self._get_sources(src_pslice), src_weights[src_pslice],
+                        tgt_center, self.nterms)
+                if ier:
+                    raise RuntimeError("formta failed")
+
+                ctr_coeffs += coeffs
+
+            local_exps[tgt_icenter] += ctr_coeffs
+
+        return local_exps
+
+    def translate_box_multipoles_to_qbx_local(self, multipole_exps):
+        local_exps = self.qbx_local_expansion_zeros()
+
+        geo_data = self.geo_data
+        qbx_center_to_target_box = geo_data.qbx_center_to_target_box()
+        qbx_centers = geo_data.centers()
+        centers = self.tree.box_centers
+
+        mploc = self.get_translation_routine("%ddmploc")
+
+        for isrc_level, ssn in enumerate(geo_data.traversal().sep_smaller_by_level):
+            source_level_start_ibox, source_mpoles_view = \
+                    self.multipole_expansions_view(multipole_exps, isrc_level)
+
+            # FIXME
+            rscale = 1
+
+            kwargs = {}
+            if self.dim == 3:
+                # FIXME Is this right?
+                kwargs["radius"] = self.tree.root_extent * 2**(-isrc_level)
+
+            for itgt_center, tgt_icenter in enumerate(geo_data.global_qbx_centers()):
+                ctr_loc = 0
+
+                icontaining_tgt_box = qbx_center_to_target_box[tgt_icenter]
+
+                tgt_center = qbx_centers[:, tgt_icenter]
+
+                for isrc_box in range(
+                        ssn.starts[icontaining_tgt_box],
+                        ssn.starts[icontaining_tgt_box+1]):
+
+                    src_ibox = ssn.lists[isrc_box]
+                    src_center = centers[:, src_ibox]
+
+                    ctr_loc = ctr_loc + mploc(
+                        self.helmholtz_k,
+                        rscale, src_center, multipole_exps[src_ibox],
+                        rscale, tgt_center, self.nterms, **kwargs)[..., 0]
+
+                local_exps[tgt_icenter] += ctr_loc
+
+        return local_exps
+
+    def translate_box_local_to_qbx_local(self, local_exps):
+        qbx_expansions = self.qbx_local_expansion_zeros()
+
+        geo_data = self.geo_data
+        if geo_data.ncenters == 0:
+            return qbx_expansions
+        trav = geo_data.traversal()
+        qbx_center_to_target_box = geo_data.qbx_center_to_target_box()
+        qbx_centers = geo_data.centers()
+
+        rscale = 1  # FIXME
+
+        locloc = self.get_translation_routine("%ddlocloc")
+
+        for isrc_level in range(geo_data.tree().nlevels):
+            local_order = self.level_orders[isrc_level]
+
+            lev_box_start, lev_box_stop = self.tree.level_start_box_nrs[
+                    isrc_level:isrc_level+2]
+            target_level_start_ibox, target_locals_view = \
+                    self.local_expansions_view(local_exps, isrc_level)
+            assert target_level_start_ibox == lev_box_start
+
+            kwargs = {}
+            if self.dim == 3:
+                # FIXME Is this right?
+                kwargs["radius"] = self.tree.root_extent * 2**(-isrc_level)
+
+            for tgt_icenter in range(geo_data.ncenters):
+                isrc_box = qbx_center_to_target_box[tgt_icenter]
+
+                tgt_center = qbx_centers[:, tgt_icenter]
+
+                # The box's expansions which we're translating here
+                # (our source) is, globally speaking, a target box.
+
+                src_ibox = trav.target_boxes[isrc_box]
+
+                # Is the box number on the level currently under
+                # consideration?
+                in_range = (lev_box_start <= src_ibox and src_ibox < lev_box_stop)
+
+                if in_range:
+                    src_center = self.tree.box_centers[:, src_ibox]
+                    tmp_loc_exp = locloc(
+                                self.helmholtz_k,
+                                rscale, src_center, local_exps[src_ibox],
+                                rscale, tgt_center, local_order, **kwargs)[..., 0]
+
+                    qbx_expansions[tgt_icenter] += tmp_loc_exp
+
+        return qbx_expansions
+
+    def eval_qbx_expansions(self, qbx_expansions):
+        output = self.full_output_zeros()
+
+        geo_data = self.geo_data
+        ctt = geo_data.center_to_tree_targets()
+        global_qbx_centers = geo_data.global_qbx_centers()
+        qbx_centers = geo_data.centers()
+
+        all_targets = geo_data.all_targets()
+
+        rscale = 1  # FIXME
+
+        taeval = self.get_expn_eval_routine("ta")
+
+        for isrc_center, src_icenter in enumerate(global_qbx_centers):
+            for icenter_tgt in range(
+                    ctt.starts[src_icenter],
+                    ctt.starts[src_icenter+1]):
+
+                center_itgt = ctt.lists[icenter_tgt]
+
+                center = qbx_centers[:, src_icenter]
+
+                pot, grad = taeval(self.helmholtz_k, rscale,
+                        center, qbx_expansions[src_icenter],
+                        all_targets[:, center_itgt])
+
+                self.add_potgrad_onto_output(output, center_itgt, pot, grad)
+
+        return output
+
+    def finalize_potential(self, potential):
+        return cl.array.to_device(self.queue, potential)
+
+    # }}}
+
+# }}}
+
+# vim: foldmethod=marker
diff --git a/pytential/qbx/interactions.py b/pytential/qbx/interactions.py
index 38271cdf96623cbd226b11835e118c52940294e6..172953f3e91297fc83a36117ccb51c938cb111f3 100644
--- a/pytential/qbx/interactions.py
+++ b/pytential/qbx/interactions.py
@@ -241,6 +241,8 @@ class L2QBXL(E2EBase):
                     <> src_ibox = target_boxes[isrc_box] \
                         {id=read_src_ibox}
 
+                    # Is the box number on the level currently under
+                    # consideration?
                     <> in_range = (target_base_ibox <= src_ibox
                             and src_ibox < target_base_ibox + nboxes)
 
diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py
index 68cceb5eb2c1c64048cca722fb370b3c559aedbc..9096998a9663c5fe43150505595c413d6d0c287a 100644
--- a/pytential/symbolic/pde/maxwell/__init__.py
+++ b/pytential/symbolic/pde/maxwell/__init__.py
@@ -37,7 +37,9 @@ class PECAugmentedMFIEOperator:
     see notes/mfie.tm
     """
 
-    def __init__(self, k):
+    def __init__(self, k=sym.var("k")):
+        from sumpy.kernel import HelmholtzKernel
+        self.kernel = HelmholtzKernel(3)
         self.k = k
 
     def j_operator(self, loc, Jt):
@@ -58,9 +60,8 @@ class PECAugmentedMFIEOperator:
 
     def scattered_boundary_field(self, Jt, rho, loc):
         Jxyz = cse(tangential_to_xyz(Jt), "Jxyz")
-        k = self.k
 
-        A = S(k, Jxyz)
+        A = S(self.kernel, Jxyz, k=self.k)
         grad_phi = grad_S(k, rho, 3)
 
         # use - n x n x v = v_tangential
@@ -72,18 +73,16 @@ class PECAugmentedMFIEOperator:
         return join_fields(E_scat, H_scat)
 
     def scattered_volume_field(self, Jt, rho):
-        Jxyz = cse(tangential_to_xyz(Jt), "Jxyz")
-
-        k = self.k
+        Jxyz = sym.cse(sym.tangential_to_xyz(Jt), "Jxyz")
 
-        A = S(k, Jxyz)
-        grad_phi = grad_S(k, rho, 3)
+        A = sym.S(self.kernel, Jxyz, k=self.k, qbx_forced_limit=None)
+        #grad_phi = grad_S(self.kernel, rho, 3, k=self.k)
 
-        E_scat = 1j*k*A - grad_phi
-        H_scat = curl_S_volume(k, Jxyz)
+        E_scat = 1j*self.k*A# - grad_phi
+        #H_scat = curl_S_volume(k, Jxyz)
 
         from pytools.obj_array import join_fields
-        return join_fields(E_scat, H_scat)
+        return E_scat #join_fields(E_scat, H_scat)
 
 # }}}
 
diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py
index b4d1c61520eeaf57b7349f4300a08359cee51ce8..c4c63b99e80a89ccbe515a2cba99f8edfd6845c4 100644
--- a/pytential/symbolic/primitives.py
+++ b/pytential/symbolic/primitives.py
@@ -894,6 +894,55 @@ def Dp(kernel, *args, **kwargs):  # noqa
 # }}}
 
 
+# {{{ geometric operations
+
+def tangential_onb(ambient_dim, dim=None, where=None):
+    if dim is None:
+        dim = ambient_dim - 1
+
+    pd_mat = cse(parametrization_derivative_matrix(ambient_dim, dim, where),
+            "pd_matrix", cse_scope.DISCRETIZATION)
+
+    # {{{ Gram-Schmidt
+
+    orth_pd_mat = np.zeros_like(pd_mat)
+    for k in range(pd_mat.shape[1]):
+        avec = pd_mat[:, k]
+        q = avec
+        for j in range(k):
+            q = q - np.dot(avec, orth_pd_mat[:, j])*orth_pd_mat[:, j]
+
+        orth_pd_mat[:, k] = cse(q/sqrt(np.sum(q**2)), "orth_pd_vec")
+
+    # }}}
+
+    return orth_pd_mat
+
+
+def xyz_to_tangential(xyz_vec, where=None):
+    ambient_dim = len(xyz_vec)
+    tonb = tangential_onb(ambient_dim, where=where)
+    return make_obj_array([
+        np.dot(tonb[:, i], xyz_vec)
+        for i in range(ambient_dim - 1)
+        ])
+
+
+def tangential_to_xyz(tangential_vec, where=None):
+    ambient_dim = len(tangential_vec) + 1
+    tonb = tangential_onb(ambient_dim, where=where)
+    return sum(
+        tonb[:, i] * tangential_vec[i]
+        for i in range(ambient_dim - 1))
+
+
+def project_to_tangential(xyz_vec, which=None):
+    return tangential_to_xyz(
+            cse(xyz_to_tangential(xyz_vec, which), which))
+
+# }}}
+
+
 def pretty(expr):
     # Doesn't quite belong here, but this is exposed to the user as
     # "pytential.sym", so in here it goes.