diff --git a/src/preprocess.jl b/src/preprocess.jl
index 9d04d36..57e1bea 100644
--- a/src/preprocess.jl
+++ b/src/preprocess.jl
@@ -21,11 +21,10 @@ mutable struct Mesh
     element_sets :: Dict{Symbol, Set{Int}}
     surface_sets :: Dict{Symbol, Vector{Tuple{Int, Symbol}}}
     surface_types :: Dict{Symbol, Symbol}
-    coloring::Union{Nothing, Vector{Vector{Int}}} # Each vector contains a list of elements that do not share nodes
 end
 
 function Mesh()
-    return Mesh(Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), nothing)
+    return Mesh(Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict())
 end
 
 """
@@ -322,15 +321,14 @@ function JuliaFEM.Problem(mesh::Mesh, ::Type{P}, name, dimension, parent_field_n
 end
 
 """
-    create_coloring!(mesh::Mesh)
+    create_coloring!(mesh::Mesh) -> Dict{Int, Int}
 
 Greedy algorithm for coloring a grid such that no two cells with the same node
 have the same color.
-This function sets the `coloring` field in `mesh` to a `Vector{Vector{Int}}` where
-each vector contains vectors of elements that do not share any nodes.
-It is therefore safe to assemble in parallel each element vector by vector.
+The returned value is a mapping between an element id and its color.
+It is safe to assemble elements with the same color in parallel
 """
-function create_coloring!(mesh::Mesh)
+function create_coloring(mesh::Mesh)
     # Contains the elements that each node contain
     cell_containing_node = Dict{Int, Set{Int}}()
     for (cellid, nodes) in mesh.elements
@@ -359,7 +357,7 @@ function create_coloring!(mesh::Mesh)
     # cell -> color of cell
     cell_colors = Dict{Int, Int}()
     # color -> list of cells
-    final_colors = Vector{Int}[]
+    final_colors = Set{Int}[]
     occupied_colors = Set{Int}()
     # Zero represents no color set yet
     for (cellid, _) in mesh.elements
@@ -389,13 +387,12 @@ function create_coloring!(mesh::Mesh)
         if free_color == 0 # no free color found, need to bump max colors
             total_colors += 1
             free_color = total_colors
-            push!(final_colors, Int[])
+            push!(final_colors, Set{Int}())
         end
 
         cell_colors[cellid] = free_color
         push!(final_colors[free_color], cellid)
     end
 
-    mesh.coloring = final_colors
-    return mesh
+    return cell_colors
 end
diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl
index 65cdc5c..e2b773d 100644
--- a/src/problems_elasticity.jl
+++ b/src/problems_elasticity.jl
@@ -73,11 +73,30 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity},
 end
 
 function assemble!(assembly::Assembly, problem::Problem{Elasticity},
-                   elements::Vector{<:Element}, time, formulation)
-    local_buffer = allocate_buffer(problem, elements)
-    assembler = FEMSparse.start_assemble(assembly.K, assembly.f)
-    for i in 1:length(elements)
-        assemble_element!(assembly, assembler, problem, elements[i], local_buffer, time, formulation)
+                   elements::Vector{T}, time, formulation) where {T <: Element}
+
+    if problem.assemble_parallel
+        # Threaded assembly
+
+        assemblers = [FEMSparse.start_assemble(assembly.K, assembly.f) for i in 1:Threads.nthreads()]
+        local_buffers = [allocate_buffer(problem, elements) for i in 1:Threads.nthreads()]
+        #TODO: We have to be a bit careful here, the index of the element is no longer
+        #
+        # should only loop over elements that exist in `elements` here
+        for (color, elements) in FEMBase.get_color_ranges(elements)
+            Threads.@threads for i in 1:length(elements)
+                element = elements[i]
+                tid = Threads.threadid()
+                assemble_element!(assembly, assemblers[tid], problem, element, local_buffers[tid], time, formulation)
+            end
+        end
+    else
+        # Normal assembly
+        local_buffer = allocate_buffer(problem, elements)
+        assembler = FEMSparse.start_assemble(assembly.K, assembly.f)
+        for i in 1:length(elements)
+            assemble_element!(assembly, assembler, problem, elements[i], local_buffer, time, formulation)
+        end
     end
 end
 
@@ -247,6 +266,8 @@ function assemble_element!(assembly::Assembly,
 
         E = element("youngs modulus", ip, time)::Float64
         nu = element("poissons ratio", ip, time)::Float64
+        #E = 200e3
+        #nu = 0.3
         la = E*nu/((1.0+nu)*(1.0-2.0*nu))
         mu = E/(2.0*(1.0+nu))
         D[1,1] = D[2,2] = D[3,3] = 2*mu + la
diff --git a/test/test_mesh_coloring.jl b/test/test_mesh_coloring.jl
index b1be331..5e43966 100644
--- a/test/test_mesh_coloring.jl
+++ b/test/test_mesh_coloring.jl
@@ -6,8 +6,8 @@ datadir = first(splitext(basename(@__FILE__)))
     fn = joinpath(datadir, "cube_tet4.inp")
     mesh = JuliaFEM.Mesh(open(parse_abaqus, fn))
 
-    JuliaFEM.create_coloring!(mesh)
-    for colors in mesh.coloring
+    coloring = JuliaFEM.create_coloring(mesh)
+    for colors in coloring
         for ele_i in colors
             for ele_j in colors
                 if ele_i == ele_j