diff --git a/examples/mesh-to-tikz.py b/examples/mesh-to-tikz.py
new file mode 100644
index 0000000000000000000000000000000000000000..2041e995d61da7f2c8ccade21250106a9ef5da6d
--- /dev/null
+++ b/examples/mesh-to-tikz.py
@@ -0,0 +1,14 @@
+from meshmode.mesh.io import generate_gmsh, FileSource
+
+h = 0.3
+order = 1
+
+mesh = generate_gmsh(
+        FileSource("../test/blob-2d.step"), 2, order=order,
+        force_ambient_dim=2,
+        other_options=[
+            "-string", "Mesh.CharacteristicLengthMax = %s;" % h]
+        )
+
+from meshmode.mesh.visualization import mesh_to_tikz
+print(mesh_to_tikz(mesh))
diff --git a/meshmode/mesh/visualization.py b/meshmode/mesh/visualization.py
index afee1b3419103909ac790b78d83d58348a104e45..abc586af42c754409612cfd47e3c2e7578c22189 100644
--- a/meshmode/mesh/visualization.py
+++ b/meshmode/mesh/visualization.py
@@ -29,6 +29,7 @@ __doc__ = """
 .. autofunction:: draw_2d_mesh
 .. autofunction:: draw_curve
 .. autofunction:: write_vertex_vtk_file
+.. autofunction:: mesh_to_tikz
 """
 
 
@@ -237,4 +238,44 @@ def write_vertex_vtk_file(mesh, file_name,
 
 # }}}
 
+
+# {{{ mesh_to_tikz
+
+def mesh_to_tikz(mesh):
+    lines = []
+
+    lines.append(r"\def\nelements{%s}" % (sum(grp.nelements for grp in mesh.groups)))
+    lines.append(r"\def\nvertices{%s}" % mesh.nvertices)
+    lines.append("")
+
+    drawel_lines = []
+    drawel_lines.append(r"\def\drawelements#1{")
+
+    for grp in mesh.groups:
+        for iel, el in enumerate(grp.vertex_indices):
+            el_nr = grp.element_nr_base+iel+1
+            elverts = mesh.vertices[:, el]
+
+            centroid = np.average(elverts, axis=1)
+            lines.append(r"\coordinate (cent%d) at (%s);"
+                    % (el_nr,
+                        ", ".join("%.5f" % (vi) for vi in centroid)))
+
+            for ivert, vert in enumerate(elverts.T):
+                lines.append(r"\coordinate (v%d-%d) at (%s);"
+                        % (el_nr, ivert+1, ", ".join("%.5f" % vi for vi in vert)))
+            drawel_lines.append(
+                    r"\draw [#1] %s -- cycle;"
+                    % " -- ".join(
+                        "(v%d-%d)" % (el_nr, vi+1)
+                        for vi in range(elverts.shape[1])))
+    drawel_lines.append("}")
+    lines.append("")
+
+    lines.extend(drawel_lines)
+
+    return "\n".join(lines)
+
+# }}}
+
 # vim: foldmethod=marker
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index e7ccf52afa6a321e2b5149ec55c725a05787930d..1ee97d94678e3f8c52fa1374b6e78de9e6541869 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -1002,6 +1002,8 @@ def test_quad_multi_element():
 # }}}
 
 
+# {{{ test_vtk_overwrite
+
 def test_vtk_overwrite(ctx_getter):
     pytest.importorskip("pyvisfile")
 
@@ -1049,6 +1051,28 @@ def test_vtk_overwrite(ctx_getter):
     _try_write_vtk(lambda x, y, **kwargs:
             write_nodal_adjacency_vtk_file(x, discr.mesh, **kwargs), discr.mesh)
 
+# }}}
+
+
+# {{{ test_mesh_to_tikz
+def test_mesh_to_tikz():
+    from meshmode.mesh.io import generate_gmsh, FileSource
+
+    h = 0.3
+    order = 1
+
+    mesh = generate_gmsh(
+            FileSource("../test/blob-2d.step"), 2, order=order,
+            force_ambient_dim=2,
+            other_options=[
+                "-string", "Mesh.CharacteristicLengthMax = %s;" % h]
+            )
+
+    from meshmode.mesh.visualization import mesh_to_tikz
+    mesh_to_tikz(mesh)
+
+# }}}
+
 
 if __name__ == "__main__":
     import sys