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/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/__init__.py b/meshmode/mesh/__init__.py
index bef25c2ef2d66f0b08ab79249d83e4056b7b12e0..c9917052ded9861366ffb63751b54db94e6a98a7 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"
 
@@ -407,18 +409,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..55a49387bf1d97c19bc77bf4f3092c143910df46 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"
 
@@ -300,9 +302,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 +342,8 @@ 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):
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..dfa46078b101d65b59459d97a19605d213171fee 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"
 
@@ -64,8 +67,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]
 
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 5febf0308759c11b2bb6b7c8596bf623b2329d56..3c8f20d29a57c42863138962cf95a903258715e1 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
 
 
@@ -264,7 +267,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 +281,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,16 +300,16 @@ 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