diff --git a/doc/conf.py b/doc/conf.py
index 66729e86760a52b858e78cba3a87b5d81b875bf1..f43c6d3596d715838132b0dfc6bfbe6d29eaa386 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -57,7 +57,7 @@ copyright = u'2014, Andreas Klöckner'
 #
 # The short X.Y version.
 ver_dic = {}
-execfile("../meshmode/version.py", ver_dic)
+exec(compile(open("../meshmode/version.py").read(), "../meshmode/version.py", 'exec'), ver_dic)
 version = ".".join(str(x) for x in ver_dic["VERSION"])
 # The full version, including alpha/beta/rc tags.
 release = ver_dic["VERSION_TEXT"]
diff --git a/examples/plot-connectivity.py b/examples/plot-connectivity.py
index 83d76b899a92e0f59d866f198322bac2feebbab0..09d795418a05f16c033f3b969862a1eda3ab07e0 100644
--- a/examples/plot-connectivity.py
+++ b/examples/plot-connectivity.py
@@ -28,7 +28,7 @@ def main():
             cl_ctx, mesh, PolynomialWarpAndBlendGroupFactory(order))
 
     from meshmode.discretization.visualization import make_visualizer
-    vis = make_visualizer(queue, discr, order)
+    vis = make_visualizer(queue, discr, order=1)
     os.remove("geometry.vtu")
     os.remove("connectivity.vtu")
     vis.write_vtk_file("geometry.vtu", [
diff --git a/meshmode/discretization/__init__.py b/meshmode/discretization/__init__.py
index e7a09c236d20f154aec31a2eb381a80f4b353424..6088753222fcdfde70e51b0cca237f2af2dd9f78 100644
--- a/meshmode/discretization/__init__.py
+++ b/meshmode/discretization/__init__.py
@@ -241,10 +241,14 @@ class Discretization(object):
                     result[d, k, i] = \
                         sum(j, resampling_mat[i, j] * nodes[d, k, j])
                     """,
-                name="nodes")
+                name="nodes",
+                default_offset=lp.auto)
 
             knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
-            return lp.tag_inames(knl, dict(k="g.0"))
+            knl = lp.tag_inames(knl, dict(k="g.0"))
+            knl = lp.tag_data_axes(knl, "result",
+                    "stride:auto,stride:auto,stride:auto")
+            return knl
 
         result = self.empty(self.real_dtype, extra_dims=(self.ambient_dim,))
 
diff --git a/meshmode/discretization/connection.py b/meshmode/discretization/connection.py
index 3fdd1db5570f515488dfa8e36ab387683ea27fbf..3f1614170e1deddb8e2aeca626c590114cae34c5 100644
--- a/meshmode/discretization/connection.py
+++ b/meshmode/discretization/connection.py
@@ -1,4 +1,8 @@
 from __future__ import division
+from __future__ import absolute_import
+import six
+from six.moves import range
+from six.moves import zip
 
 __copyright__ = "Copyright (C) 2014 Andreas Kloeckner"
 
@@ -224,7 +228,7 @@ def _build_boundary_connection(queue, vol_discr, bdry_discr, connection_data):
         connection_batches = []
         mgrp = vol_grp.mesh_el_group
 
-        for face_id in xrange(len(mgrp.face_vertex_indices())):
+        for face_id in range(len(mgrp.face_vertex_indices())):
             data = connection_data[igrp, face_id]
 
             bdry_unit_nodes_01 = (bdry_grp.unit_nodes + 1)*0.5
@@ -268,7 +272,7 @@ def make_boundary_restriction(queue, discr, group_factory):
     for igrp, mgrp in enumerate(discr.mesh.groups):
         grp_face_vertex_indices = mgrp.face_vertex_indices()
 
-        for iel_grp in xrange(mgrp.nelements):
+        for iel_grp in range(mgrp.nelements):
             for fid, loc_face_vertices in enumerate(grp_face_vertex_indices):
                 face_vertices = frozenset(
                         mgrp.vertex_indices[iel_grp, fvi]
@@ -283,11 +287,11 @@ def make_boundary_restriction(queue, discr, group_factory):
 
     boundary_faces = [
             face_ids[0]
-            for face_vertices, face_ids in face_map.iteritems()
+            for face_vertices, face_ids in six.iteritems(face_map)
             if len(face_ids) == 1]
 
     from pytools import flatten
-    bdry_vertex_vol_nrs = sorted(set(flatten(face_map.iterkeys())))
+    bdry_vertex_vol_nrs = sorted(set(flatten(six.iterkeys(face_map))))
 
     vol_to_bdry_vertices = np.empty(
             discr.mesh.vertices.shape[-1],
@@ -338,7 +342,7 @@ def make_boundary_restriction(queue, discr, group_factory):
 
         batch_base = 0
 
-        for face_id in xrange(len(grp_face_vertex_indices)):
+        for face_id in range(len(grp_face_vertex_indices)):
             batch_boundary_el_numbers_in_grp = np.array(
                     [
                         ibface_el
diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py
index a2b251ef1c5d3b32c29ff0aa7df14edc4510b04b..3c2565fdb5a7ac206d76b0c2fb5ed44853939e54 100644
--- a/meshmode/discretization/visualization.py
+++ b/meshmode/discretization/visualization.py
@@ -1,4 +1,6 @@
 from __future__ import division
+from __future__ import absolute_import
+from six.moves import range
 
 __copyright__ = "Copyright (C) 2014 Andreas Kloeckner"
 
@@ -270,7 +272,7 @@ def write_mesh_connectivity_vtk_file(file_name, mesh,  compressor=None,):
     connections = np.empty((nconnections, 2), dtype=np.int32)
 
     nb_starts = cnx.neighbors_starts
-    for iel in xrange(mesh.nelements):
+    for iel in range(mesh.nelements):
         connections[nb_starts[iel]:nb_starts[iel+1], 0] = iel
 
     connections[:, 1] = cnx.neighbors
diff --git a/meshmode/mesh.py b/meshmode/mesh.py
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py
index a2e90d9f65f240a872aee35c799f2772ccdcd2da..2347a87a7a287f0fc4d23dc60a167799b1b34eb5 100644
--- a/meshmode/mesh/__init__.py
+++ b/meshmode/mesh/__init__.py
@@ -1,4 +1,6 @@
 from __future__ import division
+from __future__ import absolute_import
+from six.moves import range
 
 __copyright__ = "Copyright (C) 2010,2012,2013 Andreas Kloeckner, Michael Tom"
 
@@ -406,18 +408,18 @@ def _compute_connectivity_from_vertices(mesh):
     # FIXME Native code would make this faster
 
     _, nvertices = mesh.vertices.shape
-    vertex_to_element = [[] for i in xrange(nvertices)]
+    vertex_to_element = [[] for i in range(nvertices)]
 
     for grp in mesh.groups:
         iel_base = grp.element_nr_base
-        for iel_grp in xrange(grp.nelements):
+        for iel_grp in range(grp.nelements):
             for ivertex in grp.vertex_indices[iel_grp]:
                 vertex_to_element[ivertex].append(iel_base + iel_grp)
 
-    element_to_element = [set() for i in xrange(mesh.nelements)]
+    element_to_element = [set() for i in range(mesh.nelements)]
     for grp in mesh.groups:
         iel_base = grp.element_nr_base
-        for iel_grp in xrange(grp.nelements):
+        for iel_grp in range(grp.nelements):
             for ivertex in grp.vertex_indices[iel_grp]:
                 element_to_element[iel_base + iel_grp].update(
                         vertex_to_element[ivertex])
diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index f5ffdee6dc7bd982806bca6baa3aa9e2da0770ff..5b925dcb62dfacc73b90b1180800a737b1ba0a4d 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -1,4 +1,6 @@
 from __future__ import division
+from __future__ import absolute_import
+from six.moves import range
 
 __copyright__ = "Copyright (C) 2013 Andreas Kloeckner"
 
@@ -51,6 +53,7 @@ Surfaces
 .. autofunction:: generate_icosahedron
 .. autofunction:: generate_icosphere
 .. autofunction:: generate_torus
+.. autofunction:: generate_regular_rect_mesh
 
 """
 
@@ -165,6 +168,8 @@ def qbx_peanut(t):
 # }}}
 
 
+# {{{ make_curve_mesh
+
 def make_curve_mesh(curve_f, element_boundaries, order):
     """
     :arg curve_f: A callable representing a parametrization for a curve,
@@ -205,6 +210,10 @@ def make_curve_mesh(curve_f, element_boundaries, order):
     return Mesh(vertices=vertices, groups=[egroup],
             element_connectivity=None)
 
+# }}}
+
+
+# {{{ make_group_from_vertices
 
 def make_group_from_vertices(vertices, vertex_indices, order):
     el_vertices = vertices[:, vertex_indices]
@@ -233,6 +242,10 @@ def make_group_from_vertices(vertices, vertex_indices, order):
             order, vertex_indices, nodes,
             unit_nodes=unit_nodes)
 
+# }}}
+
+
+# {{{ generate_icosahedron
 
 def generate_icosahedron(r, order):
     # http://en.wikipedia.org/w/index.php?title=Icosahedron&oldid=387737307
@@ -270,6 +283,10 @@ def generate_icosahedron(r, order):
     return Mesh(vertices, [grp],
             element_connectivity=None)
 
+# }}}
+
+
+# {{{ generate_icosphere
 
 def generate_icosphere(r, order):
     mesh = generate_icosahedron(r, order)
@@ -283,6 +300,10 @@ def generate_icosphere(r, order):
     return Mesh(mesh.vertices, [grp],
             element_connectivity=None)
 
+# }}}
+
+
+# {{{ generate_torus_and_cycle_vertices
 
 def generate_torus_and_cycle_vertices(r_outer, r_inner,
         n_outer=20, n_inner=10, order=1):
@@ -300,9 +321,9 @@ def generate_torus_and_cycle_vertices(r_outer, r_inner,
     def idx(i, j):
         return (i % n_outer) + (j % n_inner) * n_outer
     vertex_indices = ([(idx(i, j), idx(i+1, j), idx(i, j+1))
-            for i in xrange(n_outer) for j in xrange(n_inner)]
+            for i in range(n_outer) for j in range(n_inner)]
             + [(idx(i+1, j), idx(i+1, j+1), idx(i, j+1))
-            for i in xrange(n_outer) for j in xrange(n_inner)])
+            for i in range(n_outer) for j in range(n_inner)])
 
     vertex_indices = np.array(vertex_indices, dtype=np.int32)
     grp = make_group_from_vertices(vertices, vertex_indices, order)
@@ -340,8 +361,10 @@ def generate_torus_and_cycle_vertices(r_outer, r_inner,
 
     from meshmode.mesh import Mesh
     return (Mesh(vertices, [grp.copy(nodes=nodes)], element_connectivity=None),
-            [idx(i, 0) for i in xrange(n_outer)],
-            [idx(0, j) for j in xrange(n_inner)])
+            [idx(i, 0) for i in range(n_outer)],
+            [idx(0, j) for j in range(n_inner)])
+
+# }}}
 
 
 def generate_torus(r_outer, r_inner, n_outer=20, n_inner=10, order=1):
@@ -349,4 +372,58 @@ def generate_torus(r_outer, r_inner, n_outer=20, n_inner=10, order=1):
             r_outer, r_inner, n_outer, n_inner, order)
     return mesh
 
+
+# {{{ generate_regular_rect_mesh
+
+def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1):
+    """Create a semi-structured rectangular mesh.
+
+    :param a: the lower left hand point of the rectangle
+    :param b: the upper right hand point of the rectangle
+    :param n: a tuple of integers indicating the total number of points
+      on [a,b].
+    """
+    if min(n) < 2:
+        raise ValueError("need at least two points in each direction")
+
+    node_dict = {}
+    vertices = []
+    points_1d = [np.linspace(a_i, b_i, n_i)
+            for a_i, b_i, n_i in zip(a, b, n)]
+
+    for j in range(n[1]):
+        for i in range(n[0]):
+            node_dict[i, j] = len(vertices)
+            vertices.append(np.array([points_1d[0][i], points_1d[1][j]]))
+
+    vertices = np.array(vertices).T.copy()
+
+    vertex_indices = []
+
+    for i in range(n[0]-1):
+        for j in range(n[1]-1):
+
+            # c--d
+            # |  |
+            # a--b
+
+            a = node_dict[i, j]
+            b = node_dict[i+1, j]
+            c = node_dict[i, j+1]
+            d = node_dict[i+1, j+1]
+
+            vertex_indices.append((a, b, c))
+            vertex_indices.append((d, c, b))
+
+    vertex_indices = np.array(vertex_indices, dtype=np.int32)
+
+    grp = make_group_from_vertices(vertices, vertex_indices, order)
+
+    from meshmode.mesh import Mesh
+    return Mesh(vertices, [grp],
+            element_connectivity=None)
+
+# }}}
+
+
 # vim: fdm=marker
diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py
index dc540a8d24704af2544a4b35fe9cc8d7d9a37c62..292ee55077a178c70694bf8fd719b2fe42eef143 100644
--- a/meshmode/mesh/io.py
+++ b/meshmode/mesh/io.py
@@ -1,4 +1,8 @@
 from __future__ import division
+from __future__ import absolute_import
+import six
+from six.moves import range
+from six.moves import zip
 
 __copyright__ = "Copyright (C) 2014 Andreas Kloeckner"
 
@@ -97,7 +101,7 @@ class GmshMeshReceiver(GmshMeshReceiverBase):
         ambient_dim = self.points.shape[-1]
 
         mesh_bulk_dim = max(
-                el_type.dimensions for el_type in el_type_hist.iterkeys())
+                el_type.dimensions for el_type in six.iterkeys(el_type_hist))
 
         # {{{ build vertex numbering
 
@@ -114,7 +118,7 @@ class GmshMeshReceiver(GmshMeshReceiverBase):
         # {{{ build vertex array
 
         gmsh_vertex_indices, my_vertex_indices = \
-                zip(*vertex_index_gmsh_to_mine.iteritems())
+                list(zip(*six.iteritems(vertex_index_gmsh_to_mine)))
         vertices = np.empty(
                 (ambient_dim, len(vertex_index_gmsh_to_mine)), dtype=np.float64)
         vertices[:, np.array(my_vertex_indices, np.intp)] = \
@@ -124,7 +128,7 @@ class GmshMeshReceiver(GmshMeshReceiverBase):
 
         from meshmode.mesh import SimplexElementGroup, Mesh
 
-        for group_el_type, ngroup_elements in el_type_hist.iteritems():
+        for group_el_type, ngroup_elements in six.iteritems(el_type_hist):
             if group_el_type.dimensions != mesh_bulk_dim:
                 continue
 
diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 7d8d1b7898996aaa1edc7a34085e142aed6d6e08..94cd8b0951f5eb937c556bc239a2ebd7110bfb04 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -1,4 +1,7 @@
 from __future__ import division
+from __future__ import absolute_import
+from six.moves import range
+from functools import reduce
 
 __copyright__ = "Copyright (C) 2014 Andreas Kloeckner"
 
@@ -31,6 +34,8 @@ __doc__ = """
 .. autofunction:: find_volume_mesh_element_orientations
 .. autofunction:: perform_flips
 .. autofunction:: find_bounding_box
+.. autofunction:: merge_dijsoint_meshes
+.. autofunction:: affine_map
 """
 
 
@@ -64,8 +69,8 @@ def find_volume_mesh_element_group_orientation(mesh, grp):
             (nspan_vectors, ambient_dim),
             dtype=np.object)
 
-    for ispan in xrange(nspan_vectors):
-        for idim in xrange(ambient_dim):
+    for ispan in range(nspan_vectors):
+        for idim in range(ambient_dim):
             spanning_object_array[ispan, idim] = \
                     spanning_vectors[idim, :, ispan]
 
@@ -212,4 +217,82 @@ def find_bounding_box(mesh):
 
 # }}}
 
+
+# {{{ merging
+
+def merge_dijsoint_meshes(meshes, skip_tests=False):
+    if not meshes:
+        raise ValueError("must pass at least one mesh")
+
+    from pytools import is_single_valued
+    if not is_single_valued(mesh.ambient_dim for mesh in meshes):
+        raise ValueError("all meshes must share the same ambient dimension")
+
+    # {{{ assemble combined vertex array
+
+    ambient_dim = meshes[0].ambient_dim
+    nvertices = sum(
+            mesh.vertices.shape[-1]
+            for mesh in meshes)
+
+    vert_dtype = np.find_common_type(
+            [mesh.vertices.dtype for mesh in meshes],
+            [])
+    vertices = np.empty(
+            (ambient_dim, nvertices), vert_dtype)
+
+    current_vert_base = 0
+    vert_bases = []
+    for mesh in meshes:
+        mesh_nvert = mesh.vertices.shape[-1]
+        vertices[:, current_vert_base:current_vert_base+mesh_nvert] = \
+                mesh.vertices
+
+        vert_bases.append(current_vert_base)
+        current_vert_base += mesh_nvert
+
+    # }}}
+
+    # {{{ assemble new groups list
+
+    new_groups = []
+
+    for mesh, vert_base in zip(meshes, vert_bases):
+        for group in mesh.groups:
+            new_vertex_indices = group.vertex_indices + vert_base
+            new_group = group.copy(vertex_indices=new_vertex_indices)
+            new_groups.append(new_group)
+
+    # }}}
+
+    from meshmode.mesh import Mesh
+    return Mesh(vertices, new_groups, skip_tests=skip_tests)
+
+# }}}
+
+
+# {{{ affine map
+
+def affine_map(mesh, A=None, b=None):
+    """Apply the affine map *f(x)=Ax+b* to the geometry of *mesh*."""
+
+    vertices = np.einsum("ij,jv->iv", A, mesh.vertices) + b[:, np.newaxis]
+
+    # {{{ assemble new groups list
+
+    new_groups = []
+
+    for group in mesh.groups:
+        nodes = (
+                np.einsum("ij,jen->ien", A, group.nodes)
+                + b[:, np.newaxis, np.newaxis])
+        new_groups.append(group.copy(nodes=nodes))
+
+    # }}}
+
+    from meshmode.mesh import Mesh
+    return Mesh(vertices, new_groups, skip_tests=True)
+
+# }}}
+
 # vim: foldmethod=marker
diff --git a/meshmode/mesh/visualization.py b/meshmode/mesh/visualization.py
new file mode 100644
index 0000000000000000000000000000000000000000..7cacd1234fa68e4c794dfe4536bf804a74fffff4
--- /dev/null
+++ b/meshmode/mesh/visualization.py
@@ -0,0 +1,73 @@
+from __future__ import division
+from __future__ import absolute_import
+from six.moves import range
+
+__copyright__ = "Copyright (C) 2014 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
+
+
+# {{{ draw_2d_mesh
+
+def draw_2d_mesh(mesh, draw_numbers=True, **kwargs):
+    assert mesh.ambient_dim == 2
+
+    import matplotlib.pyplot as pt
+    import matplotlib.patches as mpatches
+    from matplotlib.path import Path
+
+    for igrp, grp in enumerate(mesh.groups):
+        for iel, el in enumerate(grp.vertex_indices):
+            elverts = mesh.vertices[:, el]
+
+            pathdata = [
+                (Path.MOVETO, (elverts[0, 0], elverts[1, 0])),
+                ]
+            for i in range(1, elverts.shape[1]):
+                pathdata.append(
+                    (Path.LINETO, (elverts[0, i], elverts[1, i]))
+                    )
+            pathdata.append(
+                (Path.CLOSEPOLY, (elverts[0, 0], elverts[1, 0])))
+
+            codes, verts = zip(*pathdata)
+            path = Path(verts, codes)
+            patch = mpatches.PathPatch(path, **kwargs)
+            pt.gca().add_patch(patch)
+
+            if draw_numbers:
+                centroid = (np.sum(elverts, axis=1)
+                        / elverts.shape[1])
+
+                if len(mesh.groups) == 1:
+                    el_label = str(iel)
+                else:
+                    el_label = "%d:%d" % (igrp, iel)
+
+                pt.text(centroid[0], centroid[1], el_label, fontsize=17,
+                        ha="center", va="center",
+                        bbox=dict(facecolor='white', alpha=0.5, lw=0))
+
+# }}}
+
+# vim: foldmethod=marker
diff --git a/setup.cfg b/setup.cfg
index eca7b43fd4a3ea8f7a21eb2a4fbdb1aab0cd4fa1..6faef2e65abe138f9ab22c94452c43be0e52c075 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,3 +1,3 @@
 [flake8]
-ignore = E126,E127,E128,E123,E226,E241,E242
+ignore = E126,E127,E128,E123,E226,E241,E242,E265
 max-line-length=85
diff --git a/setup.py b/setup.py
index 130908d0a7b2f97e6637e9adb97f69e1a3148d9b..0b15c94320e2073c270f68b247422c0cbbb1b844 100644
--- a/setup.py
+++ b/setup.py
@@ -3,7 +3,7 @@
 
 
 def main():
-    from setuptools import setup
+    from setuptools import setup, find_packages
 
     version_dict = {}
     init_filename = "meshmode/version.py"
@@ -39,11 +39,7 @@ def main():
               'Topic :: Utilities',
               ],
 
-          packages=[
-              "meshmode",
-              "meshmode.mesh",
-              "meshmode.discretization",
-              ],
+          packages=find_packages(),
           install_requires=[
               "numpy",
               "modepy",
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 5febf0308759c11b2bb6b7c8596bf623b2329d56..34521057e135ae0e9eae634685f1a4191d44d9ce 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -1,4 +1,7 @@
 from __future__ import division
+from __future__ import absolute_import
+from __future__ import print_function
+from six.moves import range
 
 __copyright__ = "Copyright (C) 2014 Andreas Kloeckner"
 
@@ -54,19 +57,19 @@ def test_boundary_interpolation(ctx_getter):
     order = 4
 
     for h in [1e-1, 3e-2, 1e-2]:
-        print "BEGIN GEN"
+        print("BEGIN GEN")
         mesh = generate_gmsh(
                 FileSource("blob-2d.step"), 2, order=order,
                 force_ambient_dim=2,
                 other_options=[
                     "-string", "Mesh.CharacteristicLengthMax = %s;" % h]
                 )
-        print "END GEN"
+        print("END GEN")
 
         vol_discr = Discretization(cl_ctx, mesh,
                 InterpolatoryQuadratureSimplexGroupFactory(order))
-        print "h=%s -> %d elements" % (
-                h, sum(mgrp.nelements for mgrp in mesh.groups))
+        print("h=%s -> %d elements" % (
+                h, sum(mgrp.nelements for mgrp in mesh.groups)))
 
         x = vol_discr.nodes()[0].with_queue(queue)
         f = 0.1*cl.clmath.sin(30*x)
@@ -81,7 +84,7 @@ def test_boundary_interpolation(ctx_getter):
         err = la.norm((bdry_f-bdry_f_2).get(), np.inf)
         eoc_rec.add_data_point(h, err)
 
-    print eoc_rec
+    print(eoc_rec)
     assert eoc_rec.order_estimate() >= order-0.5
 
 
@@ -114,6 +117,37 @@ def test_element_orientation():
     assert ((mesh_orient < 0) == (flippy > 0)).all()
 
 
+def test_merge_and_map(ctx_getter, visualize=False):
+    from meshmode.mesh.io import generate_gmsh, FileSource
+
+    mesh_order = 3
+
+    mesh = generate_gmsh(
+            FileSource("blob-2d.step"), 2, order=mesh_order,
+            force_ambient_dim=2,
+            other_options=["-string", "Mesh.CharacteristicLengthMax = 0.02;"]
+            )
+
+    from meshmode.mesh.processing import merge_dijsoint_meshes, affine_map
+    mesh2 = affine_map(mesh, A=np.eye(2), b=np.array([5, 0]))
+
+    mesh3 = merge_dijsoint_meshes((mesh2, mesh))
+
+    if visualize:
+        from meshmode.discretization import Discretization
+        from meshmode.discretization.poly_element import \
+                PolynomialWarpAndBlendGroupFactory
+        cl_ctx = ctx_getter()
+        queue = cl.CommandQueue(cl_ctx)
+
+        discr = Discretization(cl_ctx, mesh3,
+                PolynomialWarpAndBlendGroupFactory(3))
+
+        from meshmode.discretization.visualization import make_visualizer
+        vis = make_visualizer(queue, discr, 1)
+        vis.write_vtk_file("merged.vtu", [])
+
+
 @pytest.mark.parametrize("dim", [2, 3])
 @pytest.mark.parametrize("order", [1, 3])
 def test_sanity_single_element(ctx_getter, dim, order, visualize=False):
@@ -220,7 +254,7 @@ def test_sanity_balls(ctx_getter, src_file, dim, mesh_order,
 
     from pytential import bind, sym
 
-    for h in [0.2, 0.15, 0.1]:
+    for h in [0.2, 0.14, 0.1]:
         from meshmode.mesh.io import generate_gmsh, FileSource
         mesh = generate_gmsh(
                 FileSource(src_file), dim, order=mesh_order,
@@ -239,7 +273,8 @@ def test_sanity_balls(ctx_getter, src_file, dim, mesh_order,
 
         from meshmode.discretization.connection import make_boundary_restriction
         bdry_mesh, bdry_discr, bdry_connection = make_boundary_restriction(
-                queue, vol_discr, InterpolatoryQuadratureSimplexGroupFactory(quad_order))
+                queue, vol_discr,
+                InterpolatoryQuadratureSimplexGroupFactory(quad_order))
 
         # }}}
 
@@ -264,7 +299,7 @@ def test_sanity_balls(ctx_getter, src_file, dim, mesh_order,
         comp_vol = integral(vol_discr, queue, vol_one)
         rel_vol_err = abs(true_vol - comp_vol) / true_vol
         vol_eoc_rec.add_data_point(h, rel_vol_err)
-        print "VOL", true_vol, comp_vol
+        print("VOL", true_vol, comp_vol)
 
         bdry_x = bdry_discr.nodes().with_queue(queue)
 
@@ -278,7 +313,7 @@ def test_sanity_balls(ctx_getter, src_file, dim, mesh_order,
         comp_surf = integral(bdry_discr, queue, bdry_one)
         rel_surf_err = abs(true_surf - comp_surf) / true_surf
         surf_eoc_rec.add_data_point(h, rel_surf_err)
-        print "SURF", true_surf, comp_surf
+        print("SURF", true_surf, comp_surf)
 
         if visualize:
             vol_vis.write_vtk_file("volume-h=%g.vtu" % h, [
@@ -297,19 +332,30 @@ def test_sanity_balls(ctx_getter, src_file, dim, mesh_order,
 
         # }}}
 
-    print "---------------------------------"
-    print "VOLUME"
-    print "---------------------------------"
-    print vol_eoc_rec
+    print("---------------------------------")
+    print("VOLUME")
+    print("---------------------------------")
+    print(vol_eoc_rec)
     assert vol_eoc_rec.order_estimate() >= mesh_order
 
-    print "---------------------------------"
-    print "SURFACE"
-    print "---------------------------------"
-    print surf_eoc_rec
+    print("---------------------------------")
+    print("SURFACE")
+    print("---------------------------------")
+    print(surf_eoc_rec)
     assert surf_eoc_rec.order_estimate() >= mesh_order
 
 
+def test_rect_mesh(do_plot=False):
+    from meshmode.mesh.generation import generate_regular_rect_mesh
+    mesh = generate_regular_rect_mesh()
+
+    if do_plot:
+        from meshmode.mesh.visualization import draw_2d_mesh
+        draw_2d_mesh(mesh, fill=None)
+        import matplotlib.pyplot as pt
+        pt.show()
+
+
 if __name__ == "__main__":
     import sys
     if len(sys.argv) > 1: