diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py
index 382860642476922378bb79c191a795b567c8851c..da0fdfe242536997dc3dc2f71c465f2eddf07560 100644
--- a/meshmode/mesh/__init__.py
+++ b/meshmode/mesh/__init__.py
@@ -490,6 +490,7 @@ class Mesh(Record):
     face_id_dtype = np.int8
 
     def __init__(self, vertices, groups, skip_tests=False,
+            node_vertex_consistency_tolerance=None,
             nodal_adjacency=False,
             facial_adjacency_groups=False,
             boundary_tags=None,
@@ -500,6 +501,9 @@ class Mesh(Record):
 
         :arg skip_tests: Skip mesh tests, in case you want to load a broken
             mesh anyhow and then fix it inside of this data structure.
+        :arg node_vertex_consistency_tolerance: If *False*, do not check
+            for consistency between vertex and nodal data. If *None*, use
+            the (small, near FP-epsilon) tolerance.
         :arg nodal_adjacency: One of three options:
             *None*, in which case this information
             will be deduced from vertex adjacency. *False*, in which case
@@ -576,7 +580,8 @@ class Mesh(Record):
                 )
 
         if not skip_tests:
-            assert _test_node_vertex_consistency(self)
+            assert _test_node_vertex_consistency(
+                    self, node_vertex_consistency_tolerance)
             for g in self.groups:
                 assert g.vertex_indices.dtype == self.vertex_id_dtype
 
@@ -691,7 +696,7 @@ class Mesh(Record):
 
 # {{{ node-vertex consistency test
 
-def _test_node_vertex_consistency_simplex(mesh, mgrp):
+def _test_node_vertex_consistency_simplex(mesh, mgrp, tol):
     if mgrp.nelements == 0:
         return True
 
@@ -710,7 +715,8 @@ def _test_node_vertex_consistency_simplex(mesh, mgrp):
             np.sum((map_vertices - grp_vertices)**2, axis=0),
             axis=-1))
 
-    tol = 1e3 * np.finfo(per_element_vertex_errors.dtype).eps
+    if tol is None:
+        tol = 1e3 * np.finfo(per_element_vertex_errors.dtype).eps
 
     from meshmode.mesh.processing import find_bounding_box
 
@@ -723,14 +729,17 @@ def _test_node_vertex_consistency_simplex(mesh, mgrp):
     return True
 
 
-def _test_node_vertex_consistency(mesh):
+def _test_node_vertex_consistency(mesh, tol):
     """Ensure that order of by-index vertices matches that of mapped
     unit vertices.
     """
 
+    if tol is False:
+        return
+
     for mgrp in mesh.groups:
         if isinstance(mgrp, SimplexElementGroup):
-            assert _test_node_vertex_consistency_simplex(mesh, mgrp)
+            assert _test_node_vertex_consistency_simplex(mesh, mgrp, tol)
         else:
             from warnings import warn
             warn("not implemented: node-vertex consistency check for '%s'"
diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index 65761586e01c8b0aa60f82d6f8b4bc00941eeed7..ce2a5f236d0ec524a411c78c87b7904abd423db6 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -175,7 +175,10 @@ def qbx_peanut(t):
 
 # {{{ make_curve_mesh
 
-def make_curve_mesh(curve_f, element_boundaries, order):
+def make_curve_mesh(curve_f, element_boundaries, order,
+        unit_nodes=None,
+        node_vertex_consistency_tolerance=None,
+        return_parametrization_points=False):
     """
     :arg curve_f: A callable representing a parametrization for a curve,
         accepting a vector of point locations and returning
@@ -183,15 +186,20 @@ def make_curve_mesh(curve_f, element_boundaries, order):
     :arg element_boundaries: a vector of element boundary locations in
         :math:`[0,1]`, in order. 0 must be the first entry, 1 the
         last one.
-    :returns: a :class:`meshmode.mesh.Mesh`
+    :arg unit_nodes: if given, the unit nodes to use. Must have shape
+        ``(dim, nnoodes)``.
+    :returns: a :class:`meshmode.mesh.Mesh`, or if *return_parametrization_points*
+        is True, a tuple ``(mesh, par_points)``, where *par_points* is an array of
+        parametrization points.
     """
 
     assert element_boundaries[0] == 0
     assert element_boundaries[-1] == 1
     nelements = len(element_boundaries) - 1
 
-    unodes = mp.warp_and_blend_nodes(1, order)
-    nodes_01 = 0.5*(unodes+1)
+    if unit_nodes is None:
+        unit_nodes = mp.warp_and_blend_nodes(1, order)
+    nodes_01 = 0.5*(unit_nodes+1)
 
     vertices = curve_f(element_boundaries)
 
@@ -200,7 +208,8 @@ def make_curve_mesh(curve_f, element_boundaries, order):
 
     # (el_nr, node_nr)
     t = el_starts[:, np.newaxis] + el_lengths[:, np.newaxis]*nodes_01
-    nodes = curve_f(t.ravel()).reshape(vertices.shape[0], nelements, -1)
+    t = t.ravel()
+    nodes = curve_f(t).reshape(vertices.shape[0], nelements, -1)
 
     from meshmode.mesh import Mesh, SimplexElementGroup
     egroup = SimplexElementGroup(
@@ -210,12 +219,18 @@ def make_curve_mesh(curve_f, element_boundaries, order):
                 np.arange(1, nelements+1, dtype=np.int32) % nelements,
                 ]).T,
             nodes=nodes,
-            unit_nodes=unodes)
+            unit_nodes=unit_nodes)
 
-    return Mesh(
+    mesh = Mesh(
             vertices=vertices, groups=[egroup],
             nodal_adjacency=None,
-            facial_adjacency_groups=None)
+            facial_adjacency_groups=None,
+            node_vertex_consistency_tolerance=node_vertex_consistency_tolerance)
+
+    if return_parametrization_points:
+        return mesh, t
+    else:
+        return mesh
 
 # }}}