From ce6ee3855d3b5eb6b3b86c7828441b727a6887ac Mon Sep 17 00:00:00 2001
From: Jukka Aho <ahojukka5@gmail.com>
Date: Sun, 13 Aug 2017 19:57:19 +0300
Subject: [PATCH] Add new type BasisInfo

BasisInfo is an optimized version to calculate all necessary stuff using
basis functions, including e.g. basis, basis derivatives wrt to spatial
coordinates, Jacobian and it's inverse and determinant. Everything is
nicely packes to `BasisInfo`.
---
 src/FEMBasis.jl   |   1 +
 src/math.jl       | 106 ++++++++++++++++++++++++++++++++++++++++++++++
 test/test_math.jl |  16 +++++++
 3 files changed, 123 insertions(+)

diff --git a/src/FEMBasis.jl b/src/FEMBasis.jl
index 357a0a4..150df92 100644
--- a/src/FEMBasis.jl
+++ b/src/FEMBasis.jl
@@ -35,5 +35,6 @@ export NSolid
 include("math.jl")
 export interpolate, interpolate!
 export grad, grad!
+export BasisInfo
 
 end
diff --git a/src/math.jl b/src/math.jl
index c43aacc..18d6969 100644
--- a/src/math.jl
+++ b/src/math.jl
@@ -132,3 +132,109 @@ function grad(B, T, X, xi)
     dTdX = sum(kron(G[:,i], T[i]') for i=1:nbasis)
     return dTdX'
 end
+
+
+"""
+Data type for fast FEM.
+"""
+type BasisInfo{B<:AbstractBasis,T}
+    N::Matrix{T}
+    dN::Matrix{T}
+    grad::Matrix{T}
+    J::Matrix{T}
+    invJ::Matrix{T}
+    detJ::T
+end
+
+"""
+Initialization of data type `BasisInfo`.
+
+# Examples
+
+```jldoctest
+
+BasisInfo(Tri3)
+
+# output
+
+FEMBasis.BasisInfo{FEMBasis.Tri3,Float64}([0.0 0.0 0.0], [0.0 0.0 0.0; 0.0 0.0 0.0], [0.0 0.0 0.0; 0.0 0.0 0.0], [0.0 0.0; 0.0 0.0], [0.0 0.0; 0.0 0.0], 0.0)
+
+```
+
+"""
+function BasisInfo{B<:AbstractBasis}(::Type{B}, T=Float64)
+    dim, nbasis = size(B)
+    N = zeros(T, 1, nbasis)
+    dN = zeros(T, dim, nbasis)
+    grad = zeros(T, dim, nbasis)
+    J = zeros(T, dim, dim)
+    invJ = zeros(T, dim, dim)
+    detJ = zero(T)
+    return BasisInfo{B,T}(N, dN, grad, J, invJ, detJ)
+end
+
+"""
+Evaluate basis, gradient and so on for some point `xi`.
+
+# Examples
+
+```jldoctest
+
+b = BasisInfo(Quad4)
+X = ((0.0,0.0), (1.0,0.0), (1.0,1.0), (0.0,1.0))
+xi = (0.0, 0.0)
+eval_basis!(b, X, xi)
+
+# output
+
+FEMBasis.BasisInfo{FEMBasis.Quad4,Float64}([0.25 0.25 0.25 0.25], [-0.25 0.25 0.25 -0.25; -0.25 -0.25 0.25 0.25], [-0.5 0.5 0.5 -0.5; -0.5 -0.5 0.5 0.5], [0.5 0.0; 0.0 0.5], [2.0 -0.0; -0.0 2.0], 0.25)
+
+```
+"""
+function eval_basis!{B}(bi::BasisInfo{B}, X, xi)
+
+    # evaluate basis and derivatives
+    eval_basis!(B, bi.N, xi)
+    eval_dbasis!(B, bi.dN, xi)
+
+    # calculate Jacobian
+    fill!(bi.J, 0.0)
+    dim, nbasis = size(B)
+    for i=1:nbasis
+        for j=1:dim
+            for k=1:dim
+                @inbounds bi.J[j,k] += bi.dN[j,i]*X[i][k]
+            end
+        end
+    end
+
+    # calculate inverse and determinant of Jacobian
+    if dim == 3
+        a, b, c, d, e, f, g, h, i = bi.J
+        bi.detJ = a*(e*i-f*h) + b*(f*g-d*i) + c*(d*h-e*g)
+        bi.invJ[1] = 1.0 / bi.detJ * (e*i - f*h)
+        bi.invJ[2] = 1.0 / bi.detJ * (c*h - b*i)
+        bi.invJ[3] = 1.0 / bi.detJ * (b*f - c*e)
+        bi.invJ[4] = 1.0 / bi.detJ * (f*g - d*i)
+        bi.invJ[5] = 1.0 / bi.detJ * (a*i - c*g)
+        bi.invJ[6] = 1.0 / bi.detJ * (c*d - a*f)
+        bi.invJ[7] = 1.0 / bi.detJ * (d*h - e*g)
+        bi.invJ[8] = 1.0 / bi.detJ * (b*g - a*h)
+        bi.invJ[9] = 1.0 / bi.detJ * (a*e - b*d)
+    elseif dim == 2
+        a, b, c, d = bi.J
+        bi.detJ = a*d - b*c
+        bi.invJ[1] = 1.0 / bi.detJ * d
+        bi.invJ[2] = 1.0 / bi.detJ * -b
+        bi.invJ[3] = 1.0 / bi.detJ * -c
+        bi.invJ[4] = 1.0 / bi.detJ * a
+    else
+        @assert dim == 1
+        bi.detJ = bi.J[1,1]
+        bi.invJ[1,1] = 1.0 / bi.detJ
+    end
+    
+    A_mul_B!(bi.grad, bi.invJ, bi.dN)
+
+    return bi
+end
diff --git a/test/test_math.jl b/test/test_math.jl
index c26b1b3..1c56378 100644
--- a/test/test_math.jl
+++ b/test/test_math.jl
@@ -50,3 +50,19 @@ end
     dudX_expected(X) = [X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
     @test isapprox(dudX, dudX_expected([0.5, 0.5]))
 end
+
+@testset "test BasisInfo" begin
+    B = BasisInfo(Quad4)
+    X = ((0.0,0.0), (1.1,0.0), (1.0,1.0), (0.0,1.0))
+    xi = (0.0,0.0)
+    eval_basis!(B, X, xi)
+    @test isapprox(B.invJ, inv(B.J))
+    @test isapprox(B.detJ, det(B.J))
+    B = BasisInfo(Hex8)
+    X = ((0.0,0.0,0.0), (1.1,0.0,0.0), (1.0,1.0,0.0), (0.0,1.0,0.0),
+         (0.0,0.0,1.0), (1.0,0.0,1.0), (1.0,1.0,1.0), (0.0,1.0,1.0))
+    xi = (0.0,0.0,0.0)
+    eval_basis!(B, X, xi)
+    @test isapprox(B.invJ, inv(B.J))
+    @test isapprox(B.detJ, det(B.J))
+end