diff --git a/AUTHORS b/AUTHORS
index c8e8e130f3..40532a43d6 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -8,6 +8,7 @@ Yohann Dudouit
 Leila Ghaffari
 Tzanio Kolev
 David Medina
+Will Pazner
 Thilina Rathnayake
 Jeremy L. Thompson
 Stan Tomov
diff --git a/julia/LibCEED.jl/.gitignore b/julia/LibCEED.jl/.gitignore
new file mode 100644
index 0000000000..3e3da309a2
--- /dev/null
+++ b/julia/LibCEED.jl/.gitignore
@@ -0,0 +1,4 @@
+Manifest.toml
+deps/build.log
+deps/deps.jl
+docs/build
diff --git a/julia/LibCEED.jl/Project.toml b/julia/LibCEED.jl/Project.toml
new file mode 100644
index 0000000000..63caac695d
--- /dev/null
+++ b/julia/LibCEED.jl/Project.toml
@@ -0,0 +1,20 @@
+name = "LibCEED"
+uuid = "2cd74e05-b976-4426-91fa-5f1011f8952b"
+authors = ["Will Pazner <will.e.p@gmail.com>"]
+version = "0.1.0"
+
+[deps]
+CEnum = "fa961155-64e5-5f13-b03f-caf6b980ea82"
+Cassette = "7057c7e9-c182-5462-911a-8362d720325c"
+Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+UnsafeArrays = "c4a57d5a-5b31-53a6-b365-19f8c011fbd6"
+libCEED_jll = "762fde13-7596-547b-826d-8223c52d51c1"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]
diff --git a/julia/LibCEED.jl/README.md b/julia/LibCEED.jl/README.md
new file mode 100644
index 0000000000..38c1032407
--- /dev/null
+++ b/julia/LibCEED.jl/README.md
@@ -0,0 +1,49 @@
+# LibCEED.jl: Julia Interface for [libCEED](https://github.com/CEED/libCEED)
+
+## Installation
+
+When the LibCEED.jl package is built, it requires the environment variable
+`JULIA_LIBCEED_LIB` to be set to the location of the compiled libCEED shared
+library.
+
+For example, the package can be installed by:
+```julia
+% JULIA_LIBCEED_LIB=/path/to/libceed.so julia
+julia> # press ] to enter package manager
+
+(@v1.5) pkg> add LibCEED
+```
+or, equivalently,
+```julia
+% julia
+
+julia> withenv("JULIA_LIBCEED_LIB" => "/path/to/libceed.so") do
+    Pkg.add("LibCEED")
+end
+```
+
+
+## Usage
+
+This package provides both a low-level and high-level interface for libCEED.
+
+### Low-Level Interface
+
+The low-level interface (provided in the `LibCEED.C` module) is in one-to-one
+correspondence with the C libCEED iterface, and is automatically generated (with
+some minor manual modifications) using the Julia package Clang.jl. The script
+used to generate bindings is available in `generate_bindings.jl`.
+
+With the low-level interface, the user is responsible for freeing all allocated
+memory (calling the appropriate `Ceed*Destroy` functions). This interface is
+not type-safe, and calling functions with the wrong arguments can cause libCEED
+to crash.
+
+### High-Level Interface
+
+The high-level interface provides a more idiomatic Julia interface to the
+libCEED library. Objects allocated using the high-level interface will
+automatically be destroyed by the garbage collector, so the user does not need
+to manually manage memory.
+
+See the documentation for more information.
diff --git a/julia/LibCEED.jl/deps/build.jl b/julia/LibCEED.jl/deps/build.jl
new file mode 100644
index 0000000000..118e024a67
--- /dev/null
+++ b/julia/LibCEED.jl/deps/build.jl
@@ -0,0 +1,18 @@
+const libceed_envvar = "JULIA_LIBCEED_LIB"
+
+if isfile("deps.jl")
+   rm("deps.jl")
+end
+
+if haskey(ENV, libceed_envvar)
+   ceedpath = ENV[libceed_envvar]
+   @info "Using libCEED library specified by $libceed_envvar at $ceedpath"
+   open("deps.jl", write=true) do f
+      println(f, "const libceed = \"$(escape_string(ceedpath))\"")
+   end
+else
+   @info "Using prebuilt libCEED binaries provided by libCEED_jll"
+   open("deps.jl", write=true) do f
+      println(f, "using libCEED_jll")
+   end
+end
diff --git a/julia/LibCEED.jl/src/Basis.jl b/julia/LibCEED.jl/src/Basis.jl
new file mode 100644
index 0000000000..fdb55b79cf
--- /dev/null
+++ b/julia/LibCEED.jl/src/Basis.jl
@@ -0,0 +1,391 @@
+abstract type AbstractBasis end
+
+"""
+    BasisCollocated()
+
+Returns the singleton object corresponding to libCEED's `CEED_BASIS_COLLOCATED`.
+"""
+struct BasisCollocated <: AbstractBasis end
+Base.getindex(::BasisCollocated) = C.CEED_BASIS_COLLOCATED[]
+
+"""
+    Basis
+
+Wraps a `CeedBasis` object, representing a finite element basis. A `Basis` object can be
+created using one of:
+
+- [`create_tensor_h1_lagrange_basis`](@ref)
+- [`create_tensor_h1_basis`](@ref)
+- [`create_h1_basis`](@ref)
+"""
+mutable struct Basis <: AbstractBasis
+    ref::RefValue{C.CeedBasis}
+    function Basis(ref)
+        obj = new(ref)
+        finalizer(obj) do x
+            # ccall(:jl_safe_printf, Cvoid, (Cstring, Cstring), "Finalizing %s.\n", repr(x))
+            destroy(x)
+        end
+        return obj
+    end
+end
+destroy(b::Basis) = C.CeedBasisDestroy(b.ref) # COV_EXCL_LINE
+Base.getindex(b::Basis) = b.ref[]
+Base.show(io::IO, ::MIME"text/plain", b::Basis) = ceed_show(io, b, C.CeedBasisView)
+
+@doc raw"""
+    create_tensor_h1_lagrange_basis(ceed, dim, ncomp, p, q, qmode)
+
+Create a tensor-product Lagrange basis.
+
+# Arguments:
+- `ceed`:  A [`Ceed`](@ref) object where the [`Basis`](@ref) will be created.
+- `dim`:   Topological dimension of element.
+- `ncomp`: Number of field components (1 for scalar fields).
+- `p`:     Number of Gauss-Lobatto nodes in one dimension.  The polynomial degree of the
+           resulting $Q_k$ element is $k=p-1$.
+- `q`:     Number of quadrature points in one dimension.
+- `qmode`: Distribution of the $q$ quadrature points (affects order of accuracy for the
+           quadrature).
+"""
+function create_tensor_h1_lagrange_basis(c::Ceed, dim, ncomp, p, q, quad_mode::QuadMode)
+    ref = Ref{C.CeedBasis}()
+    C.CeedBasisCreateTensorH1Lagrange(c[], dim, ncomp, p, q, quad_mode, ref)
+    Basis(ref)
+end
+
+@doc raw"""
+    create_tensor_h1_basis(c::Ceed, dim, ncomp, p, q, interp1d, grad1d, qref1d, qweight1d)
+
+Create a tensor-product basis for $H^1$ discretizations.
+
+# Arguments:
+- `ceed`:      A [`Ceed`](@ref) object where the [`Basis`](@ref) will be created.
+- `dim`:       Topological dimension.
+- `ncomp`:     Number of field components (1 for scalar fields).
+- `p`:         Number of nodes in one dimension.
+- `q`:         Number of quadrature points in one dimension
+- `interp1d`:  Matrix of size `(q, p)` expressing the values of nodal basis functions at
+               quadrature points.
+- `grad1d`:    Matrix of size `(p, q)` expressing derivatives of nodal basis functions at
+               quadrature points.
+- `qref1d`:    Array of length `q` holding the locations of quadrature points on the 1D
+               reference element $[-1, 1]$.
+- `qweight1d`: Array of length `q` holding the quadrature weights on the reference element.
+"""
+function create_tensor_h1_basis(
+    c::Ceed,
+    dim,
+    ncomp,
+    p,
+    q,
+    interp1d,
+    grad1d,
+    qref1d,
+    qweight1d,
+)
+    @assert size(interp1d) == (q, p)
+    @assert size(grad1d) == (q, p)
+    @assert length(qref1d) == q
+    @assert length(qweight1d) == q
+
+    # Convert from Julia matrices (column-major) to row-major format
+    interp1d_rowmajor = collect(interp1d')
+    grad1d_rowmajor = collect(grad1d')
+
+    ref = Ref{C.CeedBasis}()
+    C.CeedBasisCreateTensorH1(
+        c[],
+        dim,
+        ncomp,
+        p,
+        q,
+        interp1d_rowmajor,
+        grad1d_rowmajor,
+        qref1d,
+        qweight1d,
+        ref,
+    )
+    Basis(ref)
+end
+
+@doc raw"""
+    create_h1_basis(c::Ceed, topo::Topology, ncomp, nnodes, nqpts, interp, grad, qref, qweight)
+
+Create a non tensor-product basis for H^1 discretizations
+
+# Arguments:
+- `ceed`:    A [`Ceed`](@ref) object where the [`Basis`](@ref) will be created.
+- `topo`:    [`Topology`](@ref) of element, e.g. hypercube, simplex, etc.
+- `ncomp`:   Number of field components (1 for scalar fields).
+- `nnodes`:  Total number of nodes.
+- `nqpts`:   Total number of quadrature points.
+- `interp`:  Matrix of size `(nqpts, nnodes)` expressing the values of nodal basis functions
+             at quadrature points.
+- `grad`:    Array of size `(dim, nqpts, nnodes)` expressing derivatives of nodal basis
+             functions at quadrature points.
+- `qref`:    Array of length `nqpts` holding the locations of quadrature points on the
+             reference element $[-1, 1]$.
+- `qweight`: Array of length `nqpts` holding the quadrature weights on the reference
+             element.
+"""
+function create_h1_basis(
+    c::Ceed,
+    topo::Topology,
+    ncomp,
+    nnodes,
+    nqpts,
+    interp,
+    grad,
+    qref,
+    qweight,
+)
+    dim = getdimension(topo)
+    @assert size(interp) == (nqpts, nnodes)
+    @assert size(grad) == (dim, nqpts, nnodes)
+    @assert length(qref) == nqpts
+    @assert length(qweight) == nqpts
+
+    # Convert from Julia matrices and tensors (column-major) to row-major format
+    interp_rowmajor = collect(interp')
+    grad_rowmajor = permutedims(grad, [3, 2, 1])
+
+    ref = Ref{C.CeedBasis}()
+    C.CeedBasisCreateH1(
+        c[],
+        topo,
+        ncomp,
+        nnodes,
+        nqpts,
+        interp_rowmajor,
+        grad_rowmajor,
+        qref,
+        qweight,
+        ref,
+    )
+    Basis(ref)
+end
+
+"""
+    apply!(b::Basis, nelem, tmode::TransposeMode, emode::EvalMode, u::AbstractCeedVector, v::AbstractCeedVector)
+
+Apply basis evaluation from nodes to quadrature points or vice versa, storing the result in
+the [`CeedVector`](@ref) `v`.
+
+`nelem` specifies the number of elements to apply the basis evaluation to; the backend will
+specify the ordering in CeedElemRestrictionCreateBlocked()
+
+Set `tmode` to `CEED_NOTRANSPOSE` to evaluate from nodes to quadrature or to
+`CEED_TRANSPOSE` to apply the transpose, mapping from quadrature points to nodes.
+
+Set the [`EvalMode`](@ref) `emode` to:
+- `CEED_EVAL_NONE` to use values directly,
+- `CEED_EVAL_INTERP` to use interpolated values,
+- `CEED_EVAL_GRAD` to use gradients,
+- `CEED_EVAL_WEIGHT` to use quadrature weights.
+"""
+function apply!(
+    b::Basis,
+    nelem,
+    tmode::TransposeMode,
+    emode::EvalMode,
+    u::AbstractCeedVector,
+    v::AbstractCeedVector,
+)
+    C.CeedBasisApply(b[], nelem, tmode, emode, u[], v[])
+end
+
+"""
+    apply(b::Basis, u::AbstractVector; nelem=1, tmode=NOTRANSPOSE, emode=EVAL_INTERP)
+
+Performs the same function as the above-defined [`apply!`](@ref apply!(b::Basis, nelem,
+tmode::TransposeMode, emode::EvalMode, u::AbstractCeedVector, v::AbstractCeedVector)), but
+automatically convert from Julia arrays to [`CeedVector`](@ref) for convenience.
+
+The result will be returned in a newly allocated array of the correct size.
+"""
+function apply(b::Basis, u::AbstractVector; nelem=1, tmode=NOTRANSPOSE, emode=EVAL_INTERP)
+    ceed_ref = Ref{C.Ceed}()
+    ccall((:CeedBasisGetCeed, C.libceed), Cint, (C.CeedBasis, Ptr{C.Ceed}), b[], ceed_ref)
+    c = Ceed(ceed_ref)
+
+    u_vec = CeedVector(c, u)
+
+    len_v = (tmode == TRANSPOSE) ? getnumnodes(b) : getnumqpts(b)
+    if emode == EVAL_GRAD
+        len_v *= getdimension(b)
+    end
+
+    v_vec = CeedVector(c, len_v)
+
+    apply!(b, nelem, tmode, emode, u_vec, v_vec)
+    Vector(v_vec)
+end
+
+"""
+    getdimension(b::Basis)
+
+Return the spatial dimension of the given [`Basis`](@ref).
+"""
+function getdimension(b::Basis)
+    dim = Ref{CeedInt}()
+    C.CeedBasisGetDimension(b[], dim)
+    dim[]
+end
+
+"""
+    getdimension(t::Topology)
+
+Return the spatial dimension of the given [`Topology`](@ref).
+"""
+function getdimension(t::Topology)
+    return Int(t) >> 16
+end
+
+"""
+    gettopology(b::Basis)
+
+Return the [`Topology`](@ref) of the given [`Basis`](@ref).
+"""
+function gettopology(b::Basis)
+    topo = Ref{Topology}()
+    C.CeedBasisGetTopology(b[], topo)
+    topo[]
+end
+
+"""
+    getnumcomponents(b::Basis)
+
+Return the number of components of the given [`Basis`](@ref).
+"""
+function getnumcomponents(b::Basis)
+    ncomp = Ref{CeedInt}()
+    C.CeedBasisGetNumComponents(b[], ncomp)
+    ncomp[]
+end
+
+"""
+    getnumnodes(b::Basis)
+
+Return the number of nodes of the given [`Basis`](@ref).
+"""
+function getnumnodes(b::Basis)
+    nnodes = Ref{CeedInt}()
+    C.CeedBasisGetNumNodes(b[], nnodes)
+    nnodes[]
+end
+
+"""
+    getnumnodes1d(b::Basis)
+
+    Return the number of 1D nodes of the given (tensor-product) [`Basis`](@ref).
+"""
+function getnumnodes1d(b::Basis)
+    nnodes1d = Ref{CeedInt}()
+    C.CeedBasisGetNumNodes1D(b[], nnodes1d)
+    nnodes1d[]
+end
+
+"""
+    getnumqpts(b::Basis)
+
+Return the number of quadrature points of the given [`Basis`](@ref).
+"""
+function getnumqpts(b::Basis)
+    nqpts = Ref{CeedInt}()
+    C.CeedBasisGetNumQuadraturePoints(b[], nqpts)
+    nqpts[]
+end
+
+"""
+    getnumqpts1d(b::Basis)
+
+Return the number of 1D quadrature points of the given (tensor-product) [`Basis`](@ref).
+"""
+function getnumqpts1d(b::Basis)
+    nqpts1d = Ref{CeedInt}()
+    C.CeedBasisGetNumQuadraturePoints1D(b[], nqpts1d)
+    nqpts1d[]
+end
+
+"""
+    getqref(b::Basis)
+
+Get the reference coordinates of quadrature points (in `dim` dimensions) of the given
+[`Basis`](@ref).
+"""
+function getqref(b::Basis)
+    ref = Ref{Ptr{CeedScalar}}()
+    C.CeedBasisGetQRef(b[], ref)
+    copy(unsafe_wrap(Array, ref[], getnumqpts(b)))
+end
+
+"""
+    getqref(b::Basis)
+
+Get the quadrature weights of quadrature points (in `dim` dimensions) of the given
+[`Basis`](@ref).
+"""
+function getqweights(b::Basis)
+    ref = Ref{Ptr{CeedScalar}}()
+    C.CeedBasisGetQWeights(b[], ref)
+    copy(unsafe_wrap(Array, ref[], getnumqpts(b)))
+end
+
+"""
+    getinterp(b::Basis)
+
+Get the interpolation matrix of the given [`Basis`](@ref). Returns a matrix of size
+`(getnumqpts(b), getnumnodes(b))`.
+"""
+function getinterp(b::Basis)
+    ref = Ref{Ptr{CeedScalar}}()
+    C.CeedBasisGetInterp(b[], ref)
+    q = getnumqpts(b)
+    p = getnumnodes(b)
+    collect(unsafe_wrap(Array, ref[], (p, q))')
+end
+
+"""
+    getinterp1d(b::Basis)
+
+Get the 1D interpolation matrix of the given [`Basis`](@ref). `b` must be a tensor-product
+basis, otherwise this function will fail. Returns a matrix of size `(getnumqpts1d(b),
+getnumnodes1d(b))`.
+"""
+function getinterp1d(b::Basis)
+    ref = Ref{Ptr{CeedScalar}}()
+    C.CeedBasisGetInterp1D(b[], ref)
+    q = getnumqpts1d(b)
+    p = getnumnodes1d(b)
+    collect(unsafe_wrap(Array, ref[], (p, q))')
+end
+
+"""
+    getgad(b::Basis)
+
+Get the gradient matrix of the given [`Basis`](@ref). Returns a tensor of size
+`(getdimension(b), getnumqpts(b), getnumnodes(b))`.
+"""
+function getgrad(b::Basis)
+    ref = Ref{Ptr{CeedScalar}}()
+    C.CeedBasisGetGrad(b[], ref)
+    dim = getdimension(b)
+    q = getnumqpts(b)
+    p = getnumnodes(b)
+    permutedims(unsafe_wrap(Array, ref[], (p, q, dim)), [3, 2, 1])
+end
+
+"""
+    getgrad1d(b::Basis)
+
+Get the 1D derivative matrix of the given [`Basis`](@ref). Returns a matrix of size
+`(getnumqpts(b), getnumnodes(b))`.
+"""
+function getgrad1d(b::Basis)
+    ref = Ref{Ptr{CeedScalar}}()
+    C.CeedBasisGetGrad1D(b[], ref)
+    q = getnumqpts1d(b)
+    p = getnumnodes1d(b)
+    collect(unsafe_wrap(Array, ref[], (p, q))')
+end
diff --git a/julia/LibCEED.jl/src/C.jl b/julia/LibCEED.jl/src/C.jl
new file mode 100644
index 0000000000..40e2114d03
--- /dev/null
+++ b/julia/LibCEED.jl/src/C.jl
@@ -0,0 +1,45 @@
+# Low-level C API for libCEED
+
+module C
+
+using CEnum, Libdl
+
+# Get the path to the libCEED dynamic library, configured during the build step
+# of the LibCEED.jl package.
+const depsfile = joinpath(@__DIR__, "..", "deps", "deps.jl")
+if !isfile(depsfile)
+    error("LibCEED.jl not properly installed. Please run Pkg.build(\"LibCEED\")")
+end
+include(depsfile)
+
+include(joinpath(@__DIR__, "generated", "libceed_common.jl"))
+include(joinpath(@__DIR__, "generated", "libceed_api.jl"))
+
+const CEED_STRIDES_BACKEND = Ref{Ptr{CeedInt}}()
+const CEED_BASIS_COLLOCATED = Ref{CeedBasis}()
+const CEED_VECTOR_ACTIVE = Ref{CeedVector}()
+const CEED_VECTOR_NONE = Ref{CeedVector}()
+const CEED_ELEMRESTRICTION_NONE = Ref{CeedElemRestriction}()
+const CEED_QFUNCTION_NONE = Ref{CeedQFunction}()
+const CEED_REQUEST_IMMEDIATE = Ref{CeedRequest}()
+const CEED_REQUEST_ORDERED = Ref{CeedRequest}()
+
+function __init__()
+    global libceed_handle = dlopen(libceed)
+    # some global variables
+    CEED_STRIDES_BACKEND[] = cglobal((:CEED_STRIDES_BACKEND, libceed))
+    CEED_BASIS_COLLOCATED[] =
+        unsafe_load(cglobal((:CEED_BASIS_COLLOCATED, libceed), CeedBasis))
+    CEED_VECTOR_ACTIVE[] = unsafe_load(cglobal((:CEED_VECTOR_ACTIVE, libceed), CeedVector))
+    CEED_VECTOR_NONE[] = unsafe_load(cglobal((:CEED_VECTOR_NONE, libceed), CeedVector))
+    CEED_ELEMRESTRICTION_NONE[] =
+        unsafe_load(cglobal((:CEED_ELEMRESTRICTION_NONE, libceed), CeedElemRestriction))
+    CEED_QFUNCTION_NONE[] =
+        unsafe_load(cglobal((:CEED_QFUNCTION_NONE, libceed), CeedQFunction))
+    CEED_REQUEST_IMMEDIATE[] =
+        unsafe_load(cglobal((:CEED_REQUEST_IMMEDIATE, libceed), Ptr{CeedRequest}))
+    CEED_REQUEST_ORDERED[] =
+        unsafe_load(cglobal((:CEED_REQUEST_ORDERED, libceed), Ptr{CeedRequest}))
+end
+
+end # module
diff --git a/julia/LibCEED.jl/src/Ceed.jl b/julia/LibCEED.jl/src/Ceed.jl
new file mode 100644
index 0000000000..d919e75812
--- /dev/null
+++ b/julia/LibCEED.jl/src/Ceed.jl
@@ -0,0 +1,113 @@
+struct CeedError <: Exception
+    fname::String
+    lineno::Int
+    func::String
+    ecode::Int
+    message::String
+end
+
+# COV_EXCL_START
+function Base.showerror(io::IO, e::CeedError)
+    println(io, "libCEED error code ", e.ecode, " in ", e.func)
+    println(io, e.fname, ':', e.lineno, '\n')
+    println(io, e.message)
+end
+# COV_EXCL_STOP
+
+function handle_ceed_error(
+    ceed::C.Ceed,
+    c_fname::Cstring,
+    lineno::Cint,
+    c_func::Cstring,
+    ecode::Cint,
+    c_format::Cstring,
+    args::Ptr{Cvoid},
+)
+    c_message = ccall(
+        (:CeedErrorFormat, C.libceed),
+        Cstring,
+        (C.Ceed, Cstring, Ptr{Cvoid}),
+        ceed,
+        c_format,
+        args,
+    )
+    fname = unsafe_string(c_fname)
+    func = unsafe_string(c_func)
+    message = unsafe_string(c_message)
+    throw(CeedError(fname, lineno, func, ecode, message))
+end
+
+mutable struct Ceed
+    ref::RefValue{C.Ceed}
+end
+
+"""
+    Ceed(spec="/cpu/self")
+
+Wraps a libCEED `Ceed` object, created with the given resource specification string.
+"""
+function Ceed(spec::AbstractString="/cpu/self")
+    obj = Ceed(Ref{C.Ceed}())
+    C.CeedInit(spec, obj.ref)
+    ehandler = @cfunction(
+        handle_ceed_error,
+        Cint,
+        (C.Ceed, Cstring, Cint, Cstring, Cint, Cstring, Ptr{Cvoid})
+    )
+    C.CeedSetErrorHandler(obj.ref[], ehandler)
+    finalizer(obj) do x
+        # ccall(:jl_safe_printf, Cvoid, (Cstring, Cstring), "Finalizing %s.\n", repr(x))
+        destroy(x)
+    end
+    return obj
+end
+destroy(c::Ceed) = C.CeedDestroy(c.ref) # COV_EXCL_LINE
+Base.getindex(c::Ceed) = c.ref[]
+
+Base.show(io::IO, ::MIME"text/plain", c::Ceed) = ceed_show(io, c, C.CeedView)
+
+"""
+    getresource(c::Ceed)
+
+Returns the resource string associated with the given [`Ceed`](@ref) object.
+"""
+function getresource(c::Ceed)
+    res = Ref{Cstring}()
+    C.CeedGetResource(c[], res)
+    unsafe_string(res[])
+end
+
+"""
+    isdeterministic(c::Ceed)
+
+Returns true if backend of the given [`Ceed`](@ref) object is deterministic, and false
+otherwise.
+"""
+function isdeterministic(c::Ceed)
+    isdet = Ref{Bool}()
+    C.CeedIsDeterministic(c[], isdet)
+    isdet[]
+end
+
+"""
+    get_preferred_memtype(c::Ceed)
+
+Returns the preferred [`MemType`](@ref) (either `MEM_HOST` or `MEM_DEVICE`) of the given
+[`Ceed`](@ref) object.
+"""
+function get_preferred_memtype(c::Ceed)
+    mtype = Ref{MemType}()
+    C.CeedGetPreferredMemType(c[], mtype)
+    mtype[]
+end
+
+"""
+    iscuda(c::Ceed)
+
+Returns true if the given [`Ceed`](@ref) object has resource `"/gpu/cuda/*"` and false
+otherwise.
+"""
+function iscuda(c::Ceed)
+    res_split = split(getresource(c), "/")
+    length(res_split) >= 3 && res_split[3] == "cuda"
+end
diff --git a/julia/LibCEED.jl/src/CeedVector.jl b/julia/LibCEED.jl/src/CeedVector.jl
new file mode 100644
index 0000000000..7d4dc0304a
--- /dev/null
+++ b/julia/LibCEED.jl/src/CeedVector.jl
@@ -0,0 +1,308 @@
+import LinearAlgebra: norm
+
+abstract type AbstractCeedVector end
+
+struct CeedVectorActive <: AbstractCeedVector end
+Base.getindex(::CeedVectorActive) = C.CEED_VECTOR_ACTIVE[]
+
+struct CeedVectorNone <: AbstractCeedVector end
+Base.getindex(::CeedVectorNone) = C.CEED_VECTOR_NONE[]
+
+mutable struct CeedVector <: AbstractCeedVector
+    ref::RefValue{C.CeedVector}
+    arr::Union{Nothing,AbstractArray}
+    CeedVector(ref::Ref{C.CeedVector}) = new(ref, nothing)
+end
+
+"""
+    CeedVector(c::Ceed, len::Integer)
+
+Creates a `CeedVector` of given length.
+"""
+function CeedVector(c::Ceed, len::Integer)
+    ref = Ref{C.CeedVector}()
+    C.CeedVectorCreate(c[], len, ref)
+    obj = CeedVector(ref)
+    finalizer(obj) do x
+        # ccall(:jl_safe_printf, Cvoid, (Cstring, Cstring), "Finalizing %s.\n", repr(x))
+        destroy(x)
+    end
+    return obj
+end
+destroy(v::CeedVector) = C.CeedVectorDestroy(v.ref) # COV_EXCL_LINE
+Base.getindex(v::CeedVector) = v.ref[]
+
+Base.summary(io::IO, v::CeedVector) = print(io, length(v), "-element CeedVector")
+function Base.show(io::IO, ::MIME"text/plain", v::CeedVector)
+    summary(io, v)
+    println(io, ":")
+    witharray_read(v, MEM_HOST) do arr
+        Base.print_array(io, arr)
+    end
+end
+Base.show(io::IO, v::CeedVector) = witharray_read(a -> show(io, a), v, MEM_HOST)
+
+function Base.length(::Type{T}, v::CeedVector) where {T}
+    len = Ref{C.CeedInt}()
+    C.CeedVectorGetLength(v[], len)
+    return T(len[])
+end
+
+Base.ndims(::CeedVector) = 1
+Base.ndims(::Type{CeedVector}) = 1
+Base.axes(v::CeedVector) = (Base.OneTo(length(v)),)
+Base.size(v::CeedVector) = (length(Int, v),)
+Base.length(v::CeedVector) = length(Int, v)
+
+"""
+    setvalue!(v::CeedVector, val::CeedScalar)
+
+Set the [`CeedVector`](@ref) to a constant value.
+"""
+setvalue!(v::CeedVector, val::CeedScalar) = C.CeedVectorSetValue(v[], val)
+"""
+    setindex!(v::CeedVector, val::CeedScalar)
+    v[] = val
+
+Set the [`CeedVector`](@ref) to a constant value, synonymous to [`setvalue!`](@ref).
+"""
+Base.setindex!(v::CeedVector, val::CeedScalar) = setvalue!(v, val)
+
+"""
+    norm(v::CeedVector, ntype::NormType)
+
+Return the norm of the given [`CeedVector`](@ref).
+
+The norm type can either be specified as one of `NORM_1`, `NORM_2`, `NORM_MAX`.
+"""
+function norm(v::CeedVector, ntype::NormType)
+    nrm = Ref{CeedScalar}()
+    C.CeedVectorNorm(v[], ntype, nrm)
+    nrm[]
+end
+
+"""
+    norm(v::CeedVector, p::Real)
+
+Return the norm of the given [`CeedVector`](@ref), see [`norm(::CeedVector,
+::NormType)`](@ref).
+
+`p` can have value 1, 2, or Inf, corresponding to `NORM_1`, `NORM_2`, and `NORM_MAX`,
+respectively.
+"""
+function norm(v::CeedVector, p::Real)
+    if p == 1
+        ntype = NORM_1
+    elseif p == 2
+        ntype = NORM_2
+    elseif isinf(p)
+        ntype = NORM_MAX
+    else
+        error("norm(v::CeedVector, p): p must be 1, 2, or Inf")
+    end
+    norm(v, ntype)
+end
+
+"""
+    reciprocal!(v::CeedVector)
+
+Set `v` to be equal to its elementwise reciprocal.
+"""
+reciprocal!(v::CeedVector) = C.CeedVectorReciprocal(v[])
+
+"""
+    setarray!(v::CeedVector, mtype::MemType, cmode::CopyMode, arr)
+
+Set the array used by a [`CeedVector`](@ref), freeing any previously allocated array if
+applicable. The backend may copy values to a different [`MemType`](@ref). See also
+[`syncarray!`](@ref) and [`takearray!`](@ref).
+
+!!! warning "Avoid OWN_POINTER CopyMode"
+    The [`CopyMode`](@ref) `OWN_POINTER` is not suitable for use with arrays that are
+    allocated by Julia, since those cannot be properly freed from libCEED.
+"""
+function setarray!(v::CeedVector, mtype::MemType, cmode::CopyMode, arr)
+    C.CeedVectorSetArray(v[], mtype, cmode, arr)
+    if cmode == USE_POINTER
+        v.arr = arr
+    end
+end
+
+"""
+    syncarray!(v::CeedVector, mtype::MemType)
+
+Sync the [`CeedVector`](@ref) to a specified [`MemType`](@ref). This function is used to
+force synchronization of arrays set with [`setarray!`](@ref). If the requested memtype is
+already synchronized, this function results in a no-op.
+"""
+syncarray!(v::CeedVector, mtype::MemType) = C.CeedVectorSyncArray(v[], mtype)
+
+"""
+    takearray!(v::CeedVector, mtype::MemType)
+
+Take ownership of the [`CeedVector`](@ref) array and remove the array from the
+[`CeedVector`](@ref). The caller is responsible for managing and freeing the array. The
+array is returns as a `Ptr{CeedScalar}`.
+"""
+function takearray!(v::CeedVector, mtype::MemType)
+    ptr = Ref{Ptr{CeedScalar}}()
+    C.CeedVectorTakeArray(v[], mtype, ptr)
+    v.arr = nothing
+    ptr[]
+end
+
+# Helper function to parse arguments of @witharray and @witharray_read
+function witharray_parse(assignment, args)
+    if !Meta.isexpr(assignment, :(=))
+        error("@witharray must have first argument of the form v_arr=v") # COV_EXCL_LINE
+    end
+    arr = assignment.args[1]
+    v = assignment.args[2]
+    mtype = MEM_HOST
+    sz = :((length($(esc(v))),))
+    body = args[end]
+    for i = 1:length(args)-1
+        a = args[i]
+        if !Meta.isexpr(a, :(=))
+            error("Incorrect call to @witharray or @witharray_read") # COV_EXCL_LINE
+        end
+        if a.args[1] == :mtype
+            mtype = a.args[2]
+        elseif a.args[1] == :size
+            sz = esc(a.args[2])
+        end
+    end
+    arr, v, sz, mtype, body
+end
+
+"""
+    @witharray(v_arr=v, [size=(dims...)], [mtype=MEM_HOST], body)
+
+Executes `body`, having extracted the contents of the [`CeedVector`](@ref) `v` as an array
+with name `v_arr`. If the [`memory type`](@ref MemType) `mtype` is not provided, `MEM_HOST`
+will be used. If the size is not specified, a flat vector will be assumed.
+
+# Examples
+Negate the contents of `CeedVector` `v`:
+```
+@witharray v_arr=v v_arr *= -1.0
+```
+"""
+macro witharray(assignment, args...)
+    arr, v, sz, mtype, body = witharray_parse(assignment, args)
+    quote
+        arr_ref = Ref{Ptr{C.CeedScalar}}()
+        C.CeedVectorGetArray($(esc(v))[], $(esc(mtype)), arr_ref)
+        try
+            $(esc(arr)) = UnsafeArray(arr_ref[], Int.($sz))
+            $(esc(body))
+        finally
+            C.CeedVectorRestoreArray($(esc(v))[], arr_ref)
+        end
+    end
+end
+
+"""
+    @witharray_read(v_arr=v, [size=(dims...)], [mtype=MEM_HOST], body)
+
+Same as [`@witharray`](@ref), but provides read-only access to the data.
+"""
+macro witharray_read(assignment, args...)
+    arr, v, sz, mtype, body = witharray_parse(assignment, args)
+    quote
+        arr_ref = Ref{Ptr{C.CeedScalar}}()
+        C.CeedVectorGetArrayRead($(esc(v))[], $(esc(mtype)), arr_ref)
+        try
+            $(esc(arr)) = UnsafeArray(arr_ref[], Int.($sz))
+            $(esc(body))
+        finally
+            C.CeedVectorRestoreArrayRead($(esc(v))[], arr_ref)
+        end
+    end
+end
+
+"""
+    setindex!(v::CeedVector, v2::AbstractArray)
+    v[] = v2
+
+Sets the values of [`CeedVector`](@ref) `v` equal to those of `v2` using broadcasting.
+"""
+Base.setindex!(v::CeedVector, v2::AbstractArray) = @witharray(a = v, a .= v2)
+
+"""
+    CeedVector(c::Ceed, v2::AbstractVector; mtype=MEM_HOST, cmode=COPY_VALUES)
+
+Creates a new [`CeedVector`](@ref) using the contents of the given vector `v2`. By default,
+the contents of `v2` will be copied to the new [`CeedVector`](@ref), but this behavior can
+be changed by specifying a different `cmode`.
+"""
+function CeedVector(c::Ceed, v2::AbstractVector; mtype=MEM_HOST, cmode=COPY_VALUES)
+    v = CeedVector(c, length(v2))
+    setarray!(v, mtype, cmode, v2)
+    v
+end
+
+"""
+    Vector(v::CeedVector)
+
+Create a new `Vector` by copying the contents of `v`.
+"""
+function Base.Vector(v::CeedVector)
+    v2 = Vector{CeedScalar}(undef, length(v))
+    @witharray_read(a = v, v2 .= a)
+end
+
+"""
+    witharray(f, v::CeedVector, mtype=MEM_HOST)
+
+Calls `f` with an array containing the data of the `CeedVector` `v`, using [`memory
+type`](@ref MemType) `mtype`.
+
+Because of performance issues involving closures, if `f` is a complex operation, it may be
+more efficient to use the macro version `@witharray` (cf. the section on "Performance of
+captured variable" in the [Julia
+documentation](https://docs.julialang.org/en/v1/manual/performance-tips) and related [GitHub
+issue](https://github.com/JuliaLang/julia/issues/15276).
+
+# Examples
+
+Return the sum of a vector:
+```
+witharray(sum, v)
+```
+"""
+function witharray(f, v::CeedVector, mtype::MemType=MEM_HOST)
+    arr_ref = Ref{Ptr{C.CeedScalar}}()
+    C.CeedVectorGetArray(v[], mtype, arr_ref)
+    arr = UnsafeArray(arr_ref[], (length(v),))
+    res = try
+        f(arr)
+    finally
+        C.CeedVectorRestoreArray(v[], arr_ref)
+    end
+    return res
+end
+
+"""
+    witharray_read(f, v::CeedVector, mtype::MemType=MEM_HOST)
+
+Same as [`witharray`](@ref), but with read-only access to the data.
+
+# Examples
+
+Display the contents of a vector:
+```
+witharray_read(display, v)
+```
+"""
+function witharray_read(f, v::CeedVector, mtype::MemType=MEM_HOST)
+    arr_ref = Ref{Ptr{C.CeedScalar}}()
+    C.CeedVectorGetArrayRead(v[], mtype, arr_ref)
+    arr = UnsafeArray(arr_ref[], (length(v),))
+    res = try
+        f(arr)
+    finally
+        C.CeedVectorRestoreArrayRead(v[], arr_ref)
+    end
+    return res
+end
diff --git a/julia/LibCEED.jl/src/Context.jl b/julia/LibCEED.jl/src/Context.jl
new file mode 100644
index 0000000000..9e3c27a502
--- /dev/null
+++ b/julia/LibCEED.jl/src/Context.jl
@@ -0,0 +1,56 @@
+mutable struct Context
+    ref::RefValue{C.CeedQFunctionContext}
+    data::Any
+    function Context(ref::Ref{C.CeedQFunctionContext})
+        obj = new(ref)
+        finalizer(obj) do x
+            # ccall(:jl_safe_printf, Cvoid, (Cstring, Cstring), "Finalizing %s.\n", repr(x))
+            C.CeedQFunctionContextDestroy(x.ref)
+        end
+        return obj
+    end
+end
+Base.getindex(ctx::Context) = ctx.ref[]
+Base.show(io::IO, ::MIME"text/plain", c::Context) =
+    ceed_show(io, c, C.CeedQFunctionContextView)
+
+"""
+    Context(ceed::Ceed, data; mtype=MEM_HOST, cmode=USE_POINTER)
+
+Create a `CeedQFunctionContext` object that allows user Q-functions to access an arbitrary
+data object. `data` should be an instance of a mutable struct. If the copy mode `cmode` is
+`USE_POINTER`, then the data will be preserved from the GC when assigned to a `QFunction`
+object using `set_context!`.
+
+Copy mode `OWN_POINTER` is not supported by this interface because Julia-allocated objects
+cannot be freed from C.
+"""
+function Context(c::Ceed, data; mtype=MEM_HOST, cmode=USE_POINTER)
+    ref = Ref{C.CeedQFunctionContext}()
+    C.CeedQFunctionContextCreate(c[], ref)
+    ctx = Context(ref)
+    set_data!(ctx, mtype, cmode, data)
+    return ctx
+end
+
+function set_data!(ctx::Context, mtype, cmode::CopyMode, data)
+    # Store a reference to the context data so that it will not be GC'd before
+    # it is accessed in the user Q-function.
+    # A reference to the context object is stored in the QFunction object, and
+    # references to the QFunctions are stored in the Operator.
+    # This means that when `apply!(op, ...)` is called, the context data is
+    # ensured to be valid.
+    if cmode == USE_POINTER
+        ctx.data = data
+    elseif cmode == OWN_POINTER
+        error("set_data!: copy mode OWN_POINTER is not supported")
+    end
+
+    C.CeedQFunctionContextSetData(
+        ctx[],
+        mtype,
+        cmode,
+        sizeof(data),
+        pointer_from_objref(data),
+    )
+end
diff --git a/julia/LibCEED.jl/src/Cuda.jl b/julia/LibCEED.jl/src/Cuda.jl
new file mode 100644
index 0000000000..372871adf5
--- /dev/null
+++ b/julia/LibCEED.jl/src/Cuda.jl
@@ -0,0 +1,155 @@
+# COV_EXCL_START
+using .CUDA, Cassette
+
+cuda_is_loaded = true
+
+#! format: off
+const cudafuns = (
+    :cos, :cospi, :sin, :sinpi, :tan,
+    :acos, :asin, :atan,
+    :cosh, :sinh, :tanh,
+    :acosh, :asinh, :atanh,
+    :log, :log10, :log1p, :log2,
+    :exp, :exp2, :exp10, :expm1, :ldexp,
+    :abs,
+    :sqrt, :cbrt,
+    :ceil, :floor,
+)
+#! format: on
+
+Cassette.@context CeedCudaContext
+
+@inline function Cassette.overdub(::CeedCudaContext, ::typeof(Core.kwfunc), f)
+    return Core.kwfunc(f)
+end
+@inline function Cassette.overdub(::CeedCudaContext, ::typeof(Core.apply_type), args...)
+    return Core.apply_type(args...)
+end
+@inline function Cassette.overdub(
+    ::CeedCudaContext,
+    ::typeof(StaticArrays.Size),
+    x::Type{<:AbstractArray{<:Any,N}},
+) where {N}
+    return StaticArrays.Size(x)
+end
+
+for f in cudafuns
+    @eval @inline function Cassette.overdub(
+        ::CeedCudaContext,
+        ::typeof(Base.$f),
+        x::Union{Float32,Float64},
+    )
+        return CUDA.$f(x)
+    end
+end
+
+function setarray!(v::CeedVector, mtype::MemType, cmode::CopyMode, arr::CuArray)
+    ptr = Ptr{CeedScalar}(UInt64(pointer(arr)))
+    C.CeedVectorSetArray(v[], mtype, cmode, ptr)
+    if cmode == USE_POINTER
+        v.arr = arr
+    end
+end
+
+struct FieldsCuda
+    inputs::NTuple{16,Int}
+    outputs::NTuple{16,Int}
+end
+
+function generate_kernel(qf_name, kf, dims_in, dims_out)
+    ninputs = length(dims_in)
+    noutputs = length(dims_out)
+
+    input_sz = prod.(dims_in)
+    output_sz = prod.(dims_out)
+
+    f_ins = [Symbol("rqi$i") for i = 1:ninputs]
+    f_outs = [Symbol("rqo$i") for i = 1:noutputs]
+
+    args = Vector{Union{Symbol,Expr}}(undef, ninputs + noutputs)
+    def_ins = Vector{Expr}(undef, ninputs)
+    f_ins_j = Vector{Union{Symbol,Expr}}(undef, ninputs)
+    for i = 1:ninputs
+        if length(dims_in[i]) == 0
+            def_ins[i] = :(local $(f_ins[i]))
+            f_ins_j[i] = f_ins[i]
+            args[i] = f_ins[i]
+        else
+            def_ins[i] =
+                :($(f_ins[i]) = LibCEED.MArray{Tuple{$(dims_in[i]...)},Float64}(undef))
+            f_ins_j[i] = :($(f_ins[i])[j])
+            args[i] = :(LibCEED.SArray{Tuple{$(dims_in[i]...)},Float64}($(f_ins[i])))
+        end
+    end
+    for i = 1:noutputs
+        args[ninputs+i] = f_outs[i]
+    end
+
+    def_outs = [
+        :($(f_outs[i]) = LibCEED.MArray{Tuple{$(dims_out[i]...)},Float64}(undef))
+        for i = 1:noutputs
+    ]
+
+    device_ptr_type = Core.LLVMPtr{CeedScalar,LibCEED.AS.Global}
+
+    read_quads_in = [
+        :(
+            for j = 1:$(input_sz[i])
+                $(f_ins_j[i]) = unsafe_load(
+                    reinterpret($device_ptr_type, fields.inputs[$i]),
+                    q + (j - 1)*Q,
+                    a,
+                )
+            end
+        ) for i = 1:ninputs
+    ]
+
+    write_quads_out = [
+        :(
+            for j = 1:$(output_sz[i])
+                unsafe_store!(
+                    reinterpret($device_ptr_type, fields.outputs[$i]),
+                    $(f_outs[i])[j],
+                    q + (j - 1)*Q,
+                    a,
+                )
+            end
+        ) for i = 1:noutputs
+    ]
+
+    qf = gensym(qf_name)
+    quote
+        function $qf(ctx_ptr, Q, fields)
+            gd = LibCEED.gridDim()
+            bi = LibCEED.blockIdx()
+            bd = LibCEED.blockDim()
+            ti = LibCEED.threadIdx()
+
+            inc = bd.x*gd.x
+
+            $(def_ins...)
+            $(def_outs...)
+
+            # Alignment for data read/write
+            a = Val($(Base.datatype_alignment(CeedScalar)))
+
+            # Cassette context for replacing intrinsics with CUDA versions
+            ctx = LibCEED.CeedCudaContext()
+
+            for q = (ti.x+(bi.x-1)*bd.x):inc:Q
+                $(read_quads_in...)
+                LibCEED.Cassette.overdub(ctx, $kf, ctx_ptr, $(args...))
+                $(write_quads_out...)
+            end
+            return
+        end
+    end
+end
+
+function mk_cufunction(ceed, def_module, qf_name, kf, dims_in, dims_out)
+    k_fn = Core.eval(def_module, generate_kernel(qf_name, kf, dims_in, dims_out))
+    tt = Tuple{Ptr{Nothing},Int32,FieldsCuda}
+    host_k = cufunction(k_fn, tt; maxregs=64)
+    return host_k.fun.handle
+end
+# COV_EXCL_STOP
diff --git a/julia/LibCEED.jl/src/ElemRestriction.jl b/julia/LibCEED.jl/src/ElemRestriction.jl
new file mode 100644
index 0000000000..8c4ab94c3f
--- /dev/null
+++ b/julia/LibCEED.jl/src/ElemRestriction.jl
@@ -0,0 +1,288 @@
+abstract type AbstractElemRestriction end
+
+"""
+    ElemRestrictionNone()
+
+Returns the singleton object corresponding to libCEED's `CEED_ELEMRESTRICTION_NONE`
+"""
+struct ElemRestrictionNone <: AbstractElemRestriction end
+Base.getindex(::ElemRestrictionNone) = C.CEED_ELEMRESTRICTION_NONE[]
+
+"""
+    ElemRestriction
+
+Wraps a `CeedElemRestriction` object, representing the restriction from local vectors to
+elements. An `ElemRestriction` object can be created using [`create_elem_restriction`](@ref)
+or [`create_elem_restriction_strided`](@ref).
+"""
+mutable struct ElemRestriction <: AbstractElemRestriction
+    ref::RefValue{C.CeedElemRestriction}
+    function ElemRestriction(ref)
+        obj = new(ref)
+        finalizer(obj) do x
+            # ccall(:jl_safe_printf, Cvoid, (Cstring, Cstring), "Finalizing %s.\n", repr(x))
+            destroy(x)
+        end
+        return obj
+    end
+end
+destroy(r::ElemRestriction) = C.CeedElemRestrictionDestroy(r.ref) # COV_EXCL_LINE
+Base.getindex(r::ElemRestriction) = r.ref[]
+Base.show(io::IO, ::MIME"text/plain", e::ElemRestriction) =
+    ceed_show(io, e, C.CeedElemRestrictionView)
+
+@doc raw"""
+    create_elem_restriction(
+        ceed::Ceed,
+        nelem,
+        elemsize,
+        ncomp,
+        compstride,
+        lsize,
+        offsets::AbstractArray{CeedInt},
+        mtype::MemType=MEM_HOST,
+        cmode::CopyMode=COPY_VALUES,
+    )
+
+Create a `CeedElemRestriction`.
+
+!!! warning "Zero-based indexing"
+    In the below notation, we are using **0-based indexing**. libCEED expects the offset
+    indices to be 0-based.
+
+# Arguments:
+- `ceed`:       The [`Ceed`](@ref) object
+- `nelem`:      Number of elements described in the `offsets` array
+- `elemsize`:   Size (number of "nodes") per element
+- `ncomp`:      Number of field components per interpolation node (1 for scalar fields)
+- `compstride`: Stride between components for the same L-vector "node". Data for node $i$,
+                component $j$, element $k$ can be found in the L-vector at index `offsets[i
+                + k*elemsize] + j*compstride`.
+- `lsize`:      The size of the L-vector. This vector may be larger than the elements and
+                fields given by this restriction.
+- `offsets`:    Array of shape `(elemsize, nelem)`. Column $i$ holds the ordered list of the
+                offsets (into the input [`CeedVector`](@ref)) for the unknowns corresponding
+                to element $i$, where $0 \leq i < \textit{nelem}$. All offsets must be in
+                the range $[0, \textit{lsize} - 1]$.
+- `mtype`:      Memory type of the `offsets` array, see [`MemType`](@ref)
+- `cmode`:      Copy mode for the `offsets` array, see [`CopyMode`](@ref)
+"""
+function create_elem_restriction(
+    c::Ceed,
+    nelem,
+    elemsize,
+    ncomp,
+    compstride,
+    lsize,
+    offsets::AbstractArray{CeedInt};
+    mtype::MemType=MEM_HOST,
+    cmode::CopyMode=COPY_VALUES,
+)
+    ref = Ref{C.CeedElemRestriction}()
+    C.CeedElemRestrictionCreate(
+        c[],
+        nelem,
+        elemsize,
+        ncomp,
+        compstride,
+        lsize,
+        mtype,
+        cmode,
+        offsets,
+        ref,
+    )
+    ElemRestriction(ref)
+end
+
+@doc raw"""
+    create_elem_restriction_strided(ceed::Ceed, nelem, elemsize, ncomp, lsize, strides)
+
+Create a strided `CeedElemRestriction`.
+
+!!! warning "Zero-based indexing"
+    In the below notation, we are using **0-based indexing**. libCEED expects the offset
+    indices to be 0-based.
+
+# Arguments:
+- `ceed`:     The [`Ceed`](@ref) object
+- `nelem`:    Number of elements described by the restriction
+- `elemsize`: Size (number of "nodes") per element
+- `ncomp`:    Number of field components per interpolation node (1 for scalar fields)
+- `lsize`:    The size of the L-vector. This vector may be larger than the elements and
+              fields given by this restriction.
+- `strides`:  Array for strides between [nodes, components, elements]. Data for node $i$,
+              component $j$, element $k$ can be found in the L-vector at index `i*strides[0]
+              + j*strides[1] + k*strides[2]`. [`STRIDES_BACKEND`](@ref) may be used with
+              vectors created by a Ceed backend.
+"""
+function create_elem_restriction_strided(c::Ceed, nelem, elemsize, ncomp, lsize, strides)
+    ref = Ref{C.CeedElemRestriction}()
+    C.CeedElemRestrictionCreateStrided(c[], nelem, elemsize, ncomp, lsize, strides, ref)
+    ElemRestriction(ref)
+end
+
+"""
+    apply!(
+        r::ElemRestriction,
+        u::CeedVector,
+        ru::CeedVector;
+        tmode=NOTRANSPOSE,
+        request=RequestImmediate(),
+    )
+
+Use the [`ElemRestriction`](@ref) to convert from L-vector to an E-vector (or apply the
+tranpose operation). The input [`CeedVector`](@ref) is `u` and the result stored in `ru`.
+
+If `tmode` is `TRANSPOSE`, then the result is added to `ru`. If `tmode` is `NOTRANSPOSE`,
+then `ru` is overwritten with the result.
+"""
+function apply!(
+    r::ElemRestriction,
+    u::CeedVector,
+    ru::CeedVector;
+    tmode=NOTRANSPOSE,
+    request=RequestImmediate(),
+)
+    C.CeedElemRestrictionApply(r[], tmode, u[], ru[], request[])
+end
+
+"""
+    apply(r::ElemRestriction, u::AbstractVector; tmode=NOTRANSPOSE)
+
+Use the [`ElemRestriction`](@ref) to convert from L-vector to an E-vector (or apply the
+tranpose operation). The input is given by `u`, and the result is returned as an array of
+type `Vector{CeedScalar}`.
+"""
+function apply(r::ElemRestriction, u::AbstractVector; tmode=NOTRANSPOSE)
+    ceed_ref = Ref{C.Ceed}()
+    ccall(
+        (:CeedElemRestrictionGetCeed, C.libceed),
+        Cint,
+        (C.CeedElemRestriction, Ptr{C.Ceed}),
+        r[],
+        ceed_ref,
+    )
+    c = Ceed(ceed_ref)
+    uv = CeedVector(c, u)
+    if tmode == NOTRANSPOSE
+        ruv = create_evector(r)
+    else
+        ruv = create_lvector(r)
+        ruv[] = 0.0
+    end
+    apply!(r, uv, ruv; tmode=tmode)
+    Vector(ruv)
+end
+
+"""
+    create_evector(r::ElemRestriction)
+
+Return a new [`CeedVector`](@ref) E-vector.
+"""
+function create_evector(r::ElemRestriction)
+    ref = Ref{C.CeedVector}()
+    C.CeedElemRestrictionCreateVector(r[], C_NULL, ref)
+    CeedVector(ref)
+end
+
+"""
+    create_lvector(r::ElemRestriction)
+
+Return a new [`CeedVector`](@ref) L-vector.
+"""
+function create_lvector(r::ElemRestriction)
+    ref = Ref{C.CeedVector}()
+    C.CeedElemRestrictionCreateVector(r[], ref, C_NULL)
+    CeedVector(ref)
+end
+
+"""
+    create_vectors(r::ElemRestriction)
+
+Return an (L-vector, E-vector) pair.
+"""
+function create_vectors(r::ElemRestriction)
+    l_ref = Ref{C.CeedVector}()
+    e_ref = Ref{C.CeedVector}()
+    C.CeedElemRestrictionCreateVector(r[], l_ref, e_ref)
+    CeedVector(l_ref), CeedVector(e_ref)
+end
+
+"""
+    getcompstride(r::ElemRestriction)
+
+Get the L-vector component stride.
+"""
+function getcompstride(r::ElemRestriction)
+    lsize = Ref{CeedInt}()
+    C.CeedElemRestrictionGetCompStride(r[], lsize)
+    lsize[]
+end
+
+"""
+    getnumelements(r::ElemRestriction)
+
+Get the total number of elements in the range of an [`ElemRestriction`](@ref).
+"""
+function getnumelements(r::ElemRestriction)
+    result = Ref{CeedInt}()
+    C.CeedElemRestrictionGetNumElements(r[], result)
+    result[]
+end
+
+"""
+    getelementsize(r::ElemRestriction)
+
+Get the size of elements in the given [`ElemRestriction`](@ref).
+"""
+function getelementsize(r::ElemRestriction)
+    result = Ref{CeedInt}()
+    C.CeedElemRestrictionGetElementSize(r[], result)
+    result[]
+end
+
+"""
+    getlvectorsize(r::ElemRestriction)
+
+Get the size of an L-vector for the given [`ElemRestriction`](@ref).
+"""
+function getlvectorsize(r::ElemRestriction)
+    result = Ref{CeedInt}()
+    C.CeedElemRestrictionGetLVectorSize(r[], result)
+    result[]
+end
+
+"""
+    getnumcomponents(r::ElemRestriction)
+
+Get the number of components in the elements of an [`ElemRestriction`](@ref).
+"""
+function getnumcomponents(r::ElemRestriction)
+    result = Ref{CeedInt}()
+    C.CeedElemRestrictionGetNumComponents(r[], result)
+    result[]
+end
+
+"""
+    getmultiplicity!(r::ElemRestriction, v::AbstractCeedVector)
+
+Get the multiplicity of nodes in an [`ElemRestriction`](@ref). The [`CeedVector`](@ref) `v`
+should be an L-vector (i.e. `length(v) == getlvectorsize(r)`, see [`create_lvector`](@ref)).
+"""
+function getmultiplicity!(r::ElemRestriction, v::AbstractCeedVector)
+    @assert length(v) == getlvectorsize(r)
+    C.CeedElemRestrictionGetMultiplicity(r[], v[])
+end
+
+"""
+    getmultiplicity(r::ElemRestriction)
+
+Convenience function to get the multiplicity of nodes in the [`ElemRestriction`](@ref),
+where the result is returned in a newly allocated Julia `Vector{CeedScalar}` (see also
+[`getmultiplicity!`](@ref)).
+"""
+function getmultiplicity(r::ElemRestriction)
+    v = create_lvector(r)
+    getmultiplicity!(r, v)
+    Vector(v)
+end
diff --git a/julia/LibCEED.jl/src/Globals.jl b/julia/LibCEED.jl/src/Globals.jl
new file mode 100644
index 0000000000..41eae89c93
--- /dev/null
+++ b/julia/LibCEED.jl/src/Globals.jl
@@ -0,0 +1,117 @@
+"""
+    CeedScalar
+
+Scalar (floating point) type. Equivalent to `Float64`.
+"""
+const CeedScalar = C.CeedScalar
+"""
+    CeedInt
+
+Integer type, used for indexing. Equivalent to `Int32`.
+"""
+const CeedInt = C.CeedInt
+
+"""
+    QuadMode
+
+One of `GAUSS` or `GAUSS_LOBATTO`.
+"""
+const QuadMode = C.CeedQuadMode
+const GAUSS = C.CEED_GAUSS
+const GAUSS_LOBATTO = C.CEED_GAUSS_LOBATTO
+
+"""
+    MemType
+
+One of `MEM_HOST` or `MEM_DEVICE`.
+"""
+const MemType = C.CeedMemType
+const MEM_HOST = C.CEED_MEM_HOST
+const MEM_DEVICE = C.CEED_MEM_DEVICE
+
+"""
+    CopyMode
+
+One of `COPY_VALUES`, `USE_POINTER` or `OWN_POINTER`.
+
+`OWN_POINTER` is not typically supported for objects created in Julia, because those must be
+destroyed by the garbage collector, and cannot be freed from C.
+"""
+const CopyMode = C.CeedCopyMode
+const COPY_VALUES = C.CEED_COPY_VALUES
+const USE_POINTER = C.CEED_USE_POINTER
+const OWN_POINTER = C.CEED_OWN_POINTER
+
+"""
+    EvalMode
+
+Evaluation mode used in the specification of input and output fields for Q-functions, e.g.
+in [`@interior_qf`](@ref).
+
+One of:
+- `EVAL_NONE`
+- `EVAL_INTERP`
+- `EVAL_GRAD`
+- `EVAL_DIV`
+- `EVAL_CURL`
+- `EVAL_WEIGHT`
+
+See the [libCEED
+documentation](https://libceed.readthedocs.io/en/latest/api/CeedBasis/?highlight=EVAL_MODE#c.CeedEvalMode)
+for further information.
+"""
+const EvalMode = C.CeedEvalMode
+const EVAL_NONE = C.CEED_EVAL_NONE
+const EVAL_INTERP = C.CEED_EVAL_INTERP
+const EVAL_GRAD = C.CEED_EVAL_GRAD
+const EVAL_DIV = C.CEED_EVAL_DIV
+const EVAL_CURL = C.CEED_EVAL_CURL
+const EVAL_WEIGHT = C.CEED_EVAL_WEIGHT
+
+"""
+    NormType
+
+Denotes type of vector norm. One of `NORM_1`, `NORM_2`, or `NORM_MAX`.
+"""
+const NormType = C.CeedNormType
+const NORM_1 = C.CEED_NORM_1
+const NORM_2 = C.CEED_NORM_2
+const NORM_MAX = C.CEED_NORM_MAX
+
+"""
+    TransposeMose
+
+Denotes whether a linear transformation or its transpose should be applied. Either
+`NOTRANSPOSE` or `TRANSPOSE`.
+"""
+const TransposeMode = C.CeedTransposeMode
+const NOTRANSPOSE = C.CEED_NOTRANSPOSE
+const TRANSPOSE = C.CEED_TRANSPOSE
+
+"""
+    Topology
+
+Type of basis shape to create non-tensor H1 element basis. One of `LINE`, `TRIANGLE`,
+`QUAD`, `TET`, `PYRAMID`, `PRISM`, or `HEX`.
+
+The dimension can be extracted with bitshift:
+
+    dim = Int(topology) >> 16
+"""
+const Topology = C.CeedElemTopology
+const LINE = C.CEED_LINE
+const TRIANGLE = C.CEED_TRIANGLE
+const QUAD = C.CEED_QUAD
+const TET = C.CEED_TET
+const PYRAMID = C.CEED_PYRAMID
+const PRISM = C.CEED_PRISM
+const HEX = C.CEED_HEX
+
+function set_globals()
+    @doc """
+        STRIDES_BACKEND
+
+    Indicate that the stride is determined by the backend.
+    """
+    global STRIDES_BACKEND = C.CEED_STRIDES_BACKEND[]
+end
diff --git a/julia/LibCEED.jl/src/LibCEED.jl b/julia/LibCEED.jl/src/LibCEED.jl
new file mode 100644
index 0000000000..6202d7c0e8
--- /dev/null
+++ b/julia/LibCEED.jl/src/LibCEED.jl
@@ -0,0 +1,150 @@
+module LibCEED
+
+using StaticArrays, UnsafeArrays, Requires
+using Base: RefValue
+
+# import low-level C interface
+include("C.jl")
+import .C
+
+export @interior_qf,
+    @witharray,
+    @witharray_read,
+    Abscissa,
+    AbscissaAndWeights,
+    Basis,
+    BasisCollocated,
+    COPY_VALUES,
+    Ceed,
+    CeedDim,
+    CeedInt,
+    CeedScalar,
+    CeedVector,
+    CeedVectorActive,
+    CeedVectorNone,
+    Context,
+    CopyMode,
+    EVAL_CURL,
+    EVAL_DIV,
+    EVAL_GRAD,
+    EVAL_INTERP,
+    EVAL_NONE,
+    EVAL_WEIGHT,
+    ElemRestriction,
+    ElemRestrictionNone,
+    EvalMode,
+    GAUSS,
+    GAUSS_LOBATTO,
+    HEX,
+    LINE,
+    MEM_DEVICE,
+    MEM_HOST,
+    MemType,
+    NORM_1,
+    NORM_2,
+    NORM_MAX,
+    NOTRANSPOSE,
+    NormType,
+    OWN_POINTER,
+    Operator,
+    PRISM,
+    PYRAMID,
+    QFunction,
+    QFunctionNone,
+    QUAD,
+    QuadMode,
+    RequestImmediate,
+    RequestOrdered,
+    STRIDES_BACKEND,
+    TET,
+    TRANSPOSE,
+    TRIANGLE,
+    Topology,
+    TransposeMode,
+    USE_POINTER,
+    UserQFunction,
+    add_input!,
+    add_output!,
+    apply!,
+    apply_add!,
+    apply,
+    assemble,
+    assemble_add_diagonal!,
+    assemble_diagonal!,
+    create_composite_operator,
+    create_elem_restriction,
+    create_elem_restriction_strided,
+    create_evector,
+    create_h1_basis,
+    create_identity_qfunction,
+    create_interior_qfunction,
+    create_lvector,
+    create_tensor_h1_basis,
+    create_tensor_h1_lagrange_basis,
+    create_vectors,
+    det,
+    extract_array,
+    extract_context,
+    gauss_quadrature,
+    get_preferred_memtype,
+    getcompstride,
+    getnumelements,
+    getelementsize,
+    getlvectorsize,
+    getmultiplicity!,
+    getmultiplicity,
+    getdimension,
+    getgrad,
+    getgrad1d,
+    getinterp,
+    getinterp1d,
+    getnumcomponents,
+    getnumnodes,
+    getnumnodes1d,
+    getnumqpts,
+    getnumqpts1d,
+    getqref,
+    getqweights,
+    getresource,
+    gettopology,
+    getvoigt!,
+    getvoigt,
+    iscuda,
+    isdeterministic,
+    lobatto_quadrature,
+    norm,
+    reciprocal!,
+    set_context!,
+    set_cufunction!,
+    set_data!,
+    set_field!,
+    setarray!,
+    setvalue!,
+    setvoigt!,
+    setvoigt,
+    syncarray!,
+    takearray!,
+    witharray,
+    witharray_read
+
+include("Globals.jl")
+include("Ceed.jl")
+include("CeedVector.jl")
+include("Basis.jl")
+include("ElemRestriction.jl")
+include("Quadrature.jl")
+include("Context.jl")
+include("UserQFunction.jl")
+include("QFunction.jl")
+include("Request.jl")
+include("Operator.jl")
+include("Misc.jl")
+
+cuda_is_loaded = false
+
+function __init__()
+    @require CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" include("Cuda.jl")
+    set_globals()
+end
+
+end # module
diff --git a/julia/LibCEED.jl/src/Misc.jl b/julia/LibCEED.jl/src/Misc.jl
new file mode 100644
index 0000000000..1ace1ad1b2
--- /dev/null
+++ b/julia/LibCEED.jl/src/Misc.jl
@@ -0,0 +1,119 @@
+import LinearAlgebra: det
+
+"""
+    CeedDim(dim)
+
+The singleton object of type `CeedDim{dim}`, used for dispatch to linear algebra operations
+specialized for small matrices (1, 2, or 3 dimensions).
+"""
+struct CeedDim{dim} end
+@inline CeedDim(dim) = CeedDim{Int(dim)}()
+
+"""
+    det(J, ::CeedDim{dim})
+
+Specialized determinant calculations for matrices of size 1, 2, or 3.
+"""
+@inline det(J, ::CeedDim{1}) = @inbounds J[1]
+@inline det(J, ::CeedDim{2}) = @inbounds J[1]*J[4] - J[3]*J[2]
+#! format: off
+@inline det(J, ::CeedDim{3}) = @inbounds (
+    J[1]*(J[5]*J[9] - J[6]*J[8]) -
+    J[2]*(J[4]*J[9] - J[6]*J[7]) +
+    J[3]*(J[4]*J[8] - J[5]*J[7])
+)
+#! format: on
+
+"""
+    setvoigt(J::StaticArray{Tuple{D,D},T,2})
+    setvoigt(J, ::CeedDim{dim})
+
+Given a symmetric matrix `J`, return a `SVector` that encodes `J` using the [Voigt
+convention](https://en.wikipedia.org/wiki/Voigt_notation).
+
+The size of the symmetric matrix `J` must be known statically, either specified using
+[`CeedDim`](@ref) or `StaticArray`.
+"""
+@inline setvoigt(J::StaticArray{Tuple{D,D}}) where {D} = setvoigt(J, CeedDim(D))
+@inline setvoigt(J, ::CeedDim{1}) = @inbounds @SVector [J[1]]
+@inline setvoigt(J, ::CeedDim{2}) = @inbounds @SVector [J[1], J[4], J[2]]
+@inline setvoigt(J, ::CeedDim{3}) = @inbounds @SVector [J[1], J[5], J[9], J[6], J[3], J[2]]
+
+@inline function setvoigt!(V, J, ::CeedDim{1})
+    @inbounds V[1] = J[1]
+end
+
+@inline function setvoigt!(V, J, ::CeedDim{2})
+    @inbounds begin
+        V[1] = J[1]
+        V[2] = J[4]
+        V[3] = J[2]
+    end
+end
+
+@inline function setvoigt!(V, J, ::CeedDim{3})
+    @inbounds begin
+        V[1] = J[1]
+        V[2] = J[5]
+        V[3] = J[9]
+        V[4] = J[6]
+        V[5] = J[3]
+        V[6] = J[2]
+    end
+end
+
+"""
+    getvoigt(V, ::CeedDim{dim})
+
+Given a vector `V` that encodes a symmetric matrix using the [Voigt
+convention](https://en.wikipedia.org/wiki/Voigt_notation), return the corresponding
+`SMatrix`.
+"""
+@inline getvoigt(V, ::CeedDim{1}) = @inbounds @SMatrix [V[1]]
+@inline getvoigt(V, ::CeedDim{2}) = @inbounds @SMatrix [V[1] V[3]; V[3] V[2]]
+@inline getvoigt(V, ::CeedDim{3}) = @inbounds @SMatrix [
+    V[1] V[6] V[5]
+    V[6] V[2] V[4]
+    V[5] V[4] V[3]
+]
+@inline getvoigt(V::StaticArray{Tuple{1}}) = getvoigt(V, CeedDim(1))
+@inline getvoigt(V::StaticArray{Tuple{3}}) = getvoigt(V, CeedDim(2))
+@inline getvoigt(V::StaticArray{Tuple{6}}) = getvoigt(V, CeedDim(3))
+
+@inline function getvoigt!(J, V, ::CeedDim{1})
+    @inbounds J[1, 1] = V[1]
+end
+
+@inline function getvoigt!(J, V, ::CeedDim{2})
+    @inbounds begin
+        #! format: off
+        J[1,1] = V[1] ; J[1,2] = V[3]
+        J[2,1] = V[3] ; J[2,2] = V[2]
+        #! format: on
+    end
+end
+
+@inline function getvoigt!(J, V, ::CeedDim{3})
+    @inbounds begin
+        #! format: off
+        J[1,1] = V[1] ; J[1,2] = V[6] ; J[1,3] = V[5]
+        J[2,1] = V[6] ; J[2,2] = V[2] ; J[2,3] = V[4]
+        J[3,1] = V[5] ; J[3,2] = V[4] ; J[3,3] = V[3]
+        #! format: on
+    end
+end
+
+function tmp_view(obj, view_fn)
+    str = mktemp() do fname, f
+        cf = Libc.FILE(f)
+        er = view_fn(obj, cf.ptr)
+        ccall(:fflush, Cint, (Ptr{Cvoid},), cf)
+        seek(f, 0)
+        read(f, String)
+    end
+    chomp(str)
+end
+
+function ceed_show(io::IO, obj, view_fn)
+    print(io, tmp_view(obj[], view_fn))
+end
diff --git a/julia/LibCEED.jl/src/Operator.jl b/julia/LibCEED.jl/src/Operator.jl
new file mode 100644
index 0000000000..74fa56297e
--- /dev/null
+++ b/julia/LibCEED.jl/src/Operator.jl
@@ -0,0 +1,149 @@
+mutable struct Operator
+    ref::RefValue{C.CeedOperator}
+    qf::AbstractQFunction
+    dqf::AbstractQFunction
+    dqfT::AbstractQFunction
+    sub_ops::Vector{Operator}
+    function Operator(ref, qf, dqf, dqfT)
+        obj = new(ref, qf, dqf, dqfT, [])
+        finalizer(obj) do x
+            # ccall(:jl_safe_printf, Cvoid, (Cstring, Cstring), "Finalizing %s.\n", repr(x))
+            destroy(x)
+        end
+        return obj
+    end
+end
+destroy(op::Operator) = C.CeedOperatorDestroy(op.ref) # COV_EXCL_LINE
+Base.getindex(op::Operator) = op.ref[]
+Base.show(io::IO, ::MIME"text/plain", op::Operator) = ceed_show(io, op, C.CeedOperatorView)
+
+"""
+    Operator(ceed::Ceed; qf, dqf=QFunctionNone(), dqfT=QFunctionNone(), fields)
+
+Creates a libCEED `CeedOperator` object using the given Q-function `qf`, and optionally its
+derivative and derivative transpose.
+
+An array of fields must be provided, where each element of the array is a tuple containing
+the name of the field (as a string or symbol), the corresponding element restriction, basis,
+and vector.
+
+# Examples
+
+Create the operator that builds the Q-data associated with the mass matrix.
+```
+build_oper = Operator(
+    ceed,
+    qf=build_qfunc,
+    fields=[
+        (:J, mesh_restr, mesh_basis, CeedVectorActive()),
+        (:w, ElemRestrictionNone(), mesh_basis, CeedVectorNone()),
+        (:qdata, sol_restr_i, BasisCollocated(), CeedVectorActive())
+    ]
+)
+```
+"""
+function Operator(c::Ceed; qf, dqf=QFunctionNone(), dqfT=QFunctionNone(), fields)
+    op = Operator(c, qf, dqf, dqfT)
+    for f ∈ fields
+        set_field!(op, String(f[1]), f[2], f[3], f[4])
+    end
+    op
+end
+
+function Operator(
+    c::Ceed,
+    qf::AbstractQFunction,
+    dqf::AbstractQFunction,
+    dqfT::AbstractQFunction,
+)
+    ref = Ref{C.CeedOperator}()
+    C.CeedOperatorCreate(c[], qf[], dqf[], dqfT[], ref)
+    Operator(ref, qf, dqf, dqfT)
+end
+
+"""
+    create_composite_operator(c::Ceed, ops)
+
+Create an [`Operator`](@ref) whose action represents the sum of the operators in the
+collection `ops`.
+"""
+function create_composite_operator(c::Ceed, ops)
+    ref = Ref{C.CeedOperator}()
+    C.CeedCompositeOperatorCreate(c[], ref)
+    comp_op = Operator(ref, QFunctionNone(), QFunctionNone(), QFunctionNone())
+    comp_op.sub_ops = ops
+    for op ∈ ops
+        C.CeedCompositeOperatorAddSub(comp_op[], op[])
+    end
+    comp_op
+end
+
+function set_field!(
+    op::Operator,
+    fieldname::AbstractString,
+    r::AbstractElemRestriction,
+    b::AbstractBasis,
+    v::AbstractCeedVector,
+)
+    C.CeedOperatorSetField(op[], fieldname, r[], b[], v[])
+end
+
+"""
+    apply!(op::Operator, vin, vout; request=RequestImmediate())
+
+Apply the action of the operator `op` to the input vector `vin`, and store the result in the
+output vector `vout`.
+
+For non-blocking application, the user can specify a request object. By default, immediate
+(synchronous) completion is requested.
+"""
+function apply!(
+    op::Operator,
+    vin::AbstractCeedVector,
+    vout::AbstractCeedVector;
+    request=RequestImmediate(),
+)
+    C.CeedOperatorApply(op[], vin[], vout[], request[])
+end
+
+"""
+    apply_add!(op::Operator, vin, vout; request=RequestImmediate())
+
+Apply the action of the operator `op` to the input vector `vin`, and add the result to the
+output vector `vout`.
+
+For non-blocking application, the user can specify a request object. By default, immediate
+(synchronous) completion is requested.
+"""
+function apply_add!(
+    op::Operator,
+    vin::AbstractCeedVector,
+    vout::AbstractCeedVector;
+    request=RequestImmediate(),
+)
+    C.CeedOperatorApplyAdd(op[], vin[], vout[], request[])
+end
+
+"""
+    assemble_diagonal!(op::Operator, diag::CeedVector; request=RequestImmediate())
+
+Overwrites a [`CeedVector`](@ref) with the diagonal of a linear [`Operator`](@ref).
+
+!!! note "Note:"
+    Currently only [`Operator`](@ref)s with a single field are supported.
+"""
+function assemble_diagonal!(op::Operator, diag::CeedVector; request=RequestImmediate())
+    C.CeedOperatorLinearAssembleDiagonal(op[], diag[], request[])
+end
+
+"""
+    assemble_diagonal!(op::Operator, diag::CeedVector; request=RequestImmediate())
+
+Adds the diagonal of a linear [`Operator`](@ref) to the given [`CeedVector`](@ref).
+
+!!! note "Note:"
+    Currently only [`Operator`](@ref)s with a single field are supported.
+"""
+function assemble_add_diagonal!(op::Operator, diag::CeedVector; request=RequestImmediate())
+    C.CeedOperatorLinearAssembleAddDiagonal(op[], diag[], request[])
+end
diff --git a/julia/LibCEED.jl/src/QFunction.jl b/julia/LibCEED.jl/src/QFunction.jl
new file mode 100644
index 0000000000..da5c769cc9
--- /dev/null
+++ b/julia/LibCEED.jl/src/QFunction.jl
@@ -0,0 +1,148 @@
+abstract type AbstractQFunction end
+
+struct QFunctionNone <: AbstractQFunction end
+Base.getindex(::QFunctionNone) = C.CEED_QFUNCTION_NONE[]
+
+"""
+    QFunction
+
+A libCEED `CeedQFunction` object, typically created using the [`@interior_qf`](@ref) macro.
+
+A `QFunction` can also be created from the "Q-function gallery" using
+[`create_interior_qfunction`](@ref). The identity Q-function can be created using
+[`create_identity_qfunction`](@ref).
+"""
+mutable struct QFunction <: AbstractQFunction
+    ref::RefValue{C.CeedQFunction}
+    user_qf::Union{Nothing,UserQFunction}
+    ctx::Union{Nothing,Context}
+    function QFunction(ref, user_qf)
+        obj = new(ref, user_qf, nothing)
+        finalizer(obj) do x
+            # ccall(:jl_safe_printf, Cvoid, (Cstring, Cstring), "Finalizing %s.\n", repr(x))
+            destroy(x)
+        end
+        return obj
+    end
+end
+QFunction(ref::Ref{C.CeedQFunction}) = QFunction(ref, nothing)
+destroy(qf::QFunction) = C.CeedQFunctionDestroy(qf.ref) # COV_EXCL_LINE
+Base.getindex(qf::QFunction) = qf.ref[]
+Base.show(io::IO, ::MIME"text/plain", qf::QFunction) =
+    ceed_show(io, qf, C.CeedQFunctionView)
+
+function create_interior_qfunction(c::Ceed, f::UserQFunction; vlength=1)
+    ref = Ref{C.CeedQFunction}()
+    # Use empty string as source location to indicate to libCEED that there is
+    # no C source for this Q-function
+    C.CeedQFunctionCreateInterior(c[], vlength, f.fptr, "", ref)
+    # COV_EXCL_START
+    if !isnothing(f.cuf)
+        C.CeedQFunctionSetCUDAUserFunction(ref[], f.cuf)
+    elseif iscuda(c) && !cuda_is_loaded
+        error(string(
+            "In order to use user Q-functions with a CUDA backend, the CUDA.jl package ",
+            "must be loaded",
+        ))
+    end
+    # COV_EXCL_STOP
+    QFunction(ref, f)
+end
+
+"""
+    create_interior_qfunction(ceed::Ceed, name::AbstractString)
+
+Create a [`QFunction`](@ref) from the Q-function gallery, using the provided name.
+
+# Examples
+
+- Build and apply the 3D mass operator
+```
+build_mass_qf = create_interior_qfunction(c, "Mass3DBuild")
+apply_mass_qf = create_interior_qfunction(c, "MassApply")
+```
+- Build and apply the 3D Poisson operator
+```
+build_poi_qf = create_interior_qfunction(c, "Poisson3DBuild")
+apply_poi_qf = create_interior_qfunction(c, "Poisson3DApply")
+```
+"""
+function create_interior_qfunction(c::Ceed, name::AbstractString)
+    ref = Ref{C.CeedQFunction}()
+    C.CeedQFunctionCreateInteriorByName(c.ref[], name, ref)
+    QFunction(ref)
+end
+
+"""
+    create_identity_qfunction(c::Ceed, size, inmode::EvalMode, outmode::EvalMode)
+
+Create an identity [`QFunction`](@ref). Inputs are written into outputs in the order given.
+This is useful for [`Operators`](@ref Operator) that can be represented with only the action
+of a [`ElemRestriction`](@ref) and [`Basis`](@ref), such as restriction and prolongation
+operators for p-multigrid. Backends may optimize `CeedOperators` with this Q-function to
+avoid the copy of input data to output fields by using the same memory location for both.
+"""
+function create_identity_qfunction(c::Ceed, size, inmode::EvalMode, outmode::EvalMode)
+    ref = Ref{C.CeedQFunction}()
+    C.CeedQFunctionCreateIdentity(c[], size, inmode, outmode, ref)
+    QFunction(ref)
+end
+
+function add_input!(qf::AbstractQFunction, name::AbstractString, size, emode)
+    C.CeedQFunctionAddInput(qf[], name, size, emode)
+end
+
+function add_output!(qf::AbstractQFunction, name::AbstractString, size, emode)
+    C.CeedQFunctionAddOutput(qf[], name, size, emode)
+end
+
+"""
+    set_context!(qf::QFunction, ctx::Context)
+
+Associate a [`Context`](@ref) object `ctx` with the given Q-function `qf`.
+"""
+function set_context!(qf::QFunction, ctx)
+    # Preserve the context data from the GC by storing a reference
+    qf.ctx = ctx
+    C.CeedQFunctionSetContext(qf[], ctx[])
+end
+
+function get_field_sizes(qf::AbstractQFunction)
+    ninputs = Ref{CeedInt}()
+    noutputs = Ref{CeedInt}()
+
+    C.CeedQFunctionGetNumArgs(qf[], ninputs, noutputs)
+
+    inputs = Ref{Ptr{C.CeedQFunctionField}}()
+    outputs = Ref{Ptr{C.CeedQFunctionField}}()
+    C.CeedQFunctionGetFields(qf[], inputs, outputs)
+
+    input_sizes = zeros(CeedInt, ninputs[])
+    output_sizes = zeros(CeedInt, noutputs[])
+
+    for i = 1:ninputs[]
+        field = unsafe_load(inputs[], i)
+        C.CeedQFunctionFieldGetSize(field, pointer(input_sizes, i))
+    end
+
+    for i = 1:noutputs[]
+        field = unsafe_load(outputs[], i)
+        C.CeedQFunctionFieldGetSize(field, pointer(output_sizes, i))
+    end
+
+    input_sizes, output_sizes
+end
+
+"""
+    apply!(qf::QFunction, Q, vin, vout)
+
+Apply the action of a [`QFunction`](@ref) to an array of input vectors, and store the result
+in an array of output vectors.
+"""
+function apply!(qf::QFunction, Q, vin, vout)
+    vins = map(x -> x[], vin)
+    vouts = map(x -> x[], vout)
+    GC.@preserve vin vout begin
+        C.CeedQFunctionApply(qf[], Q, pointer(vins), pointer(vouts))
+    end
+end
diff --git a/julia/LibCEED.jl/src/Quadrature.jl b/julia/LibCEED.jl/src/Quadrature.jl
new file mode 100644
index 0000000000..610fcbbe64
--- /dev/null
+++ b/julia/LibCEED.jl/src/Quadrature.jl
@@ -0,0 +1,37 @@
+@doc raw"""
+    gauss_quadrature(q)
+
+Return the Gauss-Legendre quadrature rule with `q` points (integrates polynomials of degree
+$2q-1$ exactly).
+
+A tuple `(x,w)` is returned.
+"""
+function gauss_quadrature(q)
+    x = zeros(CeedScalar, q)
+    w = zeros(CeedScalar, q)
+    C.CeedGaussQuadrature(q, x, w)
+    x, w
+end
+
+struct QuadratureMode{T} end
+const Abscissa = QuadratureMode{:Abscissa}()
+const AbscissaAndWeights = QuadratureMode{:AbscissaAndWeights}()
+
+@doc raw"""
+    lobatto_quadrature(q, mode::Mode=Abscissa)
+
+Return the Gauss-Lobatto quadrature rule with `q` points (integrates polynomials of degree
+$2q-3$ exactly).
+
+If `mode` is `AbscissaAndWeights`, then both the weights and abscissa are returned as a
+tuple `(x,w)`.
+
+Otherwise, (if `mode` is `Abscissa`), then only the abscissa `x` are returned.
+"""
+function lobatto_quadrature(q, mode::Mode=Abscissa) where {Mode}
+    return_weights = (mode == AbscissaAndWeights)
+    x = zeros(CeedScalar, q)
+    w = (return_weights) ? zeros(CeedScalar, q) : C_NULL
+    C.CeedLobattoQuadrature(q, x, w)
+    return_weights ? (x, w) : x
+end
diff --git a/julia/LibCEED.jl/src/Request.jl b/julia/LibCEED.jl/src/Request.jl
new file mode 100644
index 0000000000..61ea58e17d
--- /dev/null
+++ b/julia/LibCEED.jl/src/Request.jl
@@ -0,0 +1,20 @@
+abstract type AbstractRequest end
+
+struct RequestImmediate <: AbstractRequest end
+Base.getindex(::RequestImmediate) = C.CEED_REQUEST_IMMEDIATE[]
+
+struct RequestOrdered <: AbstractRequest end
+Base.getindex(::RequestOrdered) = C.CEED_REQUEST_ORDERED[]
+
+#=
+# CeedRequest is not fully implemented in libCEED. When it is implemented, the
+# following can be used as a starting point for the Julia interface.
+
+struct Request <: AbstractRequest
+    ref::RefValue{C.CeedRequest}
+end
+
+Request() = Request(Ref{C.CeedRequest}())
+
+Base.wait(req::Request) = C.CeedRequestWait(req[])
+=#
diff --git a/julia/LibCEED.jl/src/UserQFunction.jl b/julia/LibCEED.jl/src/UserQFunction.jl
new file mode 100644
index 0000000000..17d9603207
--- /dev/null
+++ b/julia/LibCEED.jl/src/UserQFunction.jl
@@ -0,0 +1,278 @@
+struct UserQFunction{F,K}
+    f::F
+    fptr::Ptr{Nothing}
+    kf::K
+    cuf::Union{Nothing,Ptr{Nothing}}
+end
+
+@inline function extract_context(ptr, ::Type{T}) where {T}
+    unsafe_load(Ptr{T}(ptr))
+end
+
+@inline function extract_array(ptr, idx, dims)
+    UnsafeArray(Ptr{CeedScalar}(unsafe_load(ptr, idx)), dims)
+end
+
+function generate_user_qfunction(
+    ceed,
+    def_module,
+    qf_name,
+    constants,
+    array_names,
+    ctx,
+    dims_in,
+    dims_out,
+    body,
+)
+    idx = gensym(:i)
+    Q = gensym(:Q)
+    ctx_ptr = gensym(:ctx_ptr)
+    in_ptr = gensym(:in_ptr)
+    out_ptr = gensym(:out_ptr)
+
+    const_assignments = Vector{Expr}(undef, length(constants))
+    for (i, c) ∈ enumerate(constants)
+        const_assignments[i] = :($(c[1]) = $(c[2]))
+    end
+
+    narrays = length(array_names)
+    arrays = Vector{Expr}(undef, narrays)
+    array_views = Vector{Expr}(undef, narrays)
+    n_in = length(dims_in)
+    for (i, arr_name) ∈ enumerate(array_names)
+        i_inout = (i <= n_in) ? i : i - n_in
+        dims = (i <= n_in) ? dims_in[i] : dims_out[i-n_in]
+        ptr = (i <= n_in) ? in_ptr : out_ptr
+        arr_name_gen = gensym(arr_name)
+        arrays[i] = :($arr_name_gen = extract_array($ptr, $i_inout, (Int($Q), $(dims...))))
+        ndims = length(dims)
+        slice = Expr(:ref, arr_name_gen, idx, (:(:) for i = 1:ndims)...)
+        if i <= n_in
+            if ndims == 0
+                array_views[i] = :($arr_name = $slice)
+            else
+                S = Tuple{dims...}
+                array_views[i] = :($arr_name = LibCEED.SArray{$S}(@view $slice))
+            end
+        else
+            array_views[i] = :($arr_name = @view $slice)
+        end
+    end
+
+    if isnothing(ctx)
+        ctx_assignment = nothing
+    else
+        ctx_assignment = :($(ctx.name) = extract_context($ctx_ptr, $(ctx.type)))
+    end
+
+    qf1 = gensym(qf_name)
+    f = Core.eval(
+        def_module,
+        quote
+            @inline function $qf1(
+                $ctx_ptr::Ptr{Cvoid},
+                $Q::CeedInt,
+                $in_ptr::Ptr{Ptr{CeedScalar}},
+                $out_ptr::Ptr{Ptr{CeedScalar}},
+            )
+                $(const_assignments...)
+                $ctx_assignment
+                $(arrays...)
+                @inbounds @simd for $idx = 1:$Q
+                    $(array_views...)
+                    $body
+                end
+                CeedInt(0)
+            end
+        end,
+    )
+    f_qn = QuoteNode(f)
+    rt = :CeedInt
+    at = :(Core.svec(Ptr{Cvoid}, CeedInt, Ptr{Ptr{CeedScalar}}, Ptr{Ptr{CeedScalar}}))
+    fptr = eval(Expr(:cfunction, Ptr{Cvoid}, f_qn, rt, at, QuoteNode(:ccall)))
+
+    # COV_EXCL_START
+    if iscuda(ceed)
+        getresource(ceed) == "/gpu/cuda/gen" && error(string(
+            "/gpu/cuda/gen is not compatible with user Q-functions defined with ",
+            "libCEED.jl.\nPlease use a different backend, for example: /gpu/cuda/shared ",
+            "or /gpu/cuda/ref",
+        ))
+        if cuda_is_loaded
+            !has_cuda() && error("No valid CUDA installation found")
+            qf2 = gensym(qf_name)
+            kf = Core.eval(
+                def_module,
+                quote
+                    @inline function $qf2($ctx_ptr::Ptr{Cvoid}, $(array_names...))
+                        $(const_assignments...)
+                        $ctx_assignment
+                        $body
+                        nothing
+                    end
+                end,
+            )
+            cuf = mk_cufunction(ceed, def_module, qf_name, kf, dims_in, dims_out)
+        else
+            error(string(
+                "User Q-functions with CUDA backends require the CUDA.jl package to be ",
+                "loaded.\nThe libCEED backend is: $(getresource(ceed))\n",
+                "Please ensure that the CUDA.jl package is installed and loaded.",
+            ))
+        end
+    else
+        kf = nothing
+        cuf = nothing
+    end
+    # COV_EXCL_STOP
+
+    UserQFunction(f, fptr, kf, cuf)
+end
+
+function meta_user_qfunction(ceed, def_module, qf, args)
+    qf_name = Meta.quot(qf)
+
+    ctx = nothing
+    constants = Expr[]
+    dims_in = Expr[]
+    dims_out = Expr[]
+    names_in = Symbol[]
+    names_out = Symbol[]
+
+    for a ∈ args[1:end-1]
+        if Meta.isexpr(a, :(=))
+            a1 = Meta.quot(a.args[1])
+            a2 = esc(a.args[2])
+            push!(constants, :(($a1, $a2)))
+        elseif Meta.isexpr(a, :tuple)
+            arr_name = a.args[1]
+            inout = a.args[2].value
+            ndim = length(a.args) - 3
+            dims = Vector{Expr}(undef, ndim)
+            for d = 1:ndim
+                dims[d] = :(Int($(a.args[d+3])))
+            end
+            dims_expr = :(Int[$(esc.(a.args[4:end])...)])
+            if inout == :in
+                push!(dims_in, dims_expr)
+                push!(names_in, arr_name)
+            elseif inout == :out
+                push!(dims_out, dims_expr)
+                push!(names_out, arr_name)
+            else
+                error("Array specification must be either :in or :out. Given $inout.")
+            end
+        elseif Meta.isexpr(a, :(::))
+            ctx = (name=a.args[1], type=a.args[2])
+        else
+            error("Bad argument to @user_qfunction")
+        end
+    end
+
+    body = Meta.quot(args[end])
+
+    return :(generate_user_qfunction(
+        $ceed,
+        $def_module,
+        $qf_name,
+        [$(constants...)],
+        $([names_in; names_out]),
+        $ctx,
+        [$(dims_in...)],
+        [$(dims_out...)],
+        $body,
+    ))
+end
+
+"""
+    @interior_qf name=def
+
+Creates a user-defined interior (volumetric) Q-function, and assigns it to a variable named
+`name`. The definition of the Q-function is given as:
+```
+@interior_qf user_qf=(
+    ceed::CEED,
+    [const1=val1, const2=val2, ...],
+    [ctx::ContextType],
+    (I1, :in, EvalMode, dims...),
+    (I2, :in, EvalMode, dims...),
+    (O1, :out, EvalMode, dims...),
+    body
+)
+```
+The definitions of form `const=val` are used for definitions which will be compile-time
+constants in the Q-function. For example, if `dim` is a variable set to the dimension of the
+problem, then `dim=dim` will make `dim` available in the body of the Q-function as a
+compile-time constant.
+
+If the user wants to provide a context struct to the Q-function, that can be achieved by
+optionally including `ctx::ContextType`, where `ContextType` is the type of the context
+struct, and `ctx` is the name to which is will be bound in the body of the Q-function.
+
+This is followed by the definition of the input and output arrays, which take the form
+`(arr_name, (:in|:out), EvalMode, dims...)`. Each array will be bound to a variable named
+`arr_name`. Input arrays should be tagged with :in, and output arrays with :out. An
+`EvalMode` should be specified, followed by the dimensions of the array. If the array
+consists of scalars (one number per Q-point) then `dims` should be omitted.
+
+# Examples
+
+- Q-function to compute the "Q-data" for the mass operator, which is given by the quadrature
+  weight times the Jacobian determinant. The mesh Jacobian (the gradient of the nodal mesh
+  points) and the quadrature weights are given as input arrays, and the Q-data is the output
+  array. `dim` is given as a compile-time constant, and so the array `J` is statically
+  sized, and therefore `det(J)` will automatically dispatch to an optimized implementation
+  for the given dimension.
+```
+@interior_qf build_qfunc = (
+    ceed, dim=dim,
+    (J, :in, EVAL_GRAD, dim, dim),
+    (w, :in, EVAL_WEIGHT),
+    (qdata, :out, EVAL_NONE),
+    qdata[] = w*det(J)
+)
+```
+"""
+macro interior_qf(args)
+    if !Meta.isexpr(args, :(=))
+        error("@interior_qf must be of form `qf = (body)`") # COV_EXCL_LINE
+    end
+
+    qf = args.args[1]
+    user_qf = esc(qf)
+    args = args.args[2].args
+    ceed = esc(args[1])
+
+    # Calculate field sizes
+    fields_in = Expr[]
+    fields_out = Expr[]
+    for a ∈ args
+        if Meta.isexpr(a, :tuple)
+            field_name = String(a.args[1])
+            inout = a.args[2].value
+            evalmode = a.args[3]
+            ndim = length(a.args) - 3
+            dims = Vector{Expr}(undef, ndim)
+            for d = 1:ndim
+                dims[d] = esc(:(Int($(a.args[d+3]))))
+            end
+            sz_expr = :(prod(($(dims...),)))
+            if inout == :in
+                push!(fields_in, :(add_input!($user_qf, $field_name, $sz_expr, $evalmode)))
+            elseif inout == :out
+                push!(
+                    fields_out,
+                    :(add_output!($user_qf, $field_name, $sz_expr, $evalmode)),
+                )
+            end
+        end
+    end
+
+    gen_user_qf = meta_user_qfunction(ceed, __module__, qf, args[2:end])
+
+    quote
+        $user_qf = create_interior_qfunction($ceed, $gen_user_qf)
+        $(fields_in...)
+        $(fields_out...)
+    end
+end
diff --git a/julia/LibCEED.jl/src/generated/generate_bindings.jl b/julia/LibCEED.jl/src/generated/generate_bindings.jl
new file mode 100644
index 0000000000..d93ca73a2d
--- /dev/null
+++ b/julia/LibCEED.jl/src/generated/generate_bindings.jl
@@ -0,0 +1,28 @@
+using Clang
+using Clang.LibClang.Clang_jll
+
+"""
+    generate_ceed_wrapper(ceed_path)
+
+Generate Julia bindings for the libCEED library, where `ceed_path` is the path
+of the libCEED directory. The generated bindings are used (with some manual
+modifications) to create the low-level C interface LibCEED.C for the LibCEED.jl
+package.
+"""
+function generate_ceed_wrapper(ceed_path)
+    headers = ["ceed.h", "ceed-cuda.h", "ceed-backend.h"]
+    ceed_include = joinpath(ceed_path, "include")
+    ceed_headers = [joinpath(ceed_include, header) for header in headers]
+
+    wc = init(;
+        headers=ceed_headers,
+        output_file=joinpath(@__DIR__, "libceed_api_gen.jl"),
+        common_file=joinpath(@__DIR__, "libceed_common_gen.jl"),
+        clang_includes=[ceed_include, CLANG_INCLUDE],
+        clang_args=map(x -> "-I"*x, find_std_headers()),
+        header_wrapped=(root, current) -> root == current,
+        header_library=x -> "libceed",
+        clang_diagnostics=true,
+    )
+    run(wc, false)
+end
diff --git a/julia/LibCEED.jl/src/generated/libceed_api.jl b/julia/LibCEED.jl/src/generated/libceed_api.jl
new file mode 100644
index 0000000000..a647d5a434
--- /dev/null
+++ b/julia/LibCEED.jl/src/generated/libceed_api.jl
@@ -0,0 +1,773 @@
+# Julia wrapper for header: ceed.h
+# Automatically generated using Clang.jl
+#! format: off
+
+function CeedInit(resource, ceed)
+    ccall((:CeedInit, libceed), Cint, (Cstring, Ptr{Ceed}), resource, ceed)
+end
+
+function CeedGetResource(ceed, resource)
+    ccall((:CeedGetResource, libceed), Cint, (Ceed, Ptr{Cstring}), ceed, resource)
+end
+
+function CeedIsDeterministic(ceed, isDeterministic)
+    ccall((:CeedIsDeterministic, libceed), Cint, (Ceed, Ptr{Bool}), ceed, isDeterministic)
+end
+
+function CeedView(ceed, stream)
+    ccall((:CeedView, libceed), Cint, (Ceed, Ptr{FILE}), ceed, stream)
+end
+
+function CeedDestroy(ceed)
+    ccall((:CeedDestroy, libceed), Cint, (Ptr{Ceed},), ceed)
+end
+
+function CeedErrorReturn(arg1, arg2, arg3, arg4, arg5, arg6, arg7)
+    ccall((:CeedErrorReturn, libceed), Cint, (Ceed, Cstring, Cint, Cstring, Cint, Cstring, Ptr{Cvoid}), arg1, arg2, arg3, arg4, arg5, arg6, arg7)
+end
+
+function CeedErrorStore(arg1, arg2, arg3, arg4, arg5, arg6, arg7)
+    ccall((:CeedErrorStore, libceed), Cint, (Ceed, Cstring, Cint, Cstring, Cint, Cstring, Ptr{Cvoid}), arg1, arg2, arg3, arg4, arg5, arg6, arg7)
+end
+
+function CeedErrorAbort(arg1, arg2, arg3, arg4, arg5, arg6, arg7)
+    ccall((:CeedErrorAbort, libceed), Cint, (Ceed, Cstring, Cint, Cstring, Cint, Cstring, Ptr{Cvoid}), arg1, arg2, arg3, arg4, arg5, arg6, arg7)
+end
+
+function CeedErrorExit(arg1, arg2, arg3, arg4, arg5, arg6, arg7)
+    ccall((:CeedErrorExit, libceed), Cint, (Ceed, Cstring, Cint, Cstring, Cint, Cstring, Ptr{Cvoid}), arg1, arg2, arg3, arg4, arg5, arg6, arg7)
+end
+
+function CeedSetErrorHandler(ceed, eh)
+    ccall((:CeedSetErrorHandler, libceed), Cint, (Ceed, Ptr{Cvoid}), ceed, eh)
+end
+
+function CeedGetErrorMessage(arg1, errmsg)
+    ccall((:CeedGetErrorMessage, libceed), Cint, (Ceed, Ptr{Cstring}), arg1, errmsg)
+end
+
+function CeedResetErrorMessage(arg1, errmsg)
+    ccall((:CeedResetErrorMessage, libceed), Cint, (Ceed, Ptr{Cstring}), arg1, errmsg)
+end
+
+function CeedGetPreferredMemType(ceed, type)
+    ccall((:CeedGetPreferredMemType, libceed), Cint, (Ceed, Ptr{CeedMemType}), ceed, type)
+end
+
+function CeedVectorCreate(ceed, len, vec)
+    ccall((:CeedVectorCreate, libceed), Cint, (Ceed, CeedInt, Ptr{CeedVector}), ceed, len, vec)
+end
+
+function CeedVectorSetArray(vec, mtype, cmode, array)
+    ccall((:CeedVectorSetArray, libceed), Cint, (CeedVector, CeedMemType, CeedCopyMode, Ptr{CeedScalar}), vec, mtype, cmode, array)
+end
+
+function CeedVectorSetValue(vec, value)
+    ccall((:CeedVectorSetValue, libceed), Cint, (CeedVector, CeedScalar), vec, value)
+end
+
+function CeedVectorSyncArray(vec, mtype)
+    ccall((:CeedVectorSyncArray, libceed), Cint, (CeedVector, CeedMemType), vec, mtype)
+end
+
+function CeedVectorTakeArray(vec, mtype, array)
+    ccall((:CeedVectorTakeArray, libceed), Cint, (CeedVector, CeedMemType, Ptr{Ptr{CeedScalar}}), vec, mtype, array)
+end
+
+function CeedVectorGetArray(vec, mtype, array)
+    ccall((:CeedVectorGetArray, libceed), Cint, (CeedVector, CeedMemType, Ptr{Ptr{CeedScalar}}), vec, mtype, array)
+end
+
+function CeedVectorGetArrayRead(vec, mtype, array)
+    ccall((:CeedVectorGetArrayRead, libceed), Cint, (CeedVector, CeedMemType, Ptr{Ptr{CeedScalar}}), vec, mtype, array)
+end
+
+function CeedVectorRestoreArray(vec, array)
+    ccall((:CeedVectorRestoreArray, libceed), Cint, (CeedVector, Ptr{Ptr{CeedScalar}}), vec, array)
+end
+
+function CeedVectorRestoreArrayRead(vec, array)
+    ccall((:CeedVectorRestoreArrayRead, libceed), Cint, (CeedVector, Ptr{Ptr{CeedScalar}}), vec, array)
+end
+
+function CeedVectorNorm(vec, type, norm)
+    ccall((:CeedVectorNorm, libceed), Cint, (CeedVector, CeedNormType, Ptr{CeedScalar}), vec, type, norm)
+end
+
+function CeedVectorReciprocal(vec)
+    ccall((:CeedVectorReciprocal, libceed), Cint, (CeedVector,), vec)
+end
+
+function CeedVectorView(vec, fpfmt, stream)
+    ccall((:CeedVectorView, libceed), Cint, (CeedVector, Cstring, Ptr{FILE}), vec, fpfmt, stream)
+end
+
+function CeedVectorGetLength(vec, length)
+    ccall((:CeedVectorGetLength, libceed), Cint, (CeedVector, Ptr{CeedInt}), vec, length)
+end
+
+function CeedVectorDestroy(vec)
+    ccall((:CeedVectorDestroy, libceed), Cint, (Ptr{CeedVector},), vec)
+end
+
+function CeedRequestWait(req)
+    ccall((:CeedRequestWait, libceed), Cint, (Ptr{CeedRequest},), req)
+end
+
+function CeedElemRestrictionCreate(ceed, nelem, elemsize, ncomp, compstride, lsize, mtype, cmode, offsets, rstr)
+    ccall((:CeedElemRestrictionCreate, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt, CeedMemType, CeedCopyMode, Ptr{CeedInt}, Ptr{CeedElemRestriction}), ceed, nelem, elemsize, ncomp, compstride, lsize, mtype, cmode, offsets, rstr)
+end
+
+function CeedElemRestrictionCreateStrided(ceed, nelem, elemsize, ncomp, lsize, strides, rstr)
+    ccall((:CeedElemRestrictionCreateStrided, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, Ptr{CeedInt}, Ptr{CeedElemRestriction}), ceed, nelem, elemsize, ncomp, lsize, strides, rstr)
+end
+
+function CeedElemRestrictionCreateBlocked(ceed, nelem, elemsize, blksize, ncomp, compstride, lsize, mtype, cmode, offsets, rstr)
+    ccall((:CeedElemRestrictionCreateBlocked, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt, CeedMemType, CeedCopyMode, Ptr{CeedInt}, Ptr{CeedElemRestriction}), ceed, nelem, elemsize, blksize, ncomp, compstride, lsize, mtype, cmode, offsets, rstr)
+end
+
+function CeedElemRestrictionCreateBlockedStrided(ceed, nelem, elemsize, blksize, ncomp, lsize, strides, rstr)
+    ccall((:CeedElemRestrictionCreateBlockedStrided, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt, Ptr{CeedInt}, Ptr{CeedElemRestriction}), ceed, nelem, elemsize, blksize, ncomp, lsize, strides, rstr)
+end
+
+function CeedElemRestrictionCreateVector(rstr, lvec, evec)
+    ccall((:CeedElemRestrictionCreateVector, libceed), Cint, (CeedElemRestriction, Ptr{CeedVector}, Ptr{CeedVector}), rstr, lvec, evec)
+end
+
+function CeedElemRestrictionApply(rstr, tmode, u, ru, request)
+    ccall((:CeedElemRestrictionApply, libceed), Cint, (CeedElemRestriction, CeedTransposeMode, CeedVector, CeedVector, Ptr{CeedRequest}), rstr, tmode, u, ru, request)
+end
+
+function CeedElemRestrictionApplyBlock(rstr, block, tmode, u, ru, request)
+    ccall((:CeedElemRestrictionApplyBlock, libceed), Cint, (CeedElemRestriction, CeedInt, CeedTransposeMode, CeedVector, CeedVector, Ptr{CeedRequest}), rstr, block, tmode, u, ru, request)
+end
+
+function CeedElemRestrictionGetCompStride(rstr, compstride)
+    ccall((:CeedElemRestrictionGetCompStride, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, compstride)
+end
+
+function CeedElemRestrictionGetNumElements(rstr, numelem)
+    ccall((:CeedElemRestrictionGetNumElements, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, numelem)
+end
+
+function CeedElemRestrictionGetElementSize(rstr, elemsize)
+    ccall((:CeedElemRestrictionGetElementSize, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, elemsize)
+end
+
+function CeedElemRestrictionGetLVectorSize(rstr, lsize)
+    ccall((:CeedElemRestrictionGetLVectorSize, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, lsize)
+end
+
+function CeedElemRestrictionGetNumComponents(rstr, numcomp)
+    ccall((:CeedElemRestrictionGetNumComponents, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, numcomp)
+end
+
+function CeedElemRestrictionGetNumBlocks(rstr, numblk)
+    ccall((:CeedElemRestrictionGetNumBlocks, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, numblk)
+end
+
+function CeedElemRestrictionGetBlockSize(rstr, blksize)
+    ccall((:CeedElemRestrictionGetBlockSize, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, blksize)
+end
+
+function CeedElemRestrictionGetMultiplicity(rstr, mult)
+    ccall((:CeedElemRestrictionGetMultiplicity, libceed), Cint, (CeedElemRestriction, CeedVector), rstr, mult)
+end
+
+function CeedElemRestrictionView(rstr, stream)
+    ccall((:CeedElemRestrictionView, libceed), Cint, (CeedElemRestriction, Ptr{FILE}), rstr, stream)
+end
+
+function CeedElemRestrictionDestroy(rstr)
+    ccall((:CeedElemRestrictionDestroy, libceed), Cint, (Ptr{CeedElemRestriction},), rstr)
+end
+
+function CeedBasisCreateTensorH1Lagrange(ceed, dim, ncomp, P, Q, qmode, basis)
+    ccall((:CeedBasisCreateTensorH1Lagrange, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, CeedQuadMode, Ptr{CeedBasis}), ceed, dim, ncomp, P, Q, qmode, basis)
+end
+
+function CeedBasisCreateTensorH1(ceed, dim, ncomp, P1d, Q1d, interp1d, grad1d, qref1d, qweight1d, basis)
+    ccall((:CeedBasisCreateTensorH1, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedBasis}), ceed, dim, ncomp, P1d, Q1d, interp1d, grad1d, qref1d, qweight1d, basis)
+end
+
+function CeedBasisCreateH1(ceed, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight, basis)
+    ccall((:CeedBasisCreateH1, libceed), Cint, (Ceed, CeedElemTopology, CeedInt, CeedInt, CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedBasis}), ceed, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight, basis)
+end
+
+function CeedBasisView(basis, stream)
+    ccall((:CeedBasisView, libceed), Cint, (CeedBasis, Ptr{FILE}), basis, stream)
+end
+
+function CeedBasisApply(basis, nelem, tmode, emode, u, v)
+    ccall((:CeedBasisApply, libceed), Cint, (CeedBasis, CeedInt, CeedTransposeMode, CeedEvalMode, CeedVector, CeedVector), basis, nelem, tmode, emode, u, v)
+end
+
+function CeedBasisGetDimension(basis, dim)
+    ccall((:CeedBasisGetDimension, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, dim)
+end
+
+function CeedBasisGetTopology(basis, topo)
+    ccall((:CeedBasisGetTopology, libceed), Cint, (CeedBasis, Ptr{CeedElemTopology}), basis, topo)
+end
+
+function CeedBasisGetNumComponents(basis, numcomp)
+    ccall((:CeedBasisGetNumComponents, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, numcomp)
+end
+
+function CeedBasisGetNumNodes(basis, P)
+    ccall((:CeedBasisGetNumNodes, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, P)
+end
+
+function CeedBasisGetNumNodes1D(basis, P1d)
+    ccall((:CeedBasisGetNumNodes1D, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, P1d)
+end
+
+function CeedBasisGetNumQuadraturePoints(basis, Q)
+    ccall((:CeedBasisGetNumQuadraturePoints, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, Q)
+end
+
+function CeedBasisGetNumQuadraturePoints1D(basis, Q1d)
+    ccall((:CeedBasisGetNumQuadraturePoints1D, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, Q1d)
+end
+
+function CeedBasisGetQRef(basis, qref)
+    ccall((:CeedBasisGetQRef, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, qref)
+end
+
+function CeedBasisGetQWeights(basis, qweight)
+    ccall((:CeedBasisGetQWeights, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, qweight)
+end
+
+function CeedBasisGetInterp(basis, interp)
+    ccall((:CeedBasisGetInterp, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, interp)
+end
+
+function CeedBasisGetInterp1D(basis, interp1d)
+    ccall((:CeedBasisGetInterp1D, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, interp1d)
+end
+
+function CeedBasisGetGrad(basis, grad)
+    ccall((:CeedBasisGetGrad, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, grad)
+end
+
+function CeedBasisGetGrad1D(basis, grad1d)
+    ccall((:CeedBasisGetGrad1D, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, grad1d)
+end
+
+function CeedBasisDestroy(basis)
+    ccall((:CeedBasisDestroy, libceed), Cint, (Ptr{CeedBasis},), basis)
+end
+
+function CeedGaussQuadrature(Q, qref1d, qweight1d)
+    ccall((:CeedGaussQuadrature, libceed), Cint, (CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}), Q, qref1d, qweight1d)
+end
+
+function CeedLobattoQuadrature(Q, qref1d, qweight1d)
+    ccall((:CeedLobattoQuadrature, libceed), Cint, (CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}), Q, qref1d, qweight1d)
+end
+
+function CeedQRFactorization(ceed, mat, tau, m, n)
+    ccall((:CeedQRFactorization, libceed), Cint, (Ceed, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedInt, CeedInt), ceed, mat, tau, m, n)
+end
+
+function CeedSymmetricSchurDecomposition(ceed, mat, lambda, n)
+    ccall((:CeedSymmetricSchurDecomposition, libceed), Cint, (Ceed, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedInt), ceed, mat, lambda, n)
+end
+
+function CeedSimultaneousDiagonalization(ceed, matA, matB, x, lambda, n)
+    ccall((:CeedSimultaneousDiagonalization, libceed), Cint, (Ceed, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedInt), ceed, matA, matB, x, lambda, n)
+end
+
+function CeedQFunctionCreateInterior(ceed, vlength, f, source, qf)
+    ccall((:CeedQFunctionCreateInterior, libceed), Cint, (Ceed, CeedInt, CeedQFunctionUser, Cstring, Ptr{CeedQFunction}), ceed, vlength, f, source, qf)
+end
+
+function CeedQFunctionCreateInteriorByName(ceed, name, qf)
+    ccall((:CeedQFunctionCreateInteriorByName, libceed), Cint, (Ceed, Cstring, Ptr{CeedQFunction}), ceed, name, qf)
+end
+
+function CeedQFunctionCreateIdentity(ceed, size, inmode, outmode, qf)
+    ccall((:CeedQFunctionCreateIdentity, libceed), Cint, (Ceed, CeedInt, CeedEvalMode, CeedEvalMode, Ptr{CeedQFunction}), ceed, size, inmode, outmode, qf)
+end
+
+function CeedQFunctionAddInput(qf, fieldname, size, emode)
+    ccall((:CeedQFunctionAddInput, libceed), Cint, (CeedQFunction, Cstring, CeedInt, CeedEvalMode), qf, fieldname, size, emode)
+end
+
+function CeedQFunctionAddOutput(qf, fieldname, size, emode)
+    ccall((:CeedQFunctionAddOutput, libceed), Cint, (CeedQFunction, Cstring, CeedInt, CeedEvalMode), qf, fieldname, size, emode)
+end
+
+function CeedQFunctionSetContext(qf, ctx)
+    ccall((:CeedQFunctionSetContext, libceed), Cint, (CeedQFunction, CeedQFunctionContext), qf, ctx)
+end
+
+function CeedQFunctionView(qf, stream)
+    ccall((:CeedQFunctionView, libceed), Cint, (CeedQFunction, Ptr{FILE}), qf, stream)
+end
+
+function CeedQFunctionApply(qf, Q, u, v)
+    ccall((:CeedQFunctionApply, libceed), Cint, (CeedQFunction, CeedInt, Ptr{CeedVector}, Ptr{CeedVector}), qf, Q, u, v)
+end
+
+function CeedQFunctionDestroy(qf)
+    ccall((:CeedQFunctionDestroy, libceed), Cint, (Ptr{CeedQFunction},), qf)
+end
+
+function CeedQFunctionContextCreate(ceed, ctx)
+    ccall((:CeedQFunctionContextCreate, libceed), Cint, (Ceed, Ptr{CeedQFunctionContext}), ceed, ctx)
+end
+
+function CeedQFunctionContextSetData(ctx, mtype, cmode, size, data)
+    ccall((:CeedQFunctionContextSetData, libceed), Cint, (CeedQFunctionContext, CeedMemType, CeedCopyMode, Csize_t, Ptr{Cvoid}), ctx, mtype, cmode, size, data)
+end
+
+function CeedQFunctionContextGetData(ctx, mtype, data)
+    ccall((:CeedQFunctionContextGetData, libceed), Cint, (CeedQFunctionContext, CeedMemType, Ptr{Cvoid}), ctx, mtype, data)
+end
+
+function CeedQFunctionContextRestoreData(ctx, data)
+    ccall((:CeedQFunctionContextRestoreData, libceed), Cint, (CeedQFunctionContext, Ptr{Cvoid}), ctx, data)
+end
+
+function CeedQFunctionContextView(ctx, stream)
+    ccall((:CeedQFunctionContextView, libceed), Cint, (CeedQFunctionContext, Ptr{FILE}), ctx, stream)
+end
+
+function CeedQFunctionContextDestroy(ctx)
+    ccall((:CeedQFunctionContextDestroy, libceed), Cint, (Ptr{CeedQFunctionContext},), ctx)
+end
+
+function CeedOperatorCreate(ceed, qf, dqf, dqfT, op)
+    ccall((:CeedOperatorCreate, libceed), Cint, (Ceed, CeedQFunction, CeedQFunction, CeedQFunction, Ptr{CeedOperator}), ceed, qf, dqf, dqfT, op)
+end
+
+function CeedCompositeOperatorCreate(ceed, op)
+    ccall((:CeedCompositeOperatorCreate, libceed), Cint, (Ceed, Ptr{CeedOperator}), ceed, op)
+end
+
+function CeedOperatorSetField(op, fieldname, r, b, v)
+    ccall((:CeedOperatorSetField, libceed), Cint, (CeedOperator, Cstring, CeedElemRestriction, CeedBasis, CeedVector), op, fieldname, r, b, v)
+end
+
+function CeedCompositeOperatorAddSub(compositeop, subop)
+    ccall((:CeedCompositeOperatorAddSub, libceed), Cint, (CeedOperator, CeedOperator), compositeop, subop)
+end
+
+function CeedOperatorLinearAssembleQFunction(op, assembled, rstr, request)
+    ccall((:CeedOperatorLinearAssembleQFunction, libceed), Cint, (CeedOperator, Ptr{CeedVector}, Ptr{CeedElemRestriction}, Ptr{CeedRequest}), op, assembled, rstr, request)
+end
+
+function CeedOperatorLinearAssembleDiagonal(op, assembled, request)
+    ccall((:CeedOperatorLinearAssembleDiagonal, libceed), Cint, (CeedOperator, CeedVector, Ptr{CeedRequest}), op, assembled, request)
+end
+
+function CeedOperatorLinearAssembleAddDiagonal(op, assembled, request)
+    ccall((:CeedOperatorLinearAssembleAddDiagonal, libceed), Cint, (CeedOperator, CeedVector, Ptr{CeedRequest}), op, assembled, request)
+end
+
+function CeedOperatorLinearAssemblePointBlockDiagonal(op, assembled, request)
+    ccall((:CeedOperatorLinearAssemblePointBlockDiagonal, libceed), Cint, (CeedOperator, CeedVector, Ptr{CeedRequest}), op, assembled, request)
+end
+
+function CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, request)
+    ccall((:CeedOperatorLinearAssembleAddPointBlockDiagonal, libceed), Cint, (CeedOperator, CeedVector, Ptr{CeedRequest}), op, assembled, request)
+end
+
+function CeedOperatorMultigridLevelCreate(opFine, PMultFine, rstrCoarse, basisCoarse, opCoarse, opProlong, opRestrict)
+    ccall((:CeedOperatorMultigridLevelCreate, libceed), Cint, (CeedOperator, CeedVector, CeedElemRestriction, CeedBasis, Ptr{CeedOperator}, Ptr{CeedOperator}, Ptr{CeedOperator}), opFine, PMultFine, rstrCoarse, basisCoarse, opCoarse, opProlong, opRestrict)
+end
+
+function CeedOperatorMultigridLevelCreateTensorH1(opFine, PMultFine, rstrCoarse, basisCoarse, interpCtoF, opCoarse, opProlong, opRestrict)
+    ccall((:CeedOperatorMultigridLevelCreateTensorH1, libceed), Cint, (CeedOperator, CeedVector, CeedElemRestriction, CeedBasis, Ptr{CeedScalar}, Ptr{CeedOperator}, Ptr{CeedOperator}, Ptr{CeedOperator}), opFine, PMultFine, rstrCoarse, basisCoarse, interpCtoF, opCoarse, opProlong, opRestrict)
+end
+
+function CeedOperatorMultigridLevelCreateH1(opFine, PMultFine, rstrCoarse, basisCoarse, interpCtoF, opCoarse, opProlong, opRestrict)
+    ccall((:CeedOperatorMultigridLevelCreateH1, libceed), Cint, (CeedOperator, CeedVector, CeedElemRestriction, CeedBasis, Ptr{CeedScalar}, Ptr{CeedOperator}, Ptr{CeedOperator}, Ptr{CeedOperator}), opFine, PMultFine, rstrCoarse, basisCoarse, interpCtoF, opCoarse, opProlong, opRestrict)
+end
+
+function CeedOperatorCreateFDMElementInverse(op, fdminv, request)
+    ccall((:CeedOperatorCreateFDMElementInverse, libceed), Cint, (CeedOperator, Ptr{CeedOperator}, Ptr{CeedRequest}), op, fdminv, request)
+end
+
+function CeedOperatorView(op, stream)
+    ccall((:CeedOperatorView, libceed), Cint, (CeedOperator, Ptr{FILE}), op, stream)
+end
+
+function CeedOperatorApply(op, in, out, request)
+    ccall((:CeedOperatorApply, libceed), Cint, (CeedOperator, CeedVector, CeedVector, Ptr{CeedRequest}), op, in, out, request)
+end
+
+function CeedOperatorApplyAdd(op, in, out, request)
+    ccall((:CeedOperatorApplyAdd, libceed), Cint, (CeedOperator, CeedVector, CeedVector, Ptr{CeedRequest}), op, in, out, request)
+end
+
+function CeedOperatorDestroy(op)
+    ccall((:CeedOperatorDestroy, libceed), Cint, (Ptr{CeedOperator},), op)
+end
+
+function CeedIntPow(base, power)
+    ccall((:CeedIntPow, libceed), CeedInt, (CeedInt, CeedInt), base, power)
+end
+
+function CeedIntMin(a, b)
+    ccall((:CeedIntMin, libceed), CeedInt, (CeedInt, CeedInt), a, b)
+end
+
+function CeedIntMax(a, b)
+    ccall((:CeedIntMax, libceed), CeedInt, (CeedInt, CeedInt), a, b)
+end
+# Julia wrapper for header: ceed-cuda.h
+# Automatically generated using Clang.jl
+
+
+function CeedQFunctionSetCUDAUserFunction(qf, f)
+    ccall((:CeedQFunctionSetCUDAUserFunction, libceed), Cint, (CeedQFunction, Cint), qf, f)
+end
+# Julia wrapper for header: ceed-backend.h
+# Automatically generated using Clang.jl
+
+
+function CeedMallocArray(n, unit, p)
+    ccall((:CeedMallocArray, libceed), Cint, (Csize_t, Csize_t, Ptr{Cvoid}), n, unit, p)
+end
+
+function CeedCallocArray(n, unit, p)
+    ccall((:CeedCallocArray, libceed), Cint, (Csize_t, Csize_t, Ptr{Cvoid}), n, unit, p)
+end
+
+function CeedReallocArray(n, unit, p)
+    ccall((:CeedReallocArray, libceed), Cint, (Csize_t, Csize_t, Ptr{Cvoid}), n, unit, p)
+end
+
+function CeedFree(p)
+    ccall((:CeedFree, libceed), Cint, (Ptr{Cvoid},), p)
+end
+
+function CeedRegister(prefix, init, priority)
+    ccall((:CeedRegister, libceed), Cint, (Cstring, Ptr{Cvoid}, UInt32), prefix, init, priority)
+end
+
+function CeedIsDebug(ceed, isDebug)
+    ccall((:CeedIsDebug, libceed), Cint, (Ceed, Ptr{Bool}), ceed, isDebug)
+end
+
+function CeedGetParent(ceed, parent)
+    ccall((:CeedGetParent, libceed), Cint, (Ceed, Ptr{Ceed}), ceed, parent)
+end
+
+function CeedGetDelegate(ceed, delegate)
+    ccall((:CeedGetDelegate, libceed), Cint, (Ceed, Ptr{Ceed}), ceed, delegate)
+end
+
+function CeedSetDelegate(ceed, delegate)
+    ccall((:CeedSetDelegate, libceed), Cint, (Ceed, Ceed), ceed, delegate)
+end
+
+function CeedGetObjectDelegate(ceed, delegate, objname)
+    ccall((:CeedGetObjectDelegate, libceed), Cint, (Ceed, Ptr{Ceed}, Cstring), ceed, delegate, objname)
+end
+
+function CeedSetObjectDelegate(ceed, delegate, objname)
+    ccall((:CeedSetObjectDelegate, libceed), Cint, (Ceed, Ceed, Cstring), ceed, delegate, objname)
+end
+
+function CeedGetOperatorFallbackResource(ceed, resource)
+    ccall((:CeedGetOperatorFallbackResource, libceed), Cint, (Ceed, Ptr{Cstring}), ceed, resource)
+end
+
+function CeedSetOperatorFallbackResource(ceed, resource)
+    ccall((:CeedSetOperatorFallbackResource, libceed), Cint, (Ceed, Cstring), ceed, resource)
+end
+
+function CeedGetOperatorFallbackParentCeed(ceed, parent)
+    ccall((:CeedGetOperatorFallbackParentCeed, libceed), Cint, (Ceed, Ptr{Ceed}), ceed, parent)
+end
+
+function CeedSetDeterministic(ceed, isDeterministic)
+    ccall((:CeedSetDeterministic, libceed), Cint, (Ceed, Bool), ceed, isDeterministic)
+end
+
+function CeedSetBackendFunction(ceed, type, object, fname, f)
+    ccall((:CeedSetBackendFunction, libceed), Cint, (Ceed, Cstring, Ptr{Cvoid}, Cstring, Ptr{Cvoid}), ceed, type, object, fname, f)
+end
+
+function CeedGetData(ceed, data)
+    ccall((:CeedGetData, libceed), Cint, (Ceed, Ptr{Cvoid}), ceed, data)
+end
+
+function CeedSetData(ceed, data)
+    ccall((:CeedSetData, libceed), Cint, (Ceed, Ptr{Cvoid}), ceed, data)
+end
+
+function CeedVectorGetCeed(vec, ceed)
+    ccall((:CeedVectorGetCeed, libceed), Cint, (CeedVector, Ptr{Ceed}), vec, ceed)
+end
+
+function CeedVectorGetState(vec, state)
+    ccall((:CeedVectorGetState, libceed), Cint, (CeedVector, Ptr{UInt64}), vec, state)
+end
+
+function CeedVectorAddReference(vec)
+    ccall((:CeedVectorAddReference, libceed), Cint, (CeedVector,), vec)
+end
+
+function CeedVectorGetData(vec, data)
+    ccall((:CeedVectorGetData, libceed), Cint, (CeedVector, Ptr{Cvoid}), vec, data)
+end
+
+function CeedVectorSetData(vec, data)
+    ccall((:CeedVectorSetData, libceed), Cint, (CeedVector, Ptr{Cvoid}), vec, data)
+end
+
+function CeedElemRestrictionGetCeed(rstr, ceed)
+    ccall((:CeedElemRestrictionGetCeed, libceed), Cint, (CeedElemRestriction, Ptr{Ceed}), rstr, ceed)
+end
+
+function CeedElemRestrictionGetStrides(rstr, strides)
+    ccall((:CeedElemRestrictionGetStrides, libceed), Cint, (CeedElemRestriction, Ptr{NTuple{3, CeedInt}}), rstr, strides)
+end
+
+function CeedElemRestrictionGetOffsets(rstr, mtype, offsets)
+    ccall((:CeedElemRestrictionGetOffsets, libceed), Cint, (CeedElemRestriction, CeedMemType, Ptr{Ptr{CeedInt}}), rstr, mtype, offsets)
+end
+
+function CeedElemRestrictionRestoreOffsets(rstr, offsets)
+    ccall((:CeedElemRestrictionRestoreOffsets, libceed), Cint, (CeedElemRestriction, Ptr{Ptr{CeedInt}}), rstr, offsets)
+end
+
+function CeedElemRestrictionIsStrided(rstr, isstrided)
+    ccall((:CeedElemRestrictionIsStrided, libceed), Cint, (CeedElemRestriction, Ptr{Bool}), rstr, isstrided)
+end
+
+function CeedElemRestrictionHasBackendStrides(rstr, hasbackendstrides)
+    ccall((:CeedElemRestrictionHasBackendStrides, libceed), Cint, (CeedElemRestriction, Ptr{Bool}), rstr, hasbackendstrides)
+end
+
+function CeedElemRestrictionGetELayout(rstr, layout)
+    ccall((:CeedElemRestrictionGetELayout, libceed), Cint, (CeedElemRestriction, Ptr{NTuple{3, CeedInt}}), rstr, layout)
+end
+
+function CeedElemRestrictionSetELayout(rstr, layout)
+    ccall((:CeedElemRestrictionSetELayout, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, layout)
+end
+
+function CeedElemRestrictionGetData(rstr, data)
+    ccall((:CeedElemRestrictionGetData, libceed), Cint, (CeedElemRestriction, Ptr{Cvoid}), rstr, data)
+end
+
+function CeedElemRestrictionSetData(rstr, data)
+    ccall((:CeedElemRestrictionSetData, libceed), Cint, (CeedElemRestriction, Ptr{Cvoid}), rstr, data)
+end
+
+function CeedBasisGetCollocatedGrad(basis, colograd1d)
+    ccall((:CeedBasisGetCollocatedGrad, libceed), Cint, (CeedBasis, Ptr{CeedScalar}), basis, colograd1d)
+end
+
+function CeedHouseholderApplyQ(A, Q, tau, tmode, m, n, k, row, col)
+    ccall((:CeedHouseholderApplyQ, libceed), Cint, (Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedTransposeMode, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt), A, Q, tau, tmode, m, n, k, row, col)
+end
+
+function CeedBasisGetCeed(basis, ceed)
+    ccall((:CeedBasisGetCeed, libceed), Cint, (CeedBasis, Ptr{Ceed}), basis, ceed)
+end
+
+function CeedBasisIsTensor(basis, istensor)
+    ccall((:CeedBasisIsTensor, libceed), Cint, (CeedBasis, Ptr{Bool}), basis, istensor)
+end
+
+function CeedBasisGetData(basis, data)
+    ccall((:CeedBasisGetData, libceed), Cint, (CeedBasis, Ptr{Cvoid}), basis, data)
+end
+
+function CeedBasisSetData(basis, data)
+    ccall((:CeedBasisSetData, libceed), Cint, (CeedBasis, Ptr{Cvoid}), basis, data)
+end
+
+function CeedBasisGetTopologyDimension(topo, dim)
+    ccall((:CeedBasisGetTopologyDimension, libceed), Cint, (CeedElemTopology, Ptr{CeedInt}), topo, dim)
+end
+
+function CeedBasisGetTensorContract(basis, contract)
+    ccall((:CeedBasisGetTensorContract, libceed), Cint, (CeedBasis, Ptr{CeedTensorContract}), basis, contract)
+end
+
+function CeedBasisSetTensorContract(basis, contract)
+    ccall((:CeedBasisSetTensorContract, libceed), Cint, (CeedBasis, Ptr{CeedTensorContract}), basis, contract)
+end
+
+function CeedTensorContractCreate(ceed, basis, contract)
+    ccall((:CeedTensorContractCreate, libceed), Cint, (Ceed, CeedBasis, Ptr{CeedTensorContract}), ceed, basis, contract)
+end
+
+function CeedTensorContractApply(contract, A, B, C, J, t, tmode, Add, u, v)
+    ccall((:CeedTensorContractApply, libceed), Cint, (CeedTensorContract, CeedInt, CeedInt, CeedInt, CeedInt, Ptr{CeedScalar}, CeedTransposeMode, CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}), contract, A, B, C, J, t, tmode, Add, u, v)
+end
+
+function CeedTensorContractGetCeed(contract, ceed)
+    ccall((:CeedTensorContractGetCeed, libceed), Cint, (CeedTensorContract, Ptr{Ceed}), contract, ceed)
+end
+
+function CeedTensorContractGetData(contract, data)
+    ccall((:CeedTensorContractGetData, libceed), Cint, (CeedTensorContract, Ptr{Cvoid}), contract, data)
+end
+
+function CeedTensorContractSetData(contract, data)
+    ccall((:CeedTensorContractSetData, libceed), Cint, (CeedTensorContract, Ptr{Cvoid}), contract, data)
+end
+
+function CeedTensorContractDestroy(contract)
+    ccall((:CeedTensorContractDestroy, libceed), Cint, (Ptr{CeedTensorContract},), contract)
+end
+
+function CeedQFunctionRegister(arg1, arg2, arg3, arg4, init)
+    ccall((:CeedQFunctionRegister, libceed), Cint, (Cstring, Cstring, CeedInt, CeedQFunctionUser, Ptr{Cvoid}), arg1, arg2, arg3, arg4, init)
+end
+
+function CeedQFunctionSetFortranStatus(qf, status)
+    ccall((:CeedQFunctionSetFortranStatus, libceed), Cint, (CeedQFunction, Bool), qf, status)
+end
+
+function CeedQFunctionGetCeed(qf, ceed)
+    ccall((:CeedQFunctionGetCeed, libceed), Cint, (CeedQFunction, Ptr{Ceed}), qf, ceed)
+end
+
+function CeedQFunctionGetVectorLength(qf, vlength)
+    ccall((:CeedQFunctionGetVectorLength, libceed), Cint, (CeedQFunction, Ptr{CeedInt}), qf, vlength)
+end
+
+function CeedQFunctionGetNumArgs(qf, numinputfields, numoutputfields)
+    ccall((:CeedQFunctionGetNumArgs, libceed), Cint, (CeedQFunction, Ptr{CeedInt}, Ptr{CeedInt}), qf, numinputfields, numoutputfields)
+end
+
+function CeedQFunctionGetSourcePath(qf, source)
+    ccall((:CeedQFunctionGetSourcePath, libceed), Cint, (CeedQFunction, Ptr{Cstring}), qf, source)
+end
+
+function CeedQFunctionGetUserFunction(qf, f)
+    ccall((:CeedQFunctionGetUserFunction, libceed), Cint, (CeedQFunction, Ptr{CeedQFunctionUser}), qf, f)
+end
+
+function CeedQFunctionGetContext(qf, ctx)
+    ccall((:CeedQFunctionGetContext, libceed), Cint, (CeedQFunction, Ptr{CeedQFunctionContext}), qf, ctx)
+end
+
+function CeedQFunctionGetInnerContext(qf, ctx)
+    ccall((:CeedQFunctionGetInnerContext, libceed), Cint, (CeedQFunction, Ptr{CeedQFunctionContext}), qf, ctx)
+end
+
+function CeedQFunctionIsIdentity(qf, isidentity)
+    ccall((:CeedQFunctionIsIdentity, libceed), Cint, (CeedQFunction, Ptr{Bool}), qf, isidentity)
+end
+
+function CeedQFunctionGetData(qf, data)
+    ccall((:CeedQFunctionGetData, libceed), Cint, (CeedQFunction, Ptr{Cvoid}), qf, data)
+end
+
+function CeedQFunctionSetData(qf, data)
+    ccall((:CeedQFunctionSetData, libceed), Cint, (CeedQFunction, Ptr{Cvoid}), qf, data)
+end
+
+function CeedQFunctionGetFields(qf, inputfields, outputfields)
+    ccall((:CeedQFunctionGetFields, libceed), Cint, (CeedQFunction, Ptr{Ptr{CeedQFunctionField}}, Ptr{Ptr{CeedQFunctionField}}), qf, inputfields, outputfields)
+end
+
+function CeedQFunctionFieldGetName(qffield, fieldname)
+    ccall((:CeedQFunctionFieldGetName, libceed), Cint, (CeedQFunctionField, Ptr{Cstring}), qffield, fieldname)
+end
+
+function CeedQFunctionFieldGetSize(qffield, size)
+    ccall((:CeedQFunctionFieldGetSize, libceed), Cint, (CeedQFunctionField, Ptr{CeedInt}), qffield, size)
+end
+
+function CeedQFunctionFieldGetEvalMode(qffield, emode)
+    ccall((:CeedQFunctionFieldGetEvalMode, libceed), Cint, (CeedQFunctionField, Ptr{CeedEvalMode}), qffield, emode)
+end
+
+function CeedQFunctionContextGetCeed(cxt, ceed)
+    ccall((:CeedQFunctionContextGetCeed, libceed), Cint, (CeedQFunctionContext, Ptr{Ceed}), cxt, ceed)
+end
+
+function CeedQFunctionContextGetState(ctx, state)
+    ccall((:CeedQFunctionContextGetState, libceed), Cint, (CeedQFunctionContext, Ptr{UInt64}), ctx, state)
+end
+
+function CeedQFunctionContextGetContextSize(ctx, ctxsize)
+    ccall((:CeedQFunctionContextGetContextSize, libceed), Cint, (CeedQFunctionContext, Ptr{Csize_t}), ctx, ctxsize)
+end
+
+function CeedQFunctionContextGetBackendData(ctx, data)
+    ccall((:CeedQFunctionContextGetBackendData, libceed), Cint, (CeedQFunctionContext, Ptr{Cvoid}), ctx, data)
+end
+
+function CeedQFunctionContextSetBackendData(ctx, data)
+    ccall((:CeedQFunctionContextSetBackendData, libceed), Cint, (CeedQFunctionContext, Ptr{Cvoid}), ctx, data)
+end
+
+function CeedOperatorGetCeed(op, ceed)
+    ccall((:CeedOperatorGetCeed, libceed), Cint, (CeedOperator, Ptr{Ceed}), op, ceed)
+end
+
+function CeedOperatorGetNumElements(op, numelem)
+    ccall((:CeedOperatorGetNumElements, libceed), Cint, (CeedOperator, Ptr{CeedInt}), op, numelem)
+end
+
+function CeedOperatorGetNumQuadraturePoints(op, numqpts)
+    ccall((:CeedOperatorGetNumQuadraturePoints, libceed), Cint, (CeedOperator, Ptr{CeedInt}), op, numqpts)
+end
+
+function CeedOperatorGetNumArgs(op, numargs)
+    ccall((:CeedOperatorGetNumArgs, libceed), Cint, (CeedOperator, Ptr{CeedInt}), op, numargs)
+end
+
+function CeedOperatorIsSetupDone(op, issetupdone)
+    ccall((:CeedOperatorIsSetupDone, libceed), Cint, (CeedOperator, Ptr{Bool}), op, issetupdone)
+end
+
+function CeedOperatorGetQFunction(op, qf)
+    ccall((:CeedOperatorGetQFunction, libceed), Cint, (CeedOperator, Ptr{CeedQFunction}), op, qf)
+end
+
+function CeedOperatorIsComposite(op, iscomposite)
+    ccall((:CeedOperatorIsComposite, libceed), Cint, (CeedOperator, Ptr{Bool}), op, iscomposite)
+end
+
+function CeedOperatorGetNumSub(op, numsub)
+    ccall((:CeedOperatorGetNumSub, libceed), Cint, (CeedOperator, Ptr{CeedInt}), op, numsub)
+end
+
+function CeedOperatorGetSubList(op, suboperators)
+    ccall((:CeedOperatorGetSubList, libceed), Cint, (CeedOperator, Ptr{Ptr{CeedOperator}}), op, suboperators)
+end
+
+function CeedOperatorGetData(op, data)
+    ccall((:CeedOperatorGetData, libceed), Cint, (CeedOperator, Ptr{Cvoid}), op, data)
+end
+
+function CeedOperatorSetData(op, data)
+    ccall((:CeedOperatorSetData, libceed), Cint, (CeedOperator, Ptr{Cvoid}), op, data)
+end
+
+function CeedOperatorSetSetupDone(op)
+    ccall((:CeedOperatorSetSetupDone, libceed), Cint, (CeedOperator,), op)
+end
+
+function CeedOperatorGetFields(op, inputfields, outputfields)
+    ccall((:CeedOperatorGetFields, libceed), Cint, (CeedOperator, Ptr{Ptr{CeedOperatorField}}, Ptr{Ptr{CeedOperatorField}}), op, inputfields, outputfields)
+end
+
+function CeedOperatorFieldGetElemRestriction(opfield, rstr)
+    ccall((:CeedOperatorFieldGetElemRestriction, libceed), Cint, (CeedOperatorField, Ptr{CeedElemRestriction}), opfield, rstr)
+end
+
+function CeedOperatorFieldGetBasis(opfield, basis)
+    ccall((:CeedOperatorFieldGetBasis, libceed), Cint, (CeedOperatorField, Ptr{CeedBasis}), opfield, basis)
+end
+
+function CeedOperatorFieldGetVector(opfield, vec)
+    ccall((:CeedOperatorFieldGetVector, libceed), Cint, (CeedOperatorField, Ptr{CeedVector}), opfield, vec)
+end
+
+function CeedMatrixMultiply(ceed, matA, matB, matC, m, n, kk)
+    ccall((:CeedMatrixMultiply, libceed), Cint, (Ceed, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedInt, CeedInt, CeedInt), ceed, matA, matB, matC, m, n, kk)
+end
diff --git a/julia/LibCEED.jl/src/generated/libceed_common.jl b/julia/LibCEED.jl/src/generated/libceed_common.jl
new file mode 100644
index 0000000000..aad84aa497
--- /dev/null
+++ b/julia/LibCEED.jl/src/generated/libceed_common.jl
@@ -0,0 +1,99 @@
+# Automatically generated using Clang.jl
+#! format: off
+
+const FILE = Cvoid
+
+# Skipping MacroDefinition: CEED_QFUNCTION ( name ) static const char name ## _loc [ ] = __FILE__ ":" # name ; static int name
+
+# Skipping MacroDefinition: CeedError ( ceed , ecode , ... ) ( CeedErrorImpl ( ( ceed ) , __FILE__ , __LINE__ , __func__ , ( ecode ) , __VA_ARGS__ ) ? : ( ecode ) )
+
+const CeedInt = Int32
+const CeedScalar = Cdouble
+const Ceed_private = Cvoid
+const Ceed = Ptr{Ceed_private}
+const CeedRequest_private = Cvoid
+const CeedRequest = Ptr{CeedRequest_private}
+const CeedVector_private = Cvoid
+const CeedVector = Ptr{CeedVector_private}
+const CeedElemRestriction_private = Cvoid
+const CeedElemRestriction = Ptr{CeedElemRestriction_private}
+const CeedBasis_private = Cvoid
+const CeedBasis = Ptr{CeedBasis_private}
+const CeedQFunction_private = Cvoid
+const CeedQFunction = Ptr{CeedQFunction_private}
+const CeedQFunctionContext_private = Cvoid
+const CeedQFunctionContext = Ptr{CeedQFunctionContext_private}
+const CeedOperator_private = Cvoid
+const CeedOperator = Ptr{CeedOperator_private}
+
+@cenum CeedMemType::UInt32 begin
+    CEED_MEM_HOST = 0
+    CEED_MEM_DEVICE = 1
+end
+
+@cenum CeedCopyMode::UInt32 begin
+    CEED_COPY_VALUES = 0
+    CEED_USE_POINTER = 1
+    CEED_OWN_POINTER = 2
+end
+
+@cenum CeedNormType::UInt32 begin
+    CEED_NORM_1 = 0
+    CEED_NORM_2 = 1
+    CEED_NORM_MAX = 2
+end
+
+@cenum CeedTransposeMode::UInt32 begin
+    CEED_NOTRANSPOSE = 0
+    CEED_TRANSPOSE = 1
+end
+
+@cenum CeedEvalMode::UInt32 begin
+    CEED_EVAL_NONE = 0
+    CEED_EVAL_INTERP = 1
+    CEED_EVAL_GRAD = 2
+    CEED_EVAL_DIV = 4
+    CEED_EVAL_CURL = 8
+    CEED_EVAL_WEIGHT = 16
+end
+
+@cenum CeedQuadMode::UInt32 begin
+    CEED_GAUSS = 0
+    CEED_GAUSS_LOBATTO = 1
+end
+
+@cenum CeedElemTopology::UInt32 begin
+    CEED_LINE = 65536
+    CEED_TRIANGLE = 131073
+    CEED_QUAD = 131074
+    CEED_TET = 196611
+    CEED_PYRAMID = 196612
+    CEED_PRISM = 196613
+    CEED_HEX = 196614
+end
+
+
+const CeedQFunctionUser = Ptr{Cvoid}
+
+# Skipping MacroDefinition: CEED_INTERN CEED_EXTERN __attribute__ ( ( visibility ( "hidden" ) ) )
+
+const CEED_MAX_RESOURCE_LEN = 1024
+const CEED_ALIGN = 64
+const CEED_COMPOSITE_MAX = 16
+const CEED_EPSILON = 1.0e-16
+const CEED_DEBUG_COLOR = 0
+
+# Skipping MacroDefinition: CeedDebug1 ( ceed , format , ... ) CeedDebugImpl ( ceed , format , ## __VA_ARGS__ )
+# Skipping MacroDefinition: CeedDebug256 ( ceed , color , ... ) CeedDebugImpl256 ( ceed , color , ## __VA_ARGS__ )
+# Skipping MacroDefinition: CeedDebug ( ... ) CeedDebug256 ( ceed , ( unsigned char ) CEED_DEBUG_COLOR , ## __VA_ARGS__ )
+# Skipping MacroDefinition: CeedChk ( ierr ) do { if ( ierr ) return ierr ; } while ( 0 )
+# Skipping MacroDefinition: CeedMalloc ( n , p ) CeedMallocArray ( ( n ) , sizeof ( * * ( p ) ) , p )
+# Skipping MacroDefinition: CeedCalloc ( n , p ) CeedCallocArray ( ( n ) , sizeof ( * * ( p ) ) , p )
+# Skipping MacroDefinition: CeedRealloc ( n , p ) CeedReallocArray ( ( n ) , sizeof ( * * ( p ) ) , p )
+
+const CeedTensorContract_private = Cvoid
+const CeedTensorContract = Ptr{CeedTensorContract_private}
+const CeedQFunctionField_private = Cvoid
+const CeedQFunctionField = Ptr{CeedQFunctionField_private}
+const CeedOperatorField_private = Cvoid
+const CeedOperatorField = Ptr{CeedOperatorField_private}