diff --git a/examples/cooks_membrane/script.jl b/examples/cooks_membrane/script.jl
index dc70f67a..f8fcb672 100644
--- a/examples/cooks_membrane/script.jl
+++ b/examples/cooks_membrane/script.jl
@@ -165,8 +165,10 @@ FEC.update_dofs!(
     (p_u.periodic_bcs, p_p.periodic_bcs)
 )
 
-pp_u = PostProcessor(mesh_u, "u.exo", u)
-pp_p = PostProcessor(mesh_p, "p.exo", p)
+# pp_u = PostProcessor(mesh_u, "u.exo", u)
+# pp_p = PostProcessor(mesh_p, "p.exo", p)
+
+
 # solver = NewtonSolver(DirectLinearSolver(asm))
 # integrator = QuasiStaticIntegrator(solver)
 
@@ -186,39 +188,49 @@ Uu = create_unknowns(asm)
 # temp = K_up .- K_pu'
 # display(K_up .- K_pu')
 
-for n in 1:10
-    FiniteElementContainers.update_time!(params)
-    FiniteElementContainers.update_bc_values!(params, asm)
-    r0 = -1e6
-    for iter in 1:10
-        assemble_vector!(asm, residual, Uu, params)
-        R = residual(asm)
+function test_func(asm, u, p)
+    FiniteElementContainers.update_time!(p)
+    FiniteElementContainers.update_bc_values!(p, asm)
+    assemble_vector!(asm, residual, u, p)
+    return nothing
+end
+
+@time test_func(asm, Uu, params)
+@time test_func(asm, Uu, params)
+
+# for n in 1:10
+#     FiniteElementContainers.update_time!(params)
+#     FiniteElementContainers.update_bc_values!(params, asm)
+#     r0 = -1e6
+#     for iter in 1:10
+#         assemble_vector!(asm, residual, Uu, params)
+#         R = residual(asm)
     
-        rnorm = norm(R)
-        if iter == 1
-            r0 = rnorm
-        end
-
-        if rnorm / r0 < 1e-8
-            break
-        end
+#         rnorm = norm(R)
+#         if iter == 1
+#             r0 = rnorm
+#         end
+
+#         if rnorm / r0 < 1e-8
+#             break
+#         end
     
-        assemble_stiffness!(asm, stiffness, Uu, params)
-        K = stiffness(asm)
+#         assemble_stiffness!(asm, stiffness, Uu, params)
+#         K = stiffness(asm)
     
-        ΔU = K \ R
-        Uu .-= ΔU
+#         ΔU = K \ R
+#         Uu .-= ΔU
 
-        println("iter = $iter, |R| = $(rnorm / r0), |ΔU| = $(norm(ΔU))")
+#         println("iter = $iter, |R| = $(rnorm / r0), |ΔU| = $(norm(ΔU))")
     
-    end
+#     end
 
-    write_times(pp_u, n + 1, params[1].times.time_current)
-    write_times(pp_p, n + 1, params[2].times.time_current)
-    write_field(pp_u, n + 1, ("displ_x", "displ_y"), params[1].field)
-    write_field(pp_p, n + 1, ("pressure",), params[2].field)
-end
-close(pp_u)
-close(pp_p)
+#     write_times(pp_u, n + 1, params[1].times.time_current)
+#     write_times(pp_p, n + 1, params[2].times.time_current)
+#     write_field(pp_u, n + 1, ("displ_x", "displ_y"), params[1].field)
+#     write_field(pp_p, n + 1, ("pressure",), params[2].field)
+# end
+# close(pp_u)
+# close(pp_p)
 
 # Δu = K \ R
diff --git a/examples/mpi-example/example.jl b/examples/mpi-example/example.jl
deleted file mode 100644
index 3a6cfc1b..00000000
--- a/examples/mpi-example/example.jl
+++ /dev/null
@@ -1,48 +0,0 @@
-import FiniteElementContainers:
-    communication_graph,
-    cpu, 
-    decompose_mesh,
-    distribute,
-    global_colorings,
-    DistributedDevice,
-    getdata,
-    num_dofs_per_rank,
-    ParVector,
-    rank_devices,
-    scatter_ghosts!
-using Exodus
-using KernelAbstractions
-using MPI
-
-backend = CPU()
-num_dofs = 1
-num_ranks = 4
-ranks = 1:4 |> collect
-mesh_file = Base.source_dir() * "/square.g"
-
-decompose_mesh(mesh_file, num_ranks)
-global_dofs_to_colors = global_colorings(mesh_file, num_dofs, num_ranks)
-
-comm = MPI.COMM_WORLD
-# comm = nothing
-
-ranks = distribute(ranks, comm)
-# n_dofs_per_rank = distribute(num_dofs_per_rank(global_dofs_to_colors), comm)
-
-# shards = shard_indices(global_dofs_to_colors, )
-comm_graphs = map(ranks) do rank
-    mesh_file = Base.source_dir() * "/square.g"
-    mesh_file = mesh_file * ".$num_ranks" * ".$(lpad(rank - 1, Exodus.exodus_pad(num_ranks |> Int32), '0'))"
-    comm_graph = communication_graph(mesh_file, global_dofs_to_colors, num_dofs, rank)
-end
-
-par_vec = ParVector(comm_graphs)
-parts = map(par_vec.parts, ranks) do part, rank
-    part .= rank
-    part
-end
-par_vec = ParVector(parts)
-par_vec_2 = scatter_ghosts!(par_vec)
-MPI.Barrier(comm)
-# par_vec_2 = par_vec_2.parts
-print("part on $(MPI.Comm_rank(comm) + 1) is $par_vec_2\n")
diff --git a/examples/mpi-example/square.g b/examples/mpi-example/square.g
deleted file mode 100644
index f41996e8..00000000
Binary files a/examples/mpi-example/square.g and /dev/null differ
diff --git a/examples/poisson/script.jl b/examples/poisson/script.jl
index 1ed8f848..07aba849 100644
--- a/examples/poisson/script.jl
+++ b/examples/poisson/script.jl
@@ -35,7 +35,7 @@ physics = Poisson(f)
 props = create_properties(physics)
 u = ScalarFunction(V, "u")
 # asm = SparseMatrixAssembler(u; use_condensed = true)
