diff --git a/Project.toml b/Project.toml
index 6886d8fd..8c73ec3e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -23,6 +23,12 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
+[weakdeps]
+SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240"
+
+[extensions]
+StartUpDGSummationByPartsOperatorsExt = "SummationByPartsOperators"
+
 [compat]
 ConstructionBase = "1"
 FillArrays = "0.13, 1"
@@ -37,6 +43,10 @@ Requires = "1"
 Setfield = "1"
 SimpleUnPack = "1"
 StaticArrays = "1"
+SummationByPartsOperators = "0.5"
 Triangulate = "2"
 WriteVTK = "1"
 julia = "1.7"
+
+[extras]
+SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240"
diff --git a/docs/src/RefElemData.md b/docs/src/RefElemData.md
index dfda3b96..ab997da2 100644
--- a/docs/src/RefElemData.md
+++ b/docs/src/RefElemData.md
@@ -72,9 +72,22 @@ By default, `RefElemData` is constructed for a nodal basis (in order to facilita
 
 ## `RefElemData` based on SBP finite differences
 
-It is also possible to construct a [`RefElemData`](@ref) based on [multi-dimensional SBP finite difference operators](https://doi.org/10.1137/15M1038360). These utilize nodes constructed by [Tianheng Chen and Chi-Wang Shu](https://doi.org/10.1016/j.jcp.2017.05.025), [Ethan Kubatko](https://sites.google.com/site/chilatosu/ethan-bio), and [Jason Hicken](https://doi.org/10.1007/s10915-020-01154-8).
+It is also possible to construct a [`RefElemData`](@ref) based on both traditional finite difference SBP operators and [multi-dimensional SBP finite difference operators](https://doi.org/10.1137/15M1038360). The traditional finite difference SBP operators are built using [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl). The multi-dimensional SBP operators utilize nodes constructed by [Tianheng Chen and Chi-Wang Shu](https://doi.org/10.1016/j.jcp.2017.05.025), [Ethan Kubatko](https://sites.google.com/site/chilatosu/ethan-bio), and [Jason Hicken](https://doi.org/10.1007/s10915-020-01154-8). 
 
-Some examples:
+Some examples of traditional finite difference SBP operators on tensor product domains:
+```julia
+using StartUpDG
+using SummationByPartsOperators # this package must be loaded to enable extension
+
+D = derivative_operator(MattssonNordström2004(), 
+                        derivative_order=1, accuracy_order=N+1,
+                        xmin=-1.0, xmax=1.0, N=20)
+
+rd = RefElemData(Line(), D) 
+rd = RefElemData(Quad(), D)
+rd = RefElemData(Hex(), D)
+```
+Some examples of multi-dimensional SBP operators:
 ```julia
 N = 3
 rd = RefElemData(Quad(), SBP(), N) # defaults to SBP{TensorProductLobatto}
diff --git a/ext/StartUpDGSummationByPartsOperatorsExt.jl b/ext/StartUpDGSummationByPartsOperatorsExt.jl
new file mode 100644
index 00000000..e036cda8
--- /dev/null
+++ b/ext/StartUpDGSummationByPartsOperatorsExt.jl
@@ -0,0 +1,413 @@
+module StartUpDGSummationByPartsOperatorsExt
+
+using LinearAlgebra: LinearAlgebra, Diagonal, diag, norm, UniformScaling
+using SparseArrays: sparse, droptol!, spzeros
+
+using StartUpDG
+
+# Required for visualization code
+if isdefined(Base, :get_extension)
+    using SummationByPartsOperators:
+        SummationByPartsOperators,
+        DerivativeOperator,
+        grid,
+        AbstractDerivativeOperator,
+        AbstractNonperiodicDerivativeOperator,
+        PeriodicDerivativeOperator,
+        AbstractPeriodicDerivativeOperator
+else
+    # Until Julia v1.9 is the minimum required version for Trixi.jl, we still support Requires.jl
+    using ..SummationByPartsOperators
+    using ..SummationByPartsOperators:
+        AbstractDerivativeOperator,
+        AbstractPeriodicDerivativeOperator,
+        AbstractNonperiodicDerivativeOperator,
+        DerivativeOperator,
+        PeriodicDerivativeOperator,
+        grid
+end
+
+function construct_1d_operators(D::AbstractDerivativeOperator, tol)
+    nodes_1d = collect(grid(D))
+    M = SummationByPartsOperators.mass_matrix(D)
+    if M isa UniformScaling
+        weights_1d = M * ones(Bool, length(nodes_1d))
+    else
+        weights_1d = diag(M)
+    end
+
+    # StartUpDG assumes nodes from -1 to +1. Thus, we need to re-scale everything.
+    # We can adjust the grid spacing as follows.
+    xmin = SummationByPartsOperators.xmin(D)
+    xmax = SummationByPartsOperators.xmax(D)
+    factor = 2 / (xmax - xmin)
+    @. nodes_1d = factor * (nodes_1d - xmin) - 1
+    @. weights_1d = factor * weights_1d
+
+    D_1d = droptol!(inv(factor) * sparse(D), tol)
+    I_1d = Diagonal(ones(Bool, length(nodes_1d)))
+
+    return nodes_1d, weights_1d, D_1d, I_1d
+end
+
+function StartUpDG.RefElemData(element_type::Line,
+    D::AbstractDerivativeOperator;
+    tol = 100 * eps(),)
+    approximation_type = D
+    N = SummationByPartsOperators.accuracy_order(D) # kind of polynomial degree
+
+    # 1D operators
+    nodes_1d, weights_1d, D_1d = construct_1d_operators(D, tol)
+
+    # volume
+    rq = r = nodes_1d
+    wq = weights_1d
+    Dr = D_1d
+    M = Diagonal(wq)
+    Pq = LinearAlgebra.I
+    Vq = LinearAlgebra.I
+
+    VDM = nothing # unused generalized Vandermonde matrix
+
+    rst = (r,)
+    rstq = (rq,)
+    Drst = (Dr,)
+
+    # face
+    face_vertices = StartUpDG.face_vertices(element_type)
+    face_mask = [1, length(nodes_1d)]
+
+    rf = [-1.0; 1.0]
+    nrJ = [-1.0; 1.0]
+    wf = [1.0; 1.0]
+    if D isa AbstractPeriodicDerivativeOperator
+        # we do not need any face stuff for periodic operators
+        Vf = spzeros(length(wf), length(wq))
+    else
+        Vf = sparse([1, 2], [1, length(nodes_1d)], [1.0, 1.0])
+    end
+    LIFT = Diagonal(wq) \ (Vf' * Diagonal(wf))
+
+    rstf = (rf,)
+    nrstJ = (nrJ,)
+
+    # low order interpolation nodes
+    r1 = StartUpDG.nodes(element_type, 1)
+    V1 = StartUpDG.vandermonde(element_type, 1, r) /
+         StartUpDG.vandermonde(element_type, 1, r1)
+
+    return RefElemData(element_type,
+        approximation_type,
+        N,
+        face_vertices,
+        V1,
+        rst,
+        VDM,
+        face_mask,
+        rst,
+        LinearAlgebra.I, # plotting
+        rstq,
+        wq,
+        Vq, # quadrature
+        rstf,
+        wf,
+        Vf,
+        nrstJ, # faces
+        M,
+        Pq,
+        Drst,
+        LIFT)
+end
+
+function StartUpDG.RefElemData(element_type::Quad,
+    D::AbstractDerivativeOperator;
+    tol = 100 * eps(),)
+    approximation_type = D
+    N = SummationByPartsOperators.accuracy_order(D) # kind of polynomial degree
+
+    # 1D operators
+    nodes_1d, weights_1d, D_1d, I_1d = construct_1d_operators(D, tol)
+
+    # volume
+    s, r = vec.(StartUpDG.NodesAndModes.meshgrid(nodes_1d)) # this is to match
+    # ordering of nrstJ
+    rq = r
+    sq = s
+    wr, ws = vec.(StartUpDG.NodesAndModes.meshgrid(weights_1d))
+    wq = wr .* ws
+    Dr = kron(I_1d, D_1d)
+    Ds = kron(D_1d, I_1d)
+    M = Diagonal(wq)
+    Pq = LinearAlgebra.I
+    Vq = LinearAlgebra.I
+
+    VDM = nothing # unused generalized Vandermonde matrix
+
+    rst = (r, s)
+    rstq = (rq, sq)
+    Drst = (Dr, Ds)
+
+    # face
+    face_vertices = StartUpDG.face_vertices(element_type)
+    face_mask = vcat(StartUpDG.find_face_nodes(element_type, r, s)...)
+
+    rf, sf, wf, nrJ, nsJ = StartUpDG.init_face_data(element_type,
+        quad_rule_face = (nodes_1d, weights_1d))
+    if D isa AbstractPeriodicDerivativeOperator
+        # we do not need any face stuff for periodic operators
+        Vf = spzeros(length(wf), length(wq))
+    else
+        Vf = sparse(eachindex(face_mask), face_mask, ones(Bool, length(face_mask)))
+    end
+    LIFT = Diagonal(wq) \ (Vf' * Diagonal(wf))
+
+    rstf = (rf, sf)
+    nrstJ = (nrJ, nsJ)
+
+    # low order interpolation nodes
+    r1, s1 = StartUpDG.nodes(element_type, 1)
+    V1 = StartUpDG.vandermonde(element_type, 1, r, s) /
+         StartUpDG.vandermonde(element_type, 1, r1, s1)
+
+    return RefElemData(element_type,
+        approximation_type,
+        N,
+        face_vertices,
+        V1,
+        rst,
+        VDM,
+        face_mask,
+        rst,
+        LinearAlgebra.I, # plotting
+        rstq,
+        wq,
+        Vq, # quadrature
+        rstf,
+        wf,
+        Vf,
+        nrstJ, # faces
+        M,
+        Pq,
+        Drst,
+        LIFT)
+end
+
+function StartUpDG.RefElemData(element_type::Hex,
+    D::AbstractDerivativeOperator;
+    tol = 100 * eps(),)
+    approximation_type = D
+    N = SummationByPartsOperators.accuracy_order(D) # kind of polynomial degree
+
+    # 1D operators
+    nodes_1d, weights_1d, D_1d, I_1d = construct_1d_operators(D, tol)
+
+    # volume
+    # to match ordering of nrstJ
+    s, r, t = vec.(StartUpDG.NodesAndModes.meshgrid(nodes_1d, nodes_1d, nodes_1d))
+    rq = r
+    sq = s
+    tq = t
+    wr, ws, wt = vec.(StartUpDG.NodesAndModes.meshgrid(weights_1d, weights_1d, weights_1d))
+    wq = wr .* ws .* wt
+    Dr = kron(I_1d, I_1d, D_1d)
+    Ds = kron(I_1d, D_1d, I_1d)
+    Dt = kron(D_1d, I_1d, I_1d)
+    M = Diagonal(wq)
+    Pq = LinearAlgebra.I
+    Vq = LinearAlgebra.I
+
+    VDM = nothing # unused generalized Vandermonde matrix
+
+    rst = (r, s, t)
+    rstq = (rq, sq, tq)
+    Drst = (Dr, Ds, Dt)
+
+    # face
+    face_vertices = StartUpDG.face_vertices(element_type)
+    face_mask = vcat(StartUpDG.find_face_nodes(element_type, r, s, t)...)
+
+    rf, sf, tf, wf, nrJ, nsJ, ntJ = let
+        rf, sf = vec.(StartUpDG.NodesAndModes.meshgrid(nodes_1d, nodes_1d))
+        wr, ws = vec.(StartUpDG.NodesAndModes.meshgrid(weights_1d, weights_1d))
+        wf = wr .* ws
+        StartUpDG.init_face_data(element_type, quad_rule_face = (rf, sf, wf))
+    end
+    Vf = sparse(eachindex(face_mask), face_mask, ones(Bool, length(face_mask)))
+    LIFT = Diagonal(wq) \ (Vf' * Diagonal(wf))
+
+    rstf = (rf, sf, tf)
+    nrstJ = (nrJ, nsJ, ntJ)
+
+    # low order interpolation nodes
+    r1, s1, t1 = StartUpDG.nodes(element_type, 1)
+    V1 = StartUpDG.vandermonde(element_type, 1, r, s, t) /
+         StartUpDG.vandermonde(element_type, 1, r1, s1, t1)
+
+    return RefElemData(element_type,
+        approximation_type,
+        N,
+        face_vertices,
+        V1,
+        rst,
+        VDM,
+        face_mask,
+        rst,
+        LinearAlgebra.I, # plotting
+        rstq,
+        wq,
+        Vq, # quadrature
+        rstf,
+        wf,
+        Vf,
+        nrstJ, # faces
+        M,
+        Pq,
+        Drst,
+        LIFT)
+end
+
+# specialized Hex constructor in 3D to reduce memory usage.
+function StartUpDG.RefElemData(element_type::Hex,
+    D::AbstractPeriodicDerivativeOperator;
+    tol = 100 * eps(),)
+    approximation_type = D
+    N = SummationByPartsOperators.accuracy_order(D) # kind of polynomial degree
+
+    # 1D operators
+    nodes_1d, weights_1d, D_1d, I_1d = construct_1d_operators(D, tol)
+
+    # volume
+    # to match ordering of nrstJ
+    s, r, t = vec.(StartUpDG.NodesAndModes.meshgrid(nodes_1d, nodes_1d, nodes_1d))
+    rq = r
+    sq = s
+    tq = t
+    wr, ws, wt = vec.(StartUpDG.NodesAndModes.meshgrid(weights_1d, weights_1d, weights_1d))
+    wq = wr .* ws .* wt
+    Dr = kron(I_1d, I_1d, D_1d)
+    Ds = kron(I_1d, D_1d, I_1d)
+    Dt = kron(D_1d, I_1d, I_1d)
+    M = Diagonal(wq)
+    Pq = LinearAlgebra.I
+    Vq = LinearAlgebra.I
+
+    VDM = nothing # unused generalized Vandermonde matrix
+
+    rst = (r, s, t)
+    rstq = (rq, sq, tq)
+    Drst = (Dr, Ds, Dt)
+
+    # face
+    # We do not need any face data for periodic operators. Thus, we just
+    # pass `nothing` to save memory.
+    face_vertices = ntuple(_ -> nothing, 3)
+    face_mask = nothing
+    wf = nothing
+    rstf = ntuple(_ -> nothing, 3)
+    nrstJ = ntuple(_ -> nothing, 3)
+    Vf = nothing
+    LIFT = nothing
+
+    # low order interpolation nodes
+    V1 = nothing # do not need to store V1, since we specialize StartUpDG.MeshData to avoid using it.
+
+    return RefElemData(element_type,
+        approximation_type,
+        N,
+        face_vertices,
+        V1,
+        rst,
+        VDM,
+        face_mask,
+        rst,
+        LinearAlgebra.I, # plotting
+        rstq,
+        wq,
+        Vq, # quadrature
+        rstf,
+        wf,
+        Vf,
+        nrstJ, # faces
+        M,
+        Pq,
+        Drst,
+        LIFT)
+end
+
+function Base.show(io::IO,
+    mime::MIME"text/plain",
+    rd::RefElemData{NDIMS, ElementType, ApproximationType}) where {
+    NDIMS,
+    ElementType <: StartUpDG.AbstractElemShape,
+    ApproximationType <: AbstractDerivativeOperator,
+}
+    @nospecialize rd
+    print(io, "RefElemData for an approximation using an ")
+    show(IOContext(io, :compact => true), rd.approximation_type)
+    print(io, " on $(rd.element_type) element")
+end
+
+function Base.show(io::IO,
+    rd::RefElemData{NDIMS, ElementType, ApproximationType}) where {
+    NDIMS,
+    ElementType <: StartUpDG.AbstractElemShape,
+    ApproximationType <: AbstractDerivativeOperator,
+}
+    @nospecialize rd
+    print(io, "RefElemData{", summary(rd.approximation_type), ", ", rd.element_type, "}")
+end
+
+function StartUpDG.inverse_trace_constant(rd::RefElemData{
+    NDIMS,
+    ElementType,
+    ApproximationType,
+}) where {
+    NDIMS,
+    ElementType <: Union{Line, Quad, Hex},
+    ApproximationType <: AbstractDerivativeOperator,
+}
+    D = rd.approximation_type
+
+    # the inverse trace constant is the maximum eigenvalue corresponding to
+    #       M_f * v = λ * M * v
+    # where M_f is the face mass matrix and M is the volume mass matrix.
+    # Since M is diagonal and since M_f is just the boundary "mask" matrix
+    # (which extracts the first and last entries of a vector), the maximum
+    # eigenvalue is the inverse of the first or last mass matrix diagonal.
+    left_weight = SummationByPartsOperators.left_boundary_weight(D)
+    right_weight = SummationByPartsOperators.right_boundary_weight(D)
+    max_eigenvalue = max(inv(left_weight), inv(right_weight))
+
+    # For tensor product elements, the trace constant for higher dimensional
+    # elements is the one-dimensional trace constant multiplied by `NDIMS`. See
+    #     "GPU-accelerated discontinuous Galerkin methods on hybrid meshes."
+    #     Chan, Jesse, et al (2016), https://doi.org/10.1016/j.jcp.2016.04.003
+    # for more details (specifically, Appendix A.1, Theorem A.4).
+    return NDIMS * max_eigenvalue
+end
+
+# This is used in `estimate_dt`. `estimate_h` uses that `Jf / J = O(h^{NDIMS-1}) / O(h^{NDIMS}) = O(1/h)`.
+# However, since we do not initialize `Jf` for periodic FDSBP operators, we specialize `estimate_h`
+# based on the reference grid provided by SummationByPartsOperators.jl and information about the domain size
+# provided by `md::MeshData``.
+function StartUpDG.estimate_h(e,
+    rd::RefElemData{NDIMS, ElementType, ApproximationType},
+    md::MeshData) where {
+    NDIMS,
+    ElementType <: StartUpDG.AbstractElemShape,
+    ApproximationType <: SummationByPartsOperators.AbstractPeriodicDerivativeOperator,
+}
+    D = rd.approximation_type
+    x = grid(D)
+
+    # we assume all SummationByPartsOperators.jl reference grids are rescaled to [-1, 1]
+    xmin = SummationByPartsOperators.xmin(D)
+    xmax = SummationByPartsOperators.xmax(D)
+    factor = 2 / (xmax - xmin)
+
+    # If the domain has size L^NDIMS, then `minimum(md.J)^(1 / NDIMS) = L`.
+    # WARNING: this is not a good estimate on anisotropic grids.
+    return minimum(diff(x)) * factor * minimum(md.J)^(1 / NDIMS)
+end
+
+end
diff --git a/src/StartUpDG.jl b/src/StartUpDG.jl
index 3e1db386..fb98506e 100644
--- a/src/StartUpDG.jl
+++ b/src/StartUpDG.jl
@@ -1,4 +1,4 @@
-module StartUpDG 
+module StartUpDG
 
 using Reexport: @reexport
 
@@ -6,7 +6,8 @@ using ConstructionBase: ConstructionBase
 using FillArrays: Ones, Zeros, Fill
 using HDF5: h5open # used to read in SBP triangular node data
 using Kronecker: kronecker # for Hex element matrix manipulations
-using LinearAlgebra: cond, diagm, eigvals, Diagonal, UniformScaling, I, mul!, norm, qr, ColumnNorm, Symmetric
+using LinearAlgebra:
+    cond, diagm, eigvals, Diagonal, UniformScaling, I, mul!, norm, qr, ColumnNorm, Symmetric
 using NodesAndModes: meshgrid, find_face_nodes, face_vertices
 @reexport using NodesAndModes # for basis functions
 using PathIntersections: PathIntersections
@@ -77,8 +78,8 @@ export uniform_mesh
 include("mesh/gmsh_utilities.jl")
 export read_Gmsh_2D # unifies v2.2.8 and v4.1 mesh reading
 export readGmsh2D, readGmsh2D_v4 # TODO: deprecate
-export read_Gmsh_2D_v2, read_Gmsh_2D_v4 
-export MeshImportOptions 
+export read_Gmsh_2D_v2, read_Gmsh_2D_v4
+export MeshImportOptions
 include("mesh/hohqmesh_utilities.jl")
 export read_HOHQMesh
 
@@ -104,10 +105,16 @@ export ck45 # LSERK 45
 using Requires: @require
 
 function __init__()
-  @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
-    include("TriangulatePlots.jl")
-  end
+    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
+        include("TriangulatePlots.jl")
+    end
+
+    # Until Julia v1.9 is the minimum required version for StartUpDG.jl, we still support Requires.jl
+    @static if !isdefined(Base, :get_extension)
+        @require SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240" begin
+            include("../ext/StartUpDGSummationByPartsOperatorsExt.jl")
+        end
+    end
 end
 
-
 end # module
diff --git a/test/Project.toml b/test/Project.toml
index 013caf76..b65cc55d 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -4,6 +4,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
+SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
@@ -12,4 +13,5 @@ HOHQMesh = "0.2"
 Plots = "1"
 StaticArrays = "1"
 StructArrays = "0.6"
+SummationByPartsOperators = "0.5"
 Suppressor = "0.2"
diff --git a/test/SummationByPartsOperatorsExt_tests.jl b/test/SummationByPartsOperatorsExt_tests.jl
new file mode 100644
index 00000000..90ad63b3
--- /dev/null
+++ b/test/SummationByPartsOperatorsExt_tests.jl
@@ -0,0 +1,113 @@
+
+@testset "SummationByPartsOperators ext tests for nonperiodic operators" begin
+    tol = 200*eps()
+    N = 3
+    f(N,r,s) = r^N + s^N # function for testing differentiation
+
+    D = derivative_operator(MattssonNordström2004(), 
+                            derivative_order=1, accuracy_order=N+1,
+                            xmin=-1.0, xmax=1.0, N=20)
+    @testset "Line" begin
+        rd = RefElemData(Line(), D)
+        @test rd.r == rd.rst[1]
+        @test rd.Np == length(rd.r)  
+        @test rd.Nq == length(rd.rq)    
+        @test abs(sum(rd.wq)) ≈ 2
+        @test abs(sum(rd.wf .* rd.nrJ)) < tol
+        @test rd.Pq*rd.Vq ≈ I
+        @test all(vec.(rd.rstf) .≈ (x->getindex(x,rd.Fmask)).(rd.rst))
+        @test StartUpDG.eigenvalue_inverse_trace_constant(rd) ≈ inverse_trace_constant(rd)    
+
+        md = MeshData(uniform_mesh(Line(), 2), rd)
+        @test estimate_h(1, rd, md) ≈ 0.5
+    end
+
+    @testset "Quad" begin
+        rd = RefElemData(Quad(), D)
+        @test rd.r == rd.rst[1]
+        @test rd.Np == length(rd.r)  
+        @test rd.Nq == length(rd.rq)    
+        @test abs(sum(rd.wq)) ≈ 4
+        @test abs(sum(rd.wf)) ≈ 8
+        @test abs(sum(rd.wf .* rd.nrJ)) + abs(sum(rd.wf .* rd.nsJ)) < tol
+        @test rd.Pq*rd.Vq ≈ I
+        @test all(vec.(rd.rstf) .≈ (x->getindex(x,rd.Fmask)).(rd.rst))
+        @test StartUpDG.eigenvalue_inverse_trace_constant(rd) ≈ inverse_trace_constant(rd)    
+    end
+    
+    @testset "Hex" begin
+        rd = RefElemData(Hex(), D)
+        @test propertynames(rd)[1] == :element_type
+        @test rd.t == rd.rst[3]
+        @test rd.tf == rd.rstf[3]    
+        @test rd.tq == rd.rstq[3]
+        @test rd.rp == rd.rstp[1]    
+        @test rd.sp == rd.rstp[2]
+        @test rd.tp == rd.rstp[3]    
+        @test rd.Np == length(rd.r)  
+        @test rd.Nq == length(rd.rq)    
+        @test abs(sum(rd.wq)) ≈ 8
+        @test abs(sum(rd.wf)) ≈ 6*4
+        @test abs(sum(rd.wf .* rd.nrJ)) < tol
+        @test abs(sum(rd.wf .* rd.nsJ)) < tol
+        @test abs(sum(rd.wf .* rd.ntJ)) < tol
+        @test rd.Pq*rd.Vq ≈ I
+        @test StartUpDG.eigenvalue_inverse_trace_constant(rd) ≈ inverse_trace_constant(rd)
+    end
+
+end
+
+@testset "SummationByPartsOperators ext tests for periodic operators" begin
+    tol = 200*eps()
+    N = 3
+    D = periodic_derivative_operator(derivative_order=1, accuracy_order=N+1,
+                                     xmin=0.0, xmax=2.0, N=20)
+
+    @testset "Line" begin
+        rd = RefElemData(Line(), D)
+        @test rd.r == rd.rst[1]
+        @test rd.Np == length(rd.r)  
+        @test rd.Nq == length(rd.rq)    
+        @test abs(sum(rd.wq)) ≈ 2
+
+        # diff matrices should be skew symmetric
+        Qrst = (A -> rd.M * A).(rd.Drst)
+        @test all((A -> norm(A+A')).(Qrst) .< tol)
+
+        @test rd.Pq*rd.Vq ≈ I
+
+    end
+
+    @testset "Quad" begin
+        rd = RefElemData(Quad(), D)
+        @test rd.r == rd.rst[1]
+        @test rd.Np == length(rd.r)  
+        @test rd.Nq == length(rd.rq)    
+        @test abs(sum(rd.wq)) ≈ 4
+
+        # diff matrices should be skew symmetric
+        Qrst = (A -> rd.M * A).(rd.Drst)
+        @test all((A -> norm(A+A')).(Qrst) .< tol)
+
+        @test rd.Pq*rd.Vq ≈ I
+    end
+    
+    @testset "Hex" begin
+        rd = RefElemData(Hex(), D)
+        @test propertynames(rd)[1] == :element_type
+        @test rd.t == rd.rst[3]
+        @test rd.tq == rd.rstq[3]
+        @test rd.rp == rd.rstp[1]    
+        @test rd.sp == rd.rstp[2]
+        @test rd.tp == rd.rstp[3]    
+        @test rd.Np == length(rd.r)  
+        @test rd.Nq == length(rd.rq)    
+        @test abs(sum(rd.wq)) ≈ 8
+        
+        # diff matrices should be skew symmetric
+        Qrst = (A -> rd.M * A).(rd.Drst)
+        @test all((A -> norm(A+A')).(Qrst) .< tol)
+
+        @test rd.Pq*rd.Vq ≈ I
+    end
+end
\ No newline at end of file
diff --git a/test/sbp_tests.jl b/test/multidim_sbp_tests.jl
similarity index 98%
rename from test/sbp_tests.jl
rename to test/multidim_sbp_tests.jl
index 164cf167..2ff50614 100644
--- a/test/sbp_tests.jl
+++ b/test/multidim_sbp_tests.jl
@@ -1,4 +1,4 @@
-@testset "SBP elements" begin
+@testset "Multidimensional SBP tests" begin
 
     tol = 200*eps()
 
@@ -64,5 +64,4 @@
         @test StartUpDG.eigenvalue_inverse_trace_constant(rd) ≈ inverse_trace_constant(rd)
     end
 
-
 end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 76e6f3a9..af4b0d34 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,13 +1,19 @@
 using Test
 using Suppressor
 using LinearAlgebra
+
+# we need to load this first before StartUpDG for Julia pre-v1.9 
+# since Requires.jl needs us to load this first.
+using SummationByPartsOperators
+
 using StartUpDG
 
 include("write_vtk_tests.jl")
 include("named_array_partition_tests.jl")
 include("triangulate_tests.jl")
 include("reference_elem_tests.jl")
-include("sbp_tests.jl")
+include("multidim_sbp_tests.jl")
+include("SummationByPartsOperatorsExt_tests.jl")
 include("MeshData_tests.jl")
 include("boundary_util_tests.jl")
 include("hybrid_mesh_tests.jl")