diff --git a/src/RefElemData_SBP.jl b/src/RefElemData_SBP.jl
index 14fba196..c3e1a848 100644
--- a/src/RefElemData_SBP.jl
+++ b/src/RefElemData_SBP.jl
@@ -26,6 +26,20 @@ function RefElemData(elementType::Union{Quad, Hex}, approxType::SBP{TensorProduc
                      Polynomial(TensorProductQuadrature(gauss_lobatto_quad(0, 0, N))), 
                      N; kwargs...)
 
+    # brute-force determine Fmask so that rd.rf = rd.r[rd.Fmask], etc.
+    new_Fmask = copy(rd.Fmask)    
+    for fid in eachindex(rd.rf)
+        rstf = SVector(getindex.(rd.rstf, fid))
+        for i in eachindex(rd.r)
+            rst = SVector(getindex.(rd.rst, i))
+            if rstf ≈ rst
+                new_Fmask[fid] = i
+                break
+            end
+        end
+    end
+    rd = @set rd.Fmask = new_Fmask;                     
+
     rd = @set rd.Vf = droptol!(sparse(rd.Vf), tol)
     rd = @set rd.LIFT = Diagonal(rd.wq) \ (rd.Vf' * Diagonal(rd.wf)) # TODO: make this more efficient with LinearMaps?
 
diff --git a/test/multidim_sbp_tests.jl b/test/multidim_sbp_tests.jl
index 2ff50614..d8b2bdee 100644
--- a/test/multidim_sbp_tests.jl
+++ b/test/multidim_sbp_tests.jl
@@ -47,6 +47,9 @@
     @testset "TensorProductLobatto Hex" begin
         rd = RefElemData(Hex(),SBP(),N)
         @test propertynames(rd)[1] == :element_type
+        @test rd.rf == rd.r[rd.Fmask]
+        @test rd.sf == rd.s[rd.Fmask]
+        @test rd.tf == rd.t[rd.Fmask]
         @test rd.t == rd.rst[3]
         @test rd.tf == rd.rstf[3]    
         @test rd.tq == rd.rstq[3]