-asm = SparseMatrixAssembler(u)
+asm = SparseMatrixAssembler(u; use_inplace_methods = true)
 
 
 # dbcs = nothing
diff --git a/src/Formulations.jl b/src/Formulations.jl
index e47a0832..b84e510a 100644
--- a/src/Formulations.jl
+++ b/src/Formulations.jl
@@ -191,6 +191,99 @@ function scatter_with_gradients_and_gradients!(
   return nothing
 end
 
+"""
+General stiffness like matrix scatter method
+"""
+function scatter_with_gradients_and_gradients!(
+  storage::AbstractVector,
+  form::AbstractElementFormulation{ND, NF},
+  e,
+  conns,
+  ∇N_X,
+  K::AbstractMatrix{T},
+) where {ND, NF, T <: Number}
+
+  @assert size(K, 1) == NF^2
+  @assert size(K, 2) == NF^2
+  @assert ND == size(∇N_X, 2)
+
+  N    = size(∇N_X, 1)
+  NDPE = N * NF
+
+  start_id = (e - 1) * NDPE^2 + 1
+  ids      = start_id:(start_id + NDPE^2 - 1)
+
+  inc = 1
+
+  for n2 in 1:N
+    for d2 in 1:NF
+      for n1 in 1:N
+        for d1 in 1:NF
+          contrib = zero(T)
+
+          for J in 1:ND
+            a = (J - 1) * ND + d1
+
+            for L in 1:ND
+              b = (L - 1) * ND + d2
+              contrib += ∇N_X[n1, J] * K[a, b] * ∇N_X[n2, L]
+            end
+          end
+
+          storage[ids[inc]] += contrib
+          inc += 1
+        end
+      end
+    end
+  end
+
+  return nothing
+end
+
+function scatter_gradient_gradient!(
+  storage::AbstractField,
+  form::AbstractElementFormulation{ND, NF},
+  e,
+  conns,
+  ∇N_X,
+  A::AbstractMatrix{T},
+  v_el::AbstractVector{T},
+) where {ND, NF, T <: Number}
+
+  @assert ND == size(∇N_X, 2)
+  @assert size(A) == (ND * NF, ND * NF)
+
+  N = size(∇N_X, 1)
+  NDPE = N * NF
+
+  @assert length(v_el) == NDPE
+
+  for n1 in 1:N
+    for d1 in 1:NF
+      global_id = NF * (conns[n1] - 1) + d1
+      contrib = zero(T)
+
+      for n2 in 1:N
+        for d2 in 1:NF
+          local_id = NF * (n2 - 1) + d2
+          for j1 in 1:ND
+            a = (j1 - 1) * NF + d1
+
+            for j2 in 1:ND
+                b = (j2 - 1) * NF + d2
+                contrib += ∇N_X[n1, j1] * A[a, b] * ∇N_X[n2, j2] * v_el[local_id]
+            end
+          end
+        end
+      end
+
+      fec_atomic_add!(storage, global_id, contrib)
+    end
+  end
+
+  return nothing
+end
+
 # implement for those that have it
 """
 Scalar equation specialization
diff --git a/test/poisson/TestPoissonCommon.jl b/test/poisson/TestPoissonCommon.jl
index 69b07d55..b48df482 100644
--- a/test/poisson/TestPoissonCommon.jl
+++ b/test/poisson/TestPoissonCommon.jl
@@ -1,6 +1,8 @@
 # since this type of physics has no state
 # just return the old state which will be empty
 
+using ReferenceFiniteElements
+
 struct Poisson{F <: Function} <: AbstractPhysics{1, 0, 0}
   func::F
 end
@@ -91,7 +93,6 @@ end
   interps = map_interpolants(interps, x_el)
   (; X_q, N, ∇N_X, JxW) = interps
   ∇u_q = interpolate_field_gradients(physics, interps, u_el)
-  # R_q = ∇u_q * ∇N_X' - N' * physics.func(X_q, 0.0)
   form = GeneralFormulation{size(X_q, 1), num_fields(physics)}()
   scatter_with_gradients!(storage, form, e, conn, ∇N_X, JxW * ∇u_q)
   scatter_with_values!(storage, form, e, conn, N, -JxW * physics.func(X_q, 0.0))