#====================================================================
  Copyright (c) 2020 Samadrita Karmakar samadritakarmakar@gmail.com

  This Source Code Form is subject to the terms of the Mozilla Public
  License, v. 2.0. If a copy of the MPL was not distributed with this
  file, You can obtain one at http://mozilla.org/MPL/2.0/.
 =====================================================================#
 
function lagrangeHexOrder1_gmsh(ξ_Array::Array{T})::Array{T} where T
    ξ::T = ξ_Array[1]
    η::T =  ξ_Array[2]
    ζ::T = ξ_Array[3]
    ξ_η_ζ::Array{T,2} = [ -1.0  -1.0  -1.0;
  1.0  -1.0  -1.0;
  1.0   1.0  -1.0;
 -1.0   1.0  -1.0;
 -1.0  -1.0   1.0;
  1.0  -1.0   1.0;
  1.0   1.0   1.0;
 -1.0   1.0   1.0]
    ϕ::Array{T}= Array{T}(undef, size(ξ_η_ζ, 1))
    for a ∈ 1:size(ξ_η_ζ, 1)
        ξₐ::T = ξ_η_ζ[a,1]
        ηₐ::T = ξ_η_ζ[a,2]
        ζₐ::T = ξ_η_ζ[a,3]
        ϕ[a] = (0.125*(1+ ξₐ*ξ)*(1+ηₐ*η)*(1+ζₐ*ζ))
    end
    return ϕ
end

function lagrangeHexOrder2_gmsh(ξ_Array::Array{T})::Array{T} where T
    ξ::T = ξ_Array[1]
    η::T =  ξ_Array[2]
    ζ::T = ξ_Array[3]
    ξ_η_ζ::Array{T,2} = [-1.0 -1.0 -1.0
    1.0 -1.0 -1.0;
    1.0 1.0 -1.0;
    -1.0 1.0 -1.0;
    -1.0 -1.0 1.0;
    1.0 -1.0 1.0;
    1.0 1.0 1.0;
    -1.0 1.0 1.0;
    0.0 -1.0 -1.0;
    -1.0 0.0 -1.0;
    -1.0 -1.0 0.0;
    1.0 0.0 -1.0;
    1.0 -1.0 0.0;
    0.0 1.0 -1.0;
    1.0 1.0 0.0;
    -1.0 1.0 0.0;
    0.0 -1.0 1.0;
    -1.0 0.0 1.0;
    1.0 0.0 1.0;
    0.0 1.0 1.0;
    0.0 0.0 -1.0;
    0.0 -1.0 0.0;
    -1.0 0.0 0.0;
    1.0 0.0 0.0;
    0.0 1.0 0.0;
    0.0 0.0 1.0;
    0.0 0.0 0.0]
    ϕ::Array{T}= Array{T}(undef, size(ξ_η_ζ, 1))
    for a ∈ 1:8
        ξₐ::T = ξ_η_ζ[a,1]
        ηₐ::T = ξ_η_ζ[a,2]
        ζₐ::T = ξ_η_ζ[a,3]
        ϕ[a] = (0.125*(1+ ξₐ*ξ)*(1+ηₐ*η)*(1+ζₐ*ζ)*ξₐ*ξ*ηₐ*η*ζₐ*ζ)
    end
    for a ∈ [9 14 17 20]
        ξₐ::T = ξ_η_ζ[a,1]
        ηₐ::T = ξ_η_ζ[a,2]
        ζₐ::T = ξ_η_ζ[a,3]
        ϕ[a] = 0.25*(1-ξ^2)*(1+ηₐ*η)*(1+ζₐ*ζ)*ηₐ*η*ζₐ*ζ
    end
    for a ∈ [10 12 18 19]
        ξₐ::T = ξ_η_ζ[a,1]
        ηₐ::T = ξ_η_ζ[a,2]
        ζₐ::T = ξ_η_ζ[a,3]
        ϕ[a] = 0.25*(1+ξₐ*ξ)*(1-η^2)*(1+ζₐ*ζ)*ξₐ*ξ*ζₐ*ζ
    end
    for a ∈ [11 13 15 16]
        ξₐ::T = ξ_η_ζ[a,1]
        ηₐ::T = ξ_η_ζ[a,2]
        ζₐ::T = ξ_η_ζ[a,3]
        ϕ[a] = 0.25*(1+ξₐ*ξ)*(1+ηₐ*η)*(1-ζ^2)*ξₐ*ξ*ηₐ*η
    end
    for a ∈ [21 26]
        ξₐ::T = ξ_η_ζ[a,1]
        ηₐ::T = ξ_η_ζ[a,2]
        ζₐ::T = ξ_η_ζ[a,3]
        ϕ[a] = 0.5*(1-ξ^2)*(1-η^2)*(1+ζₐ*ζ)*ζₐ*ζ
    end
    for a ∈ [22 25]
        ξₐ::T = ξ_η_ζ[a,1]
        ηₐ::T = ξ_η_ζ[a,2]
        ζₐ::T = ξ_η_ζ[a,3]
        ϕ[a] = 0.5*(1-ξ^2)*(1+ηₐ*η)*(1-ζ^2)*ηₐ*η
    end
    for a ∈ [23 24]
        ξₐ::T = ξ_η_ζ[a,1]
        ηₐ::T = ξ_η_ζ[a,2]
        ζₐ::T = ξ_η_ζ[a,3]
        ϕ[a] = 0.5*(1+ξₐ*ξ)*(1-η^2)*(1-ζ^2)*ξₐ*ξ
    end
    ϕ[27] = (1-ξ^2)*(1-η^2)*(1-ζ^2)
    return ϕ
