export Tria
mutable struct Tria <: AbstractElement
    id::String
    hid::Int

    node1::Node
    node2::Node
    node3::Node
    material::Material
    t::Float64
    t2::Float64
    t3::Float64

    elm_type::String
    mass_type::String

    center::Array{Float64}
    A::Float64
    T::Matrix{Float64}
end

function Tria(id,hid,node1,node2,node3,material,t,t2=0,t3=0,elm_type="TMT",mass_type="concentrate")
    #Initialize local CSys,could be optimized by using a MSE plane
    o=(node1.loc+node2.loc+node3.loc)/3
    pt1 = (node1.loc+node2.loc)/2
    pt2 = (node2.loc+node3.loc)/2
    csys = CSys(o, pt1, pt2)

    #area is considered as the average of trangles generated by splitting the quand with diagonals
    A=det([[1;1;1] node2.loc-node1.loc node3.loc-node1.loc]')/2

    E=material.E
    ν=material.ν

    #3D to local 2D
    T=zeros(18,18)
    for i in 1:6
        T[3i-2:3i,3i-2:3i]=csys.T
    end

    Tria(id,hid,node1,node2,node3,material,t,t2,t3,elm_type,mass_type,o,A,T)
end

for (root,dirs,files) in walkdir(joinpath(@__DIR__,"trias"))
    for file in files
        if file[end-2:end]==".jl"
            include(joinpath(@__DIR__,"trias",file))
        end
    end
end

function integrateK(elm::Tria)
    if elm.elm_type=="GT9"
        return K_GT9(elm)
    elseif elm.elm_type=="TMT"
        return K_TMT(elm)
    # elseif elm.elm_type=="DKT"
    #     K=K_DKT(elm)
    # elseif elm.elm_type=="DKGT"
    #     K=K_DKGT(elm)
    elseif elm.elm_type=="TMGT"
        return K_TMGT(elm)
    else
        throw("Quad element type error!")
    end
end

#WIP
function integrateKσ(elm::Tria,σ::Vector{Float64})
    E₀,ν₀=elm.material.E,elm.material.ν
    center=elm.center
    t=elm.t
    T=elm.T[1:3,1:3]
    x₁,y₁,z₁=T*(elm.node1.loc-center)
    x₂,y₂,z₂=T*(elm.node2.loc-center)
    x₃,y₃,z₃=T*(elm.node3.loc-center)
    K=Matrix{Float64}(undef,8,8)
    J=Matrix{Float64}(undef,2,2)

    o=(node1.loc+node2.loc+node3.loc+node4.loc)/4
    n1=node1.loc+[u[1];u[2];0]
    n2=node2.loc+[u[3];u[4];0]
    n3=node3.loc+[u[5];u[6];0]
    pt1 = n1+n2
    pt2 = n2+n3
    csys = CSys(o, pt1, pt2)
    V=csys.T
    T̄ᵉ=zeros(8,8)
    T̄ᵉ[1:2,1:2]=V
    T̄ᵉ[3:4,3:4]=V
    T̄ᵉ[5:6,5:6]=V
    T̄ᵉ[7:8,7:8]=V
end

#WIP
function integrateKu(elm::Tria,u)
    E₀,ν₀=elm.material.E,elm.material.ν
    center=elm.center
    t=elm.t
    T=elm.T[1:3,1:3]
    x₁,y₁,z₁=T*(elm.node1.loc-center)
    x₂,y₂,z₂=T*(elm.node2.loc-center)
    x₃,y₃,z₃=T*(elm.node3.loc-center)
    x₄,y₄,z₄=T*(elm.node4.loc-center)
    K=Matrix{Float64}(undef,8,8)
    J=Matrix{Float64}(undef,2,2)

    o=(node1.loc+node2.loc+node3.loc+node4.loc)/4
    n1=node1.loc+[u[1];u[2];0]
    n2=node2.loc+[u[3];u[4];0]
    n3=node3.loc+[u[5];u[6];0]
    n4=node4.loc+[u[7];u[8];0]
    pt1 = n1+n2
    pt2 = n2+n3
    csys = CSys(o, pt1, pt2)
    V=csys.T
    T̄ᵉ=zeros(8,8)
    T̄ᵉ[1:2,1:2]=V
    T̄ᵉ[3:4,3:4]=V
    T̄ᵉ[5:6,5:6]=V
    T̄ᵉ[7:8,7:8]=V
end

function integrateM!(elm::Tria)
    membrane,bending=elm.membrane,elm.bending
    ρ=elm.material.ρ
    center=elm.center
    t=elm.t
    A=elm.A
    T=elm.T[1:3,1:3]
    x₁,y₁,z₁=T*(elm.node1.loc-center)
    x₂,y₂,z₂=T*(elm.node2.loc-center)
    x₃,y₃,z₃=T*(elm.node3.loc-center)
    # M=Matrix{Float64}(ρ*t*A/18,18,18)
    M=spzeros(18,18)
    elm.Mᵉ=M
end

function integrateP(elm::Tria,elm_force)
    N=elm.N
    B=elm.B
    D=elm.D
    J=elm.J

    f,s,σ₀,ϵ₀=elmforce.f,elmforce.s,elmforce.σ₀,elmforce.ϵ₀
    a=[-1 -1]
    b=[1 1]

    function f1(x)
        N'*f*det(J)
    end

    function f2(x)
        N'*s*det(J)
    end

    function f3(x)
        x->B'*σ₀*det(J)
    end

    function f4(x)
        x->B'*D*ϵ₀*det(J)
    end

    Pᵉ_f=hquadrature(f1,a,b)
    Pᵉ_s=hquadrature(f2,[-1],[1])#ξ=-1
    Pᵉ_σ₀=hquadrature(f3,a,b)
    Pᵉ_ϵ₀=hquadrature(f4,a,b)

    Pᵉ_f,Pᵉ_s,Pᵉ_σ₀,Pᵉ_ϵ₀=calc_force(elm,elm_force) #Pᵉ=Pᵉf+Pᵉs+Pᵉσ₀+Pᵉϵ₀
    Pᵉ=Pᵉ_f+Pᵉ_s+Pᵉ_σ₀+Pᵉ_ϵ₀
end

# end
