diff --git a/assets/quad.geo b/assets/quad.geo
new file mode 100644
index 00000000..ff5ad2f6
--- /dev/null
+++ b/assets/quad.geo
@@ -0,0 +1,3 @@
+//+
+SetFactory("OpenCASCADE");
+Rectangle(1) = {0, 0, 0, 1, 1, 0};
diff --git a/assets/quad.msh b/assets/quad.msh
new file mode 100644
index 00000000..d9938461
--- /dev/null
+++ b/assets/quad.msh
@@ -0,0 +1,125 @@
+$MeshFormat
+4.1 0 8
+$EndMeshFormat
+$Entities
+4 4 1 0
+1 0 0 0 0 
+2 1 0 0 0 
+3 1 1 0 0 
+4 0 1 0 0 
+1 -9.999999994736442e-08 -1e-07 -1e-07 1.0000001 1e-07 1e-07 0 2 1 -2 
+2 0.9999999000000001 -9.999999994736442e-08 -1e-07 1.0000001 1.0000001 1e-07 0 2 2 -3 
+3 -9.999999994736442e-08 0.9999999000000001 -1e-07 1.0000001 1.0000001 1e-07 0 2 3 -4 
+4 -1e-07 -9.999999994736442e-08 -1e-07 1e-07 1.0000001 1e-07 0 2 4 -1 
+1 -9.999999994736442e-08 -9.999999994736442e-08 -1e-07 1.0000001 1.0000001 1e-07 0 4 1 2 3 4 
+$EndEntities
+$Nodes
+9 25 1 25
+0 1 0 1
+1
+0 0 0
+0 2 0 1
+2
+1 0 0
+0 3 0 1
+3
+1 1 0
+0 4 0 1
+4
+0 1 0
+1 1 0 3
+5
+6
+7
+0.2499999999999998 0 0
+0.4999999999999988 0 0
+0.7499999999999991 0 0
+1 2 0 3
+8
+9
+10
+1 0.2499999999999998 0
+1 0.4999999999999988 0
+1 0.7499999999999991 0
+1 3 0 3
+11
+12
+13
+0.7500000000000002 1 0
+0.5000000000000011 1 0
+0.2500000000000009 1 0
+1 4 0 3
+14
+15
+16
+0 0.7500000000000002 0
+0 0.5000000000000011 0
+0 0.2500000000000009 0
+2 1 0 9
+17
+18
+19
+20
+21
+22
+23
+24
+25
+0.7499999999999997 0.2499999999999997 0
+0.2500000000000001 0.5000000000000002 0
+0.5000000000000004 0.7499999999999999 0
+0.4999999999999997 0.2499999999999998 0
+0.2500000000000003 0.7500000000000001 0
+0.2499999999999998 0.2500000000000002 0
+0.75 0.4999999999999995 0
+0.7500000000000002 0.7499999999999996 0
+0.5000000000000001 0.4999999999999998 0
+$EndNodes
+$Elements
+9 36 1 36
+0 1 15 1
+1 1 
+0 2 15 1
+2 2 
+0 3 15 1
+3 3 
+0 4 15 1
+4 4 
+1 1 1 4
+5 1 5 
+6 5 6 
+7 6 7 
+8 7 2 
+1 2 1 4
+9 2 8 
+10 8 9 
+11 9 10 
+12 10 3 
+1 3 1 4
+13 3 11 
+14 11 12 
+15 12 13 
+16 13 4 
+1 4 1 4
+17 4 14 
+18 14 15 
+19 15 16 
+20 16 1 
+2 1 3 16
+21 20 25 18 22 
+22 25 19 21 18 
+23 25 23 24 19 
+24 2 8 17 7 
+25 24 23 9 10 
+26 21 19 12 13 
+27 22 18 15 16 
+28 7 17 20 6 
+29 17 8 9 23 
+30 23 25 20 17 
+31 14 15 18 21 
+32 16 1 5 22 
+33 13 4 14 21 
+34 10 3 11 24 
+35 22 5 6 20 
+36 24 11 12 19 
+$EndElements
diff --git a/src/geometry.jl b/src/geometry.jl
index 67d2ac64..22fe1172 100644
--- a/src/geometry.jl
+++ b/src/geometry.jl
@@ -325,7 +325,13 @@ end
 function node_coordinates(fe::AbstractLagrangeMeshFace)
     @assert fe |> geometry |> is_unitary
     mexps = monomial_exponents(fe)
-    node_coordinates_from_monomials_exponents(mexps,fe.order_per_dir,fe.geometry |> real_type)
+    lib_node_to_coords = node_coordinates_from_monomials_exponents(mexps,fe.order_per_dir,fe.geometry |> real_type)
+    if length(fe.lib_to_user_nodes) == 0
+        return lib_node_to_coords
+    end
+    user_node_to_coords = similar(lib_node_to_coords)
+    user_node_to_coords[fe.lib_to_user_nodes] = lib_node_to_coords
+    user_node_to_coords
 end
 
 function primal_basis(fe::AbstractLagrangeMeshFace)
diff --git a/test/geometry_tests.jl b/test/geometry_tests.jl
index 440017ee..21908667 100644
--- a/test/geometry_tests.jl
+++ b/test/geometry_tests.jl
@@ -67,6 +67,18 @@ x = gk.node_coordinates(fe)
 A = tabulator(gk.value,x)
 
 outdir = mkpath(joinpath(@__DIR__,"..","output"))
+msh =  joinpath(@__DIR__,"..","assets","quad.msh")
+mesh = gk.mesh_from_gmsh(msh)
+vtk_grid(joinpath(outdir,"quad"),gk.vtk_args(mesh)...) do vtk
+    gk.vtk_physical_faces!(vtk,mesh)
+    gk.vtk_physical_nodes!(vtk,mesh)
+end
+for d in 0:gk.num_dims(mesh)
+    vtk_grid(joinpath(outdir,"quad_$d"),gk.vtk_args(mesh,d)...) do vtk
+        gk.vtk_physical_faces!(vtk,mesh,d)
+        gk.vtk_physical_nodes!(vtk,mesh,d)
+    end
+end
 
 msh =  joinpath(@__DIR__,"..","assets","demo.msh")
 mesh = gk.mesh_from_gmsh(msh;complexify=false)