Special handling of constant metric Tet10

Mass matrix can be analytically solved if Tet10 metric is constant, i.e.
the midnodes are in the midpoint of corner nodes. This should increase
assembling speed of mass matrix.

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