From 2bb66a9568f3940600a69251fe3686aeebd47925 Mon Sep 17 00:00:00 2001
From: Jukka Aho <ahojukka5@gmail.com>
Date: Tue, 22 Aug 2017 08:15:54 +0300
Subject: [PATCH] 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 | 76 +++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 76 insertions(+)

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