diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index 394f0cb592..c2733f48b5 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -2979,14 +2979,14 @@ def refine_marked_elements(self, mark, netgen_flags=None):
                     netgen_flags=netgen_flags)
 
     @PETSc.Log.EventDecorator()
-    def curve_field(self, order, permutation_tol=1e-8, cg_field=None):
+    def curve_field(self, order, permutation_tol=None, cg_field=None):
         '''Return a function containing the curved coordinates of the mesh.
 
         This method requires that the mesh has been constructed from a
         netgen mesh.
 
         :arg order: the order of the curved mesh.
-        :arg permutation_tol: tolerance used to construct the permutation of the reference element.
+        :arg permutation_tol: ignored.
         :arg cg_field: return a CG function field representing the mesh, as opposed to a DG field.
             Defaults to the continuity of the coordinates of the original mesh.
 
@@ -2998,6 +2998,12 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None):
 
         if not hasattr(self, "netgen_mesh"):
             raise ValueError("Cannot curve a mesh that has not been generated by netgen.")
+        if permutation_tol is not None:
+            warnings.warn(
+                "permutation_tol is no longer required to obtain the curved coordinates. "
+                "This kwarg will be removed in a future release.",
+                FutureWarning,
+            )
 
         if cg_field is None:
             cg_field = not self.coordinates.function_space().finat_element.is_dg()
@@ -3020,10 +3026,13 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None):
         fiat_element = new_coordinates.function_space().finat_element.fiat_equivalent
         nodes = fiat_element.dual_basis()
         ref_pts = []
-        for node in nodes:
-            # Assert singleton point for each node.
-            pt, = node.get_point_dict().keys()
-            ref_pts.append(pt)
+        entity_ids = fiat_element.entity_dofs()
+        for dim in sorted(entity_ids):
+            for entity in sorted(entity_ids[dim]):
+                for i in entity_ids[dim][entity]:
+                    # Assert singleton point for each node.
+                    pt, = nodes[i].get_point_dict().keys()
+                    ref_pts.append(pt)
         reference_points = np.array(ref_pts)
 
         # Construct numpy arrays for physical domain data
@@ -3033,8 +3042,8 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None):
         curved_points = np.zeros(
             (ng_dimension, reference_points.shape[0], self.geometric_dimension)
         )
+        self.netgen_mesh.Curve(1)
         self.netgen_mesh.CalcElementMapping(reference_points, physical_points)
-        # NOTE: This will segfault for MeshHierarchy on a netgen CSG geometry
         self.netgen_mesh.Curve(order)
         self.netgen_mesh.CalcElementMapping(reference_points, curved_points)
         curved = ng_element.NumPy()["curved"]
@@ -3059,7 +3068,6 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None):
         permutation = find_permutation(
             own_physical_points,
             new_coordinates.dat.data_ro_with_halos[broken_indices].real,
-            tol=permutation_tol,
         )
         self.comm.Barrier()
         # Apply the permutation to each cell in turn
@@ -3425,7 +3433,7 @@ def Mesh(meshfile, **kwargs):
         # Curve the mesh, if requested
         degree = netgen_flags.get("degree", 1)
         if degree != 1:
-            permutation_tol = netgen_flags.get("permutation_tol", 1e-8)
+            permutation_tol = netgen_flags.get("permutation_tol", None)
             cg = netgen_flags.get("cg", None)
             coordinates = mesh.curve_field(
                 order=degree,
diff --git a/firedrake/mg/netgen.py b/firedrake/mg/netgen.py
index 83f98bbbef..f6088528ed 100644
--- a/firedrake/mg/netgen.py
+++ b/firedrake/mg/netgen.py
@@ -242,7 +242,7 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None):
     order = flags.get("degree", 1)
     if isinstance(order, int):
         order = [order]*(levs+1)
-    permutation_tol = flags.get("permutation_tol", 1e-8)
+    permutation_tol = flags.get("permutation_tol", None)
     refType = flags.get("refinement_type", "uniform")
     optMoves = flags.get("optimisation_moves", False)
     snap = flags.get("snap_to", "geometry")
diff --git a/firedrake/netgen.py b/firedrake/netgen.py
index 57b442b83a..e44af74140 100644
--- a/firedrake/netgen.py
+++ b/firedrake/netgen.py
@@ -80,8 +80,7 @@ def netgen_distribute(V: firedrake.functionspaceimpl.WithGeometryBase,
 
 
 @PETSc.Log.EventDecorator()
-def find_permutation(points_a: np.ndarray, points_b: np.ndarray,
-                     tol: float = 1e-5):
+def find_permutation(points_a: np.ndarray, points_b: np.ndarray):
     """ Find all permutations between a list of two sets of points.
 
     Given two numpy arrays of shape (ncells, npoints, dim) containing
@@ -95,7 +94,20 @@ def find_permutation(points_a: np.ndarray, points_b: np.ndarray,
     if points_a.shape != points_b.shape:
         raise ValueError("`points_a` and `points_b` must have the same shape.")
 
-    p = [np.where(cdist(a, b).T < tol)[1] for a, b in zip(points_a, points_b)]
+    # Match reference points instead of physical points to ensure scale invariance
+    dim = points_a.shape[-1]
+    vids = list(range(dim+1))
+    # Infer the affine mapping (A, b) from the image of the vertices (first dim+1 dofs)
+    bs = points_a[:, vids[:1], :]
+    As = points_a[:, vids[1:], :]
+    As -= bs
+    Ainvs = np.linalg.inv(As)
+    # x_phys = A * x_ref + b <==> x_ref = inv(A) * (x_phys - b)
+    # Multiply inv(A) from the right, since the data is row-major
+    ref_points_a = np.matmul(points_a - bs, Ainvs)
+    ref_points_b = np.matmul(points_b - bs, Ainvs)
+
+    p = [np.argmin(cdist(a, b), axis=0) for a, b in zip(ref_points_a, ref_points_b)]
 
     if len(p) == 0:
         return p
diff --git a/tests/firedrake/regression/test_netgen.py b/tests/firedrake/regression/test_netgen.py
index 8524b52748..64abf21c31 100644
--- a/tests/firedrake/regression/test_netgen.py
+++ b/tests/firedrake/regression/test_netgen.py
@@ -34,14 +34,34 @@ def test_netgen_csg_mesh_high_order():
     assert mesh3.coordinates.function_space().ufl_element().degree() == order
 
 
-def square_geometry(h):
+def square_geometry(h, L=np.pi):
     from netgen.geom2d import SplineGeometry
     geo = SplineGeometry()
-    geo.AddRectangle((0, 0), (np.pi, np.pi), bc="rect")
+    geo.AddRectangle((0, 0), (L, L), bc="rect")
     ngmesh = geo.GenerateMesh(maxh=h)
     return ngmesh
 
 
+def circle_geometry(h, R=1.0):
+    from netgen.geom2d import SplineGeometry
+    geo = SplineGeometry()
+    geo.AddCircle((0, 0), R, bc="circ")
+    ngmesh = geo.GenerateMesh(maxh=h)
+    return ngmesh
+
+
+@pytest.mark.parametrize("scale", (1E-5, 1E5))
+def test_high_order(scale):
+    # Test scale independence of high-order geometry
+    expected = np.pi * scale * scale
+    ngmesh = circle_geometry(h=scale/4, R=scale)
+
+    degree = 3
+    msh = Mesh(ngmesh, netgen_flags={"degree": degree})
+    assert msh.coordinates.function_space().ufl_element().degree() == degree
+    assert np.isclose(assemble(1*dx(domain=msh)), expected)
+
+
 def poisson(h, degree=2):
     import netgen
     comm = COMM_WORLD