end

function lagrangeHexOrder3_gmsh(ξ_Array::Array{T})::Array{T} where T
    ξ::T = ξ_Array[1]
    η::T =  ξ_Array[2]
    ζ::T = ξ_Array[3]
    ξ_η_ζ::Array{T,2} = [-1.0		-1.0		-1.0;
1.0		-1.0		-1.0;
1.0		1.0		    -1.0;
-1.0	1.0	       	-1.0;
-1.0	-1.0		1.0;
1.0		-1.0		1.0;
1.0		1.0		     1.0;
-1.0	1.0		      1.0;
-1.0/3.0 -1.0		-1.0;
1.0/3.0	-1.0		-1.0;
-1.0	-1.0/3.0	-1.0;
-1.0	1.0/3.0		-1.0;
-1.0	-1.0		-1.0/3.0;
-1.0	-1.0		1.0/3.0;
1.0		-1.0/3.0	-1.0;
1.0		1.0/3.0		-1.0;
1.0		-1.0		-1.0/3.0;
1.0		-1.0		1.0/3.0;
1.0/3.0		1.0		-1.0;
-1.0/3.0	1.0		-1.0;
1.0		   1.0		-1.0/3.0;
1.0		   1.0		1.0/3.0;
-1.0		1.0		-1.0/3.0;
-1.0		1.0		1.0/3.0;
-1.0/3.0	-1.0		1.0;
1.0/3.0		-1.0		1.0;
-1.0		-1.0/3.0	1.0;
-1.0		1.0/3.0		1.0;
1.0		   -1.0/3.0		1.0;
1.0		   1.0/3.0		1.0;
1.0/3.0		1.0		    1.0;
-1.0/3.0	1.0		     1.0;
-1.0/3.0	-1.0/3.0	-1.0;
-1.0/3.0	1.0/3.0		-1.0;
1.0/3.0		1.0/3.0		-1.0;
1.0/3.0		-1.0/3.0	-1.0;
-1.0/3.0		-1.0	-1.0/3.0;
1.0/3.0		-1.0		-1.0/3.0;
1.0/3.0		-1.0		1.0/3.0;
-1.0/3.0		-1.0	1.0/3.0;
-1.0		-1.0/3.0	-1.0/3.0;
-1.0		-1.0/3.0	1.0/3.0;
-1.0		1.0/3.0		1.0/3.0;
-1.0		1.0/3.0		-1.0/3.0;
1.0		-1.0/3.0		-1.0/3.0;
1.0		1.0/3.0		     -1.0/3.0;
1.0		1.0/3.0		     1.0/3.0;
1.0		-1.0/3.0		 1.0/3.0;
1.0/3.0		1.0		     -1.0/3.0;
-1.0/3.0	1.0		      -1.0/3.0;
-1.0/3.0	1.0		      1.0/3.0;
1.0/3.0		1.0		      1.0/3.0;
-1.0/3.0	-1.0/3.0		1.0;
1.0/3.0		-1.0/3.0		1.0;
1.0/3.0		 1.0/3.0		1.0;
-1.0/3.0	   1.0/3.0		1.0;
-1.0/3.0	-1.0/3.0	  -1.0/3.0;
1.0/3.0		-1.0/3.0	  -1.0/3.0;
1.0/3.0		1.0/3.0		  -1.0/3.0;
-1.0/3.0	1.0/3.0		  -1.0/3.0;
-1.0/3.0	-1.0/3.0	   1.0/3.0;
1.0/3.0		-1.0/3.0	   1.0/3.0;
1.0/3.0		1.0/3.0		   1.0/3.0;
-1.0/3.0	1.0/3.0		   1.0/3.0]
    ϕ::Array{T}= Array{T}(undef, size(ξ_η_ζ, 1))
    ϕ = zeros(size(ξ_η_ζ, 1))
    for a ∈ 1:8
        ξₐ::T = ξ_η_ζ[a,1]
        ηₐ::T = ξ_η_ζ[a,2]
        ζₐ::T = ξ_η_ζ[a,3]
        ϕ[a] = -(729.0/4096.0*(1.0+ξₐ*ξ)*(1.0+ηₐ*η)*(1.0/9.0-ξ^2)*(1.0/9.0-η^2)*(1.0+ζₐ*ζ)*(1.0/9.0-ζ^2))
    end
    for a ∈ [9 10 19 20 25 26 31 32]
        ξₐ::T = ξ_η_ζ[a,1]
        ηₐ::T = ξ_η_ζ[a,2]
        ζₐ::T = ξ_η_ζ[a,3]
        ϕ[a] = (2187/4096*(1.0/3.0+3.0*ξₐ*ξ)*(1.0-ξ^2)*(1.0+ηₐ*η)*(1.0/9.0-η^2)*(1.0+ζₐ*ζ)*(1.0/9.0-ζ^2))
    end
    for a ∈ [11 12 15 16 27 28 29 30]
        ξₐ::T = ξ_η_ζ[a,1]
        ηₐ::T = ξ_η_ζ[a,2]
        ζₐ::T = ξ_η_ζ[a,3]
        ϕ[a] = (2187/4096*(1.0+ξₐ*ξ)*(1.0/9.0-ξ^2)*(1.0/3.0+3.0*ηₐ*η)*(1.0-η^2)*(1.0+ζₐ*ζ)*(1.0/9.0-ζ^2))
    end
    for a ∈ [13 14 17 18 21 22 23 24]
        ξₐ::T = ξ_η_ζ[a,1]
        ηₐ::T = ξ_η_ζ[a,2]
        ζₐ::T = ξ_η_ζ[a,3]
        ϕ[a] = (2187/4096*(1.0+ξₐ*ξ)*(1.0/9.0-ξ^2)*(1.0+ηₐ*η)*(1.0/9.0-η^2)*(1.0/3.0+3.0*ζₐ*ζ)*(1.0-ζ^2))
    end
    for a ∈ [33 34 35 36 53 54 55 56]#[33 34 35 36 49 50 51 52]
        ξₐ::T = ξ_η_ζ[a,1]
        ηₐ::T = ξ_η_ζ[a,2]
        ζₐ::T = ξ_η_ζ[a,3]
        ϕ[a] = -(6561/4096*(1.0/3.0+3.0*ξₐ*ξ)*(1.0-ξ^2)*(1.0/3.0+3.0*ηₐ*η)*(1.0-η^2)*(1.0+ζₐ*ζ)*(1.0/9.0-ζ^2))
    end
    for a ∈ [37 38 39 40 49 50 51 52]#[37 38 39 40 53 54 55 56]
        ξₐ::T = ξ_η_ζ[a,1]
        ηₐ::T = ξ_η_ζ[a,2]
        ζₐ::T = ξ_η_ζ[a,3]
        ϕ[a] = -(6561/4096*(1.0/3.0+3.0*ξₐ*ξ)*(1.0-ξ^2)*(1.0+ηₐ*η)*(1.0/9.0-η^2)*(1.0/3.0+3.0*ζₐ*ζ)*(1.0-ζ^2))
    end
    for a ∈ 41:48
        ξₐ::T = ξ_η_ζ[a,1]
        ηₐ::T = ξ_η_ζ[a,2]
        ζₐ::T = ξ_η_ζ[a,3]
        ϕ[a] = -(6561/4096*(1.0+ξₐ*ξ)*(1.0/9.0-ξ^2)*(1.0/3.0+3.0*ηₐ*η)*(1.0-η^2)*(1.0/3.0+3.0*ζₐ*ζ)*(1.0-ζ^2))
    end
    for a ∈ 57:64
        ξₐ::T = ξ_η_ζ[a,1]
        ηₐ::T = ξ_η_ζ[a,2]
        ζₐ::T = ξ_η_ζ[a,3]
        ϕ[a] = (19683/4096*(1.0/3.0+3.0*ξₐ*ξ)*(1.0-ξ^2)*(1.0/3.0+3.0*ηₐ*η)*(1.0-η^2)*(1.0/3.0+3.0*ζₐ*ζ)*(1.0-ζ^2))
    end
    return ϕ
