diff --git a/Project.toml b/Project.toml
index 19aa4a1..29b9ba7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "FiniteElementContainers"
 uuid = "d08262e4-672f-4e7f-a976-f2cea5767631"
 authors = ["Craig M. Hamel <cmhamel32@gmail.com> and contributors"]
-version = "0.10.1"
+version = "0.10.2"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -27,16 +27,16 @@ MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
 
 [extensions]
-AbaqusReaderExt = "AbaqusReader"
 AMDGPUExt = "AMDGPU"
+AbaqusReaderExt = "AbaqusReader"
 BlockArraysExt = "BlockArrays"
 CUDAExt = "CUDA"
 MPIExt = "MPI"
 PartitionedArraysExt = "PartitionedArrays"
 
 [compat]
-AbaqusReader = "0.2.7"
 AMDGPU = "2"
+AbaqusReader = "0.2.7"
 Adapt = "4"
 Aqua = "0.8"
 Atomix = "1"
@@ -60,7 +60,6 @@ julia = "1"
 
 [extras]
 AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e"
-Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
diff --git a/examples/moose-tutorials-clone/ex01/mug.e b/examples/moose-tutorials-clone/ex01/mug.e
new file mode 100644
index 0000000..1b461db
Binary files /dev/null and b/examples/moose-tutorials-clone/ex01/mug.e differ
diff --git a/examples/moose-tutorials-clone/ex01/script.jl b/examples/moose-tutorials-clone/ex01/script.jl
new file mode 100644
index 0000000..95cc486
--- /dev/null
+++ b/examples/moose-tutorials-clone/ex01/script.jl
@@ -0,0 +1,76 @@
+using FiniteElementContainers
+using LinearAlgebra
+
+# physics code we actually need to implement
+struct Poisson{F <: Function} <: AbstractPhysics{1, 0, 0}
+    func::F
+end
+
+@inline function FiniteElementContainers.residual(
+    physics::Poisson, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    ∇u_q = interpolate_field_gradients(physics, interps, u_el)
+    ∇u_q = unpack_field(∇u_q, 1)
+    R_q = ∇N_X * ∇u_q - N * physics.func(X_q, 0.0)
+    return JxW * R_q[:]
+end
+  
+@inline function FiniteElementContainers.stiffness(
+    physics::Poisson, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    K_q = ∇N_X * ∇N_X'
+    return JxW * K_q
+end
+
+# define some helper funcs upfront
+f(_, _) = 0.0
+one_func(_, _) = 1.0
+zero_func(_, _) = 0.0
+
+function simulate()
+    # load a mesh
+    mesh = UnstructuredMesh("mug.e")
+    V = FunctionSpace(mesh, H1Field, Lagrange)
+    physics = Poisson(f)
+    props = create_properties(physics)
+    u = ScalarFunction(V, :u)
+    asm = SparseMatrixAssembler(u; use_condensed=true)
+    dbcs = [
+        DirichletBC(:u, :bottom, one_func)
+        DirichletBC(:u, :top, zero_func)
+    ]
+    U = create_field(asm)
+    p = create_parameters(mesh, asm, physics, props; dirichlet_bcs = dbcs)
+
+    # solver = NewtonSolver(DirectLinearSolver(asm))
+    # integrator = QuasiStaticIntegrator(solver)
+    # evolve!(integrator, p)
+
+    FiniteElementContainers.update_bc_values!(p)
+    residual_int = VectorIntegral(asm, residual)
+    stiffness_int = MatrixIntegral(asm, stiffness)
+
+    K = stiffness_int(U, p)
+    remove_fixed_dofs!(stiffness_int)
+
+    for n in 1:2
+        R = residual_int(U, p)
+        remove_fixed_dofs!(residual_int)
+        ΔU = -K \ R.data
+        U.data .+= ΔU
+        @show n norm(R) norm(ΔU)
+    end
+
+    U = p.h1_field
+
+    pp = PostProcessor(mesh, "output.e", u)
+    write_times(pp, 1, 0.0)
+    write_field(pp, 1, ("u",), U)
+    close(pp)
+end
+
+simulate()
diff --git a/examples/moose-tutorials-clone/ex02/mug.e b/examples/moose-tutorials-clone/ex02/mug.e
new file mode 100644
index 0000000..1b461db
Binary files /dev/null and b/examples/moose-tutorials-clone/ex02/mug.e differ
diff --git a/examples/moose-tutorials-clone/ex02/script.jl b/examples/moose-tutorials-clone/ex02/script.jl
new file mode 100644
index 0000000..16f63fb
--- /dev/null
+++ b/examples/moose-tutorials-clone/ex02/script.jl
@@ -0,0 +1,78 @@
+using FiniteElementContainers
+using LinearAlgebra
+using StaticArrays
+
+# physics code we actually need to implement
+struct Advection{N} <: AbstractPhysics{1, 0, 0}
+    v::SVector{N, Float64}
+end
+
+@inline function FiniteElementContainers.residual(
+    physics::Advection, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    ∇u_q = interpolate_field_gradients(physics, interps, u_el)
+    ∇u_q = unpack_field(∇u_q, 1)
+    term = dot(∇u_q, physics.v)
+    R_q = ∇N_X * ∇u_q + term * N
+    return JxW * R_q[:]
+end
+  
+@inline function FiniteElementContainers.stiffness(
+    physics::Advection, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    term = ∇N_X * physics.v
+    K_q = ∇N_X * ∇N_X' + N * term'
+    return JxW * K_q
+end
+
+# define some helper funcs upfront
+one_func(_, _) = 1.0
+zero_func(_, _) = 0.0
+
+function simulate()
+    # load a mesh
+    mesh = UnstructuredMesh("mug.e")
+    V = FunctionSpace(mesh, H1Field, Lagrange)
+    physics = Advection(SVector{3, Float64}(0., 0., 1.))
+    props = create_properties(physics)
+    u = ScalarFunction(V, :u)
+    asm = SparseMatrixAssembler(u; use_condensed=true)
+    dbcs = [
+        DirichletBC(:u, :bottom, one_func)
+        DirichletBC(:u, :top, zero_func)
+    ]
+    U = create_field(asm)
+    p = create_parameters(mesh, asm, physics, props; dirichlet_bcs = dbcs)
+
+    # solver = NewtonSolver(DirectLinearSolver(asm))
+    # integrator = QuasiStaticIntegrator(solver)
+    # evolve!(integrator, p)
+
+    FiniteElementContainers.update_bc_values!(p)
+    residual_int = VectorIntegral(asm, residual)
+    stiffness_int = MatrixIntegral(asm, stiffness)
+
+    K = stiffness_int(U, p)
+    remove_fixed_dofs!(stiffness_int)
+
+    for n in 1:2
+        R = residual_int(U, p)
+        remove_fixed_dofs!(residual_int)
+        ΔU = -K \ R.data
+        U.data .+= ΔU
+        @show n norm(R) norm(ΔU)
+    end
+
+    U = p.h1_field
+
+    pp = PostProcessor(mesh, "output.e", u)
+    write_times(pp, 1, 0.0)
+    write_field(pp, 1, ("u",), U)
+    close(pp)
+end
+
+simulate()
diff --git a/examples/moose-tutorials-clone/ex03/mug.e b/examples/moose-tutorials-clone/ex03/mug.e
new file mode 100644
index 0000000..1b461db
Binary files /dev/null and b/examples/moose-tutorials-clone/ex03/mug.e differ
diff --git a/examples/moose-tutorials-clone/ex03/script.jl b/examples/moose-tutorials-clone/ex03/script.jl
new file mode 100644
index 0000000..16d584f
--- /dev/null
+++ b/examples/moose-tutorials-clone/ex03/script.jl
@@ -0,0 +1,107 @@
+using FiniteElementContainers
+using ForwardDiff
+using LinearAlgebra
+using StaticArrays
+
+struct CoupledPhysics <: AbstractPhysics{2, 0, 0}
+end
+
+@inline function FiniteElementContainers.residual(
+    physics::CoupledPhysics, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    ∇u_q = interpolate_field_gradients(physics, interps, u_el)
+    ∇u = unpack_field(∇u_q, 1)
+    ∇v = unpack_field(∇u_q, 2)
+    term = dot(∇u, ∇v)
+    R_u = ∇N_X * ∇u + term * N
+    R_v = ∇N_X * ∇v
+    R = vcat(R_u', R_v')
+    return JxW * R[:]
+end
+  
+@inline function FiniteElementContainers.stiffness(
+    physics::CoupledPhysics, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    # interps = map_interpolants(interps, x_el)
+    # (; X_q, N, ∇N_X, JxW) = interps
+    # ∇u_q = interpolate_field_gradients(physics, interps, u_el)
+    # ∇u = unpack_field(∇u_q, 1)
+    # ∇v = unpack_field(∇u_q, 2)
+
+    # # term = ∇N_X * physics.v
+    # term_1 = ∇N_X * ∇v
+    # term_2 = ∇N_X * ∇v
+    # K_uu = ∇N_X * ∇N_X' + N * term_1'
+    # K_uv = N * term_2'
+    # K_vu = zeros(typeof(K_uv))
+    # K_vv = ∇N_X * ∇N_X'
+
+    # # K_q = ∇N_X * ∇N_X' + N * term'
+    # K = vcat(
+    #     hcat(K_uu, K_uv),
+    #     hcat(K_vu, K_vv)
+    # )
+    # return JxW * K
+    ForwardDiff.jacobian(
+        z -> FiniteElementContainers.residual(
+            physics, interps, x_el, t, dt, z, u_el_old, state_old_q, state_new_q, props_el
+        ), u_el
+    )
+end
+
+# define some helper funcs upfront
+one_func(_, _) = 1.0
+two_func(_, _) = 2.0
+zero_func(_, _) = 0.0
+
+function simulate()
+    # load a mesh
+    mesh = UnstructuredMesh("mug.e")
+    V = FunctionSpace(mesh, H1Field, Lagrange)
+    physics = CoupledPhysics()
+    props = create_properties(physics)
+    u = FiniteElementContainers.GeneralFunction(
+        ScalarFunction(V, :u),
+        ScalarFunction(V, :v)
+    )
+
+    asm = SparseMatrixAssembler(u; use_condensed=true)
+    dbcs = [
+        DirichletBC(:u, :bottom, two_func)
+        DirichletBC(:u, :top, zero_func)
+        DirichletBC(:v, :bottom, one_func)
+        DirichletBC(:v, :top, zero_func)
+    ]
+    U = create_field(asm)
+    p = create_parameters(mesh, asm, physics, props; dirichlet_bcs = dbcs)
+
+    solver = NewtonSolver(DirectLinearSolver(asm))
+    integrator = QuasiStaticIntegrator(solver)
+    evolve!(integrator, p)
+
+    # FiniteElementContainers.update_bc_values!(p)
+    # residual_int = VectorIntegral(asm, residual)
+    # stiffness_int = MatrixIntegral(asm, stiffness)
+
+    # K = stiffness_int(U, p)
+    # remove_fixed_dofs!(stiffness_int)
+
+    # for n in 1:4
+    #     R = residual_int(U, p)
+    #     remove_fixed_dofs!(residual_int)
+    #     ΔU = -K \ R.data
+    #     U.data .+= ΔU
+    #     @show n norm(R) norm(ΔU)
+    # end
+
+    U = p.h1_field
+
+    pp = PostProcessor(mesh, "output.e", u)
+    write_times(pp, 1, 0.0)
+    write_field(pp, 1, ("u", "v"), U)
+    close(pp)
+end
+
+simulate()
\ No newline at end of file
diff --git a/examples/moose-tutorials-clone/ex06/cyl-tet.e b/examples/moose-tutorials-clone/ex06/cyl-tet.e
new file mode 100644
index 0000000..23e8639
Binary files /dev/null and b/examples/moose-tutorials-clone/ex06/cyl-tet.e differ
diff --git a/examples/moose-tutorials-clone/ex06/script.jl b/examples/moose-tutorials-clone/ex06/script.jl
new file mode 100644
index 0000000..9b57c5d
--- /dev/null
+++ b/examples/moose-tutorials-clone/ex06/script.jl
@@ -0,0 +1,69 @@
+using FiniteElementContainers
+using LinearAlgebra
+
+# physics code we actually need to implement
+struct Transient <: AbstractPhysics{1, 0, 0}
+end
+
+@inline function FiniteElementContainers.residual(
+    physics::Transient, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    u_q, ∇u_q = interpolate_field_values_and_gradients(physics, interps, u_el)
+    u_q_old = interpolate_field_values(physics, interps, u_el_old)
+    ∇u = unpack_field(∇u_q, 1)
+    dudt = 20. * (u_q[1] - u_q_old[1]) / dt
+    R = dudt * N + ∇N_X * ∇u
+    return JxW * R[:]
+end
+  
+@inline function FiniteElementContainers.stiffness(
+    physics::Transient, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    K_q = (20. / dt) * N * N' + ∇N_X * ∇N_X'
+    return JxW * K_q
+end
+
+# define some helper funcs upfront
+one_func(_, _) = 1.0
+zero_func(_, _) = 0.0
+
+function simulate()
+    # load a mesh
+    mesh = UnstructuredMesh("cyl-tet.e")
+    times = TimeStepper(0., 75., 75)
+    V = FunctionSpace(mesh, H1Field, Lagrange)
+    physics = Transient()
+    props = create_properties(physics)
+    u = ScalarFunction(V, :u)
+    asm = SparseMatrixAssembler(u; use_condensed=true)
+    dbcs = [
+        DirichletBC(:u, :bottom, zero_func)
+        DirichletBC(:u, :top, one_func)
+    ]
+    U = create_field(asm)
+    p = create_parameters(
+        mesh, asm, physics, props; 
+        dirichlet_bcs = dbcs,
+        times = times
+    )
+    solver = NewtonSolver(DirectLinearSolver(asm))
+    integrator = QuasiStaticIntegrator(solver)
+    pp = PostProcessor(mesh, "output.e", u)
+
+    n = 1
+    while times.time_current[1] < 75.0
+        evolve!(integrator, p)
+        U = p.h1_field
+        write_times(pp, n, times.time_current[1])
+        write_field(pp, n, ("u",), U)
+        n = n + 1
+    end
+
+    close(pp)
+end
+
+simulate()
diff --git a/ext/AMDGPUExt.jl b/ext/AMDGPUExt.jl
index 68e1ee2..2f23bbd 100644
--- a/ext/AMDGPUExt.jl
+++ b/ext/AMDGPUExt.jl
@@ -8,6 +8,23 @@ using KernelAbstractions
 # need to double check if it's ROCArray
 FiniteElementContainers.rocm(x) = Adapt.adapt_structure(ROCArray, x)
 
+function AMDGPU.rocSPARSE.ROCSparseMatrixCOO(asm::SparseMatrixAssembler)
+  @assert typeof(get_backend(asm)) <: ROCBackend
+
+  if FiniteElementContainers._is_condensed(asm.dof)
+    n_dofs = length(asm.dof)
+  else
+    n_dofs = length(asm.dof.unknown_dofs)
+  end
+  rows, cols = asm.matrix_pattern.Is, asm.matrix_pattern.Js
+  vals = asm.stiffness_storage[asm.matrix_pattern.unknown_dofs]
+  perm = asm.matrix_pattern.permutation
+  return AMDGPU.rocSPARSE.ROCSparseMatrixCOO(
+    rows[perm], cols[perm], vals[perm],
+    (n_dofs, n_dofs), length(asm.matrix_pattern.Is)
+  )
+end
+
 # this method need to have the assembler initialized first
 # if the stored values in asm.pattern.cscnzval or zero
 # AMDGPU will error out
@@ -36,7 +53,10 @@ function AMDGPU.rocSPARSE.ROCSparseMatrixCSC(asm::SparseMatrixAssembler)
 end
 
 function FiniteElementContainers._stiffness(asm::SparseMatrixAssembler, ::ROCBackend)
-  K = AMDGPU.rocSPARSE.ROCSparseMatrixCSC(asm)
+  K = AMDGPU.rocSPARSE.ROCSparseMatrixCSC(
+    AMDGPU.rocSPARSE.ROCSparseMatrixCOO(asm)
+  )
+  # K = AMDGPU.rocSPARSE.ROCSparseMatrixCSC(asm)
 
   if FiniteElementContainers._is_condensed(asm.dof)
     FiniteElementContainers._adjust_matrix_entries_for_constraints!(
diff --git a/ext/CUDAExt.jl b/ext/CUDAExt.jl
index 3444e1c..af3216c 100644
--- a/ext/CUDAExt.jl
+++ b/ext/CUDAExt.jl
@@ -7,6 +7,23 @@ using KernelAbstractions
 
 FiniteElementContainers.cuda(x) = adapt(CuArray, x)
 
+function CUDA.CUSPARSE.CuSparseMatrixCOO(asm::SparseMatrixAssembler)
+  @assert typeof(get_backend(asm)) <: CUDABackend
+
+  if FiniteElementContainers._is_condensed(asm.dof)
+    n_dofs = length(asm.dof)
+  else
+    n_dofs = length(asm.dof.unknown_dofs)
+  end
+  rows, cols = asm.matrix_pattern.Is, asm.matrix_pattern.Js
+  vals = asm.stiffness_storage[asm.matrix_pattern.unknown_dofs]
+  perm = asm.matrix_pattern.permutation
+  return CUDA.CUSPARSE.CuSparseMatrixCOO(
+    rows[perm], cols[perm], vals[perm],
+    (n_dofs, n_dofs), length(asm.matrix_pattern.Is)
+  )
+end
+
 # this method need to have the assembler initialized first
 # if the stored values in asm.pattern.cscnzval or zero
 # CUDA will error out
@@ -31,7 +48,9 @@ function CUDA.CUSPARSE.CuSparseMatrixCSC(asm::SparseMatrixAssembler)
 end
 
 function FiniteElementContainers._stiffness(asm::SparseMatrixAssembler, ::CUDABackend)
-  K = CUDA.CUSPARSE.CuSparseMatrixCSC(asm)
+  K = CUDA.CUSPARSE.CuSparseMatrixCSC(
+    CUDA.CUSPARSE.CuSparseMatrixCOO(asm)
+  )
   if FiniteElementContainers._is_condensed(asm.dof)
     FiniteElementContainers._adjust_matrix_entries_for_constraints!(
       K, asm.constraint_storage, get_backend(asm)
diff --git a/src/FiniteElementContainers.jl b/src/FiniteElementContainers.jl
index 9b17232..87bb075 100644
--- a/src/FiniteElementContainers.jl
+++ b/src/FiniteElementContainers.jl
@@ -71,6 +71,13 @@ export InitialConditions
 export update_field_ics!
 export update_ic_values!
 
+# Integrals
+export MatrixIntegral
+export ScalarIntegral
+export VectorIntegral
+export integrate
+export remove_fixed_dofs!
+
 # Integrators
 export AbstractIntegrator
 export QuasiStaticIntegrator
@@ -115,6 +122,7 @@ export map_interpolants
 export num_properties
 export num_states
 export reshape_element_level_field
+export unpack_field
 
 export energy
 export hvp
diff --git a/src/Parameters.jl b/src/Parameters.jl
index 76fc7ba..c52ed69 100644
--- a/src/Parameters.jl
+++ b/src/Parameters.jl
@@ -139,14 +139,13 @@ function Parameters(
   )
 
   update_dofs!(assembler, p)
-  Uu = create_unknowns(assembler.dof)
-
-  # assemble the stiffness at least once for 
-  # making easier to use on GPU
-  # assemble!(assembler, Uu, p, Val{:stiffness}(), H1Field)
-  # TODO should we also assemble mass if necessary?
-  assemble_stiffness!(assembler, stiffness, Uu, p)
-  K = stiffness(assembler)
+  # Uu = create_unknowns(assembler.dof)
+
+  # # assemble the stiffness at least once for 
+  # # making easier to use on GPU
+  # # TODO should we also assemble mass if necessary?
+  # assemble_stiffness!(assembler, stiffness, Uu, p)
+  # K = stiffness(assembler)
 
   return p
 end
diff --git a/src/assemblers/SparsityPatterns.jl b/src/assemblers/SparsityPatterns.jl
index 1eb9e6c..0fc62d3 100644
--- a/src/assemblers/SparsityPatterns.jl
+++ b/src/assemblers/SparsityPatterns.jl
@@ -23,6 +23,8 @@ struct SparseMatrixPattern{
   csccolptr::I
   cscrowval::I
   cscnzval::R
+  # trying something new
+  permutation::I
 end
 
 # TODO won't work for H(div) or H(curl) yet
@@ -91,15 +93,25 @@ function SparseMatrixPattern(dof::DofManager)
   cscrowval = Vector{Int64}(undef, 0)
   cscnzval  = Vector{Float64}(undef, 0)
 
-  return SparseMatrixPattern(
+  # set permutation
+  ks = map((i, j) -> (i, j), Is, Js)
+  permutation = sortperm(ks)
+
+  pattern = SparseMatrixPattern(
     Is, Js, 
     unknown_dofs, 
     block_start_indices, block_el_level_sizes, 
     # cache arrays
     klasttouch, csrrowptr, csrcolval, csrnzval,
     # additional cache arrays
-    csccolptr, cscrowval, cscnzval
+    csccolptr, cscrowval, cscnzval,
+    permutation
   )
+
+  # # assemble this once here with a dummy vector
+  # Vs_temp = zeros(length(Is))
+  # SparseArrays.sparse!(pattern, Vs_temp)
+  return pattern
 end
 
 function Adapt.adapt_structure(to, asm::SparseMatrixPattern)
@@ -117,11 +129,12 @@ function Adapt.adapt_structure(to, asm::SparseMatrixPattern)
   csccolptr = adapt(to, asm.csccolptr)
   cscrowval = adapt(to, asm.cscrowval)
   cscnzval = adapt(to, asm.cscnzval)
+  perm = adapt(to, asm.permutation)
   return SparseMatrixPattern(
     Is, Js,
     unknown_dofs, asm.block_start_indices, asm.block_el_level_sizes,
     klasttouch, csrrowptr, csrcolval, csrnzval,
-    csccolptr, cscrowval, cscnzval
+    csccolptr, cscrowval, cscnzval, perm
   )
 end
 
@@ -181,6 +194,11 @@ function _update_dofs!(pattern::SparseMatrixPattern, dof, dirichlet_dofs)
   resize!(pattern.csrcolval, length(pattern.Is))
   resize!(pattern.csrnzval, length(pattern.Is))
 
+  ks = map((i, j) -> (i, j), pattern.Is, pattern.Js)
+  permutation = sortperm(ks)
+  resize!(pattern.permutation, length(permutation))
+  pattern.permutation .= permutation
+
   return nothing
 end
 
diff --git a/src/integrals/Integrals.jl b/src/integrals/Integrals.jl
index d09bbae..bc7eb69 100644
--- a/src/integrals/Integrals.jl
+++ b/src/integrals/Integrals.jl
@@ -1,26 +1,33 @@
 abstract type AbstractIntegral{
-    A         <: AbstractAssembler,
-    Integrand <: Function
+    A <: AbstractAssembler,
+    C,
+    I <: Function
 } end
 
+function Adapt.adapt_storage(to, int::AbstractIntegral)
+    return eval(typeof(int).name.name)(
+        adapt(to, int.assembler),
+        adapt(to, int.cache),
+        int.integrand
+    )
+end
+
+KA.get_backend(int::AbstractIntegral) = KA.get_backend(int.cache)
+
 function integrate end
 
-struct ScalarIntegral{A, I, C <: NamedTuple} <: AbstractIntegral{A, I}
+function (integral::AbstractIntegral)(U, p)
+    return integrate(integral, U, p)
+end
+
+struct ScalarIntegral{A, C <: NamedTuple, I} <: AbstractIntegral{A, C, I}
     assembler::A
     cache::C
     integrand::I
 end
 
 function ScalarIntegral(asm, integrand)
-    fspace = function_space(asm.dof)
-    cache = Matrix{Float64}[]
-    for (key, val) in pairs(fspace.ref_fes)
-        NQ = ReferenceFiniteElements.num_quadrature_points(val)
-        NE = size(getfield(fspace.elem_conns, key), 2)
-        field = L2ElementField(zeros(Float64, NQ, NE))
-        push!(cache, field)
-    end
-    cache = NamedTuple{keys(fspace.ref_fes)}(tuple(cache...))
+    cache = create_assembler_cache(asm, AssembledScalar())
     return ScalarIntegral(asm, cache, integrand)
 end
 
@@ -36,24 +43,76 @@ end
 function integrate(integral::ScalarIntegral, U, p)
     cache, dof = integral.cache, integral.assembler.dof
     func = integral.integrand
-    assemble_quadrature_quantity!(cache, dof, func, U, p)
+    assemble_quadrature_quantity!(cache, dof, nothing, func, U, p)
     return mapreduce(sum, sum, values(cache))
 end
 
-struct VectorIntegral{A, I, C <: AbstractField} <: AbstractIntegral{A, I}
+struct MatrixIntegral{A, C <: AbstractVector, I} <: AbstractIntegral{A, C, I}
+    assembler::A
+    cache::C
+    integrand::I
+end
+
+function MatrixIntegral(asm, integrand)
+    cache = create_assembler_cache(asm, AssembledMatrix())
+    return MatrixIntegral(asm, cache, integrand)
+end
+
+function integrate(integral::MatrixIntegral, U, p)
+    assemble_matrix!(
+        integral.cache, 
+        integral.assembler.matrix_pattern,
+        integral.assembler.dof,
+        integral.integrand,
+        U, p
+    )
+    return SparseArrays.sparse!(
+        integral.assembler.matrix_pattern,
+        integral.cache
+    )
+end
+
+function remove_fixed_dofs!(integral::MatrixIntegral)
+    backend = KA.get_backend(integral)
+
+    # TODO change API so we don't need to
+    # bring this guy to life
+    A = SparseArrays.sparse!(
+        integral.assembler.matrix_pattern,
+        integral.cache
+    )
+    _adjust_matrix_entries_for_constraints!(
+        A, integral.assembler.constraint_storage, backend
+    )
+    return nothing
+end
+
+struct VectorIntegral{A, C <: AbstractField, I} <: AbstractIntegral{A, C, I}
     assembler::A
     cache::C
     integrand::I
 end
 
 function VectorIntegral(asm, integrand)
-    cache = create_field(asm)
+    cache = create_assembler_cache(asm, AssembledVector())
     return VectorIntegral(asm, cache, integrand)
 end
 
 function integrate(integral::VectorIntegral, U, p)
-    cache, dof = integral.cache, integral.assembler.dof
-    func = integral.integrand
-    assemble_vector!(cache, dof, func, U, p)
-    return cache
+    assemble_vector!(
+        integral.cache, 
+        integral.assembler.vector_pattern,
+        integral.assembler.dof,
+        integral.integrand,
+        U, p
+    )
+    return integral.cache
+end
+
+function remove_fixed_dofs!(integral::VectorIntegral)
+    backend = KA.get_backend(integral)
+    _adjust_vector_entries_for_constraints!(
+        integral.cache, integral.assembler.constraint_storage, backend
+    )
+    return nothing
 end
diff --git a/src/integrators/QuasiStaticIntegrator.jl b/src/integrators/QuasiStaticIntegrator.jl
index 194ab04..59dec04 100644
--- a/src/integrators/QuasiStaticIntegrator.jl
+++ b/src/integrators/QuasiStaticIntegrator.jl
@@ -18,5 +18,9 @@ function evolve!(integrator::QuasiStaticIntegrator, p)
   update_bc_values!(p)
   solve!(integrator.solver, integrator.solution, p)
   KA.synchronize(KA.get_backend(integrator))
+  # the call below will ensure the return fields have bcs properly enfroced
+  # before being saved as the old solution for the next step.
+  _update_for_assembly!(p, integrator.solver.linear_solver.assembler.dof, integrator.solution)
+  p.h1_field_old.data .= p.h1_field.data
   return nothing
 end
diff --git a/src/meshes/Meshes.jl b/src/meshes/Meshes.jl
index 58cbada..017e713 100644
--- a/src/meshes/Meshes.jl
+++ b/src/meshes/Meshes.jl
@@ -8,6 +8,7 @@ const elem_type_map = Dict{String, Type{<:ReferenceFiniteElements.AbstractElemen
   "TRI3"    => Tri3,
   "TRI6"    => Tri6,
   "TET"     => Tet4,
+  "TETRA"   => Tet4,
   "TETRA4"  => Tet4,
   "TETRA10" => Tet10
 )
diff --git a/src/physics/Physics.jl b/src/physics/Physics.jl
index e38f00d..697bf18 100644
--- a/src/physics/Physics.jl
+++ b/src/physics/Physics.jl
@@ -72,6 +72,42 @@ $(TYPEDSIGNATURES)
   return SMatrix{NDof, N, T, NxNDof}(u_el)
 end
 
+"""
+$(TYPEDSIGNATURES)
+Unpacks a single fields values from a SMatrix.
+This is useful for extracting specific components
+from an interpolated field gradient at a
+quadrature point ∇u_q.
+
+Return a SVector
+"""
+@inline function unpack_field(field::SMatrix{M, N, T, L}, dof::Int) where {M, N, T <: Number, L}
+  ids = SVector{N, Int}(1:N)
+  return SVector{N, T}(field[dof, ids])
+end
+
+"""
+$(TYPEDSIGNATURES)
+Unpacks a range of fields values from a SMatrix.
+This is useful for extracting specific components
+from an interpolated field gradient at a
+quadrature point ∇u_q.
+
+returns a SMatrix.
+
+Note the Val{D} that is a necessary input. This is
+crucial for performance with ```StaticArrays```.
+"""
+@inline function unpack_field(
+  field::SMatrix{M, N, T, L}, dof_start::Int, dof_end::Int,
+  ::Val{D}
+) where {M, N, T <: Number, L, D}
+  @assert D == dof_end - dof_start + 1 "Number of expected dofs wrapped in Val is not equal to dof_end - dof_start"
+  ids = SVector{N, Int}(1:N)
+  dofs = SVector{D, Int}(dof_start:dof_end)
+  return SMatrix{D, N, T, D * N}(field[dofs, ids])
+end
+
 # Can we make something that makes interfacing with kernels easier?
 # How can we make something like this work nicely with AD?
 # struct PhysicsQuadratureState{T, ND, NN, NF, NP, NS, NDxNF, NNxND, NNxNNxND}