diff --git a/src/MeshData.jl b/src/MeshData.jl
index 0e368151..a37097c6 100644
--- a/src/MeshData.jl
+++ b/src/MeshData.jl
@@ -357,7 +357,9 @@ function MeshData(VX, VY, VZ, EToV, rd::RefElemData{3}; is_periodic=(false, fals
                   periodicity)
 
     if any(is_periodic)
-        md = make_periodic(md, is_periodic)
+        # loosen the tolerance if N >> 1
+        tol = length(rd.r) * 100 * eps() 
+        md = make_periodic(md, is_periodic; tol)
     end
 
     return md
diff --git a/src/RefElemData_polynomial.jl b/src/RefElemData_polynomial.jl
index a2a1fdfa..748f4a03 100644
--- a/src/RefElemData_polynomial.jl
+++ b/src/RefElemData_polynomial.jl
@@ -129,7 +129,8 @@ function RefElemData(elem::Union{Tet, Hex},
 
     # Construct matrices on reference elements
     r, s, t = nodes(elem, N)
-    Fmask = hcat(find_face_nodes(elem, r, s, t)...)
+    tol = 1e2 * eps() * length(r) # loosen the tolerance if N >> 1
+    Fmask = hcat(find_face_nodes(elem, r, s, t, tol)...)
     VDM, Vr, Vs, Vt = basis(elem, N, r, s, t)
     Dr, Ds, Dt = (A -> A / VDM).((Vr, Vs, Vt))
 
diff --git a/test/MeshData_tests.jl b/test/MeshData_tests.jl
index d939f90d..653dd1dc 100644
--- a/test/MeshData_tests.jl
+++ b/test/MeshData_tests.jl
@@ -289,3 +289,26 @@ approx_elem_types_to_test = [(Polynomial(), Hex()),
         
     end
 end
+
+@testset "Very high orders" begin
+    @testset "Tet" begin 
+        N = 12
+        rd = RefElemData(Tet(), N)
+        md = MeshData(uniform_mesh(Tet(), 2), rd; is_periodic=true)
+
+        (; x, y, z) = md
+        u = @. sin(pi * x) * sin(pi * y) * sin(pi * z)
+        dudx_exact = @. pi * cos(pi * x) * sin(pi * y) * sin(pi * z)
+        dudy_exact = @. pi * sin(pi * x) * cos(pi * y) * sin(pi * z)
+        dudz_exact = @. pi * sin(pi * x) * sin(pi * y) * cos(pi * z)
+
+        dudr, duds, dudt = (D -> D * u).(rd.Drst)
+        dudx = @. (md.rxJ * dudr + md.sxJ * duds + md.txJ * dudt) / md.J
+        dudy = @. (md.ryJ * dudr + md.syJ * duds + md.tyJ * dudt) / md.J
+        dudz = @. (md.rzJ * dudr + md.szJ * duds + md.tzJ * dudt) / md.J
+
+        @test norm(dudx - dudx_exact, Inf) < 1e-3
+        @test norm(dudy - dudy_exact, Inf) < 1e-3
+        @test norm(dudz - dudz_exact, Inf) < 1e-3
+    end
+end