end

function checkLagrangeHex()
    ξ_η_ζ::Array{T,2} = [-1.0 -1.0 -1.0;
    1.0 -1.0 -1.0;
    1.0 1.0 -1.0;
    -1.0 1.0 -1.0;
    -1.0 -1.0 1.0;
    1.0 -1.0 1.0;
    1.0 1.0 1.0;
    -1.0 1.0 1.0;
    -1.0/3.0 -1.0 -1.0;
    1.0/3.0 -1.0 -1.0;
    -1.0 -1.0/3.0 -1.0;
    -1.0 1.0/3.0 -1.0;
    -1.0 -1.0 -1.0/3.0;
    -1.0 -1.0 1.0/3.0;
    1.0 -1.0/3.0 -1.0;
    1.0 1.0/3.0 -1.0;
    1.0 -1.0 -1.0/3.0;
    1.0 -1.0 1.0/3.0;
    1.0/3.0 1.0 -1.0;
    -1.0/3.0 1.0 -1.0;
    1.0 1.0 -1.0/3.0;
    1.0 1.0 1.0/3.0;
    -1.0 1.0 -1.0/3.0;
    -1.0 1.0 1.0/3.0;
    -1.0/3.0 -1.0 1.0;
    1.0/3.0 -1.0 1.0;
    -1.0 -1.0/3.0 1.0;
    -1.0 1.0/3.0 1.0;
    1.0 -1.0/3.0 1.0;
    1.0 1.0/3.0 1.0;
    1.0/3.0 1.0 1.0;
    -1.0/3.0 1.0 1.0;
    -1.0/3.0 -1.0/3.0 -1.0;
    -1.0/3.0 1.0/3.0 -1.0;
    1.0/3.0 1.0/3.0 -1.0;
    1.0/3.0 -1.0/3.0 -1.0;
    -1.0/3.0 -1.0 -1.0/3.0;
    1.0/3.0 -1.0 -1.0/3.0;
    1.0/3.0 -1.0 1.0/3.0;
    -1.0/3.0 -1.0 1.0/3.0;
    -1.0 -1.0/3.0 -1.0/3.0;
    -1.0 -1.0/3.0 1.0/3.0;
    -1.0 1.0/3.0 1.0/3.0;
    -1.0 1.0/3.0 -1.0/3.0;
    1.0 -1.0/3.0 -1.0/3.0;
    1.0 1.0/3.0 -1.0/3.0;
    1.0 1.0/3.0 1.0/3.0;
    1.0 -1.0/3.0 1.0/3.0;
    1.0/3.0 1.0 -1.0/3.0;
    -1.0/3.0 1.0 -1.0/3.0;
    -1.0/3.0 1.0 1.0/3.0;
    1.0/3.0 1.0 1.0/3.0;
    -1.0/3.0 -1.0/3.0 1.0;
    1.0/3.0 -1.0/3.0 1.0;
    1.0/3.0 1.0/3.0 1.0;
    -1.0/3.0 1.0/3.0 1.0;
    -1.0/3.0 -1.0/3.0 -1.0/3.0;
    1.0/3.0 -1.0/3.0 -1.0/3.0;
    1.0/3.0 1.0/3.0 -1.0/3.0;
    -1.0/3.0 1.0/3.0 -1.0/3.0;
    -1.0/3.0 -1.0/3.0 1.0/3.0;
    1.0/3.0 -1.0/3.0 1.0/3.0;
    1.0/3.0 1.0/3.0 1.0/3.0;
    -1.0/3.0 1.0/3.0 1.0/3.0]
    for pnt ∈ 1:size(ξ_η_ζ,1)
        ξ=ξ_η_ζ[pnt,1]
        η=ξ_η_ζ[pnt,2]
        ζ = ξ_η_ζ[pnt,3]
        ϕ = lagrangeHexOrder3_gmsh([ξ, η, ζ])
        println(pnt)
        println("ξ= ", ξ, " η = ", η, " ζ = ", ζ)
        for i ∈ 1:length(ϕ)
            print(ϕ[i], " ")
        end
        println()
        println()
    end
end
