diff --git a/src/preprocess.jl b/src/preprocess.jl
index 2698986..9d04d36 100644
--- a/src/preprocess.jl
+++ b/src/preprocess.jl
@@ -21,10 +21,11 @@ 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())
+    return Mesh(Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), nothing)
 end
 
 """
@@ -319,3 +320,82 @@ function JuliaFEM.Problem(mesh::Mesh, ::Type{P}, name, dimension, parent_field_n
     problem.elements = create_elements(mesh, name)
     return problem
 end
+
+"""
+    create_coloring!(mesh::Mesh)
+
+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.
+"""
+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
+        for v in nodes
+            if !haskey(cell_containing_node, v)
+                cell_containing_node[v] = Set{Int}()
+            end
+            push!(cell_containing_node[v], cellid)
+        end
+    end
+
+    I, J, V = Int[], Int[], Bool[]
+    for (node, cells) in cell_containing_node
+        for cell1 in cells # All these cells have a neighboring node
+            for cell2 in cells
+                if cell1 != cell2
+                    push!(I, cell1)
+                    push!(J, cell2)
+                    push!(V, true)
+                end
+            end
+        end
+    end
+
+    incidence_matrix = sparse(I, J, V)
+    # cell -> color of cell
+    cell_colors = Dict{Int, Int}()
+    # color -> list of cells
+    final_colors = Vector{Int}[]
+    occupied_colors = Set{Int}()
+    # Zero represents no color set yet
+    for (cellid, _) in mesh.elements
+        cell_colors[cellid] = 0
+    end
+    total_colors = 0
+    for (cellid, _) in mesh.elements
+        empty!(occupied_colors)
+        # loop over neighbors
+        for r in nzrange(incidence_matrix, cellid)
+            cell_neighbour = incidence_matrix.rowval[r]
+            color = cell_colors[cell_neighbour]
+            if color != 0
+                push!(occupied_colors, color)
+            end
+        end
+
+        # occupied colors now contains all the colors we are not allowed to use
+        free_color = 0
+        for attempt_color in 1:total_colors
+            if attempt_color ∉ occupied_colors
+                free_color = attempt_color
+                break
+            end
+        end
+
+        if free_color == 0 # no free color found, need to bump max colors
+            total_colors += 1
+            free_color = total_colors
+            push!(final_colors, Int[])
+        end
+
+        cell_colors[cellid] = free_color
+        push!(final_colors[free_color], cellid)
+    end
+
+    mesh.coloring = final_colors
+    return mesh
+end
diff --git a/test/test_mesh_coloring.jl b/test/test_mesh_coloring.jl
new file mode 100644
index 0000000..b1be331
--- /dev/null
+++ b/test/test_mesh_coloring.jl
@@ -0,0 +1,20 @@
+using JuliaFEM, LinearAlgebra, Test
+
+datadir = first(splitext(basename(@__FILE__)))
+
+@testset "Test that mesh coloring provides a good coloring"
+    fn = joinpath(datadir, "cube_tet4.inp")
+    mesh = JuliaFEM.Mesh(open(parse_abaqus, fn))
+
+    JuliaFEM.create_coloring!(mesh)
+    for colors in mesh.coloring
+        for ele_i in colors
+            for ele_j in colors
+                if ele_i == ele_j
+                    continue
+                end
+                @test intersect(Set(mesh.elements[ele_i]), Set(mesh.elements[ele_j])) |> isempty
+            end
+        end
+    end
+end
\ No newline at end of file