diff --git a/src/assembly.jl b/src/assembly.jl
index 05a04bc..4756283 100644
--- a/src/assembly.jl
+++ b/src/assembly.jl
@@ -95,3 +95,79 @@ function assemble_mass_matrix!{Basis}(problem::Problem, elements::Vector{Element
     return
 end
 
+"""
+    assemble_mass_matrix!(problem, elements::Vector{Element{Tet10}}, time)
+
+Assemble Tet10 mass matrices using special method. If Tet10 has constant metric
+if can be integrated analytically to gain performance.
+"""
+function assemble_mass_matrix!(problem::Problem, elements::Vector{Element{Tet10}}, time)
+    nnodes = length(Tet10)
+    dim = get_unknown_field_dimension(problem)
+    M = zeros(nnodes, nnodes)
+    N = zeros(1, nnodes)
+    NtN = zeros(nnodes, nnodes)
+    ldofs = zeros(Int, nnodes)
+
+    M_CM = 1.0/2520.0 * [
+         6  1  1  1 -4 -6 -4 -4 -6 -6
+         1  6  1  1 -4 -4 -6 -6 -4 -6
+         1  1  6  1 -6 -4 -4 -6 -6 -4
+         1  1  1  6 -6 -6 -6 -4 -4 -4
+        -4 -4 -6 -6 32 16 16 16 16  8
+        -6 -4 -4 -6 16 32 16  8 16 16
+        -4 -6 -4 -6 16 16 32 16  8 16
+        -4 -6 -6 -4 16  8 16 32 16 16
+        -6 -4 -6 -4 16 16  8 16 32 16
+        -6 -6 -4 -4  8 16 16 16 16 32]
+
+    function is_CM(element::Element{Tet10}, X; rtol=1.0e-6)
+        isapprox(X[5],  1/2*(X[1]+X[2]); rtol=rtol) || return false
+        isapprox(X[6],  1/2*(X[2]+X[3]); rtol=rtol) || return false
+        isapprox(X[7],  1/2*(X[3]+X[1]); rtol=rtol) || return false
+        isapprox(X[8],  1/2*(X[1]+X[4]); rtol=rtol) || return false
+        isapprox(X[9],  1/2*(X[2]+X[4]); rtol=rtol) || return false
+        isapprox(X[10], 1/2*(X[3]+X[4]); rtol=rtol) || return false
+        return true
+    end
+
+
+    n_CM = 0
+    for element in elements
+        for (i, j) in enumerate(get_connectivity(element))
+            @inbounds ldofs[i] = (j-1)*dim
+        end
+
+        X = element("geometry", time)
+        rho = element("density", time)
+        if is_CM(element, X) && length(rho) == 1
+            ip = (1.0/3.0, 1.0/3.0, 1.0/3.0)
+            detJ = element(ip, time, Val{:detJ})
+            rho = element("density", ip, time)
+            CM_s = detJ*rho
+            n_CM += 1
+            for i=1:dim
+                add!(problem.assembly.M, ldofs+i, ldofs+i, CM_s * M_CM)
+            end
+        else
+            fill!(M, 0.0)
+            for ip in get_integration_points(element, 2)
+                detJ = element(ip, time, Val{:detJ})
+                rho = element("density", ip, time)
+                w = ip.weight*rho*detJ
+                eval_basis!(Tet10, N, ip)
+                N = element(ip, time)
+                At_mul_B!(NtN, N, N)
+                scale!(NtN, w)
+                for i=1:nnodes^2
+                    M[i] += NtN[i]
+                end
+            end
+            for i=1:dim
+                add!(problem.assembly.M, ldofs+i, ldofs+i, M)
+            end
+        end
+    end
+    info("$n_CM of $(length(elements)) was constant metric.")
+    return
+end