diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index 71d5ef5..0254349 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -108,6 +108,8 @@ module JuliaFEM
 using SparseArrays, LinearAlgebra, Statistics
 using Reexport, ForwardDiff, LightXML, HDF5, Parameters
 
+import FEMSparse
+
 @reexport using FEMBase
 import FEMBase: get_unknown_field_name, get_unknown_field_dimension,
                 assemble!, update!, initialize!
diff --git a/src/problems_dirichlet.jl b/src/problems_dirichlet.jl
index 441e7d9..183a924 100644
--- a/src/problems_dirichlet.jl
+++ b/src/problems_dirichlet.jl
@@ -88,9 +88,9 @@ function assemble!(problem::Problem{Dirichlet}, time::Float64=0.0;
             end
         end
         for (k, v) in field_vals
-            add!(problem.assembly.C1, k, k, 1.0)
-            add!(problem.assembly.C2, k, k, 1.0)
-            add!(problem.assembly.g, k, 1, v)
+            FEMBase.add!(problem.assembly.C1, k, k, 1.0)
+            FEMBase.add!(problem.assembly.C2, k, k, 1.0)
+            FEMBase.add!(problem.assembly.g, k, 1, v)
         end
     end
 
diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl
index 4993479..9730fc8 100644
--- a/src/problems_elasticity.jl
+++ b/src/problems_elasticity.jl
@@ -73,7 +73,7 @@ end
 function assemble!(assembly::Assembly, problem::Problem{Elasticity},
                    elements::Vector{<:Element}, time, formulation)
     local_buffer = allocate_buffer(problem, elements)
-    assembler = FEMBase.start_assemble(assembly.K)
+    assembler = FEMSparse.start_assemble(assembly.K, assembly.f)
     for element in elements
         assemble_element!(assembly, assembler, problem, element, local_buffer, time, formulation)
     end
@@ -168,7 +168,7 @@ end
 
 """ Assemble 3d continuum elements in general solid mechanics problem. """
 function assemble_element!(assembly::Assembly,
-                   assembler::FEMBase.AssemblerSparsityPattern,
+                   assembler::FEMSparse.AssemblerSparsityPattern,
                    problem::Problem{Elasticity},
                    element::Element{El},
                    local_buffer::Elasticity3DLocalBuffers,
@@ -341,12 +341,9 @@ function assemble_element!(assembly::Assembly,
         # internal load
         mul!(Bt_mul_S, transpose(BL), stress_vec)
         rmul!(Bt_mul_S, w)
-        for i=1:ndofs
-            @inbounds f_int[i] += Bt_mul_S[i]
-        end
+        f_int .+= Bt_mul_S
 
         # external load start
-
         if haskey(element, "displacement load")
             T = element("displacement load", ip, time)
             f_ext += w*vec(T*N)
@@ -358,22 +355,21 @@ function assemble_element!(assembly::Assembly,
                 f_ext[i:dim:end] += w*vec(b*N)
             end
         end
-
         # external load end
-
     end
 
     gdofs = get_gdofs(problem, element)
 
+    # Update f_ext in place to be f_ext - f_int
+    f_ext .-= f_int
+
     # add contributions to K, Kg, f
-    FEMBase.assemble_local_matrix!(assembler, gdofs, Km)
+    FEMSparse.assemble_local!(assembler, gdofs, Km, f_ext)
 
     if props.geometric_stiffness
-        add!(assembly.Kg, gdofs, gdofs, Kg)
+        FEMSparse.assemble_local_matrix!(assembler, gdofs, Kg)
     end
 
-    add!(assembly.f, gdofs, f_ext - f_int)
-
     return nothing
 end
 
@@ -424,8 +420,7 @@ function assemble!(assembly::Assembly,
         end
 
         gdofs = get_gdofs(problem, element)
-        add!(assembly.f, gdofs, f)
-
+        FEMSparse.assemble_local_vector!(assembly.f, gdofs, f)
     end
 end
 
diff --git a/src/solvers.jl b/src/solvers.jl
index 98399ba..552feb7 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -58,17 +58,17 @@ function get_field_assembly(solver::Solver)
     problems = get_field_problems(solver)
 
     problem = problems[1]
-    M = SparseMatrixCOO()
-    K = spzeros(size(problems[1].assembly.K)...)
-    Kg = SparseMatrixCOO()
-    f = SparseMatrixCOO()
-    fg = SparseMatrixCOO()
+    M = problem.assembly.M
+    K = problem.assembly.K
+    f = problem.assembly.f
+    Kg = problem.assembly.Kg
+    fg = problem.assembly.fg
 
-    for problem in problems
+    for problem in problems[2:end]
         append!(M, problem.assembly.M)
         K += problem.assembly.K
         append!(Kg, problem.assembly.Kg)
-        append!(f, problem.assembly.f)
+        f += problem.assembly.f
         append!(fg, problem.assembly.fg)
     end
 
@@ -80,7 +80,6 @@ function get_field_assembly(solver::Solver)
               "pushed to problem and formulation is correct.")
     end
     Kg = sparse(Kg, N, N)
-    f = sparse(f, N, 1)
     fg = sparse(fg, N, 1)
 
     return M, K, Kg, f, fg
@@ -146,7 +145,7 @@ function get_boundary_assembly(solver::Solver, N)
         C1_ = sparse(assembly.C1, N, N)
         C2_ = sparse(assembly.C2, N, N)
         D_ = sparse(assembly.D, N, N)
-        f_ = sparse(assembly.f, N, 1)
+        # f_ = assembly.f
         g_ = sparse(assembly.g, N, 1)
         for dof in assembly.removed_dofs
             @info("$(problem.name): removing dof $dof from assembly")
@@ -171,7 +170,7 @@ function get_boundary_assembly(solver::Solver, N)
         C1 += C1_
         C2 += C2_
         D += D_
-        f += f_
+        #f += f_
         g += g_
     end
     return K, C1, C2, D, f, g